Merge pull request #887 from chleh/sparsity-pattern

Multi-component sparsity pattern

AssemblerLib/ComputeSparsityPattern.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "MeshLib/NodeAdjacencyTable.h"
+#include "LocalToGlobalIndexMap.h"
+
+#include "ComputeSparsityPattern.h"
+
+namespace AssemblerLib
+{
+
+SparsityPattern
+computeSparsityPattern(LocalToGlobalIndexMap const& dof_table,
+        MeshLib::Mesh const& mesh
+        )
+{
+    MeshLib::NodeAdjacencyTable node_adjacency_table;
+    node_adjacency_table.createTable(mesh.getNodes());
+
+    // A mapping   mesh node id -> global indices
+    // It acts as a cache for dof table queries.
+    std::vector<std::vector<GlobalIndexType> > global_idcs;
+
+    global_idcs.reserve(mesh.getNNodes());
+    for (std::size_t n=0; n<mesh.getNNodes(); ++n)
+    {
+        MeshLib::Location l(mesh.getID(), MeshLib::MeshItemType::Node, n);
+        global_idcs.push_back(dof_table.getGlobalIndices(l));
+    }
+
+    SparsityPattern sparsity_pattern(dof_table.dofSize());
+
+    // Map adjacent mesh nodes to "adjacent global indices".
+    for (std::size_t n=0; n<mesh.getNNodes(); ++n)
+    {
+        auto const& node_ids = node_adjacency_table.getAdjacentNodes(n);
+        for (auto an : node_ids) {
+            auto const& row_ids = global_idcs[an];
+            auto const num_components = row_ids.size();
+            for (auto r : row_ids) {
+                // Each component leads to an entry in the row.
+                // For the sparsity pattern only the number of entries are needed.
+                sparsity_pattern[r] += num_components;
+            }
+        }
+    }
+
+    return sparsity_pattern;
+}
+
+}

AssemblerLib/ComputeSparsityPattern.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef ASSEMBLERLIB_COMPUTESPARSITYPATTERN_H
+#define ASSEMBLERLIB_COMPUTESPARSITYPATTERN_H
+
+#include <vector>
+
+#include "ProcessLib/NumericsConfig.h"
+
+namespace AssemblerLib
+{
+
+class LocalToGlobalIndexMap;
+
+/// A vector telling how many nonzeros there are in each global matrix row.
+using SparsityPattern = std::vector<GlobalIndexType>;
+
+/**
+ * @brief Computes a sparsity pattern for the given inputs.
+ *
+ * @param dof_table            maps mesh nodes to global indices
+ * @param mesh                 mesh for which the two parameters above are defined
+ *
+ * @return The computed sparsity pattern.
+ */
+SparsityPattern
+computeSparsityPattern(
+        LocalToGlobalIndexMap const& dof_table,
+        MeshLib::Mesh const& mesh
+        );
+}
+
+#endif // ASSEMBLERLIB_COMPUTESPARSITYPATTERN_H
+

AssemblerLib/LocalToGlobalIndexMap.h

@@ -88,6 +88,12 @@ public:
         return _mesh_component_map.getGlobalIndex(l, c);
     }
 
+    /// Forwards the respective method from MeshComponentMap.
+    std::vector<GlobalIndexType> getGlobalIndices(const MeshLib::Location &l) const
+    {
+        return _mesh_component_map.getGlobalIndices(l);
+    }
+
 private:
     /// Private constructor used by internally created local-to-global index
     /// maps. The mesh_component_map is passed as argument instead of being

MathLib/LinAlg/Eigen/EigenMatrix.h

@@ -197,13 +197,10 @@ void operator()(EigenMatrix &matrix, SPARSITY_PATTERN const& sparsity_pattern)
     static_assert(EigenMatrix::RawMatrixType::IsRowMajor,
                   "Set matrix sparsity relies on the EigenMatrix to be in "
                   "row-major storage order.");
-    auto const n_rows = matrix.getNRows();
-    Eigen::VectorXi row_sizes(n_rows);
 
-    for (auto i = decltype(n_rows){0}; i < n_rows; i++)
-        row_sizes[i] = sparsity_pattern.getNodeDegree(i);
+    assert(matrix.getNRows() == sparsity_pattern.size());
 
-    matrix.getRawMatrix().reserve(row_sizes);
+    matrix.getRawMatrix().reserve(sparsity_pattern);
 }
 };
 

MathLib/LinAlg/Lis/LisMatrix.h

@@ -185,8 +185,7 @@ void operator()(LisMatrix &matrix, SPARSITY_PATTERN const& sparsity_pattern)
     row_sizes.reserve(n_rows);
 
     // LIS needs 1 more entry, otherewise it starts reallocating arrays.
-    for (auto i = decltype(n_rows){0}; i < n_rows; i++)
-        row_sizes.push_back(sparsity_pattern.getNodeDegree(i) + 1);
+    for (auto i : sparsity_pattern) row_sizes.push_back(i+1);
 
     int ierr = lis_matrix_malloc(matrix._AA, 0, row_sizes.data());
     checkLisError(ierr);

MeshLib/NodeAdjacencyTable.h

@@ -54,6 +54,11 @@ public:
         return _data[node_id].size();
     }
 
+    std::vector<std::size_t> const& getAdjacentNodes(std::size_t const node_id) const
+    {
+        return _data[node_id];
+    }
+
     void createTable(std::vector<Node*> const& nodes)
     {
         if (_data.size() != nodes.size())

ProcessLib/GroundwaterFlowProcess.h

@@ -28,6 +28,7 @@
 #include "AssemblerLib/VectorMatrixAssembler.h"
 #include "AssemblerLib/LocalDataInitializer.h"
 #include "AssemblerLib/LocalToGlobalIndexMap.h"
+#include "AssemblerLib/ComputeSparsityPattern.h"
 
 #include "FileIO/VtkIO/VtuInterface.h"
 
@@ -36,7 +37,6 @@
 
 #include "MeshLib/MeshSubset.h"
 #include "MeshLib/MeshSubsets.h"
-#include "MeshLib/NodeAdjacencyTable.h"
 #include "MeshGeoToolsLib/MeshNodeSearcher.h"
 
 #include "UniformDirichletBoundaryCondition.h"
@@ -53,7 +53,7 @@ namespace MeshLib
     class Element;
     class Mesh;
     template <typename PROP_VAL_TYPE> class PropertyVector;
-};
+}
 
 namespace ProcessLib
 {
@@ -238,7 +238,9 @@ public:
         _A.reset(_global_setup.createMatrix(num_unknowns, mat_opt));
 #else
         DBUG("Compute sparsity pattern");
-        _node_adjacency_table.createTable(_mesh.getNodes());
+        _sparsity_pattern = std::move(
+            AssemblerLib::computeSparsityPattern(
+                *_local_to_global_index_map, _mesh));
 
         DBUG("Allocate global matrix, vectors, and linear solver.");
         const std::size_t num_unknowns = _local_to_global_index_map->dofSize();
@@ -282,7 +284,7 @@ public:
         DBUG("Solve GroundwaterFlowProcess.");
 
         _A->setZero();
-        MathLib::setMatrixSparsity(*_A, _node_adjacency_table);
+        MathLib::setMatrixSparsity(*_A, _sparsity_pattern);
         *_rhs = 0;   // This resets the whole vector.
 
         // Call global assembler for each local assembly item.
@@ -430,7 +432,7 @@ private:
 
     std::vector<NeumannBc<GlobalSetup>*> _neumann_bcs;
 
-    MeshLib::NodeAdjacencyTable _node_adjacency_table;
+    AssemblerLib::SparsityPattern _sparsity_pattern;
 };
 
 }   // namespace ProcessLib