Feature eigen

Context

This MR generalizes the use of Eigen in waves.

Tasks

Following MR !50 (merged), here is the updated to-do list:

• replace gmm::matrix<double> by Eigen::Matrix<double>
• replace std::vector<double> by Eigen::(Row)Vector<double>
• refactor elementary matrices computation
• change method return type, ie void build(Matrix const &A, Matrix &K) should become Matrix K build(Matrix const &A)
• replace tbox::Pt by Eigen::VectorXd
• investigate Pardiso interface
• investigate Eigen multithreading feature
• clean the thermo-mechanical matrices

Notes

• I added a general interface tbox::LinearSolver from which derive tbox::Mumps, tbox::Pardiso, tbox::Gmres and tbox::SparseLu. The first class is a direct interface to the MUMPS solver, while the other three merely call Eigen interface or implementation of these solvers. The linear solvers are interfaced in python and the nonlinear solvers must now be initialized with a problem and a linear solver.
• I did not clean the thermo-mechanical matrices computation in tbox::Tetra4 and tbox::Hex8. This could be further addressed in the branch clean_mirrors, which brings LFS support for large geometry files and improved tests.
• Regarding the multi-threading in Eigen, it does not make sense to use it if the top level functions in waves are multi-threaded, which should be the case. However, multi-threading can be used for linear solvers. For now, only Pardiso can use multi-threading (GMRES and SparseLU are not multi-threaded and I did not look into MUMPS).

Additional changes

• I replaced the std::vector holding the shape functions by Eigen::Matrices, and consequently updated the Cache, ShapeFunction and Gauss classes. Additionally, the members of Cache and Gauss are now private and can only be accessed by getters (which I inlined).
• I replaced the std::vector by std::deque to hold the list of Eigen::Triplet for the global matrices assembly, as it proved to be faster.
• I changed the assembly procedure of periodic BCs (tbox::MshDeform) and Kutta condition (flow::Newton and flow::Picard) so that the contributions are directly assembled on the right rows, thus avoiding the costly manipulation of the rows after the matrix has been built.

Tests

Passed on ubuntu18.04 (python 2.7.15), debian4.9 (python 2.7.13), msys2 (python 3.8.1).

mrstlnos/src/Stress_computation.h

@@ -20,7 +20,7 @@
 void strain_at_gp(size_t n_n, Mem &mem, Cache &cache, size_t a, const Eigen::MatrixXd &u_e, Eigen::MatrixXd &Strain)
 {
     Eigen::MatrixXd const &J = mem.getJinv(a);
-    Eigen::MatrixXd const &dff = cache.getDsf(a);
+    Eigen::MatrixXd const &dffi = cache.getDsf(a);
     Eigen::Matrix3d du_d_Phi = Eigen::Matrix3d::Zero();
     for (size_t i = 0; i < n_n; ++i)
     {
@@ -56,7 +56,7 @@ void strain_stress_x_y_z(Element *elem, Eigen::MatrixXd &Strain, Eigen::MatrixXd
     Cache &cache = elem->getVCache();
     Mem &mem = elem->getVMem();
 
-    size_t n_gp = cache.dff.size();
+    size_t n_gp = mem.getGPn();
     size_t n_n = elem->nodes.size();
 
     Eigen::MatrixXd u_e(3, n_n);
@@ -83,7 +83,7 @@ void strain_stress_x_y_z_at_node(Element *elem, std::vector<Eigen::MatrixXd> &St
     Cache &cache = elem->getVCache();
     Mem &mem = elem->getVMem();
 
-    size_t n_gp = cache.dff.size();
+    size_t n_gp = mem.getGPn();
     size_t n_n = elem->nodes.size();
 
     std::vector<Eigen::MatrixXd> Strain_at_gp(n_gp);