Merge pull request #1240 from chleh/move-traits-back-to-mathlib

Move MatrixVectorTraits back to MathLib

MathLib/LinAlg/GlobalMatrixVectorTypes.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 MATHLIB_LINALG_GLOBALMATRIXVECTORTYPES_H
+#define MATHLIB_LINALG_GLOBALMATRIXVECTORTYPES_H
+
+#include "SparsityPattern.h"
+
+//
+// Global vector/matrix types and linear solver.
+//
+#if defined(OGS_USE_EIGENLIS)
+
+    #include "MathLib/LinAlg/Eigen/EigenMatrix.h"
+    #include "MathLib/LinAlg/Eigen/EigenVector.h"
+    #include "MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h"
+
+    namespace detail
+    {
+    using GlobalVectorType = MathLib::EigenVector;
+    using GlobalMatrixType = MathLib::EigenMatrix;
+
+    using LinearSolverType = MathLib::EigenLisLinearSolver;
+    }
+
+#elif defined(USE_PETSC)
+    #include "MathLib/LinAlg/PETSc/PETScVector.h"
+    #include "MathLib/LinAlg/PETSc/PETScMatrix.h"
+    #include "MathLib/LinAlg/PETSc/PETScLinearSolver.h"
+namespace detail
+{
+    using GlobalVectorType = MathLib::PETScVector;
+    using GlobalMatrixType = MathLib::PETScMatrix;
+
+    using LinearSolverType = MathLib::PETScLinearSolver;
+}
+
+#else
+#ifdef OGS_USE_EIGEN
+    #include "MathLib/LinAlg/Eigen/EigenVector.h"
+    #include "MathLib/LinAlg/Eigen/EigenMatrix.h"
+    #include "MathLib/LinAlg/Eigen/EigenLinearSolver.h"
+namespace detail
+{
+    using GlobalVectorType = MathLib::EigenVector;
+    using GlobalMatrixType = MathLib::EigenMatrix;
+
+    using LinearSolverType = MathLib::EigenLinearSolver;
+}
+#else   // OGS_USE_EIGEN
+    #include "MathLib/LinAlg/Dense/DenseVector.h"
+    #include "MathLib/LinAlg/Dense/GlobalDenseMatrix.h"
+    #include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
+namespace detail
+{
+    using GlobalVectorType = MathLib::DenseVector<double>;
+    using GlobalMatrixType = MathLib::GlobalDenseMatrix<double>;
+
+    using LinearSolverType =
+        MathLib::GaussAlgorithm<GlobalMatrixType, GlobalVectorType>;
+}
+
+#endif    // USE_LIS
+#endif    // OGS_USE_EIGEN
+
+
+/// A type used for indexing of global vectors and matrices. It is equal to the
+/// GlobalMatrixType::IndexType and the GlobalVectorType::IndexType.
+using GlobalIndexType = detail::GlobalMatrixType::IndexType;
+using GlobalSparsityPattern = MathLib::SparsityPattern<GlobalIndexType>;
+
+#endif // MATHLIB_LINALG_GLOBALMATRIXVECTORTYPES_H

MathLib/LinAlg/MatrixSpecifications.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 MATHLIB_LINALG_MATRIXSPECIFICATIONS_H
+#define MATHLIB_LINALG_MATRIXSPECIFICATIONS_H
+
+#include "GlobalMatrixVectorTypes.h"
+
+namespace MathLib
+{
+
+struct MatrixSpecifications
+{
+    MatrixSpecifications(std::size_t const nrows_, std::size_t const ncols_,
+                         std::vector<GlobalIndexType> const*const ghost_indices_,
+                         GlobalSparsityPattern const*const sparsity_pattern_)
+        : nrows(nrows_), ncols(ncols_), ghost_indices(ghost_indices_)
+        , sparsity_pattern(sparsity_pattern_)
+    {
+    }
+
+    std::size_t const nrows;
+    std::size_t const ncols;
+    std::vector<GlobalIndexType> const*const ghost_indices;
+    GlobalSparsityPattern const*const sparsity_pattern;
+};
+
+} // namespace MathLib
+
+
+#endif // MATHLIB_LINALG_MATRIXSPECIFICATIONS_H

NumLib/DOF/MatrixVectorTraits.cpp → MathLib/LinAlg/MatrixVectorTraits.cpp

@@ -7,9 +7,8 @@
  *
  */
 
-#include "LocalToGlobalIndexMap.h"
-#include "MatrixProviderUser.h"
 #include "MatrixVectorTraits.h"
+#include "MatrixSpecifications.h"
 
 #ifdef OGS_USE_EIGEN
 
@@ -34,10 +33,8 @@ std::unique_ptr<Eigen::MatrixXd>
 MatrixVectorTraits<Eigen::MatrixXd>::
 newInstance(MatrixSpecifications const& spec)
 {
-    auto const nrows = spec.dof_table ? spec.dof_table->dofSizeWithGhosts() : spec.nrows;
-    auto const ncols = spec.dof_table ? nrows : spec.ncols;
-
-    return std::unique_ptr<Eigen::MatrixXd>(new Eigen::MatrixXd(nrows, ncols));
+    return std::unique_ptr<Eigen::MatrixXd>(
+        new Eigen::MatrixXd(spec.nrows, spec.ncols));
 }
 
 std::unique_ptr<Eigen::VectorXd>
@@ -58,9 +55,7 @@ std::unique_ptr<Eigen::VectorXd>
 MatrixVectorTraits<Eigen::VectorXd>::
 newInstance(MatrixSpecifications const& spec)
 {
-    auto const nrows = spec.dof_table ? spec.dof_table->dofSizeWithGhosts() : spec.nrows;
-
-    return std::unique_ptr<Eigen::VectorXd>(new Eigen::VectorXd(nrows));
+    return std::unique_ptr<Eigen::VectorXd>(new Eigen::VectorXd(spec.nrows));
 }
 
 } // namespace MathLib
@@ -93,22 +88,15 @@ std::unique_ptr<PETScMatrix>
 MatrixVectorTraits<PETScMatrix>::
 newInstance(MatrixSpecifications const& spec)
 {
-    auto const nrows =
-        spec.dof_table ? spec.dof_table->dofSizeWithoutGhosts() : spec.nrows;
-    auto const ncols = spec.dof_table ? nrows : spec.ncols;
+    auto const nrows = spec.nrows;
+    auto const ncols = spec.ncols;
 
-    // TODO I guess it is not hard to make NumLib::computeSparsityPattern()
-    //      work also for NodePartitionedMesh'es.
-    assert(((!spec.sparsity_pattern) || spec.sparsity_pattern->size() == 0u) &&
-           "The sparsity pattern is not used with PETSc so I rather crash than"
-           " silently accept a precomputed sparsity pattern.");
-
-    if (spec.mesh)
+    if (spec.sparsity_pattern)
     {
-        assert(dynamic_cast<MeshLib::NodePartitionedMesh const*>(spec.mesh));
-        auto const& mesh = *static_cast<MeshLib::NodePartitionedMesh const*>(spec.mesh);
+        // Assert that the misuse of the sparsity pattern is consistent.
+        assert(spec.sparsity_pattern->size() == 1);
 
-        auto const max_nonzeroes = mesh.getMaximumNConnectedNodesToNode();
+        auto const max_nonzeroes = spec.sparsity_pattern->front();
 
         PETScMatrixOption mat_opt;
         mat_opt.d_nz = max_nonzeroes;
@@ -142,10 +130,9 @@ newInstance(MatrixSpecifications const& spec)
 {
     auto const is_global_size = false;
 
-    if (spec.dof_table) {
-        auto const& dt = *spec.dof_table;
-        return std::unique_ptr<PETScVector>(new PETScVector(
-            dt.dofSizeWithoutGhosts(), dt.getGhostIndices(), is_global_size));
+    if (spec.ghost_indices != nullptr) {
+        return std::unique_ptr<PETScVector>(
+            new PETScVector(spec.nrows, *spec.ghost_indices, is_global_size));
     } else {
         return std::unique_ptr<PETScVector>(
             new PETScVector(spec.nrows, is_global_size));
@@ -178,9 +165,7 @@ std::unique_ptr<EigenMatrix>
 MatrixVectorTraits<EigenMatrix>::
 newInstance(MatrixSpecifications const& spec)
 {
-    auto const nrows = spec.dof_table ? spec.dof_table->dofSizeWithGhosts() : spec.nrows;
-
-    auto A = std::unique_ptr<EigenMatrix>(new EigenMatrix(nrows));
+    auto A = std::unique_ptr<EigenMatrix>(new EigenMatrix(spec.nrows));
 
     if (spec.sparsity_pattern)
         setMatrixSparsity(*A, *spec.sparsity_pattern);
@@ -206,9 +191,7 @@ std::unique_ptr<EigenVector>
 MatrixVectorTraits<EigenVector>::
 newInstance(MatrixSpecifications const& spec)
 {
-    auto const nrows = spec.dof_table ? spec.dof_table->dofSizeWithGhosts() : spec.nrows;
-
-    return std::unique_ptr<EigenVector>(new EigenVector(nrows));
+    return std::unique_ptr<EigenVector>(new EigenVector(spec.nrows));
 }
 
 } // namespace MathLib

NumLib/DOF/UnifiedMatrixSetters.cpp → MathLib/LinAlg/UnifiedMatrixSetters.cpp

@@ -82,6 +82,7 @@ void setVector(Eigen::VectorXd& v, MatrixVectorTraits<Eigen::VectorXd>::Index co
 
 // Global PETScMatrix/PETScVector //////////////////////////////////////////
 
+#include <numeric>
 #include "MathLib/LinAlg/PETSc/PETScVector.h"
 #include "MathLib/LinAlg/PETSc/PETScMatrix.h"
 

NumLib/DOF/ComputeSparsityPattern.cpp

@@ -12,20 +12,36 @@
 #include "LocalToGlobalIndexMap.h"
 #include "MeshLib/NodeAdjacencyTable.h"
 
-namespace NumLib
+#ifdef USE_PETSC
+#include "MeshLib/NodePartitionedMesh.h"
+
+GlobalSparsityPattern computeSparsityPatternPETSc(
+    NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
+    MeshLib::Mesh const& mesh)
 {
-GlobalSparsityPattern computeSparsityPattern(
-    LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
+    assert(dynamic_cast<MeshLib::NodePartitionedMesh const*>(&mesh));
+    auto const& npmesh =
+        *static_cast<MeshLib::NodePartitionedMesh const*>(&mesh);
+
+    auto const max_nonzeroes = npmesh.getMaximumNConnectedNodesToNode();
+
+    // The sparsity pattern is misused here in the sense that it will only
+    // contain a single value.
+    return GlobalSparsityPattern(1, max_nonzeroes);
+}
+#else
+GlobalSparsityPattern computeSparsityPatternNonPETSc(
+    NumLib::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;
+    std::vector<std::vector<GlobalIndexType>> global_idcs;
 
     global_idcs.reserve(mesh.getNNodes());
-    for (std::size_t n=0; n<mesh.getNNodes(); ++n)
+    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));
@@ -34,15 +50,18 @@ GlobalSparsityPattern computeSparsityPattern(
     GlobalSparsityPattern sparsity_pattern(dof_table.dofSizeWithGhosts());
 
     // Map adjacent mesh nodes to "adjacent global indices".
-    for (std::size_t n=0; n<mesh.getNNodes(); ++n)
+    for (std::size_t n = 0; n < mesh.getNNodes(); ++n)
     {
         auto const& node_ids = node_adjacency_table.getAdjacentNodes(n);
-        for (auto an : node_ids) {
+        for (auto an : node_ids)
+        {
             auto const& row_ids = global_idcs[an];
             auto const num_components = row_ids.size();
-            for (auto r : row_ids) {
+            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.
+                // For the sparsity pattern only the number of entries are
+                // needed.
                 sparsity_pattern[r] += num_components;
             }
         }
@@ -50,5 +69,18 @@ GlobalSparsityPattern computeSparsityPattern(
 
     return sparsity_pattern;
 }
+#endif
+
+namespace NumLib
+{
+GlobalSparsityPattern computeSparsityPattern(
+    LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
+{
+#ifdef USE_PETSC
+    return computeSparsityPatternPETSc(dof_table, mesh);
+#else
+    return computeSparsityPatternNonPETSc(dof_table, mesh);
+#endif
+}
 
 }

NumLib/DOF/MatrixProviderUser.h

@@ -12,31 +12,11 @@
 
 #include <cstddef>
 
-#include "NumLib/NumericsConfig.h"
-
-namespace NumLib { class LocalToGlobalIndexMap; }
-namespace MeshLib { class Mesh; }
+#include "MathLib/LinAlg/MatrixSpecifications.h"
 
 namespace MathLib
 {
 
-struct MatrixSpecifications
-{
-    MatrixSpecifications(std::size_t const nrows_, std::size_t const ncols_,
-                         GlobalSparsityPattern const*const sparsity_pattern_,
-                         NumLib::LocalToGlobalIndexMap const*const dof_table_,
-                         MeshLib::Mesh const*const mesh_)
-        : nrows(nrows_), ncols(ncols_), sparsity_pattern(sparsity_pattern_)
-        , dof_table(dof_table_), mesh(mesh_)
-    {}
-
-    std::size_t const nrows;
-    std::size_t const ncols;
-    GlobalSparsityPattern const*const sparsity_pattern;
-    NumLib::LocalToGlobalIndexMap const*const dof_table;
-    MeshLib::Mesh const*const mesh;
-};
-
 class MatrixSpecificationsProvider
 {
 public:

NumLib/DOF/MeshComponentMap.cpp

@@ -141,6 +141,7 @@ MeshComponentMap::MeshComponentMap(
         }
         comp_id++;
     }
+    _num_local_dof = _dict.size();
 
     if (order == ComponentOrder::BY_LOCATION)
         renumberByLocation();

NumLib/DOF/SimpleMatrixVectorProvider-impl.h

@@ -11,7 +11,7 @@
 #include <logog/include/logog.hpp>
 
 #include "MathLib/LinAlg/BLAS.h"
-#include "MatrixVectorTraits.h"
+#include "MathLib/LinAlg/MatrixVectorTraits.h"
 #include "SimpleMatrixVectorProvider.h"
 
 namespace detail

NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h

@@ -14,7 +14,7 @@
 
 #include "MathLib/LinAlg/BLAS.h"
 #include "NumLib/Assembler/SerialExecutor.h"
-#include "NumLib/DOF/MatrixVectorTraits.h"
+#include "MathLib/LinAlg/MatrixVectorTraits.h"
 #include "NumLib/Function/Interpolation.h"
 #include "LocalLinearLeastSquaresExtrapolator.h"
 

NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h

@@ -49,15 +49,20 @@ public:
      * dof_table for one single component variable.
      */
     explicit LocalLinearLeastSquaresExtrapolator(
-            MathLib::MatrixSpecifications const& matrix_specs)
-        : _nodal_values(MathLib::GlobalVectorProvider<GlobalVector>::provider.
-                        getVector(matrix_specs))
+        MathLib::MatrixSpecifications const& matrix_specs,
+        NumLib::LocalToGlobalIndexMap const& dof_table)
+        : _nodal_values(
+              MathLib::GlobalVectorProvider<GlobalVector>::provider.getVector(
+                  matrix_specs))
 #ifndef USE_PETSC
-        , _residuals(matrix_specs.dof_table->size())
+          ,
+          _residuals(dof_table.size())
 #else
-        , _residuals(matrix_specs.dof_table->size(), false)
+          ,
+          _residuals(dof_table.size(), false)
 #endif
-        , _local_to_global(*matrix_specs.dof_table)
+          ,
+          _local_to_global(dof_table)
     {}
 
 

NumLib/NumericsConfig.h

@@ -11,8 +11,7 @@
 #define APPLICATIONS_NUMERICSCONFIG_H_
 
 #include <type_traits>
-
-#include "MathLib/LinAlg/SparsityPattern.h"
+#include "MathLib/LinAlg/GlobalMatrixVectorTypes.h"
 
 /**
  * This file provides a configuration of the global matrix/vector and
@@ -22,64 +21,6 @@
  * The existence of the GlobalSetupType is checked at the end of the file.
  */
 
-//
-// Global vector/matrix types and linear solver.
-//
-#if defined(OGS_USE_EIGENLIS)
-
-    #include "MathLib/LinAlg/Eigen/EigenMatrix.h"
-    #include "MathLib/LinAlg/Eigen/EigenVector.h"
-    #include "MathLib/LinAlg/EigenLis/EigenLisLinearSolver.h"
-
-    namespace detail
-    {
-    using GlobalVectorType = MathLib::EigenVector;
-    using GlobalMatrixType = MathLib::EigenMatrix;
-
-    using LinearSolverType = MathLib::EigenLisLinearSolver;
-    }
-
-#elif defined(USE_PETSC)
-    #include "MathLib/LinAlg/PETSc/PETScVector.h"
-    #include "MathLib/LinAlg/PETSc/PETScMatrix.h"
-    #include "MathLib/LinAlg/PETSc/PETScLinearSolver.h"
-namespace detail
-{
-    using GlobalVectorType = MathLib::PETScVector;
-    using GlobalMatrixType = MathLib::PETScMatrix;
-
-    using LinearSolverType = MathLib::PETScLinearSolver;
-}
-
-#else
-#ifdef OGS_USE_EIGEN
-    #include "MathLib/LinAlg/Eigen/EigenVector.h"
-    #include "MathLib/LinAlg/Eigen/EigenMatrix.h"
-    #include "MathLib/LinAlg/Eigen/EigenLinearSolver.h"
-namespace detail
-{
-    using GlobalVectorType = MathLib::EigenVector;
-    using GlobalMatrixType = MathLib::EigenMatrix;
-
-    using LinearSolverType = MathLib::EigenLinearSolver;
-}
-#else   // OGS_USE_EIGEN
-    #include "MathLib/LinAlg/Dense/DenseVector.h"
-    #include "MathLib/LinAlg/Dense/GlobalDenseMatrix.h"
-    #include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
-namespace detail
-{
-    using GlobalVectorType = MathLib::DenseVector<double>;
-    using GlobalMatrixType = MathLib::GlobalDenseMatrix<double>;
-
-    using LinearSolverType =
-        MathLib::GaussAlgorithm<GlobalMatrixType, GlobalVectorType>;
-}
-
-#endif    // USE_LIS
-#endif    // OGS_USE_EIGEN
-
-
 //
 // Global vector/matrix builder.
 //
@@ -127,9 +68,4 @@ static_assert(std::is_same<detail::GlobalMatrixType::IndexType,
               "The global matrix and vector index types do not match.");
 // Both types are integral types and equal, define a single GlobalIndexType.
 
-/// A type used for indexing of global vectors and matrices. It is equal to the
-/// GlobalMatrixType::IndexType and the GlobalVectorType::IndexType.
-using GlobalIndexType = detail::GlobalMatrixType::IndexType;
-using GlobalSparsityPattern = MathLib::SparsityPattern<GlobalIndexType>;
-
 #endif  // APPLICATIONS_NUMERICSCONFIG_H_

NumLib/ODESolver/ODESystem.h

@@ -11,7 +11,7 @@
 #define NUMLIB_ODESYSTEM_H
 
 #include "NumLib/IndexValueVector.h"
-#include "NumLib/DOF/MatrixVectorTraits.h"
+#include "MathLib/LinAlg/MatrixVectorTraits.h"
 
 #include "Types.h"
 #include "EquationSystem.h"

NumLib/ODESolver/TimeDiscretizedODESystem.h

@@ -13,8 +13,8 @@
 #include <memory>
 
 #include "MathLib/LinAlg/ApplyKnownSolution.h"
+#include "MathLib/LinAlg/UnifiedMatrixSetters.h"
 #include "NumLib/IndexValueVector.h"
-#include "NumLib/DOF/UnifiedMatrixSetters.h"
 
 #include "ODESystem.h"
 #include "NonlinearSystem.h"

ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h

@@ -102,7 +102,8 @@ private:
 
         // TOOD Later on the DOF table can change during the simulation!
         _extrapolator.reset(
-            new ExtrapolatorImplementation(Base::getMatrixSpecifications()));
+            new ExtrapolatorImplementation(Base::getMatrixSpecifications(),
+                                           *Base::_local_to_global_index_map));
 
         Base::_secondary_variables.addSecondaryVariable(
             "darcy_velocity_x", 1,

ProcessLib/Process.h

@@ -87,10 +87,8 @@ public:
         DBUG("Construct dof mappings.");
         constructDofTable();
 
-#ifndef USE_PETSC
         DBUG("Compute sparsity pattern");
         computeSparsityPattern();
-#endif
 
         initializeConcreteProcess(*_local_to_global_index_map, _mesh,
                                   _integration_order);
@@ -127,10 +125,9 @@ public:
 
     MathLib::MatrixSpecifications getMatrixSpecifications() const override final
     {
-        return { 0u, 0u,
-                 &_sparsity_pattern,
-                 _local_to_global_index_map.get(),
-                 &_mesh };
+        auto const& l = *_local_to_global_index_map;
+        return {l.dofSizeWithoutGhosts(), l.dofSizeWithoutGhosts(),
+                &l.getGhostIndices(), &_sparsity_pattern};
     }
 
     void assemble(const double t, GlobalVector const& x,

ProcessLib/TES/TESProcess.cpp

@@ -192,10 +192,15 @@ void TESProcess<GlobalSetup>::initializeConcreteProcess(
             // by location order is needed for output
             NumLib::ComponentOrder::BY_LOCATION));
 
-    _extrapolator.reset(
-        new ExtrapolatorImplementation(MathLib::MatrixSpecifications(
-            0u, 0u, nullptr, _local_to_global_index_map_single_component.get(),
-            &mesh)));
+    {
+        auto const& l = *_local_to_global_index_map_single_component;
+        _extrapolator.reset(new ExtrapolatorImplementation(
+            MathLib::MatrixSpecifications(l.dofSizeWithoutGhosts(),
+                                          l.dofSizeWithoutGhosts(),
+                                          &l.getGhostIndices(),
+                                          nullptr),
+            l));
+    }
 
     // secondary variables
     auto add2nd = [&](std::string const& var_name, unsigned const n_components,
@@ -343,9 +348,11 @@ TESProcess<GlobalSetup>::computeVapourPartialPressure(
 {
     assert(&dof_table == this->_local_to_global_index_map.get());
 
+    auto const& dof_table_single = *_local_to_global_index_map_single_component;
     result_cache = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
-        {0, 0, nullptr, _local_to_global_index_map_single_component.get(),
-         nullptr});
+        {dof_table_single.dofSizeWithoutGhosts(),
+         dof_table_single.dofSizeWithoutGhosts(),
+         &dof_table_single.getGhostIndices(), nullptr});
 
     GlobalIndexType nnodes = this->_mesh.getNNodes();
 
@@ -375,9 +382,11 @@ TESProcess<GlobalSetup>::computeRelativeHumidity(
 {
     assert(&dof_table == this->_local_to_global_index_map.get());
 
+    auto const& dof_table_single = *_local_to_global_index_map_single_component;
     result_cache = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
-        {0, 0, nullptr, _local_to_global_index_map_single_component.get(),
-         nullptr});
+        {dof_table_single.dofSizeWithoutGhosts(),
+         dof_table_single.dofSizeWithoutGhosts(),
+         &dof_table_single.getGhostIndices(), nullptr});
 
     GlobalIndexType nnodes = this->_mesh.getNNodes();
 
@@ -412,9 +421,11 @@ TESProcess<GlobalSetup>::computeEquilibriumLoading(
 {
     assert(&dof_table == this->_local_to_global_index_map.get());
 
+    auto const& dof_table_single = *_local_to_global_index_map_single_component;
     result_cache = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(
-        {0, 0, nullptr, _local_to_global_index_map_single_component.get(),
-         nullptr});
+        {dof_table_single.dofSizeWithoutGhosts(),
+         dof_table_single.dofSizeWithoutGhosts(),
+         &dof_table_single.getGhostIndices(), nullptr});
 
     GlobalIndexType nnodes = this->_mesh.getNNodes();
 

Tests/NumLib/ODEs.h

@@ -10,8 +10,8 @@
 #ifndef TESTS_NUMLIB_ODES_H
 #define TESTS_NUMLIB_ODES_H
 
-#include "NumLib/DOF/UnifiedMatrixSetters.h"
 #include "MathLib/LinAlg/BLAS.h"
+#include "MathLib/LinAlg/UnifiedMatrixSetters.h"
 #include "NumLib/ODESolver/ODESystem.h"
 
 // debug
@@ -55,7 +55,7 @@ public:
 
     MathLib::MatrixSpecifications getMatrixSpecifications() const override
     {
-        return { N, N, nullptr, nullptr, nullptr };
+        return { N, N, nullptr, nullptr };
     }
 
     bool isLinear() const override
@@ -133,7 +133,7 @@ public:
 
     MathLib::MatrixSpecifications getMatrixSpecifications() const override
     {
-        return { N, N, nullptr, nullptr, nullptr };
+        return { N, N, nullptr, nullptr };
     }
 
     bool isLinear() const override
@@ -271,7 +271,7 @@ public:
 
     MathLib::MatrixSpecifications getMatrixSpecifications() const override
     {
-        return { N, N, nullptr, nullptr, nullptr };
+        return { N, N, nullptr, nullptr };
     }
 
     bool isLinear() const override

Tests/NumLib/TestExtrapolation.cpp

@@ -11,8 +11,8 @@
 #include <gtest/gtest.h>
 
 #include "NumLib/DOF/MatrixProviderUser.h"
-#include "NumLib/DOF/MatrixVectorTraits.h"
-#include "NumLib/DOF/UnifiedMatrixSetters.h"
+#include "MathLib/LinAlg/MatrixVectorTraits.h"
+#include "MathLib/LinAlg/UnifiedMatrixSetters.h"
 #include "NumLib/Assembler/VectorMatrixAssembler.h"
 
 #include "MathLib/LinAlg/BLAS.h"
@@ -182,7 +182,11 @@ public:
         // Passing _dof_table works, because this process has only one variable
         // and the variable has exactly one component.
         _extrapolator.reset(new ExtrapolatorImplementation(
-            MathLib::MatrixSpecifications{0u, 0u, nullptr, _dof_table.get(), nullptr}));
+            MathLib::MatrixSpecifications{_dof_table->dofSizeWithoutGhosts(),
+                                          _dof_table->dofSizeWithoutGhosts(),
+                                          &_dof_table->getGhostIndices(),
+                                          nullptr},
+            *_dof_table));
 
         createAssemblers(mesh);
     }
@@ -333,7 +337,7 @@ TEST(NumLib, DISABLED_Extrapolation)
         TestProcess<GlobalSetup> pcs(*mesh, integration_order);
 
         // generate random nodal values
-        MathLib::MatrixSpecifications spec{nnodes, nnodes, nullptr, nullptr, nullptr};
+        MathLib::MatrixSpecifications spec{nnodes, nnodes, nullptr, nullptr};
         auto x = MathLib::MatrixVectorTraits<GlobalVector>::newInstance(spec);
 
         fillVectorRandomly(*x);