Merge pull request #16 from geomechanics/feature/fast_matrix_assembly

[Feature] Fast matrix assembly

Author: bodhinandach

include/linear_solvers/assemblers/assembler_eigen_implicit_linear.tcc

@@ -17,9 +17,9 @@ bool mpm::AssemblerEigenImplicitLinear<Tdim>::assemble_stiffness_matrix() {
     stiffness_matrix_.resize(active_dof_ * Tdim, active_dof_ * Tdim);
     stiffness_matrix_.setZero();
 
-    // Reserve storage for sparse matrix
-    stiffness_matrix_.reserve(
-        Eigen::VectorXi::Constant(active_dof_ * Tdim, sparse_row_size_ * Tdim));
+    // Triplets and reserve storage for sparse matrix
+    std::vector<Eigen::Triplet<double>> tripletList;
+    tripletList.reserve(active_dof_ * Tdim * sparse_row_size_ * Tdim);
 
     // Cell pointer
     const auto& cells = mesh_->cells();
@@ -43,17 +43,23 @@ bool mpm::AssemblerEigenImplicitLinear<Tdim>::assemble_stiffness_matrix() {
           for (unsigned j = 0; j < nids.size(); ++j) {
             for (unsigned k = 0; k < Tdim; ++k) {
               for (unsigned l = 0; l < Tdim; ++l) {
-                stiffness_matrix_.coeffRef(
-                    nactive_node * k + global_node_indices_.at(cid)(i),
-                    nactive_node * l + global_node_indices_.at(cid)(j)) +=
-                    cell_stiffness(Tdim * i + k, Tdim * j + l);
+                if (std::abs(cell_stiffness(Tdim * i + k, Tdim * j + l)) >
+                    std::numeric_limits<double>::epsilon())
+                  tripletList.push_back(Eigen::Triplet<double>(
+                      nactive_node * k + global_node_indices_.at(cid)(i),
+                      nactive_node * l + global_node_indices_.at(cid)(j),
+                      cell_stiffness(Tdim * i + k, Tdim * j + l)));
               }
             }
           }
         }
         ++cid;
       }
     }
+
+    // Fast assembly from triplets
+    stiffness_matrix_.setFromTriplets(tripletList.begin(), tripletList.end());
+
   } catch (std::exception& exception) {
     console_->error("{} #{}: {}\n", __FILE__, __LINE__, exception.what());
     status = false;

include/solvers/mpm_implicit_linear.tcc

@@ -345,16 +345,16 @@ bool mpm::MPMImplicitLinear<Tdim>::compute_equilibrium_equation() {
         Tdim * assembler_->global_active_dof());
 
     // Prepare rank global mapper
-    std::vector<int> predictor_rgm;
+    std::vector<int> matrix_rgm;
     for (unsigned dir = 0; dir < Tdim; ++dir) {
       auto dir_rgm = assembler_->rank_global_mapper();
       std::for_each(dir_rgm.begin(), dir_rgm.end(),
                     [size = assembler_->global_active_dof(),
                      dir = dir](int& rgm) { rgm += dir * size; });
-      predictor_rgm.insert(predictor_rgm.end(), dir_rgm.begin(), dir_rgm.end());
+      matrix_rgm.insert(matrix_rgm.end(), dir_rgm.begin(), dir_rgm.end());
     }
     // Assign rank global mapper to solver
-    linear_solver_["displacement"]->assign_rank_global_mapper(predictor_rgm);
+    linear_solver_["displacement"]->assign_rank_global_mapper(matrix_rgm);
 #endif
 
     // Solve matrix equation and assign solution to assembler