Reorganize material models with Eigen (#318)

* Consolidate get_pk2cc and get_pk2cc_dev(). Use only mat_models.cpp

* Update material unit tests

* Add voigt_to_cc, fix undefined Ja bug, fix no pressure in HO-ma bug, update HO-ma reference solutions, remove mat_models_carray.h includes

* Fix Ja calculation bug leading to ustruct failures, remove get_pk2cc_dev

* Fix indexing bug for struct viscosity tangent terms

* Remove mat_models_carray.h and commented code

* trying out eigen conversion

* Change from Matrix to Tensor

* Adding/moving templated Eigen matrix/tensor function implementations to header; using fixed size Matrix and Tensors, instead of just Tensors. Compiles and passes integration tests/unit tests that use NeoHookean

* Update eigen includes path

* Consolidate volumetric stress and elasticity contributions

* General function used by multiple models to project S_bar and CC_bar to S_iso and CC_iso; Add several required tensor operations for Eigen Matrices/Tensors

* Bug fixes. 1) For some reason, code runs but is incorrect when using auto instead of EigenTensor<nsd> for PP. 2) fix typo for ten_symm_prod_eigen()

* Replace some mat_fun loops with built in Eigen operations

* Add Eigen implementation for Mooney-Rivlin model

* Adding Linear, StVK, and modified StVK. Compared to original implementation, only modifying to use Eigen tensor functions.

* Add HGO model. Passes ustruct/tensile_adventitia_HGO

* Add Guccione model Eigen implementation. Passes Guccione tests

* Add Holzapfel-Ogden model. Passes HO tests

* Add HO-ma model. Passes all struct/ustruct tests

* Replacing auto with EigenMatrix<nsd>, since using auto is not recommended for Eigen

* Replacing Eigen operations with for loops for some tensor functions, faster for some reason

* Update .gitignore genBC files to reflect renaming of svFSIplus to svMultiPhysics

* Add -march=native flag to CMAKE_C_FLAGS and CMAKE_CXX_FLAGS for optimization

* Update CMakeLists.txt to clarify -march=native flag usage for Eigen performance

* Remove temporary mat_models_Dave_array.cpp

* Removing mat_fun_carray.h/cpp. Restructuring material unit tests to use Arrays instead of C-arrays.

* Comment out -march=native flag in CMakeLists.txt for C and C++ compiler flags to try to fix Ubuntu timeout issue

* Cleaning up and renaming to address @ktbolt and @mrp089 comments

* Replace dyadic_product and symmetric_dyadic_product with algebraic expressions. Avoid overloaded function and improve readability

* Addressing @ktbolt and @mrp089 comments round 2. Rename EigenMatrix and EigenTensor aliases to Matrix and Tensor. Consolidate double_dot_product functions for tensors. Rename get_() functions to compute_(). Fix throw runtime_error bug for invalid material models.

---------

Co-authored-by: Marisa Bazzi <[email protected]>

.gitignore

@@ -23,7 +23,7 @@ Documentation/html/
 
 # genBC files
 **/genBC/obj/
-**/genBC_svFSIplus/obj/
+**/genBC_svMultiPhysics/obj/
 genBC.exe
 AllData
 GenBC.int

CMakeLists.txt

@@ -61,6 +61,10 @@ message(STATUS "SV_TOP: ${SV_TOP}")
 # CMake module path inside of true simvascular source
 set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/Code/CMake"
 	"${CMAKE_MODULE_PATH}")
+
+# Add -march=native to CMAKE_C_FLAGS and CMAKE_CXX_FLAGS for Eigen performance
+#set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -march=native")
+#set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -march=native")
 #-----------------------------------------------------------------------------
 
 #-----------------------------------------------------------------------------

Code/Source/solver/CMakeLists.txt

@@ -186,7 +186,6 @@ set(CSRCS
   ls.h ls.cpp
   main.cpp
   mat_fun.h mat_fun.cpp
-  mat_fun_carray.h mat_fun_carray.cpp
   mat_models.h mat_models.cpp
   mesh.h mesh.cpp
   nn.h nn.cpp

Code/Source/solver/mat_fun.cpp

@@ -30,12 +30,8 @@
 
 #include "mat_fun.h"
 
-/*
-#include <eigen3/Eigen/Core>
-#include <eigen3/Eigen/Dense>
-*/
-
 #include "consts.h"
+
 #include "utils.h"
 
 #include <math.h>

Code/Source/solver/mat_fun.h

@@ -30,6 +30,10 @@
 
 #ifndef MAT_FUN_H 
 #define MAT_FUN_H 
+#include "eigen3/Eigen/Core"
+#include "eigen3/Eigen/Dense"
+#include "eigen3/unsupported/Eigen/CXX11/Tensor"
+#include <stdexcept>
 
 #include "Array.h"
 #include "Tensor4.h"
@@ -43,7 +47,59 @@
 /// \todo [TODO:DaveP] this should just be a namespace?
 //
 namespace mat_fun {
+    // Define templated type aliases for Eigen matrices and tensors for convenience
+    template<size_t nsd>
+    using Matrix = Eigen::Matrix<double, nsd, nsd>;
+
+    template<size_t nsd>
+    using Tensor = Eigen::TensorFixedSize<double, Eigen::Sizes<nsd, nsd, nsd, nsd>>;
+
+    // Function to convert Array<double> to Eigen::Matrix
+    template <typename MatrixType>
+    MatrixType convert_to_eigen_matrix(const Array<double>& src) {
+        MatrixType mat;
+        for (int i = 0; i < mat.rows(); ++i)
+            for (int j = 0; j < mat.cols(); ++j)
+                mat(i, j) = src(i, j);
+        return mat;
+    }
+
+    // Function to convert Eigen::Matrix to Array<double>
+    template <typename MatrixType>
+    void convert_to_array(const MatrixType& mat, Array<double>& dest) {
+        for (int i = 0; i < mat.rows(); ++i)
+            for (int j = 0; j < mat.cols(); ++j)
+                dest(i, j) = mat(i, j);
+    }
+
+    // Function to convert a higher-dimensional array like Dm
+    template <typename MatrixType>
+    void copy_Dm(const MatrixType& mat, Array<double>& dest, int rows, int cols) {
+        for (int i = 0; i < rows; ++i)
+            for (int j = 0; j < cols; ++j)
+                dest(i, j) = mat(i, j);
+    }
+
+    template <int nsd>
+    Eigen::Matrix<double, nsd, 1> cross_product(const Eigen::Matrix<double, nsd, 1>& u, const Eigen::Matrix<double, nsd, 1>& v) {
+        if constexpr (nsd == 2) {
+            return Eigen::Matrix<double, 2, 1>(v(1), - v(0));
+        }
+        else if constexpr (nsd == 3) {
+            return u.cross(v);
+        }
+        else {
+            throw std::runtime_error("[cross_product] Invalid number of spatial dimensions '" + std::to_string(nsd) + "'. Valid dimensions are 2 or 3.");
+        }
+    }
+
     double mat_ddot(const Array<double>& A, const Array<double>& B, const int nd);
+
+    template <int nsd>
+    double double_dot_product(const Matrix<nsd>& A, const Matrix<nsd>& B) {
+        return A.cwiseProduct(B).sum();
+    }
+
     double mat_det(const Array<double>& A, const int nd);
     Array<double> mat_dev(const Array<double>& A, const int nd);
 
@@ -61,20 +117,157 @@ namespace mat_fun {
 
     Array<double> mat_symm(const Array<double>& A, const int nd);
     Array<double> mat_symm_prod(const Vector<double>& u, const Vector<double>& v, const int nd);
+
     double mat_trace(const Array<double>& A, const int nd);
 
     Tensor4<double> ten_asym_prod12(const Array<double>& A, const Array<double>& B, const int nd);
     Tensor4<double> ten_ddot(const Tensor4<double>& A, const Tensor4<double>& B, const int nd);
     Tensor4<double> ten_ddot_2412(const Tensor4<double>& A, const Tensor4<double>& B, const int nd);
     Tensor4<double> ten_ddot_3424(const Tensor4<double>& A, const Tensor4<double>& B, const int nd);
 
+    /**
+     * @brief Contracts two 4th order tensors A and B over two dimensions, 
+     * 
+     */
+    template <int nsd>
+    Tensor<nsd>
+    double_dot_product(const Tensor<nsd>& A, const std::array<int, 2>& dimsA, 
+                        const Tensor<nsd>& B, const std::array<int, 2>& dimsB) {
+
+        // Define the contraction dimensions
+        Eigen::array<Eigen::IndexPair<int>, 2> contractionDims = {
+            Eigen::IndexPair<int>(dimsA[0], dimsB[0]), // Contract A's dimsA[0] with B's dimsB[0]
+            Eigen::IndexPair<int>(dimsA[1], dimsB[1])  // Contract A's dimsA[1] with B's dimsB[1]
+        };
+
+        // Return the double dot product
+        return A.contract(B, contractionDims);
+
+        // For some reason, in this case the Eigen::Tensor contract function is
+        // faster than a for loop implementation.
+    }
+
     Tensor4<double> ten_dyad_prod(const Array<double>& A, const Array<double>& B, const int nd);
+
+    /**
+     * @brief Compute the dyadic product of two 2nd order tensors A and B, C_ijkl = A_ij * B_kl
+     * 
+     * @tparam nsd, the number of spatial dimensions
+     * @param A, the first 2nd order tensor
+     * @param B, the second 2nd order tensor
+     * @return Tensor<nsd>
+     */
+    template <int nsd>
+    Tensor<nsd> 
+    dyadic_product(const Matrix<nsd>& A, const Matrix<nsd>& B) {
+        // Initialize the result tensor
+        Tensor<nsd> C;
+
+        // Compute the dyadic product: C_ijkl = A_ij * B_kl
+        for (int i = 0; i < nsd; ++i) {
+            for (int j = 0; j < nsd; ++j) {
+                for (int k = 0; k < nsd; ++k) {
+                    for (int l = 0; l < nsd; ++l) {
+                        C(i,j,k,l) = A(i,j) * B(k,l);
+                    }
+                }
+            }
+        }
+        // For some reason, in this case the Eigen::Tensor contract function is 
+        // slower than the for loop implementation
+
+        return C;
+    }
+
     Tensor4<double> ten_ids(const int nd);
+
+    /**
+     * @brief Create a 4th order identity tensor:
+     * I_ijkl = 0.5 * (δ_ik * δ_jl + δ_il * δ_jk)
+     * 
+     * @tparam nsd, the number of spatial dimensions
+     * @return Tensor<nsd> 
+     */
+    template <int nsd>
+    Tensor<nsd>
+    fourth_order_identity() {
+        // Initialize as zero
+        Tensor<nsd> I;
+        I.setZero();
+
+        // Set only non-zero entries
+        for (int i = 0; i < nsd; ++i) {
+            for (int j = 0; j < nsd; ++j) {
+                I(i,j,i,j) += 0.5;
+                I(i,j,j,i) += 0.5;
+            }
+        }
+
+        return I;
+    }
+
     Array<double> ten_mddot(const Tensor4<double>& A, const Array<double>& B, const int nd);
 
     Tensor4<double> ten_symm_prod(const Array<double>& A, const Array<double>& B, const int nd);
+
+    /// @brief Create a 4th order tensor from symmetric outer product of two matrices: C_ijkl = 0.5 * (A_ik * B_jl + A_il * B_jk)
+    ///
+    /// Reproduces 'FUNCTION TEN_SYMMPROD(A, B, nd) RESULT(C)'.
+    //
+    template <int nsd>
+    Tensor<nsd>
+    symmetric_dyadic_product(const Matrix<nsd>& A, const Matrix<nsd>& B) {
+
+        // Initialize the result tensor
+        Tensor<nsd> C;
+
+        // Compute the symmetric product: C_ijkl = 0.5 * (A_ik * B_jl + A_il * B_jk)
+        for (int i = 0; i < nsd; ++i) {
+            for (int j = 0; j < nsd; ++j) {
+                for (int k = 0; k < nsd; ++k) {
+                    for (int l = 0; l < nsd; ++l) {
+                        C(i,j,k,l) = 0.5 * (A(i,k) * B(j,l) + A(i,l) * B(j,k));
+                    }
+                }
+            }
+        }
+        // For some reason, in this case the for loop implementation is faster 
+        // than the Eigen::Tensor contract method
+
+        // Return the symmetric product
+        return C;
+    }
+
     Tensor4<double> ten_transpose(const Tensor4<double>& A, const int nd);
 
+    /**
+     * @brief Performs a tensor transpose operation on a 4th order tensor A, B_ijkl = A_klij
+     * 
+     * @tparam nsd, the number of spatial dimensions
+     * @param A, the input 4th order tensor
+     * @return Tensor<nsd>
+     */
+    template <int nsd>
+    Tensor<nsd>
+    transpose(const Tensor<nsd>& A) {
+
+        // Initialize the result tensor
+        Tensor<nsd> B;
+
+        // Permute the tensor indices to perform the transpose operation
+        for (int i = 0; i < nsd; ++i) {
+            for (int j = 0; j < nsd; ++j) {
+                for (int k = 0; k < nsd; ++k) {
+                    for (int l = 0; l < nsd; ++l) {
+                        B(i,j,k,l) = A(k,l,i,j);
+                    }
+                }
+            }
+        }
+
+        return B;
+    }
+
     Array<double> transpose(const Array<double>& A);
 
     void ten_init(const int nd);