diff --git a/AssemblerLib/ComputeSparsityPattern.cpp b/AssemblerLib/ComputeSparsityPattern.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c45c3b3231666c54b80f642e8f81e5d3b26240f4
--- /dev/null
+++ b/AssemblerLib/ComputeSparsityPattern.cpp
@@ -0,0 +1,57 @@
+/**
+ * \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;
+}
+
+}
diff --git a/AssemblerLib/ComputeSparsityPattern.h b/AssemblerLib/ComputeSparsityPattern.h
new file mode 100644
index 0000000000000000000000000000000000000000..766e2c05568624f20bf1545c3dee8d0104e0e24d
--- /dev/null
+++ b/AssemblerLib/ComputeSparsityPattern.h
@@ -0,0 +1,41 @@
+/**
+ * \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
+
diff --git a/AssemblerLib/LocalToGlobalIndexMap.h b/AssemblerLib/LocalToGlobalIndexMap.h
index 4b21e4612038b650ec48a465b7c8ab0018a05e52..ccac55bf7242cc4df6b1e40a05090fef1f7a9985 100644
--- a/AssemblerLib/LocalToGlobalIndexMap.h
+++ b/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
diff --git a/MathLib/LinAlg/Eigen/EigenMatrix.h b/MathLib/LinAlg/Eigen/EigenMatrix.h
index 44e15aff2d1604572bff397e5ee5845fbd5656f9..24396601ebcab395705d619484b8c6a9d2119de6 100644
--- a/MathLib/LinAlg/Eigen/EigenMatrix.h
+++ b/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);
}
};
diff --git a/MathLib/LinAlg/Lis/LisMatrix.h b/MathLib/LinAlg/Lis/LisMatrix.h
index 1756accbde367ad0ce5d07c490cf7a5b934f5c2e..63257364956eb4cff86c250b5b440f82f6148e76 100644
--- a/MathLib/LinAlg/Lis/LisMatrix.h
+++ b/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);
diff --git a/MeshLib/NodeAdjacencyTable.h b/MeshLib/NodeAdjacencyTable.h
index f10da9af1a2ade97e756f8509c29d4ce355359b2..e54d4dddaacb394b98ebf0a9c74193530178f368 100644
--- a/MeshLib/NodeAdjacencyTable.h
+++ b/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())
diff --git a/ProcessLib/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlowProcess.h
index 39623c1b88c78a2a994c53e9e861c12345a1f09b..99a5026f47e1fab044f7a90c81a1065405ba64db 100644
--- a/ProcessLib/GroundwaterFlowProcess.h
+++ b/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