diff --git a/Applications/ApplicationsLib/LinearSolverLibrarySetup.h b/Applications/ApplicationsLib/LinearSolverLibrarySetup.h
index 65e83055063d25e9691c47ecb08f1ffbf08bf38c..555d1a2a17547ada4e0bea6b811b124342aa1b2f 100644
--- a/Applications/ApplicationsLib/LinearSolverLibrarySetup.h
+++ b/Applications/ApplicationsLib/LinearSolverLibrarySetup.h
@@ -19,7 +19,7 @@
 /// The default implementation is empty providing polymorphic behaviour when
 /// using this class.
 
-#include "MathLib/LinAlg/GlobalMatrixProviders.h"
+#include "NumLib/DOF/GlobalMatrixProviders.h"
 
 #if defined(USE_PETSC)
 #include <petsc.h>
diff --git a/Applications/CLI/CMakeLists.txt b/Applications/CLI/CMakeLists.txt
index 6329b9e1a711d7eb8fe741dbf5280ff6890ba468..085155e363ad21aab906d779aad0feefc0f4ba10 100644
--- a/Applications/CLI/CMakeLists.txt
+++ b/Applications/CLI/CMakeLists.txt
@@ -5,7 +5,6 @@ add_executable(ogs
 
 target_link_libraries(ogs
     ApplicationsLib
-    AssemblerLib
     MeshGeoToolsLib
     ProcessLib
     NumLib
diff --git a/Applications/CLI/ogs.cpp b/Applications/CLI/ogs.cpp
index 16d95f119db3aaad82e8136402d1354178edb9a0..6fa28db83adf2fc61f4d83ba0323818396500dac 100644
--- a/Applications/CLI/ogs.cpp
+++ b/Applications/CLI/ogs.cpp
@@ -24,7 +24,7 @@
 #include "Applications/ApplicationsLib/ProjectData.h"
 #include "Applications/ApplicationsLib/UncoupledProcessesTimeLoop.h"
 
-#include "ProcessLib/NumericsConfig.h"
+#include "NumLib/NumericsConfig.h"
 
 
 int main(int argc, char *argv[])
diff --git a/Applications/DataExplorer/DataExplorer.cmake b/Applications/DataExplorer/DataExplorer.cmake
index 074582ca08d07acc2bd2172257dffa446ed14860..cbb440bf5ae83ac19a587912fe348e7b233d1405 100644
--- a/Applications/DataExplorer/DataExplorer.cmake
+++ b/Applications/DataExplorer/DataExplorer.cmake
@@ -57,7 +57,6 @@ target_link_libraries(DataExplorer
     ApplicationsFileIO
     DataHolderLib
     FileIO
-    InSituLib
     QtDataView
     QtDiagramView
     QtStratView
diff --git a/Applications/DataExplorer/DataView/ElementTreeModel.cpp b/Applications/DataExplorer/DataView/ElementTreeModel.cpp
index 840c2be8ddf34799e59a0ecdea0279498df358c6..4ea042eabb9c5151bf65bf771732648df3134d64 100644
--- a/Applications/DataExplorer/DataView/ElementTreeModel.cpp
+++ b/Applications/DataExplorer/DataView/ElementTreeModel.cpp
@@ -21,7 +21,7 @@
 
 #include "GeoLib/AABB.h"
 
-#include "InSituLib/VtkMappedMeshSource.h"
+#include "MeshLib/Vtk/VtkMappedMeshSource.h"
 
 #include "TreeItem.h"
 
@@ -46,8 +46,8 @@ void ElementTreeModel::setElement(vtkUnstructuredGridAlgorithm const*const grid,
     _mesh_source = grid;
     this->clearView();
 
-    InSituLib::VtkMappedMeshSource const*const source =
-        dynamic_cast<InSituLib::VtkMappedMeshSource const*const>(grid);
+    MeshLib::VtkMappedMeshSource const*const source =
+        dynamic_cast<MeshLib::VtkMappedMeshSource const*const>(grid);
 
     if (!source)
         return;
diff --git a/Applications/DataExplorer/DataView/MshItem.cpp b/Applications/DataExplorer/DataView/MshItem.cpp
index 4bb692bbbaca120b01b96f38f73e7fc833a12367..663662fee77d85d3786f5c16a3053bbe9502ed05 100644
--- a/Applications/DataExplorer/DataView/MshItem.cpp
+++ b/Applications/DataExplorer/DataView/MshItem.cpp
@@ -13,7 +13,7 @@
  */
 
 #include "MshItem.h"
-#include "InSituLib/VtkMappedMeshSource.h"
+#include "MeshLib/Vtk/VtkMappedMeshSource.h"
 
 /**
  * Constructor.
@@ -24,8 +24,8 @@
 MshItem::MshItem(const QList<QVariant> &data, TreeItem* parent, const MeshLib::Mesh* mesh)
     : TreeItem(data, parent)
 {
-    _mesh_source = InSituLib::VtkMappedMeshSource::New();
-    static_cast<InSituLib::VtkMappedMeshSource*>(_mesh_source)->SetMesh(mesh);
+    _mesh_source = MeshLib::VtkMappedMeshSource::New();
+    static_cast<MeshLib::VtkMappedMeshSource*>(_mesh_source)->SetMesh(mesh);
 }
 
 MshItem::~MshItem()
diff --git a/Applications/DataExplorer/DataView/MshItem.h b/Applications/DataExplorer/DataView/MshItem.h
index 2c9debd854a8150a9b892251f141c67cbd6fee06..5ad4f366d3e48409b6d6ac200f7747a05114dac3 100644
--- a/Applications/DataExplorer/DataView/MshItem.h
+++ b/Applications/DataExplorer/DataView/MshItem.h
@@ -17,7 +17,7 @@
 
 #include "TreeItem.h"
 
-#include "InSituLib/VtkMappedMeshSource.h"
+#include "MeshLib/Vtk/VtkMappedMeshSource.h"
 
 namespace MeshLib {
     class Mesh;
@@ -37,10 +37,10 @@ public:
     /// Returns the mesh.
     MeshLib::Mesh const* getMesh() const { return _mesh_source->GetMesh(); }
     /// Returns the VTK object.
-    InSituLib::VtkMappedMeshSource*  vtkSource() const { return _mesh_source; }
+    MeshLib::VtkMappedMeshSource*  vtkSource() const { return _mesh_source; }
 
 private:
-    InSituLib::VtkMappedMeshSource * _mesh_source;
+    MeshLib::VtkMappedMeshSource * _mesh_source;
 };
 
 #endif //MSHITEM_H
diff --git a/Applications/DataExplorer/DataView/MshView.h b/Applications/DataExplorer/DataView/MshView.h
index 4ca6e4079520c64fab91a81c54c8f4123798346f..c1117ffb8ec2309e05dafa792481e19133a2499e 100644
--- a/Applications/DataExplorer/DataView/MshView.h
+++ b/Applications/DataExplorer/DataView/MshView.h
@@ -25,7 +25,7 @@ namespace MeshLib {
     class Mesh;
 }
 
-namespace InSituLib {
+namespace MeshLib {
     class VtkMappedMeshSource;
 }
 
@@ -101,7 +101,7 @@ signals:
     void enableRemoveButton(bool);
     void meshSelected(MeshLib::Mesh const&);
     void openMeshFile(int);
-    void qualityCheckRequested(InSituLib::VtkMappedMeshSource*);
+    void qualityCheckRequested(MeshLib::VtkMappedMeshSource*);
     void removeSelectedMeshComponent();
     void requestCondSetupDialog(const std::string&, const GeoLib::GEOTYPE, const std::size_t, bool on_points);
     void requestMeshRemoval(const QModelIndex&);
diff --git a/Applications/DataExplorer/VtkVis/VtkVisPipeline.cpp b/Applications/DataExplorer/VtkVis/VtkVisPipeline.cpp
index f8a115b53e0f9ba64e6df036fa081a3f0ed7a2f6..799f55cfd9756c8af67709f3a7feeb664b987504 100644
--- a/Applications/DataExplorer/VtkVis/VtkVisPipeline.cpp
+++ b/Applications/DataExplorer/VtkVis/VtkVisPipeline.cpp
@@ -36,7 +36,7 @@
 #include "VtkVisPipelineItem.h"
 #include "VtkVisPointSetItem.h"
 
-#include "InSituLib/VtkMappedMeshSource.h"
+#include "MeshLib/Vtk/VtkMappedMeshSource.h"
 
 #include <vtkAlgorithm.h>
 #include <vtkCamera.h>
@@ -317,8 +317,8 @@ QModelIndex VtkVisPipeline::addPipelineItem( vtkAlgorithm* source, QModelIndex p
         vtkXMLReader* new_reader = dynamic_cast<vtkXMLReader*>(source);
         vtkImageReader2* image_reader = dynamic_cast<vtkImageReader2*>(source);
         VtkAlgorithmProperties* props = dynamic_cast<VtkAlgorithmProperties*>(source);
-        InSituLib::VtkMappedMeshSource* meshSource =
-            dynamic_cast<InSituLib::VtkMappedMeshSource*>(source);
+        MeshLib::VtkMappedMeshSource* meshSource =
+            dynamic_cast<MeshLib::VtkMappedMeshSource*>(source);
         if (old_reader)
             itemName = old_reader->GetFileName();
         else if (new_reader)
@@ -432,7 +432,7 @@ void VtkVisPipeline::listArrays(vtkDataSet* dataSet)
 }
 
 void VtkVisPipeline::showMeshElementQuality(
-    InSituLib::VtkMappedMeshSource* source,
+    MeshLib::VtkMappedMeshSource* source,
     MeshLib::MeshQualityType t, std::vector<double> const& quality)
 {
     if (!source || quality.empty())
diff --git a/Applications/DataExplorer/VtkVis/VtkVisPipeline.h b/Applications/DataExplorer/VtkVis/VtkVisPipeline.h
index cdf6b0f1f51c496cabcdf6970e075dadc3d11a2c..7623720650e87c3b3b3fedc4aef363941670bad7 100644
--- a/Applications/DataExplorer/VtkVis/VtkVisPipeline.h
+++ b/Applications/DataExplorer/VtkVis/VtkVisPipeline.h
@@ -46,7 +46,7 @@ class StationTreeModel;
 class TreeModel;
 class VtkVisPipelineItem;
 
-namespace InSituLib
+namespace MeshLib
 {
     class VtkMappedMeshSource;
 }
@@ -101,7 +101,7 @@ public:
     void setGlobalBackfaceCulling(bool enable) const;
 
     /// Checks the quality of mesh elements and adds a filter to highlight deformed elements.
-    void showMeshElementQuality(InSituLib::VtkMappedMeshSource* mesh, MeshLib::MeshQualityType t, std::vector<double> const& quality);
+    void showMeshElementQuality(MeshLib::VtkMappedMeshSource* mesh, MeshLib::MeshQualityType t, std::vector<double> const& quality);
 
 public slots:
     /// \brief Adds the given Model to the pipeline.
diff --git a/Applications/DataExplorer/VtkVis/VtkVisPointSetItem.cpp b/Applications/DataExplorer/VtkVis/VtkVisPointSetItem.cpp
index 98fa56fc5bbb96d1e2dcc02cc026ff1193ac5ef1..927bc53d48375cfb975a2fd397a317f394d85bfa 100644
--- a/Applications/DataExplorer/VtkVis/VtkVisPointSetItem.cpp
+++ b/Applications/DataExplorer/VtkVis/VtkVisPointSetItem.cpp
@@ -31,7 +31,7 @@
 #include "VtkCompositeFilter.h"
 #include "VtkCompositeContourFilter.h"
 #include "VtkCompositeThresholdFilter.h"
-#include "InSituLib/VtkMappedMeshSource.h"
+#include "MeshLib/Vtk/VtkMappedMeshSource.h"
 
 #include "QVtkDataSetMapper.h"
 #include <vtkActor.h>
@@ -185,7 +185,7 @@ void VtkVisPointSetItem::Initialize(vtkRenderer* renderer)
     }
 
     // Show edges on meshes
-    if (dynamic_cast<InSituLib::VtkMappedMeshSource*>(this->_algorithm))
+    if (dynamic_cast<MeshLib::VtkMappedMeshSource*>(this->_algorithm))
         _vtkProps->GetProperties()->SetEdgeVisibility(1);
 }
 
diff --git a/Applications/DataExplorer/mainwindow.cpp b/Applications/DataExplorer/mainwindow.cpp
index 7413056141793b57165fc7556a123914ac50fd5e..ffea4aff4af38b478104b1ac93ae4f6a1fb32639 100644
--- a/Applications/DataExplorer/mainwindow.cpp
+++ b/Applications/DataExplorer/mainwindow.cpp
@@ -57,7 +57,7 @@
 #include "VtkRaster.h"
 #include "VtkVisPipelineItem.h"
 
-#include "InSituLib/VtkMappedMeshSource.h"
+#include "MeshLib/Vtk/VtkMappedMeshSource.h"
 
 // FileIO includes
 #include "GMSInterface.h"
@@ -170,9 +170,9 @@ MainWindow::MainWindow(QWidget* parent /* = 0*/)
     connect(mshTabWidget->treeView, SIGNAL(requestMeshRemoval(const QModelIndex &)),
             _elementModel.get(), SLOT(clearView()));
     connect(mshTabWidget->treeView,
-        SIGNAL(qualityCheckRequested(InSituLib::VtkMappedMeshSource*)),
+        SIGNAL(qualityCheckRequested(MeshLib::VtkMappedMeshSource*)),
             this,
-            SLOT(showMeshQualitySelectionDialog(InSituLib::VtkMappedMeshSource*)));
+            SLOT(showMeshQualitySelectionDialog(MeshLib::VtkMappedMeshSource*)));
     connect(mshTabWidget->treeView, SIGNAL(requestMeshToGeometryConversion(const MeshLib::Mesh*)),
             this, SLOT(convertMeshToGeometry(const MeshLib::Mesh*)));
     connect(mshTabWidget->treeView, SIGNAL(elementSelected(vtkUnstructuredGridAlgorithm const*const, unsigned, bool)),
@@ -1055,7 +1055,7 @@ void MainWindow::showMergeGeometriesDialog()
         OGSError::box("Points are missing for\n at least one geometry.");
 }
 
-void MainWindow::showMeshQualitySelectionDialog(InSituLib::VtkMappedMeshSource* mshSource)
+void MainWindow::showMeshQualitySelectionDialog(MeshLib::VtkMappedMeshSource* mshSource)
 {
     if (mshSource == nullptr)
         return;
diff --git a/Applications/DataExplorer/mainwindow.h b/Applications/DataExplorer/mainwindow.h
index 24b1c4da557ec8da13f6ffc85a7d08c259467898..81436b36df9376798b09cb92da553bd56d473c82 100644
--- a/Applications/DataExplorer/mainwindow.h
+++ b/Applications/DataExplorer/mainwindow.h
@@ -33,7 +33,7 @@
 class TreeModel;
 class ProcessModel;
 
-namespace InSituLib
+namespace MeshLib
 {
     class VtkMappedMeshSource;
 }
@@ -104,7 +104,7 @@ protected slots:
     void showGMSHPrefsDialog();
     void showMergeGeometriesDialog();
     void showMeshAnalysisDialog();
-    void showMeshQualitySelectionDialog(InSituLib::VtkMappedMeshSource* mshSource);
+    void showMeshQualitySelectionDialog(MeshLib::VtkMappedMeshSource* mshSource);
     void showVisalizationPrefsDialog();
     void updateDataViews();
     void writeGeometryToFile(QString listName, QString fileName);
diff --git a/AssemblerLib/CMakeLists.txt b/AssemblerLib/CMakeLists.txt
deleted file mode 100644
index 8bb4b5be726702cce9fb160e006e9971538cdc75..0000000000000000000000000000000000000000
--- a/AssemblerLib/CMakeLists.txt
+++ /dev/null
@@ -1,19 +0,0 @@
-include(${PROJECT_SOURCE_DIR}/scripts/cmake/OGSEnabledElements.cmake)
-
-#Source files grouped by a directory
-GET_SOURCE_FILES(SOURCES_ASSEMBLERLIB)
-set(SOURCES ${SOURCES_ASSEMBLERLIB})
-
-# Create the library
-add_library(AssemblerLib STATIC ${SOURCES})
-
-target_link_libraries(AssemblerLib
-    MeshLib
-)
-
-if(TARGET Eigen)
-    add_dependencies(AssemblerLib Eigen)
-endif()
-if(TARGET Boost)
-    add_dependencies(AssemblerLib Boost)
-endif()
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1e7efe9521c49e7337e713fe809691426b254bea..f74323313f8ac79c16998e383b0149aae2ab9529 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -188,7 +188,6 @@ include_directories( SYSTEM ${CMAKE_CURRENT_SOURCE_DIR}/ThirdParty/tclap/include
 include_directories(SYSTEM ${CMAKE_CURRENT_SOURCE_DIR}/ThirdParty/vtkGUISupportQt)
 
 add_subdirectory( Applications )
-add_subdirectory( AssemblerLib )
 add_subdirectory( BaseLib )
 # TODO This is a hack but we have to make sure that Boost is built first
 if(TARGET Boost)
@@ -197,7 +196,6 @@ endif()
 add_subdirectory( DataHolderLib )
 add_subdirectory( FileIO )
 add_subdirectory( GeoLib )
-add_subdirectory( InSituLib )
 add_subdirectory( MaterialsLib )
 add_subdirectory( MathLib )
 add_subdirectory( MeshLib )
diff --git a/GeoLib/CMakeLists.txt b/GeoLib/CMakeLists.txt
index 64bb24741980bc1d0059f19c8f4fdf68661f364e..f0bb369a1e716ccaab2bf64c473e4bad885525a7 100644
--- a/GeoLib/CMakeLists.txt
+++ b/GeoLib/CMakeLists.txt
@@ -22,7 +22,7 @@ GET_SOURCE_FILES(SOURCES_IO_GMSHIO IO/Gmsh)
 set(SOURCES ${SOURCES} ${SOURCES_IO_GMSHIO})
 
 # Create the library
-add_library(GeoLib STATIC ${SOURCES}
+add_library(GeoLib ${SOURCES}
     ${CMAKE_CURRENT_SOURCE_DIR}/../ThirdParty/tetgen/predicates.cxx
 )
 
diff --git a/MathLib/CMakeLists.txt b/MathLib/CMakeLists.txt
index dd6028af94d6df408ae4e2d1dd041ca4c1252a3d..af30101b392911215d91a5894171b3f7625f21e2 100644
--- a/MathLib/CMakeLists.txt
+++ b/MathLib/CMakeLists.txt
@@ -45,12 +45,12 @@ set(SOURCES ${SOURCES} ${SOURCES_NONLINEAR})
 
 
 # Create the library
-add_library(MathLib STATIC ${SOURCES})
+add_library(MathLib ${SOURCES})
 
 set_target_properties(MathLib PROPERTIES LINKER_LANGUAGE CXX)
 
 target_link_libraries(MathLib
-    AssemblerLib
+    BaseLib
 )
 
 if (CVODE_FOUND)
diff --git a/MathLib/LinAlg/LinearSolverOptions.h b/MathLib/LinAlg/LinearSolverOptions.h
index 9a3156b05c8bd890150267cb2b0ce5d96a23e4fd..af99422535936e1786f87afc8390b21f1fed9701 100644
--- a/MathLib/LinAlg/LinearSolverOptions.h
+++ b/MathLib/LinAlg/LinearSolverOptions.h
@@ -1,4 +1,14 @@
-#pragma once
+/**
+ * \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_LINEAR_SOLVER_OPTIONS_H
+#define MATHLIB_LINEAR_SOLVER_OPTIONS_H
 
 #include "BaseLib/ConfigTree.h"
 
@@ -27,3 +37,5 @@ void ignoreOtherLinearSolvers(BaseLib::ConfigTree const& config,
                               std::string const& solver_name);
 
 }
+
+#endif
diff --git a/MathLib/LinAlg/SparsityPattern.h b/MathLib/LinAlg/SparsityPattern.h
index 5fc669c324018f8a3ef659d1a8b7ed560661ba4e..678a2339c493b37b3011e6d5aec13dab4b3c970b 100644
--- a/MathLib/LinAlg/SparsityPattern.h
+++ b/MathLib/LinAlg/SparsityPattern.h
@@ -12,12 +12,11 @@
 
 #include <vector>
 
-#include "ProcessLib/NumericsConfig.h"
-
 namespace MathLib
 {
 /// A vector telling how many nonzeros there are in each global matrix row.
-using SparsityPattern = std::vector<GlobalIndexType>;
+template <typename IndexType>
+using SparsityPattern = std::vector<IndexType>;
 }
 
 #endif // MATHLIB_LINALG_SPARSITYPATTERN_H
diff --git a/MeshLib/CMakeLists.txt b/MeshLib/CMakeLists.txt
index 1b23db93238608f23a042117ff9019c0061cca93..7ea38ca1e3f1019bf20acfa58c811192eca534b8 100644
--- a/MeshLib/CMakeLists.txt
+++ b/MeshLib/CMakeLists.txt
@@ -13,12 +13,12 @@ GET_SOURCE_FILES(SOURCES_SEARCH MeshSearch)
 GET_SOURCE_FILES(SOURCES_IO IO)
 GET_SOURCE_FILES(SOURCES_IO_LEGACY IO/Legacy)
 GET_SOURCE_FILES(SOURCES_IO_VTKIO IO/VtkIO)
-
 GET_SOURCE_FILES(SOURCES_QUALITY MeshQuality)
+GET_SOURCE_FILES(SOURCES_VTK Vtk)
 
 set(SOURCES ${SOURCES_MESHLIB} ${SOURCES_ELEMENTS} ${SOURCES_EDITING}
     ${SOURCES_GENERATORS} ${SOURCES_QUALITY} ${SOURCES_SEARCH}
-    ${SOURCES_IO} ${SOURCES_IO_LEGACY}  ${SOURCES_IO_VTKIO})
+    ${SOURCES_IO} ${SOURCES_IO_LEGACY}  ${SOURCES_IO_VTKIO} ${SOURCES_VTK})
 
 if(QT4_FOUND)
     GET_SOURCE_FILES(SOURCES_IO_FEFLOW IO/FEFLOW)
@@ -32,13 +32,12 @@ if(OGS_USE_PETSC)
 endif()
 
 # Create the library
-add_library(MeshLib STATIC ${SOURCES})
+add_library(MeshLib ${SOURCES})
 
 target_link_libraries(MeshLib
     BaseLib
     GeoLib
     MathLib
-    InSituLib # VtkMappedMeshSource
     ${VTK_LIBRARIES}
 )
 
diff --git a/MeshLib/IO/VtkIO/VtuInterface-impl.h b/MeshLib/IO/VtkIO/VtuInterface-impl.h
index 72a76eb8a0d717909ae2d302de579a89713480c1..87fa7faa950e9ae2c4c5b57f0dceef225c4147c0 100644
--- a/MeshLib/IO/VtkIO/VtuInterface-impl.h
+++ b/MeshLib/IO/VtkIO/VtuInterface-impl.h
@@ -24,7 +24,7 @@
 
 #include <logog/include/logog.hpp>
 
-#include "InSituLib/VtkMappedMeshSource.h"
+#include "MeshLib/Vtk/VtkMappedMeshSource.h"
 
 namespace MeshLib
 {
@@ -44,7 +44,7 @@ bool VtuInterface::writeVTU(std::string const &file_name, const int num_partitio
     if(_data_mode == vtkXMLWriter::Appended)
         WARN("Appended data mode is currently not supported, written file is not valid!");
 
-    vtkNew<InSituLib::VtkMappedMeshSource> vtkSource;
+    vtkNew<MeshLib::VtkMappedMeshSource> vtkSource;
     vtkSource->SetMesh(_mesh);
 
     vtkSmartPointer<UnstructuredGridWriter> vtuWriter =
diff --git a/MeshLib/IO/VtkIO/VtuInterface.cpp b/MeshLib/IO/VtkIO/VtuInterface.cpp
index eeb68c27bc960ea857eeab1c96f3a5b6a8599370..abb2d49d3c12a9ae147265e1c5f94aa1ea0ce809 100644
--- a/MeshLib/IO/VtkIO/VtuInterface.cpp
+++ b/MeshLib/IO/VtkIO/VtuInterface.cpp
@@ -31,9 +31,9 @@
 #endif
 
 #include "BaseLib/FileTools.h"
-#include "InSituLib/VtkMappedMeshSource.h"
 #include "MeshLib/Mesh.h"
 #include "MeshLib/MeshGenerators/VtkMeshConverter.h"
+#include "MeshLib/Vtk/VtkMappedMeshSource.h"
 
 namespace MeshLib
 {
diff --git a/InSituLib/VtkMappedMesh.cpp b/MeshLib/Vtk/VtkMappedMesh.cpp
similarity index 99%
rename from InSituLib/VtkMappedMesh.cpp
rename to MeshLib/Vtk/VtkMappedMesh.cpp
index c73f3dcd0df2f7129c8e60d14f3c3f81ffd9db8f..a5998250d1953f7d369c057411b79a9d2f584d95 100644
--- a/InSituLib/VtkMappedMesh.cpp
+++ b/MeshLib/Vtk/VtkMappedMesh.cpp
@@ -30,7 +30,7 @@
 #include "MeshLib/Node.h"
 #include "MeshLib/VtkOGSEnum.h"
 
-namespace InSituLib {
+namespace MeshLib {
 
 vtkStandardNewMacro(VtkMappedMesh)
 vtkStandardNewMacro(VtkMappedMeshImpl)
diff --git a/InSituLib/VtkMappedMesh.h b/MeshLib/Vtk/VtkMappedMesh.h
similarity index 99%
rename from InSituLib/VtkMappedMesh.h
rename to MeshLib/Vtk/VtkMappedMesh.h
index ee5fda3d02ca05971f277973aaf11e7e51df23e4..a80b8d809634eb9ee82425e43828f53d195ba1e7 100644
--- a/InSituLib/VtkMappedMesh.h
+++ b/MeshLib/Vtk/VtkMappedMesh.h
@@ -26,7 +26,7 @@ namespace MeshLib {
     class Node;
 }
 
-namespace InSituLib
+namespace MeshLib
 {
 
 class VtkMappedMeshImpl : public vtkObject
diff --git a/InSituLib/VtkMappedMeshSource.cpp b/MeshLib/Vtk/VtkMappedMeshSource.cpp
similarity index 98%
rename from InSituLib/VtkMappedMeshSource.cpp
rename to MeshLib/Vtk/VtkMappedMeshSource.cpp
index 47ed8c52b8e4a69e1c9f16e0654340a0e560a9ce..ee1fab6643aa532c68ac44b1dfc22c5c3309a81b 100644
--- a/InSituLib/VtkMappedMeshSource.cpp
+++ b/MeshLib/Vtk/VtkMappedMeshSource.cpp
@@ -22,7 +22,7 @@
 #include "VtkMappedMesh.h"
 #include "VtkMeshNodalCoordinatesTemplate.h"
 
-namespace InSituLib {
+namespace MeshLib {
 
 vtkStandardNewMacro(VtkMappedMeshSource)
 
@@ -116,4 +116,4 @@ int VtkMappedMeshSource::RequestInformation(
     return 1;
 }
 
-} // Namespace InSituLib
+} // Namespace MeshLib
diff --git a/InSituLib/VtkMappedMeshSource.h b/MeshLib/Vtk/VtkMappedMeshSource.h
similarity index 93%
rename from InSituLib/VtkMappedMeshSource.h
rename to MeshLib/Vtk/VtkMappedMeshSource.h
index 98116bdcd8b649473077d06b7103550884fc0f7d..4460d1bc58cb84719a804bf45fa0ad2b8893a435 100644
--- a/InSituLib/VtkMappedMeshSource.h
+++ b/MeshLib/Vtk/VtkMappedMeshSource.h
@@ -6,7 +6,7 @@
  *         vtkUnstructuredGrids.
  * Usage:
  * \code
- * vtkNew<InSituLib::VtkMappedMeshSource> vtkSource;
+ * vtkNew<MeshLib::VtkMappedMeshSource> vtkSource;
  * vtkSource->SetMesh(mesh);
  * vtkSource->Update();
  * vtkUnstructuredGrid* output = vtkSource->GetOutput();
@@ -34,10 +34,11 @@
 #include <vtkUnstructuredGrid.h>
 #include <vtkUnstructuredGridAlgorithm.h>
 
-#include "InSituLib/VtkMappedPropertyVectorTemplate.h"
 #include "MeshLib/Properties.h"
 #include "MeshLib/PropertyVector.h"
 
+#include "VtkMappedPropertyVectorTemplate.h"
+
 class vtkCellData;
 class vtkDataArrayCollection;
 class vtkPointData;
@@ -47,7 +48,7 @@ namespace MeshLib {
     class Properties;
 }
 
-namespace InSituLib {
+namespace MeshLib {
 
 /// Adapter which maps a MeshLib::Mesh to a vtkUnstructuredGridAlgorithm.
 /// Allows for zero-copy access of the mesh from the visualization side.
@@ -79,7 +80,7 @@ private:
     VtkMappedMeshSource(const VtkMappedMeshSource &); // Not implemented.
     void operator=(const VtkMappedMeshSource &);      // Not implemented.
 
-    /// Adds a zero-copy array (InSituLib::VtkMappedPropertyVectorTemplate) as
+    /// Adds a zero-copy array (MeshLib::VtkMappedPropertyVectorTemplate) as
     /// either point or cell data to the mesh.
     /// \param properties MeshLib::Properties object
     /// \param prop_name The name of the property vector to be mapped from
@@ -115,6 +116,6 @@ private:
     vtkNew<vtkCellData> CellData;
 };
 
-} // Namespace InSituLib
+} // Namespace MeshLib
 
 #endif //_VTKMAPPEDMESHSOURCE
diff --git a/InSituLib/VtkMappedPropertyVectorTemplate-impl.h b/MeshLib/Vtk/VtkMappedPropertyVectorTemplate-impl.h
similarity index 99%
rename from InSituLib/VtkMappedPropertyVectorTemplate-impl.h
rename to MeshLib/Vtk/VtkMappedPropertyVectorTemplate-impl.h
index 54fcf31ff7906d5c227bb702a5ec9ffc97e462da..b3592b68c1dd182798d2da5e2b64c744525ec19c 100644
--- a/InSituLib/VtkMappedPropertyVectorTemplate-impl.h
+++ b/MeshLib/Vtk/VtkMappedPropertyVectorTemplate-impl.h
@@ -19,7 +19,7 @@
 #include <vtkVariant.h>
 #include <vtkVariantCast.h>
 
-namespace InSituLib {
+namespace MeshLib {
 
 // Can't use vtkStandardNewMacro on a templated class.
 template <class Scalar> VtkMappedPropertyVectorTemplate<Scalar> *
@@ -471,4 +471,4 @@ template <class Scalar> vtkIdType VtkMappedPropertyVectorTemplate<Scalar>
     return -1;
 }
 
-} // end namespace InSituLib
+} // end namespace MeshLib
diff --git a/InSituLib/VtkMappedPropertyVectorTemplate.h b/MeshLib/Vtk/VtkMappedPropertyVectorTemplate.h
similarity index 98%
rename from InSituLib/VtkMappedPropertyVectorTemplate.h
rename to MeshLib/Vtk/VtkMappedPropertyVectorTemplate.h
index 77fbc650a81df7252262ad020d5447a3bb0384e2..2cef7f9c4a9f0be85823995d153b7cea662c340b 100644
--- a/InSituLib/VtkMappedPropertyVectorTemplate.h
+++ b/MeshLib/Vtk/VtkMappedPropertyVectorTemplate.h
@@ -22,7 +22,7 @@
 
 #include "MeshLib/Elements/Element.h"
 
-namespace InSituLib {
+namespace MeshLib {
 
 template <class Scalar>
 class VtkMappedPropertyVectorTemplate :
@@ -112,7 +112,7 @@ private:
     bool Save;
 };
 
-} // end namespace InSituLib
+} // end namespace MeshLib
 
 #include "VtkMappedPropertyVectorTemplate-impl.h"
 
diff --git a/InSituLib/VtkMeshNodalCoordinatesTemplate-impl.h b/MeshLib/Vtk/VtkMeshNodalCoordinatesTemplate-impl.h
similarity index 99%
rename from InSituLib/VtkMeshNodalCoordinatesTemplate-impl.h
rename to MeshLib/Vtk/VtkMeshNodalCoordinatesTemplate-impl.h
index f308af44c82d99667eff608743e69e9882ccf073..f21973a241ea104a6ca1ea68c193bbcfad66f433 100644
--- a/InSituLib/VtkMeshNodalCoordinatesTemplate-impl.h
+++ b/MeshLib/Vtk/VtkMeshNodalCoordinatesTemplate-impl.h
@@ -21,7 +21,7 @@
 
 #include "MeshLib/Node.h"
 
-namespace InSituLib {
+namespace MeshLib {
 
 // Can't use vtkStandardNewMacro with a template.
 template <class Scalar> VtkMeshNodalCoordinatesTemplate<Scalar> *
diff --git a/InSituLib/VtkMeshNodalCoordinatesTemplate.h b/MeshLib/Vtk/VtkMeshNodalCoordinatesTemplate.h
similarity index 99%
rename from InSituLib/VtkMeshNodalCoordinatesTemplate.h
rename to MeshLib/Vtk/VtkMeshNodalCoordinatesTemplate.h
index eff4b21ef90330e990b604389f2d6a11cbe2af0b..c9a2b58149a98ba618afd344e6a97c2acddca81a 100644
--- a/InSituLib/VtkMeshNodalCoordinatesTemplate.h
+++ b/MeshLib/Vtk/VtkMeshNodalCoordinatesTemplate.h
@@ -25,7 +25,7 @@ namespace MeshLib
     class Node;
 }
 
-namespace InSituLib
+namespace MeshLib
 {
 
 template <class Scalar>
diff --git a/AssemblerLib/SerialExecutor.h b/NumLib/Assembler/SerialExecutor.h
similarity index 95%
rename from AssemblerLib/SerialExecutor.h
rename to NumLib/Assembler/SerialExecutor.h
index 37dd1bd1d7917a1effdc968d8c5bfd3ad9f822a8..205ec3513db735f731d229fc301c5f593df1d6f8 100644
--- a/AssemblerLib/SerialExecutor.h
+++ b/NumLib/Assembler/SerialExecutor.h
@@ -10,10 +10,10 @@
  *
  */
 
-#ifndef ASSEMBLERLIB_SERIALEXECUTOR_H_H
-#define ASSEMBLERLIB_SERIALEXECUTOR_H_H
+#ifndef NUMLIB_SERIALEXECUTOR_H_H
+#define NUMLIB_SERIALEXECUTOR_H_H
 
-namespace AssemblerLib
+namespace NumLib
 {
 
 struct SerialExecutor
@@ -107,6 +107,6 @@ struct SerialExecutor
     }
 };
 
-}   // namespace AssemblerLib
+}   // namespace NumLib
 
-#endif  // ASSEMBLERLIB_SERIALEXECUTOR_H_H
+#endif  // NUMLIB_SERIALEXECUTOR_H_H
diff --git a/AssemblerLib/VectorMatrixAssembler.h b/NumLib/Assembler/VectorMatrixAssembler.h
similarity index 93%
rename from AssemblerLib/VectorMatrixAssembler.h
rename to NumLib/Assembler/VectorMatrixAssembler.h
index bc920630a96f5daf36a04a370486445abbc575d1..1a3a4e063a9e54e64e37de32d39be2f79f776a19 100644
--- a/AssemblerLib/VectorMatrixAssembler.h
+++ b/NumLib/Assembler/VectorMatrixAssembler.h
@@ -7,18 +7,17 @@
  *
  */
 
-#ifndef ASSEMBLERLIB_VECTORMATRIXASSEMBLER_H_
-#define ASSEMBLERLIB_VECTORMATRIXASSEMBLER_H_
-
-#include "LocalToGlobalIndexMap.h"
+#ifndef NUMLIB_VECTORMATRIXASSEMBLER_H_
+#define NUMLIB_VECTORMATRIXASSEMBLER_H_
 
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
 #include "NumLib/ODESolver/Types.h"
 
 namespace
 {
-inline AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices
+inline NumLib::LocalToGlobalIndexMap::RowColumnIndices
 getRowColumnIndices(std::size_t const id,
-                    AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+                    NumLib::LocalToGlobalIndexMap const& dof_table,
                     std::vector<GlobalIndexType>& indices)
 {
     assert(dof_table.size() > id);
@@ -33,13 +32,13 @@ getRowColumnIndices(std::size_t const id,
         indices.insert(indices.end(), idcs.begin(), idcs.end());
     }
 
-    return AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices(indices,
+    return NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices,
                                                                  indices);
 }
 
 template <typename Callback, typename GlobalVector, typename... Args>
 void passLocalVector_(Callback& cb, std::size_t const id,
-                      AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+                      NumLib::LocalToGlobalIndexMap const& dof_table,
                       GlobalVector const& x, Args&&... args)
 {
     std::vector<GlobalIndexType> indices;
@@ -57,7 +56,7 @@ void passLocalVector_(Callback& cb, std::size_t const id,
 }
 }
 
-namespace AssemblerLib
+namespace NumLib
 {
 
 /*! Calls the local assemblers of FEM processes and assembles
@@ -215,6 +214,6 @@ private:
     LocalToGlobalIndexMap const& _data_pos;
 };
 
-}   // namespace AssemblerLib
+}   // namespace NumLib
 
-#endif  // ASSEMBLERLIB_VECTORMATRIXASSEMBLER_H_
+#endif  // NUMLIB_VECTORMATRIXASSEMBLER_H_
diff --git a/AssemblerLib/VectorMatrixBuilder.h b/NumLib/Assembler/VectorMatrixBuilder.h
similarity index 85%
rename from AssemblerLib/VectorMatrixBuilder.h
rename to NumLib/Assembler/VectorMatrixBuilder.h
index a2886f0be900e26c26a8bea43d57c1e4025cbd75..fbba5d69459e353b7db87e12d93097d92a1c7f55 100644
--- a/AssemblerLib/VectorMatrixBuilder.h
+++ b/NumLib/Assembler/VectorMatrixBuilder.h
@@ -10,10 +10,10 @@
  *
  */
 
-#ifndef ASSEMBLERLIB_VECTORMATRIXBUILDER_H_
-#define ASSEMBLERLIB_VECTORMATRIXBUILDER_H_
+#ifndef NUMLIB_VECTORMATRIXBUILDER_H_
+#define NUMLIB_VECTORMATRIXBUILDER_H_
 
-namespace AssemblerLib
+namespace NumLib
 {
 
 template <typename MatrixType_, typename VectorType_>
@@ -44,6 +44,6 @@ public:
 
 };
 
-}   // namespace AssemblerLib
+}   // namespace NumLib
 
-#endif  // ASSEMBLERLIB_VECTORMATRIXBUILDER_H_
+#endif  // NUMLIB_VECTORMATRIXBUILDER_H_
diff --git a/NumLib/CMakeLists.txt b/NumLib/CMakeLists.txt
index 7d5a00bb7cf2e6d9777d5443730ab8ab6e9f808a..b23206f9fd9053f15e8aea5b150bd607265fee34 100644
--- a/NumLib/CMakeLists.txt
+++ b/NumLib/CMakeLists.txt
@@ -2,6 +2,11 @@
 GET_SOURCE_FILES(SOURCES_NUMLIB)
 set(SOURCES ${SOURCES_NUMLIB})
 
+GET_SOURCE_FILES(SOURCES_ASSEMBLER Assembler)
+set(SOURCES ${SOURCES} ${SOURCES_ASSEMBLER})
+GET_SOURCE_FILES(SOURCES_DOF DOF)
+set(SOURCES ${SOURCES} ${SOURCES_DOF})
+
 GET_SOURCE_FILES(SOURCES_FEM Fem)
 set(SOURCES ${SOURCES} ${SOURCES_FEM})
 GET_SOURCE_FILES(SOURCES_FEM_COORDINATESMAPPING Fem/CoordinatesMapping)
diff --git a/AssemblerLib/ComponentGlobalIndexDict.h b/NumLib/DOF/ComponentGlobalIndexDict.h
similarity index 92%
rename from AssemblerLib/ComponentGlobalIndexDict.h
rename to NumLib/DOF/ComponentGlobalIndexDict.h
index 21ff6dd89aaf1f62ec0e465d2184edffd3bc27bd..6a7e0eba8c30546384acec308e89f3362bf0bc86 100644
--- a/AssemblerLib/ComponentGlobalIndexDict.h
+++ b/NumLib/DOF/ComponentGlobalIndexDict.h
@@ -10,8 +10,8 @@
  *
  */
 
-#ifndef ASSEMBLERLIB_COMPONENTGLOBALINDEXDICT_H_
-#define ASSEMBLERLIB_COMPONENTGLOBALINDEXDICT_H_
+#ifndef NUMLIB_COMPONENTGLOBALINDEXDICT_H_
+#define NUMLIB_COMPONENTGLOBALINDEXDICT_H_
 
 #include <limits>
 
@@ -21,9 +21,9 @@
 #include <boost/multi_index_container.hpp>
 
 #include "MeshLib/Location.h"
-#include "ProcessLib/NumericsConfig.h"
+#include "NumLib/NumericsConfig.h"
 
-namespace AssemblerLib
+namespace NumLib
 {
 
 /// \internal
@@ -118,6 +118,6 @@ typedef boost::multi_index::multi_index_container<
     > ComponentGlobalIndexDict;
 
 }    // namespace detail
-}    // namespace AssemblerLib
+}    // namespace NumLib
 
-#endif  // ASSEMBLERLIB_COMPONENTGLOBALINDEXDICT_H_
+#endif  // NUMLIB_COMPONENTGLOBALINDEXDICT_H_
diff --git a/AssemblerLib/ComputeSparsityPattern.cpp b/NumLib/DOF/ComputeSparsityPattern.cpp
similarity index 91%
rename from AssemblerLib/ComputeSparsityPattern.cpp
rename to NumLib/DOF/ComputeSparsityPattern.cpp
index b357ac594e7ba786f8633effd052df496719b0e7..c6085cb116cb5ebc4caaf6b417fef2848eba4deb 100644
--- a/AssemblerLib/ComputeSparsityPattern.cpp
+++ b/NumLib/DOF/ComputeSparsityPattern.cpp
@@ -12,9 +12,9 @@
 #include "LocalToGlobalIndexMap.h"
 #include "MeshLib/NodeAdjacencyTable.h"
 
-namespace AssemblerLib
+namespace NumLib
 {
-MathLib::SparsityPattern computeSparsityPattern(
+GlobalSparsityPattern computeSparsityPattern(
     LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh)
 {
     MeshLib::NodeAdjacencyTable node_adjacency_table;
@@ -31,7 +31,7 @@ MathLib::SparsityPattern computeSparsityPattern(
         global_idcs.push_back(dof_table.getGlobalIndices(l));
     }
 
-    MathLib::SparsityPattern sparsity_pattern(dof_table.dofSizeWithGhosts());
+    GlobalSparsityPattern sparsity_pattern(dof_table.dofSizeWithGhosts());
 
     // Map adjacent mesh nodes to "adjacent global indices".
     for (std::size_t n=0; n<mesh.getNNodes(); ++n)
diff --git a/AssemblerLib/ComputeSparsityPattern.h b/NumLib/DOF/ComputeSparsityPattern.h
similarity index 70%
rename from AssemblerLib/ComputeSparsityPattern.h
rename to NumLib/DOF/ComputeSparsityPattern.h
index f4364c526c44fc387d69a19ef07f0bee08927702..568d4b231896c2ecdef5ed20e0d56f5d43a62b0b 100644
--- a/AssemblerLib/ComputeSparsityPattern.h
+++ b/NumLib/DOF/ComputeSparsityPattern.h
@@ -7,20 +7,19 @@
  *
  */
 
-#ifndef ASSEMBLERLIB_COMPUTESPARSITYPATTERN_H
-#define ASSEMBLERLIB_COMPUTESPARSITYPATTERN_H
+#ifndef NUMLIB_COMPUTESPARSITYPATTERN_H
+#define NUMLIB_COMPUTESPARSITYPATTERN_H
 
 #include <vector>
 
-#include "MathLib/LinAlg/SparsityPattern.h"
-#include "ProcessLib/NumericsConfig.h"
+#include "NumLib/NumericsConfig.h"
 
 namespace MeshLib
 {
 class Mesh;
 }
 
-namespace AssemblerLib
+namespace NumLib
 {
 class LocalToGlobalIndexMap;
 
@@ -32,9 +31,9 @@ class LocalToGlobalIndexMap;
  *
  * @return The computed sparsity pattern.
  */
-MathLib::SparsityPattern computeSparsityPattern(
+GlobalSparsityPattern computeSparsityPattern(
     LocalToGlobalIndexMap const& dof_table, MeshLib::Mesh const& mesh);
 }
 
-#endif // ASSEMBLERLIB_COMPUTESPARSITYPATTERN_H
+#endif // NUMLIB_COMPUTESPARSITYPATTERN_H
 
diff --git a/MathLib/LinAlg/GlobalMatrixProviders.cpp b/NumLib/DOF/GlobalMatrixProviders.cpp
similarity index 96%
rename from MathLib/LinAlg/GlobalMatrixProviders.cpp
rename to NumLib/DOF/GlobalMatrixProviders.cpp
index ff847a7b000ac4c517fafb026fea833208c25f2f..0d420288509d653626e72da155e90d158a9150d1 100644
--- a/MathLib/LinAlg/GlobalMatrixProviders.cpp
+++ b/NumLib/DOF/GlobalMatrixProviders.cpp
@@ -7,9 +7,9 @@
  *
  */
 
-#include<memory>
+#include <memory>
 
-#include "ProcessLib/NumericsConfig.h"
+#include "NumLib/NumericsConfig.h"
 
 #include "GlobalMatrixProviders.h"
 #include "SimpleMatrixVectorProvider.h"
diff --git a/MathLib/LinAlg/GlobalMatrixProviders.h b/NumLib/DOF/GlobalMatrixProviders.h
similarity index 100%
rename from MathLib/LinAlg/GlobalMatrixProviders.h
rename to NumLib/DOF/GlobalMatrixProviders.h
diff --git a/AssemblerLib/LocalToGlobalIndexMap.cpp b/NumLib/DOF/LocalToGlobalIndexMap.cpp
similarity index 97%
rename from AssemblerLib/LocalToGlobalIndexMap.cpp
rename to NumLib/DOF/LocalToGlobalIndexMap.cpp
index e1c55ac7ec89686083799f397fa1bf6dab33fcf8..27432950335095b3ea46762e9debad64452f96ff 100644
--- a/AssemblerLib/LocalToGlobalIndexMap.cpp
+++ b/NumLib/DOF/LocalToGlobalIndexMap.cpp
@@ -9,7 +9,7 @@
 
 #include "LocalToGlobalIndexMap.h"
 
-namespace AssemblerLib
+namespace NumLib
 {
 
 template <typename ElementIterator>
@@ -46,7 +46,7 @@ LocalToGlobalIndexMap::findGlobalIndices(
 
 LocalToGlobalIndexMap::LocalToGlobalIndexMap(
     std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
-    AssemblerLib::ComponentOrder const order)
+    NumLib::ComponentOrder const order)
     : _mesh_subsets(std::move(mesh_subsets)),
       _mesh_component_map(_mesh_subsets, order)
 {
@@ -71,7 +71,7 @@ LocalToGlobalIndexMap::LocalToGlobalIndexMap(
     std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
     std::size_t const component_id,
     std::vector<MeshLib::Element*> const& elements,
-    AssemblerLib::MeshComponentMap&& mesh_component_map)
+    NumLib::MeshComponentMap&& mesh_component_map)
     : _mesh_subsets(std::move(mesh_subsets)),
       _mesh_component_map(std::move(mesh_component_map))
 {
@@ -176,4 +176,4 @@ std::ostream& operator<<(std::ostream& os, LocalToGlobalIndexMap const& map)
 }
 #endif  // NDEBUG
 
-}   // namespace AssemblerLib
+}   // namespace NumLib
diff --git a/AssemblerLib/LocalToGlobalIndexMap.h b/NumLib/DOF/LocalToGlobalIndexMap.h
similarity index 93%
rename from AssemblerLib/LocalToGlobalIndexMap.h
rename to NumLib/DOF/LocalToGlobalIndexMap.h
index 657bf8bcc6a3a3ab5018d9f7fdb2bbeb7d11dd7f..9932f447094452c59577fdbd6acf1df5a505d946 100644
--- a/AssemblerLib/LocalToGlobalIndexMap.h
+++ b/NumLib/DOF/LocalToGlobalIndexMap.h
@@ -7,8 +7,8 @@
  *
  */
 
-#ifndef ASSEMBLERLIB_LOCALTOGLOBALINDEXMAP_H_
-#define ASSEMBLERLIB_LOCALTOGLOBALINDEXMAP_H_
+#ifndef NUMLIB_LOCALTOGLOBALINDEXMAP_H_
+#define NUMLIB_LOCALTOGLOBALINDEXMAP_H_
 
 #ifndef NDEBUG
 #include <iosfwd>
@@ -18,11 +18,12 @@
 
 #include <Eigen/Dense>
 
-#include "AssemblerLib/MeshComponentMap.h"
 #include "MathLib/LinAlg/RowColumnIndices.h"
 #include "MeshLib/MeshSubsets.h"
 
-namespace AssemblerLib
+#include "MeshComponentMap.h"
+
+namespace NumLib
 {
 
 /// Row and column indices in global linear algebra objects for each mesh item.
@@ -45,7 +46,7 @@ public:
     /// each mesh element of the given mesh_subsets.
     explicit LocalToGlobalIndexMap(
         std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
-        AssemblerLib::ComponentOrder const order);
+        NumLib::ComponentOrder const order);
 
     /// Derive a LocalToGlobalIndexMap constrained to a set of mesh subsets and
     /// elements. A new mesh component map will be constructed using the passed
@@ -124,7 +125,7 @@ private:
         std::vector<std::unique_ptr<MeshLib::MeshSubsets>>&& mesh_subsets,
         std::size_t const component_id,
         std::vector<MeshLib::Element*> const& elements,
-        AssemblerLib::MeshComponentMap&& mesh_component_map);
+        NumLib::MeshComponentMap&& mesh_component_map);
 
     template <typename ElementIterator>
     void
@@ -135,7 +136,7 @@ private:
 private:
     /// A vector of mesh subsets for each process variables' components.
     std::vector<std::unique_ptr<MeshLib::MeshSubsets>> const _mesh_subsets;
-    AssemblerLib::MeshComponentMap _mesh_component_map;
+    NumLib::MeshComponentMap _mesh_component_map;
 
     using Table = Eigen::Matrix<LineIndex, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
 
@@ -153,6 +154,6 @@ private:
 
 };
 
-}   // namespace AssemblerLib
+}   // namespace NumLib
 
-#endif  // ASSEMBLERLIB_LOCALTOGLOBALINDEXMAP_H_
+#endif  // NUMLIB_LOCALTOGLOBALINDEXMAP_H_
diff --git a/MathLib/LinAlg/MatrixProviderUser.h b/NumLib/DOF/MatrixProviderUser.h
similarity index 92%
rename from MathLib/LinAlg/MatrixProviderUser.h
rename to NumLib/DOF/MatrixProviderUser.h
index b6d16e1c631704b267227d279d47194072081646..82cfe2d8f6a22e8fd1345f587356d905ab67b09c 100644
--- a/MathLib/LinAlg/MatrixProviderUser.h
+++ b/NumLib/DOF/MatrixProviderUser.h
@@ -12,9 +12,9 @@
 
 #include <cstddef>
 
-#include "SparsityPattern.h"
+#include "NumLib/NumericsConfig.h"
 
-namespace AssemblerLib { class LocalToGlobalIndexMap; }
+namespace NumLib { class LocalToGlobalIndexMap; }
 namespace MeshLib { class Mesh; }
 
 namespace MathLib
@@ -23,8 +23,8 @@ namespace MathLib
 struct MatrixSpecifications
 {
     MatrixSpecifications(std::size_t const nrows_, std::size_t const ncols_,
-                         SparsityPattern const*const sparsity_pattern_,
-                         AssemblerLib::LocalToGlobalIndexMap const*const dof_table_,
+                         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_)
@@ -32,8 +32,8 @@ struct MatrixSpecifications
 
     std::size_t const nrows;
     std::size_t const ncols;
-    SparsityPattern const*const sparsity_pattern;
-    AssemblerLib::LocalToGlobalIndexMap const*const dof_table;
+    GlobalSparsityPattern const*const sparsity_pattern;
+    NumLib::LocalToGlobalIndexMap const*const dof_table;
     MeshLib::Mesh const*const mesh;
 };
 
diff --git a/MathLib/LinAlg/MatrixVectorTraits.cpp b/NumLib/DOF/MatrixVectorTraits.cpp
similarity index 97%
rename from MathLib/LinAlg/MatrixVectorTraits.cpp
rename to NumLib/DOF/MatrixVectorTraits.cpp
index d112e9791a2966b2b6441c9ed80d9a18468bbf94..2cdd4d74e1c16d122ccc071f44424ba8cc1ac55f 100644
--- a/MathLib/LinAlg/MatrixVectorTraits.cpp
+++ b/NumLib/DOF/MatrixVectorTraits.cpp
@@ -7,7 +7,7 @@
  *
  */
 
-#include "AssemblerLib/LocalToGlobalIndexMap.h"
+#include "LocalToGlobalIndexMap.h"
 #include "MatrixProviderUser.h"
 #include "MatrixVectorTraits.h"
 
@@ -97,7 +97,7 @@ newInstance(MatrixSpecifications const& spec)
         spec.dof_table ? spec.dof_table->dofSizeWithoutGhosts() : spec.nrows;
     auto const ncols = spec.dof_table ? nrows : spec.ncols;
 
-    // TODO I guess it is not hard to make AssemblerLib::computeSparsityPattern()
+    // 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"
diff --git a/MathLib/LinAlg/MatrixVectorTraits.h b/NumLib/DOF/MatrixVectorTraits.h
similarity index 100%
rename from MathLib/LinAlg/MatrixVectorTraits.h
rename to NumLib/DOF/MatrixVectorTraits.h
diff --git a/AssemblerLib/MeshComponentMap.cpp b/NumLib/DOF/MeshComponentMap.cpp
similarity index 99%
rename from AssemblerLib/MeshComponentMap.cpp
rename to NumLib/DOF/MeshComponentMap.cpp
index 0c6dccd591315e137115c8f5a63c64767a29d26b..cb6ab849ca672e6183f3bc50186fef5406955785 100644
--- a/AssemblerLib/MeshComponentMap.cpp
+++ b/NumLib/DOF/MeshComponentMap.cpp
@@ -18,7 +18,7 @@
 #include "MeshLib/NodePartitionedMesh.h"
 #endif
 
-namespace AssemblerLib
+namespace NumLib
 {
 
 using namespace detail;
@@ -321,4 +321,4 @@ GlobalIndexType MeshComponentMap::getLocalIndex(
 #endif
 }
 
-}   // namespace AssemblerLib
+}   // namespace NumLib
diff --git a/AssemblerLib/MeshComponentMap.h b/NumLib/DOF/MeshComponentMap.h
similarity index 97%
rename from AssemblerLib/MeshComponentMap.h
rename to NumLib/DOF/MeshComponentMap.h
index 221337754d8d190b635e0004cac8e871c50b68cb..4e6e2e71d6d21911cfc007ba5c1e66a52e8a8fd0 100644
--- a/AssemblerLib/MeshComponentMap.h
+++ b/NumLib/DOF/MeshComponentMap.h
@@ -10,8 +10,8 @@
  *
  */
 
-#ifndef ASSEMBLERLIB_MESHCOMPONENTMAP_H_
-#define ASSEMBLERLIB_MESHCOMPONENTMAP_H_
+#ifndef NUMLIB_MESHCOMPONENTMAP_H_
+#define NUMLIB_MESHCOMPONENTMAP_H_
 
 #include "ComponentGlobalIndexDict.h"
 
@@ -20,7 +20,7 @@ namespace MeshLib
     class MeshSubsets;
 }
 
-namespace AssemblerLib
+namespace NumLib
 {
 
 /// Ordering of components in global matrix/vector.
@@ -190,6 +190,6 @@ private:
     std::vector<GlobalIndexType> _ghosts_indices;
 };
 
-}   // namespace AssemblerLib
+}   // namespace NumLib
 
-#endif  // ASSEMBLERLIB_MESHCOMPONENTMAP_H_
+#endif  // NUMLIB_MESHCOMPONENTMAP_H_
diff --git a/MathLib/LinAlg/SimpleMatrixVectorProvider-impl.h b/NumLib/DOF/SimpleMatrixVectorProvider-impl.h
similarity index 99%
rename from MathLib/LinAlg/SimpleMatrixVectorProvider-impl.h
rename to NumLib/DOF/SimpleMatrixVectorProvider-impl.h
index 8e1d83690b1b199d7094512839a21b70551a93de..aa7d09b989c8e1328419c6ef00389585411b65d8 100644
--- a/MathLib/LinAlg/SimpleMatrixVectorProvider-impl.h
+++ b/NumLib/DOF/SimpleMatrixVectorProvider-impl.h
@@ -10,7 +10,7 @@
 #include <cassert>
 #include <logog/include/logog.hpp>
 
-#include "BLAS.h"
+#include "MathLib/LinAlg/BLAS.h"
 #include "MatrixVectorTraits.h"
 #include "SimpleMatrixVectorProvider.h"
 
diff --git a/MathLib/LinAlg/SimpleMatrixVectorProvider.h b/NumLib/DOF/SimpleMatrixVectorProvider.h
similarity index 100%
rename from MathLib/LinAlg/SimpleMatrixVectorProvider.h
rename to NumLib/DOF/SimpleMatrixVectorProvider.h
diff --git a/MathLib/LinAlg/UnifiedMatrixSetters.cpp b/NumLib/DOF/UnifiedMatrixSetters.cpp
similarity index 100%
rename from MathLib/LinAlg/UnifiedMatrixSetters.cpp
rename to NumLib/DOF/UnifiedMatrixSetters.cpp
diff --git a/MathLib/LinAlg/UnifiedMatrixSetters.h b/NumLib/DOF/UnifiedMatrixSetters.h
similarity index 100%
rename from MathLib/LinAlg/UnifiedMatrixSetters.h
rename to NumLib/DOF/UnifiedMatrixSetters.h
diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h
index c9dcf84211573bbc7bc8e5f1b0d6cb95737c577f..c8c35472679456a133338d3d011ef9e534e102f8 100644
--- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h
+++ b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator-impl.h
@@ -12,9 +12,9 @@
 #include <logog/include/logog.hpp>
 #include <Eigen/Core>
 
-#include "AssemblerLib/SerialExecutor.h"
 #include "MathLib/LinAlg/BLAS.h"
-#include "MathLib/LinAlg/MatrixVectorTraits.h"
+#include "NumLib/Assembler/SerialExecutor.h"
+#include "NumLib/DOF/MatrixVectorTraits.h"
 #include "NumLib/Function/Interpolation.h"
 #include "LocalLinearLeastSquaresExtrapolator.h"
 
@@ -37,7 +37,7 @@ extrapolate(LocalAssemblers const& local_assemblers, PropertyTag const property)
     using Self = LocalLinearLeastSquaresExtrapolator<
         GlobalVector, PropertyTag, LocalAssembler>;
 
-    AssemblerLib::SerialExecutor::executeMemberDereferenced(
+    NumLib::SerialExecutor::executeMemberDereferenced(
         *this, &Self::extrapolateElement, local_assemblers, property, *counts);
 
     MathLib::BLAS::componentwiseDivide(_nodal_values, _nodal_values, *counts);
@@ -54,7 +54,7 @@ calculateResiduals(LocalAssemblers const& local_assemblers,
     using Self = LocalLinearLeastSquaresExtrapolator<
         GlobalVector, PropertyTag, LocalAssembler>;
 
-    AssemblerLib::SerialExecutor::executeMemberDereferenced(
+    NumLib::SerialExecutor::executeMemberDereferenced(
         *this, &Self::calculateResiudalElement, local_assemblers, property);
 }
 
diff --git a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
index b177164304962a752a6c66014b7570c7c60a295c..b15c46ab4294d40104c54b5acd4ba5cdf332841d 100644
--- a/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
+++ b/NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h
@@ -10,7 +10,7 @@
 #ifndef NUMLIB_LOCAL_LLSQ_EXTRAPOLATOR_H
 #define NUMLIB_LOCAL_LLSQ_EXTRAPOLATOR_H
 
-#include "MathLib/LinAlg/GlobalMatrixProviders.h"
+#include "NumLib/DOF/GlobalMatrixProviders.h"
 #include "Extrapolator.h"
 
 namespace NumLib
@@ -100,7 +100,7 @@ private:
     GlobalVector  _residuals;    //!< extrapolation residuals
 
     //! DOF table used for writing to global vectors.
-    AssemblerLib::LocalToGlobalIndexMap const& _local_to_global;
+    NumLib::LocalToGlobalIndexMap const& _local_to_global;
 
     //! Avoids frequent reallocations.
     Eigen::MatrixXd _local_matrix_cache;
diff --git a/ProcessLib/GlobalSetup.h b/NumLib/GlobalSetup.h
similarity index 91%
rename from ProcessLib/GlobalSetup.h
rename to NumLib/GlobalSetup.h
index 7e30d81769d3e6f47600dfe571ce26350b323ca0..27bed2596c65753701e9d77d91d4938bc8c0af2b 100644
--- a/ProcessLib/GlobalSetup.h
+++ b/NumLib/GlobalSetup.h
@@ -7,12 +7,12 @@
  *
  */
 
-#ifndef PROCESS_LIB_GLOBAL_SETUP_H_
-#define PROCESS_LIB_GLOBAL_SETUP_H_
+#ifndef NUMLIB_GLOBAL_SETUP_H_
+#define NUMLIB_GLOBAL_SETUP_H_
 
 #include <functional>
 
-namespace ProcessLib
+namespace NumLib
 {
 
 /// The GlobalSetup collects vector and matrix builder and corresponding global
@@ -64,6 +64,6 @@ struct GlobalSetup
     GlobalSetup() = delete;
 };
 
-}   // namespace ProcessLib
+}   // namespace NumLib
 
-#endif  // PROCESS_LIB_GLOBAL_SETUP_H_
+#endif  // NUMLIB_GLOBAL_SETUP_H_
diff --git a/NumLib/IndexValueVector.h b/NumLib/IndexValueVector.h
new file mode 100644
index 0000000000000000000000000000000000000000..a6bdc0cb18dbadc174ee69f24c477d71c398afda
--- /dev/null
+++ b/NumLib/IndexValueVector.h
@@ -0,0 +1,27 @@
+/**
+ * \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 NUM_LIB_INDEXVALUEVECTOR_H
+#define NUM_LIB_INDEXVALUEVECTOR_H
+
+#include <vector>
+
+namespace NumLib
+{
+
+template <typename IndexType>
+struct IndexValueVector final
+{
+    std::vector<IndexType> ids;
+    std::vector<double> values;
+};
+
+}
+
+#endif  // NUM_LIB_INDEXVALUEVECTOR_H
diff --git a/ProcessLib/NumericsConfig.h b/NumLib/NumericsConfig.h
similarity index 92%
rename from ProcessLib/NumericsConfig.h
rename to NumLib/NumericsConfig.h
index dcdba10ec651473f5082158303a86f85461b699a..43d22fa4734764d3edd9b7980e43007306aaedb9 100644
--- a/ProcessLib/NumericsConfig.h
+++ b/NumLib/NumericsConfig.h
@@ -12,6 +12,8 @@
 
 #include <type_traits>
 
+#include "MathLib/LinAlg/SparsityPattern.h"
+
 /**
  * This file provides a configuration of the global matrix/vector and
  * corresponding linear solver, and the global executer types.
@@ -82,11 +84,11 @@ namespace detail
 // Global vector/matrix builder.
 //
 
-#include "AssemblerLib/VectorMatrixBuilder.h"
+#include "NumLib/Assembler/VectorMatrixBuilder.h"
 namespace detail
 {
 using GlobalVectorMatrixBuilderType =
-        AssemblerLib::VectorMatrixBuilder<
+        NumLib::VectorMatrixBuilder<
             GlobalMatrixType,
             GlobalVectorType>;
 }
@@ -94,10 +96,10 @@ using GlobalVectorMatrixBuilderType =
 //
 // Global executor
 //
-#include "AssemblerLib/SerialExecutor.h"
+#include "NumLib/Assembler/SerialExecutor.h"
 namespace detail
 {
-using GlobalExecutorType = AssemblerLib::SerialExecutor;
+using GlobalExecutorType = NumLib::SerialExecutor;
 }
 
 ///
@@ -105,7 +107,7 @@ using GlobalExecutorType = AssemblerLib::SerialExecutor;
 ///
 #include "GlobalSetup.h"
 using GlobalSetupType =
-    ProcessLib::GlobalSetup<
+    NumLib::GlobalSetup<
         detail::GlobalVectorMatrixBuilderType,
         detail::GlobalExecutorType,
         detail::LinearSolverType>;
@@ -128,5 +130,6 @@ static_assert(std::is_same<detail::GlobalMatrixType::IndexType,
 /// 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_
diff --git a/NumLib/ODESolver/EquationSystem.h b/NumLib/ODESolver/EquationSystem.h
index 42635717116e861d0478a5f71dcc51a8de7fb33f..1856bd17abc08fe704932160b70f71ee001ffbbb 100644
--- a/NumLib/ODESolver/EquationSystem.h
+++ b/NumLib/ODESolver/EquationSystem.h
@@ -10,7 +10,7 @@
 #ifndef NUMLIB_EQUATIONSYSTEM_H
 #define NUMLIB_EQUATIONSYSTEM_H
 
-#include "MathLib/LinAlg/MatrixProviderUser.h"
+#include "NumLib/DOF/MatrixProviderUser.h"
 
 namespace NumLib
 {
diff --git a/NumLib/ODESolver/NonlinearSolver-impl.h b/NumLib/ODESolver/NonlinearSolver-impl.h
index 5876ff11b79793b4b3c5b4b170d8b6f742c6aee8..1d356d06e89c9dd8f43bc1abfed7f8da3e83b8ba 100644
--- a/NumLib/ODESolver/NonlinearSolver-impl.h
+++ b/NumLib/ODESolver/NonlinearSolver-impl.h
@@ -14,11 +14,11 @@
 
 #include "BaseLib/ConfigTree.h"
 #include "MathLib/LinAlg/BLAS.h"
-#include "MathLib/LinAlg/GlobalMatrixProviders.h"
+#include "MathLib/LinAlg/VectorNorms.h"
+#include "NumLib/DOF/GlobalMatrixProviders.h"
 
 #include "NonlinearSolver.h"
 
-#include "MathLib/LinAlg/VectorNorms.h"
 
 namespace NumLib
 {
diff --git a/NumLib/ODESolver/ODESystem.h b/NumLib/ODESolver/ODESystem.h
index 1183c5bf1125661882bc9657fcfc4a27503c9682..697e198dc1f5cc44db78f130de0474c42605d6d9 100644
--- a/NumLib/ODESolver/ODESystem.h
+++ b/NumLib/ODESolver/ODESystem.h
@@ -10,18 +10,12 @@
 #ifndef NUMLIB_ODESYSTEM_H
 #define NUMLIB_ODESYSTEM_H
 
-#include "MathLib/LinAlg/MatrixVectorTraits.h"
+#include "NumLib/IndexValueVector.h"
+#include "NumLib/DOF/MatrixVectorTraits.h"
 
 #include "Types.h"
 #include "EquationSystem.h"
 
-// TODO move to other namespace
-namespace ProcessLib
-{
-template <typename IndexType>
-struct DirichletBc;
-}
-
 namespace NumLib
 {
 
@@ -67,7 +61,7 @@ public:
 
     //! Provides known solutions (Dirichlet boundary conditions) vector for
     //! the ode system at the given time \c t.
-    virtual std::vector<ProcessLib::DirichletBc<Index>> const*
+    virtual std::vector<NumLib::IndexValueVector<Index>> const*
     getKnownSolutions(double const t) const
     {
         (void)t;
diff --git a/NumLib/ODESolver/TimeDiscretization.h b/NumLib/ODESolver/TimeDiscretization.h
index 63fa996c982cb9fe86786c0b3641c9f839a67afa..e60165cb2604c7312e6ec498756ca015c866961e 100644
--- a/NumLib/ODESolver/TimeDiscretization.h
+++ b/NumLib/ODESolver/TimeDiscretization.h
@@ -13,7 +13,7 @@
 #include <vector>
 
 #include "MathLib/LinAlg/BLAS.h"
-#include "MathLib/LinAlg/GlobalMatrixProviders.h"
+#include "NumLib/DOF/GlobalMatrixProviders.h"
 #include "Types.h"
 
 
diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
index 916851ed999217a52ea60114bf10fed1589d8597..da9b2e81af28c9979ecb53babed226ef7f7b3881 100644
--- a/NumLib/ODESolver/TimeDiscretizedODESystem.h
+++ b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -13,8 +13,8 @@
 #include <memory>
 
 #include "MathLib/LinAlg/ApplyKnownSolution.h"
-#include "MathLib/LinAlg/UnifiedMatrixSetters.h"
-#include "ProcessLib/DirichletBc.h"
+#include "NumLib/IndexValueVector.h"
+#include "NumLib/DOF/UnifiedMatrixSetters.h"
 
 #include "ODESystem.h"
 #include "NonlinearSystem.h"
@@ -31,9 +31,9 @@ void applyKnownSolutions(std::vector<Solutions> const*const known_solutions,
     if (!known_solutions) return;
 
     for (auto const& bc : *known_solutions) {
-        for (std::size_t i=0; i<bc.global_ids.size(); ++i) {
+        for (std::size_t i=0; i<bc.ids.size(); ++i) {
             // TODO that might have bad performance for some Vector types, e.g., PETSc.
-            MathLib::setVector(x, bc.global_ids[i], bc.values[i]);
+            MathLib::setVector(x, bc.ids[i], bc.values[i]);
         }
     }
 }
@@ -212,7 +212,7 @@ public:
                 values.resize(bc.values.size(), 0.0);
 
                 // TODO maybe it would be faster to apply all at once
-                MathLib::applyKnownSolution(Jac, res, minus_delta_x, bc.global_ids, values);
+                MathLib::applyKnownSolution(Jac, res, minus_delta_x, bc.ids, values);
             }
         }
     }
@@ -361,7 +361,7 @@ public:
         if (known_solutions) {
             for (auto const& bc : *known_solutions) {
                 // TODO maybe it would be faster to apply all at once
-                MathLib::applyKnownSolution(A, rhs, x, bc.global_ids, bc.values);
+                MathLib::applyKnownSolution(A, rhs, x, bc.ids, bc.values);
             }
         }
     }
diff --git a/NumLib/ODESolver/TimeLoopSingleODE.h b/NumLib/ODESolver/TimeLoopSingleODE.h
index ff8f29d8e086bda2b65dc8b9b7c5fbfaff4f6a36..f7c45ad92ce0e61b2263cebbb11efca5ad1690fb 100644
--- a/NumLib/ODESolver/TimeLoopSingleODE.h
+++ b/NumLib/ODESolver/TimeLoopSingleODE.h
@@ -10,7 +10,7 @@
 #ifndef NUMLIB_TIMELOOP_H
 #define NUMLIB_TIMELOOP_H
 
-#include "MathLib/LinAlg/GlobalMatrixProviders.h"
+#include "NumLib/DOF/GlobalMatrixProviders.h"
 #include "TimeDiscretizedODESystem.h"
 #include "NonlinearSolver.h"
 
diff --git a/ProcessLib/CMakeLists.txt b/ProcessLib/CMakeLists.txt
index f60e5716d6a120f0229e5f1dd282e6eb513f68cf..280aac21594086dac29af104c64def5214621b96 100644
--- a/ProcessLib/CMakeLists.txt
+++ b/ProcessLib/CMakeLists.txt
@@ -12,7 +12,6 @@ APPEND_SOURCE_FILES(SOURCES TES)
 add_library(ProcessLib ${SOURCES})
 
 target_link_libraries(ProcessLib
-    AssemblerLib
     MaterialsLib
     MeshGeoToolsLib
     NumLib # for shape matrices
diff --git a/ProcessLib/DirichletBc.h b/ProcessLib/DirichletBc.h
index 1e608dbaa9ea199c55057e42f8b96ce9bc044f0f..924fb1a2d2e1b5bd953546974041fb20dc2bbcd8 100644
--- a/ProcessLib/DirichletBc.h
+++ b/ProcessLib/DirichletBc.h
@@ -12,7 +12,7 @@
 
 #include <vector>
 
-#include "MeshLib/MeshSearch/NodeSearch.h"
+#include "NumLib/IndexValueVector.h"
 
 namespace ProcessLib
 {
@@ -20,11 +20,7 @@ namespace ProcessLib
 /// A dirichlet boundary condition is represented by a list of global indices
 /// with corresponding values.
 template <typename IndexType>
-struct DirichletBc final
-{
-    std::vector<IndexType> global_ids;
-    std::vector<double> values;
-};
+using DirichletBc = NumLib::IndexValueVector<IndexType>;
 
 }  // namespace ProcessLib
 
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
index 331bad5a960439e03dc550aaaad0f2ccf4e9caf6..96fb30f1a55c646945a3a03ff6c15bc9d6ece4c0 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
@@ -105,7 +105,7 @@ public:
         }
     }
 
-    void addToGlobal(AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+    void addToGlobal(NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
         GlobalMatrix& /*M*/, GlobalMatrix& K, GlobalVector& b)
         const override
     {
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess-fwd.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess-fwd.h
index 57c6516da9fb9009409dd88ca5dafafb4437b869..d5e33c7c6a8367f00c0a4bd0e124dbada30cb63f 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess-fwd.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess-fwd.h
@@ -11,7 +11,7 @@
 #define PROCESS_LIB_GROUNDWATERFLOWPROCESS_FWD_H_
 
 #include "GroundwaterFlowProcess.h"
-#include "ProcessLib/NumericsConfig.h"
+#include "NumLib/NumericsConfig.h"
 
 extern template class ProcessLib::GroundwaterFlow::GroundwaterFlowProcess<GlobalSetupType>;
 
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
index 603c4978c4db0de2bf8d155ec1964a643a1f041c..8f6bc793428ae11088e5b588b3158156ffdfb54e 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowProcess.h
@@ -12,7 +12,7 @@
 
 #include <cassert>
 
-#include "AssemblerLib/VectorMatrixAssembler.h"
+#include "NumLib/Assembler/VectorMatrixAssembler.h"
 #include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
 #include "ProcessLib/Process.h"
 #include "ProcessLib/Utils/CreateLocalAssemblers.h"
@@ -78,7 +78,7 @@ private:
     using LocalAssemblerInterface =
         GroundwaterFlowLocalAssemblerInterface<GlobalMatrix, GlobalVector>;
 
-    using GlobalAssembler = AssemblerLib::VectorMatrixAssembler<
+    using GlobalAssembler = NumLib::VectorMatrixAssembler<
             GlobalMatrix, GlobalVector, LocalAssemblerInterface,
             NumLib::ODESystemTag::FirstOrderImplicitQuasilinear>;
 
@@ -88,7 +88,7 @@ private:
         GlobalVector, IntegrationPointValue, LocalAssemblerInterface>;
 
     void initializeConcreteProcess(
-            AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+            NumLib::LocalToGlobalIndexMap const& dof_table,
             MeshLib::Mesh const& mesh,
             unsigned const integration_order) override
     {
diff --git a/ProcessLib/LocalAssemblerInterface.h b/ProcessLib/LocalAssemblerInterface.h
index 063002f80b34d0b29b17752e4aaa46f1a4957794..788ba9fbab39c4cb48b2e632a659e7c0e448b17a 100644
--- a/ProcessLib/LocalAssemblerInterface.h
+++ b/ProcessLib/LocalAssemblerInterface.h
@@ -27,7 +27,7 @@ public:
     virtual void assemble(double const t, std::vector<double> const& local_x) = 0;
 
     virtual void addToGlobal(
-        AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const&,
+        NumLib::LocalToGlobalIndexMap::RowColumnIndices const&,
         GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) const = 0;
 
     virtual void assembleJacobian(double const /*t*/,
@@ -39,7 +39,7 @@ public:
         std::abort();
     }
 
-    virtual void addJacobianToGlobal(AssemblerLib::LocalToGlobalIndexMap::
+    virtual void addJacobianToGlobal(NumLib::LocalToGlobalIndexMap::
                                          RowColumnIndices const& /*indices*/,
                                      GlobalMatrix& /*Jac*/) const
     {
diff --git a/ProcessLib/NeumannBc.h b/ProcessLib/NeumannBc.h
index 4d841002a2647b69b92dcefddcbe474c9aaa88c3..f52ad1fd403a05b32e37b4f7dea8bb0f6198bf12 100644
--- a/ProcessLib/NeumannBc.h
+++ b/ProcessLib/NeumannBc.h
@@ -16,7 +16,7 @@
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 
-#include "AssemblerLib/VectorMatrixAssembler.h"
+#include "NumLib/Assembler/VectorMatrixAssembler.h"
 #include "MeshLib/MeshSubset.h"
 #include "MeshLib/MeshSearch/NodeSearch.h"
 
@@ -54,7 +54,7 @@ public:
     NeumannBc(
         NeumannBcConfig const& bc,
         unsigned const integration_order,
-        AssemblerLib::LocalToGlobalIndexMap const& local_to_global_index_map,
+        NumLib::LocalToGlobalIndexMap const& local_to_global_index_map,
         std::size_t const variable_id,
         std::size_t const component_id)
         : _function(*bc.getFunction()),
@@ -136,7 +136,7 @@ private:
 
     /// Local dof table, a subset of the global one restricted to the
     /// participating #_elements of the boundary condition.
-    std::unique_ptr<AssemblerLib::LocalToGlobalIndexMap> _local_to_global_index_map;
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap> _local_to_global_index_map;
 
     /// Integration order for integration over the lower-dimensional elements of
     /// the #_function.
@@ -146,7 +146,7 @@ private:
         GlobalMatrix, GlobalVector>;
 
     using GlobalAssembler =
-        AssemblerLib::VectorMatrixAssembler<
+        NumLib::VectorMatrixAssembler<
             GlobalMatrix, GlobalVector, LocalAssembler,
             NumLib::ODESystemTag::NeumannBC>;
 
diff --git a/ProcessLib/NeumannBcAssembler.h b/ProcessLib/NeumannBcAssembler.h
index fe124ef9996aa3659552fd23af1c3e0eddf672c1..3d77a029e6bf7eb3ab9a4bdde6cd7da8173264af 100644
--- a/ProcessLib/NeumannBcAssembler.h
+++ b/ProcessLib/NeumannBcAssembler.h
@@ -27,7 +27,7 @@ public:
 
     virtual void assemble(const double t) = 0;
 
-    virtual void addToGlobal(AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const&,
+    virtual void addToGlobal(NumLib::LocalToGlobalIndexMap::RowColumnIndices const&,
                              GlobalVector& b) const = 0;
 };
 
@@ -95,7 +95,7 @@ public:
         }
     }
 
-    void addToGlobal(AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+    void addToGlobal(NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
                      GlobalVector& b) const override
     {
         b.add(indices.rows, *_localRhs);
diff --git a/ProcessLib/Process.h b/ProcessLib/Process.h
index 8e142fe0f85c1f8b7333adaad487b715f05ee918..f543196a2f873020854f3295f1ce6971d8c13b67 100644
--- a/ProcessLib/Process.h
+++ b/ProcessLib/Process.h
@@ -10,8 +10,8 @@
 #ifndef PROCESS_LIB_PROCESS_H_
 #define PROCESS_LIB_PROCESS_H_
 
-#include "AssemblerLib/ComputeSparsityPattern.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
+#include "NumLib/DOF/ComputeSparsityPattern.h"
 #include "NumLib/ODESolver/ODESystem.h"
 #include "NumLib/ODESolver/TimeDiscretization.h"
 #include "NumLib/ODESolver/NonlinearSolver.h"
@@ -190,7 +190,7 @@ public:
 private:
     /// Process specific initialization called by initialize().
     virtual void initializeConcreteProcess(
-        AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
         MeshLib::Mesh const& mesh,
         unsigned const integration_order) = 0;
 
@@ -232,9 +232,9 @@ private:
         }
 
         _local_to_global_index_map.reset(
-            new AssemblerLib::LocalToGlobalIndexMap(
+            new NumLib::LocalToGlobalIndexMap(
                 std::move(all_mesh_subsets),
-                AssemblerLib::ComponentOrder::BY_LOCATION));
+                NumLib::ComponentOrder::BY_LOCATION));
     }
 
     /// Sets the initial condition values in the solution vector x for a given
@@ -304,7 +304,7 @@ private:
     /// DOF-table.
     void computeSparsityPattern()
     {
-        _sparsity_pattern = std::move(AssemblerLib::computeSparsityPattern(
+        _sparsity_pattern = std::move(NumLib::computeSparsityPattern(
             *_local_to_global_index_map, _mesh));
     }
 
@@ -312,7 +312,7 @@ protected:
     MeshLib::Mesh& _mesh;
     std::unique_ptr<MeshLib::MeshSubset const> _mesh_subset_all_nodes;
 
-    std::unique_ptr<AssemblerLib::LocalToGlobalIndexMap>
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap>
         _local_to_global_index_map;
 
     SecondaryVariableCollection<GlobalVector> _secondary_variables;
@@ -320,7 +320,7 @@ protected:
 
 private:
     unsigned const _integration_order = 2;
-    MathLib::SparsityPattern _sparsity_pattern;
+    GlobalSparsityPattern _sparsity_pattern;
 
     std::vector<DirichletBc<GlobalIndexType>> _dirichlet_bcs;
     std::vector<std::unique_ptr<NeumannBc<GlobalSetup>>> _neumann_bcs;
diff --git a/ProcessLib/ProcessOutput.h b/ProcessLib/ProcessOutput.h
index 0e29556db7bf46f6eb3a1c63e900f9ea6032e07c..9757f6b6e487675d5c687594a6e8e5a837a3f5b1 100644
--- a/ProcessLib/ProcessOutput.h
+++ b/ProcessLib/ProcessOutput.h
@@ -88,7 +88,7 @@ void doProcessOutput(
         std::string const& file_name,
         GlobalVector const& x,
         MeshLib::Mesh& mesh,
-        AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
         std::vector<std::reference_wrapper<ProcessVariable>> const&
         process_variables,
         SecondaryVariableCollection<GlobalVector> secondary_variables,
diff --git a/ProcessLib/ProcessVariable.h b/ProcessLib/ProcessVariable.h
index 6386cbb396e45e6c4042136803e6d08255b8ca65..3a735d3c780aac8b7e3cdb0332d456381b5779e2 100644
--- a/ProcessLib/ProcessVariable.h
+++ b/ProcessLib/ProcessVariable.h
@@ -64,7 +64,7 @@ public:
     void initializeDirichletBCs(
         OutputIterator output_bcs,
         MeshGeoToolsLib::MeshNodeSearcher& searcher,
-        const AssemblerLib::LocalToGlobalIndexMap& dof_table,
+        const NumLib::LocalToGlobalIndexMap& dof_table,
         const unsigned nodal_dof_idx)
     {
         for (auto& bc_config : _dirichlet_bc_configs)
diff --git a/ProcessLib/SecondaryVariable.h b/ProcessLib/SecondaryVariable.h
index 7e3cbfbef0c2ccd7968baa9a8444e5a44f53a172..9983763728d8d3d491c3278e45b60925f9d36279 100644
--- a/ProcessLib/SecondaryVariable.h
+++ b/ProcessLib/SecondaryVariable.h
@@ -14,7 +14,7 @@
 #include "BaseLib/uniqueInsert.h"
 #include "NumLib/Extrapolation/Extrapolator.h"
 
-namespace AssemblerLib { class LocalToGlobalIndexMap; }
+namespace NumLib { class LocalToGlobalIndexMap; }
 
 namespace ProcessLib
 {
@@ -36,7 +36,7 @@ struct SecondaryVariableFunctions final
      */
     using Function = std::function<GlobalVector const&(
         GlobalVector const& x,
-        AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
         std::unique_ptr<GlobalVector>& result_cache)>;
 
     SecondaryVariableFunctions() = default;
@@ -49,7 +49,7 @@ struct SecondaryVariableFunctions final
         // Used to detect nasty implicit conversions.
         static_assert(std::is_same<GlobalVector const&,
             typename std::result_of<F1(
-                GlobalVector const&, AssemblerLib::LocalToGlobalIndexMap const&,
+                GlobalVector const&, NumLib::LocalToGlobalIndexMap const&,
                 std::unique_ptr<GlobalVector>&
                 )>::type>::value,
             "The function eval_field_ does not return a const reference"
@@ -57,7 +57,7 @@ struct SecondaryVariableFunctions final
 
         static_assert(std::is_same<GlobalVector const&,
             typename std::result_of<F2(
-                GlobalVector const&, AssemblerLib::LocalToGlobalIndexMap const&,
+                GlobalVector const&, NumLib::LocalToGlobalIndexMap const&,
                 std::unique_ptr<GlobalVector>&
             )>::type>::value,
             "The function eval_residuals_ does not return a const reference"
@@ -71,7 +71,7 @@ struct SecondaryVariableFunctions final
         // Used to detect nasty implicit conversions.
         static_assert(std::is_same<GlobalVector const&,
             typename std::result_of<F1(
-                GlobalVector const&, AssemblerLib::LocalToGlobalIndexMap const&,
+                GlobalVector const&, NumLib::LocalToGlobalIndexMap const&,
                 std::unique_ptr<GlobalVector>&
                 )>::type>::value,
             "The function eval_field_ does not return a const reference"
@@ -212,7 +212,7 @@ makeExtrapolator(PropertyEnum const property,
 
     auto const eval_field = [property, &extrapolator, &local_assemblers](
             GlobalVector const& /*x*/,
-            AssemblerLib::LocalToGlobalIndexMap const& /*dof_table*/,
+            NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
             std::unique_ptr<GlobalVector>& /*result_cache*/
             ) -> GlobalVector const&
     {
@@ -222,7 +222,7 @@ makeExtrapolator(PropertyEnum const property,
 
     auto const eval_residuals = [property, &extrapolator, &local_assemblers](
             GlobalVector const& /*x*/,
-            AssemblerLib::LocalToGlobalIndexMap const& /*dof_table*/,
+            NumLib::LocalToGlobalIndexMap const& /*dof_table*/,
             std::unique_ptr<GlobalVector>& /*result_cache*/
             ) -> GlobalVector const&
     {
diff --git a/ProcessLib/TES/TESLocalAssembler-impl.h b/ProcessLib/TES/TESLocalAssembler-impl.h
index 3b9ad99e5c296fe86dda8f522bf79ba0e02e9acb..c5894a8c1e593196514d5067ee1e67d7b8b24231 100644
--- a/ProcessLib/TES/TESLocalAssembler-impl.h
+++ b/ProcessLib/TES/TESLocalAssembler-impl.h
@@ -190,7 +190,7 @@ template <typename ShapeFunction_, typename IntegrationMethod_,
 void TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalMatrix,
                        GlobalVector, GlobalDim>::
     addToGlobal(
-        AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+        NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
         GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) const
 {
     M.add(indices, _local_M);
diff --git a/ProcessLib/TES/TESLocalAssembler.h b/ProcessLib/TES/TESLocalAssembler.h
index 68da0989069180c12cc5ee63299b0e82e76068cd..cc80fc33e7fce097d3e79eef1e7db13b97b2e54c 100644
--- a/ProcessLib/TES/TESLocalAssembler.h
+++ b/ProcessLib/TES/TESLocalAssembler.h
@@ -33,7 +33,7 @@ public:
                           std::vector<double> const& local_x) = 0;
 
     virtual void addToGlobal(
-        AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const&,
+        NumLib::LocalToGlobalIndexMap::RowColumnIndices const&,
         GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) const = 0;
 
     virtual bool checkBounds(std::vector<double> const& local_x,
@@ -58,7 +58,7 @@ public:
     void assemble(double const t, std::vector<double> const& local_x) override;
 
     void addToGlobal(
-        AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+        NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
         GlobalMatrix& M, GlobalMatrix& K, GlobalVector& b) const override;
 
     Eigen::Map<const Eigen::VectorXd> getShapeMatrix(
diff --git a/ProcessLib/TES/TESProcess.cpp b/ProcessLib/TES/TESProcess.cpp
index 930bc4d4c925d4210b03fd900449f8598ff249ae..48b887a70609bfafec803d358367e618891728e5 100644
--- a/ProcessLib/TES/TESProcess.cpp
+++ b/ProcessLib/TES/TESProcess.cpp
@@ -10,9 +10,9 @@
 #include "TESProcess.h"
 
 // TODO Copied from VectorMatrixAssembler. Could be provided by the DOF table.
-inline AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices
+inline NumLib::LocalToGlobalIndexMap::RowColumnIndices
 getRowColumnIndices_(std::size_t const id,
-                     AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+                     NumLib::LocalToGlobalIndexMap const& dof_table,
                      std::vector<GlobalIndexType>& indices)
 {
     assert(dof_table.size() > id);
@@ -27,14 +27,14 @@ getRowColumnIndices_(std::size_t const id,
         indices.insert(indices.end(), idcs.begin(), idcs.end());
     }
 
-    return AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices(indices,
+    return NumLib::LocalToGlobalIndexMap::RowColumnIndices(indices,
                                                                  indices);
 }
 
 template <typename GlobalVector>
 void getVectorValues(
     GlobalVector const& x,
-    AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const& r_c_indices,
+    NumLib::LocalToGlobalIndexMap::RowColumnIndices const& r_c_indices,
     std::vector<double>& local_x)
 {
     auto const& indices = r_c_indices.rows;
@@ -50,7 +50,7 @@ void getVectorValues(
 // TODO that essentially duplicates code which is also present in ProcessOutput.
 template <typename GlobalVector>
 double getNodalValue(GlobalVector const& x, MeshLib::Mesh const& mesh,
-                     AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+                     NumLib::LocalToGlobalIndexMap const& dof_table,
                      std::size_t const node_id,
                      std::size_t const global_component_id)
 {
@@ -170,7 +170,7 @@ TESProcess<GlobalSetup>::TESProcess(
 
 template <typename GlobalSetup>
 void TESProcess<GlobalSetup>::initializeConcreteProcess(
-    AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+    NumLib::LocalToGlobalIndexMap const& dof_table,
     MeshLib::Mesh const& mesh, unsigned const integration_order)
 {
     DBUG("Create global assembler.");
@@ -187,10 +187,10 @@ void TESProcess<GlobalSetup>::initializeConcreteProcess(
     all_mesh_subsets_single_component.emplace_back(
         new MeshLib::MeshSubsets(this->_mesh_subset_all_nodes.get()));
     _local_to_global_index_map_single_component.reset(
-        new AssemblerLib::LocalToGlobalIndexMap(
+        new NumLib::LocalToGlobalIndexMap(
             std::move(all_mesh_subsets_single_component),
             // by location order is needed for output
-            AssemblerLib::ComponentOrder::BY_LOCATION));
+            NumLib::ComponentOrder::BY_LOCATION));
 
     _extrapolator.reset(
         new ExtrapolatorImplementation(MathLib::MatrixSpecifications(
@@ -338,7 +338,7 @@ template <typename GlobalSetup>
 typename TESProcess<GlobalSetup>::GlobalVector const&
 TESProcess<GlobalSetup>::computeVapourPartialPressure(
     typename TESProcess::GlobalVector const& x,
-    AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+    NumLib::LocalToGlobalIndexMap const& dof_table,
     std::unique_ptr<GlobalVector>& result_cache)
 {
     assert(&dof_table == this->_local_to_global_index_map.get());
@@ -370,7 +370,7 @@ template <typename GlobalSetup>
 typename TESProcess<GlobalSetup>::GlobalVector const&
 TESProcess<GlobalSetup>::computeRelativeHumidity(
     typename TESProcess::GlobalVector const& x,
-    AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+    NumLib::LocalToGlobalIndexMap const& dof_table,
     std::unique_ptr<GlobalVector>& result_cache)
 {
     assert(&dof_table == this->_local_to_global_index_map.get());
@@ -407,7 +407,7 @@ template <typename GlobalSetup>
 typename TESProcess<GlobalSetup>::GlobalVector const&
 TESProcess<GlobalSetup>::computeEquilibriumLoading(
     typename TESProcess::GlobalVector const& x,
-    AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+    NumLib::LocalToGlobalIndexMap const& dof_table,
     std::unique_ptr<GlobalVector>& result_cache)
 {
     assert(&dof_table == this->_local_to_global_index_map.get());
diff --git a/ProcessLib/TES/TESProcess.h b/ProcessLib/TES/TESProcess.h
index 5149b329af517daccfd47661aae521601d78a2ba..5034eff13eb271a70b0f9a3443a480f0959ce340 100644
--- a/ProcessLib/TES/TESProcess.h
+++ b/ProcessLib/TES/TESProcess.h
@@ -10,7 +10,7 @@
 #ifndef PROCESS_LIB_TESPROCESS_H_
 #define PROCESS_LIB_TESPROCESS_H_
 
-#include "AssemblerLib/LocalToGlobalIndexMap.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
 #include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
 #include "ProcessLib/Process.h"
 
@@ -59,7 +59,7 @@ private:
     using LocalAssembler =
         TESLocalAssemblerInterface<GlobalMatrix, GlobalVector>;
 
-    using GlobalAssembler = AssemblerLib::VectorMatrixAssembler<
+    using GlobalAssembler = NumLib::VectorMatrixAssembler<
         GlobalMatrix, GlobalVector, LocalAssembler,
         NumLib::ODESystemTag::FirstOrderImplicitQuasilinear>;
 
@@ -70,7 +70,7 @@ private:
             GlobalVector, TESIntPtVariables, LocalAssembler>;
 
     void initializeConcreteProcess(
-        AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
         MeshLib::Mesh const& mesh, unsigned const integration_order) override;
 
     void assembleConcreteProcess(const double t, GlobalVector const& x,
@@ -79,17 +79,17 @@ private:
 
     GlobalVector const& computeVapourPartialPressure(
         GlobalVector const& x,
-        AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
         std::unique_ptr<GlobalVector>& result_cache);
 
     GlobalVector const& computeRelativeHumidity(
         GlobalVector const& x,
-        AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
         std::unique_ptr<GlobalVector>& result_cache);
 
     GlobalVector const& computeEquilibriumLoading(
         GlobalVector const& x,
-        AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
         std::unique_ptr<GlobalVector>& result_cache);
 
     std::unique_ptr<GlobalAssembler> _global_assembler;
@@ -97,7 +97,7 @@ private:
 
     AssemblyParams _assembly_params;
 
-    std::unique_ptr<AssemblerLib::LocalToGlobalIndexMap>
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap>
         _local_to_global_index_map_single_component;
 
     std::unique_ptr<ExtrapolatorInterface> _extrapolator;
diff --git a/ProcessLib/UniformDirichletBoundaryCondition.h b/ProcessLib/UniformDirichletBoundaryCondition.h
index 80627072198f72420be56692d51040a9860a8468..f614c40ca426eebfbe025429061d4a83302e7bda 100644
--- a/ProcessLib/UniformDirichletBoundaryCondition.h
+++ b/ProcessLib/UniformDirichletBoundaryCondition.h
@@ -13,13 +13,12 @@
 #include <algorithm>
 #include <vector>
 
-#include "logog/include/logog.hpp"
-
-#include "NumericsConfig.h" // for GlobalIndexType
+#include <logog/include/logog.hpp>
 
 #include "BaseLib/ConfigTree.h"
-#include "AssemblerLib/LocalToGlobalIndexMap.h"
 #include "MeshGeoToolsLib/MeshNodeSearcher.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/NumericsConfig.h" // for GlobalIndexType
 
 #include "DirichletBc.h"
 
@@ -56,7 +55,7 @@ public:
     /// object.
     void initialize(
             MeshGeoToolsLib::MeshNodeSearcher& searcher,
-            AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+            NumLib::LocalToGlobalIndexMap const& dof_table,
             std::size_t component_id,
             DirichletBc<GlobalIndexType>& bc)
     {
@@ -65,7 +64,7 @@ public:
         std::vector<std::size_t> ids = searcher.getMeshNodeIDs(*_geometry);
 
         // convert mesh node ids to global index for the given component
-        bc.global_ids.reserve(bc.global_ids.size() + ids.size());
+        bc.ids.reserve(bc.ids.size() + ids.size());
         bc.values.reserve(bc.values.size() + ids.size());
         for (auto& id : ids)
         {
@@ -82,7 +81,7 @@ public:
             // PETSc routines. Therefore, the following if-condition is applied.
             if (g_idx >= 0)
             {
-                bc.global_ids.emplace_back(g_idx);
+                bc.ids.emplace_back(g_idx);
                 bc.values.emplace_back(_value);
             }
         }
diff --git a/ProcessLib/Utils/CreateLocalAssemblers.h b/ProcessLib/Utils/CreateLocalAssemblers.h
index 8d61deb840129a95b5fd2fcc33ff15282f8a1e3f..081478cee2aea9f9a4d802e4585822c79501a66a 100644
--- a/ProcessLib/Utils/CreateLocalAssemblers.h
+++ b/ProcessLib/Utils/CreateLocalAssemblers.h
@@ -12,7 +12,7 @@
 #include <vector>
 #include <logog/include/logog.hpp>
 
-#include "AssemblerLib/LocalToGlobalIndexMap.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
 
 #include "LocalDataInitializer.h"
 
@@ -29,7 +29,7 @@ template<unsigned GlobalDim, typename GlobalSetup,
          typename LocalAssemblerInterface,
          typename... ExtraCtorArgs>
 void createLocalAssemblers(
-        AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
         std::vector<MeshLib::Element*> const& mesh_elements,
         unsigned const integration_order,
         std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
@@ -83,7 +83,7 @@ template<typename GlobalSetup,
 void createLocalAssemblers(
         const unsigned dimension,
         std::vector<MeshLib::Element*> const& mesh_elements,
-        AssemblerLib::LocalToGlobalIndexMap const& dof_table,
+        NumLib::LocalToGlobalIndexMap const& dof_table,
         unsigned const integration_order,
         std::vector<std::unique_ptr<LocalAssemblerInterface>>& local_assemblers,
         ExtraCtorArgs&&... extra_ctor_args
diff --git a/ProcessLib/Utils/LocalDataInitializer.h b/ProcessLib/Utils/LocalDataInitializer.h
index 9b464d1bd687d1b8008d2d83d6b7bf45361d815c..0172936204529718d967f800412afd9bc5e638c0 100644
--- a/ProcessLib/Utils/LocalDataInitializer.h
+++ b/ProcessLib/Utils/LocalDataInitializer.h
@@ -16,8 +16,8 @@
 #include <typeinfo>
 #include <unordered_map>
 
-#include "AssemblerLib/LocalToGlobalIndexMap.h"
 #include "MeshLib/Elements/Elements.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
 #include "NumLib/Fem/Integration/GaussIntegrationPolicy.h"
 
 
@@ -133,7 +133,7 @@ public:
     using LADataIntfPtr = std::unique_ptr<LocalAssemblerInterface>;
 
     explicit LocalDataInitializer(
-        AssemblerLib::LocalToGlobalIndexMap const& dof_table)
+        NumLib::LocalToGlobalIndexMap const& dof_table)
         : _dof_table(dof_table)
     {
         // /// Lines and points ///////////////////////////////////
@@ -308,7 +308,7 @@ private:
     /// Mapping of element types to local assembler constructors.
     std::unordered_map<std::type_index, LADataBuilder> _builder;
 
-    AssemblerLib::LocalToGlobalIndexMap const& _dof_table;
+    NumLib::LocalToGlobalIndexMap const& _dof_table;
 };
 
 }   // namespace ProcessLib
diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt
index 6c3ca22f10f333e19c4efc3a02ba20393510693a..86c9fffdb30c1aed0e345c9ffecd736f69734fd1 100644
--- a/Tests/CMakeLists.txt
+++ b/Tests/CMakeLists.txt
@@ -8,12 +8,10 @@ endif()
 
 APPEND_SOURCE_FILES(TEST_SOURCES)
 APPEND_SOURCE_FILES(TEST_SOURCES ApplicationsLib)
-APPEND_SOURCE_FILES(TEST_SOURCES AssemblerLib)
 APPEND_SOURCE_FILES(TEST_SOURCES BaseLib)
 APPEND_SOURCE_FILES(TEST_SOURCES FileIO)
 APPEND_SOURCE_FILES(TEST_SOURCES GeoLib)
 APPEND_SOURCE_FILES(TEST_SOURCES GeoLib/IO)
-APPEND_SOURCE_FILES(TEST_SOURCES InSituLib)
 APPEND_SOURCE_FILES(TEST_SOURCES MathLib)
 APPEND_SOURCE_FILES(TEST_SOURCES MeshLib)
 APPEND_SOURCE_FILES(TEST_SOURCES MeshGeoToolsLib)
@@ -24,7 +22,7 @@ if(QT4_FOUND)
 endif()
 
 if(OGS_USE_PETSC OR OGS_USE_MPI)
-    list(REMOVE_ITEM TEST_SOURCES AssemblerLib/TestSerialLinearSolver.cpp)
+    list(REMOVE_ITEM TEST_SOURCES NumLib/TestSerialLinearSolver.cpp)
 endif()
 
 add_executable(testrunner ${TEST_SOURCES})
@@ -33,10 +31,8 @@ set_target_properties(testrunner PROPERTIES FOLDER Testing)
 target_link_libraries(testrunner
     ApplicationsLib
     ApplicationsFileIO
-    AssemblerLib
     FileIO
     GTest
-    InSituLib
     MeshGeoToolsLib
     MeshLib
     NumLib
diff --git a/Tests/MathLib/TestGlobalMatrixInterface.cpp b/Tests/MathLib/TestGlobalMatrixInterface.cpp
index e211df26547ed764dd5882ff75c54c9a9ca507f0..5a1b1efd200d242ac904bc60388df23d863bf651 100644
--- a/Tests/MathLib/TestGlobalMatrixInterface.cpp
+++ b/Tests/MathLib/TestGlobalMatrixInterface.cpp
@@ -28,7 +28,7 @@
 #include "MathLib/LinAlg/Dense/DenseMatrix.h"
 #include "MathLib/LinAlg/FinalizeMatrixAssembly.h"
 
-#include "ProcessLib/NumericsConfig.h"
+#include "NumLib/NumericsConfig.h"
 
 namespace
 {
diff --git a/Tests/MathLib/TestGlobalVectorInterface.cpp b/Tests/MathLib/TestGlobalVectorInterface.cpp
index 15df0229e9018d5788954894e0f1018fdda62c94..aa7a0c7491bbe7453882ef0182226d1a3ab7b906 100644
--- a/Tests/MathLib/TestGlobalVectorInterface.cpp
+++ b/Tests/MathLib/TestGlobalVectorInterface.cpp
@@ -27,7 +27,7 @@
 #endif
 
 #include "MathLib/LinAlg/FinalizeVectorAssembly.h"
-#include "ProcessLib/NumericsConfig.h"
+#include "NumLib/NumericsConfig.h"
 
 namespace
 {
diff --git a/Tests/InSituLib/TestVtkMappedMeshSource.cpp b/Tests/MeshLib/TestVtkMappedMeshSource.cpp
similarity index 98%
rename from Tests/InSituLib/TestVtkMappedMeshSource.cpp
rename to Tests/MeshLib/TestVtkMappedMeshSource.cpp
index 08b0ad22bfc471a8250e542f448bdab6876e37f2..ea8e520db7f8378158ceb53323873e6568392e6c 100644
--- a/Tests/InSituLib/TestVtkMappedMeshSource.cpp
+++ b/Tests/MeshLib/TestVtkMappedMeshSource.cpp
@@ -18,14 +18,12 @@
 #include <numeric>
 
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
-
-#include "InSituLib/VtkMappedMesh.h"
-#include "InSituLib/VtkMappedMeshSource.h"
-
 #include "MeshLib/Elements/Element.h"
 #include "MeshLib/Mesh.h"
 #include "MeshLib/MeshGenerators/MeshGenerator.h"
 #include "MeshLib/MeshGenerators/VtkMeshConverter.h"
+#include "MeshLib/Vtk/VtkMappedMesh.h"
+#include "MeshLib/Vtk/VtkMappedMeshSource.h"
 
 #include "gtest/gtest.h"
 
@@ -127,7 +125,7 @@ TEST_F(InSituMesh, MappedMesh)
 {
     ASSERT_TRUE(mesh != nullptr);
 
-    vtkNew<InSituLib::VtkMappedMesh> vtkMesh;
+    vtkNew<MeshLib::VtkMappedMesh> vtkMesh;
     vtkMesh->GetImplementation()->SetNodes(mesh->getNodes());
     vtkMesh->GetImplementation()->SetElements(mesh->getElements());
 
@@ -155,7 +153,7 @@ TEST_F(InSituMesh, DISABLED_MappedMeshSourceRoundtrip)
     std::string test_data_file(BaseLib::BuildInfo::tests_tmp_path + "/MappedMeshSourceRoundtrip.vtu");
 
     // -- Test VtkMappedMeshSource, i.e. OGS mesh to VTK mesh
-    vtkNew<InSituLib::VtkMappedMeshSource> vtkSource;
+    vtkNew<MeshLib::VtkMappedMeshSource> vtkSource;
     vtkSource->SetMesh(mesh);
     vtkSource->Update();
     vtkUnstructuredGrid* output = vtkSource->GetOutput();
diff --git a/Tests/InSituLib/TestVtkMappedPropertyVector.cpp b/Tests/MeshLib/TestVtkMappedPropertyVector.cpp
similarity index 90%
rename from Tests/InSituLib/TestVtkMappedPropertyVector.cpp
rename to Tests/MeshLib/TestVtkMappedPropertyVector.cpp
index c7099f7247a8b640d513bd4ce25feca0fce1eaf0..4760ac9365c144bac270e565af3940fbc329bdef 100644
--- a/Tests/InSituLib/TestVtkMappedPropertyVector.cpp
+++ b/Tests/MeshLib/TestVtkMappedPropertyVector.cpp
@@ -18,12 +18,12 @@
 
 #include "gtest/gtest.h"
 
-#include "InSituLib/VtkMappedPropertyVectorTemplate.h"
 #include "MeshLib/Mesh.h"
 #include "MeshLib/MeshGenerators/MeshGenerator.h"
+#include "MeshLib/Vtk/VtkMappedPropertyVectorTemplate.h"
 
 // Creates a PropertyVector<double> and maps it into a vtkDataArray-equivalent
-TEST(InSituLibMappedPropertyVector, Double)
+TEST(MeshLibMappedPropertyVector, Double)
 {
     const std::size_t mesh_size = 5;
     const double length = 1.0;
@@ -40,7 +40,7 @@ TEST(InSituLibMappedPropertyVector, Double)
     (*double_properties).resize(number_of_tuples);
     std::iota((*double_properties).begin(), (*double_properties).end(), 1);
 
-    vtkNew<InSituLib::VtkMappedPropertyVectorTemplate<double> > dataArray;
+    vtkNew<MeshLib::VtkMappedPropertyVectorTemplate<double> > dataArray;
     dataArray->SetPropertyVector(*double_properties);
 
     ASSERT_EQ(dataArray->GetNumberOfComponents(), 1);
@@ -55,7 +55,7 @@ TEST(InSituLibMappedPropertyVector, Double)
 }
 
 // Creates a PropertyVector<int> and maps it into a vtkDataArray-equivalent
-TEST(InSituLibMappedPropertyVector, Int)
+TEST(MeshLibMappedPropertyVector, Int)
 {
     const std::size_t mesh_size = 5;
     const double length = 1.0;
@@ -72,7 +72,7 @@ TEST(InSituLibMappedPropertyVector, Int)
     (*properties).resize(number_of_tuples);
     std::iota((*properties).begin(), (*properties).end(), -5);
 
-    vtkNew<InSituLib::VtkMappedPropertyVectorTemplate<int> > dataArray;
+    vtkNew<MeshLib::VtkMappedPropertyVectorTemplate<int> > dataArray;
     dataArray->SetPropertyVector(*properties);
 
     ASSERT_EQ(dataArray->GetNumberOfComponents(), 1);
@@ -87,7 +87,7 @@ TEST(InSituLibMappedPropertyVector, Int)
 }
 
 // Creates a PropertyVector<unsigned> and maps it into a vtkDataArray-equivalent
-TEST(InSituLibMappedPropertyVector, Unsigned)
+TEST(MeshLibMappedPropertyVector, Unsigned)
 {
     const std::size_t mesh_size = 5;
     const double length = 1.0;
@@ -104,7 +104,7 @@ TEST(InSituLibMappedPropertyVector, Unsigned)
     (*properties).resize(number_of_tuples);
     std::iota((*properties).begin(), (*properties).end(), 0);
 
-    vtkNew<InSituLib::VtkMappedPropertyVectorTemplate<unsigned> > dataArray;
+    vtkNew<MeshLib::VtkMappedPropertyVectorTemplate<unsigned> > dataArray;
     dataArray->SetPropertyVector(*properties);
 
     ASSERT_EQ(dataArray->GetNumberOfComponents(), 1);
diff --git a/Tests/InSituLib/TestVtkMeshNodalCoordinatesTemplate.cpp b/Tests/MeshLib/TestVtkMeshNodalCoordinatesTemplate.cpp
similarity index 91%
rename from Tests/InSituLib/TestVtkMeshNodalCoordinatesTemplate.cpp
rename to Tests/MeshLib/TestVtkMeshNodalCoordinatesTemplate.cpp
index 18e2eae5e0e6bbdf89606923bc6826dfc5391c46..2cd9721fbb2355c91949f4129e1efe2909590af0 100644
--- a/Tests/InSituLib/TestVtkMeshNodalCoordinatesTemplate.cpp
+++ b/Tests/MeshLib/TestVtkMeshNodalCoordinatesTemplate.cpp
@@ -16,11 +16,11 @@
 
 #include "gtest/gtest.h"
 
-#include "InSituLib/VtkMeshNodalCoordinatesTemplate.h"
 #include "MeshLib/Mesh.h"
 #include "MeshLib/MeshGenerators/MeshGenerator.h"
+#include "MeshLib/Vtk/VtkMeshNodalCoordinatesTemplate.h"
 
-TEST(InSituLibNodalCoordinates, Init)
+TEST(MeshLibdalCoordinates, Init)
 {
     const std::size_t subdivisions = 99;
     const double length = 10.0;
@@ -28,7 +28,7 @@ TEST(InSituLibNodalCoordinates, Init)
 
     MeshLib::Mesh* mesh = MeshLib::MeshGenerator::generateRegularQuadMesh(length, subdivisions);
 
-    vtkNew<InSituLib::VtkMeshNodalCoordinatesTemplate<double> > nodeCoords;
+    vtkNew<MeshLib::VtkMeshNodalCoordinatesTemplate<double> > nodeCoords;
     nodeCoords->SetNodes(mesh->getNodes());
     //nodeCoords->PrintSelf(std::cout, vtkIndent());
 
diff --git a/Tests/AssemblerLib/LocalToGlobalIndexMap.cpp b/Tests/NumLib/LocalToGlobalIndexMap.cpp
similarity index 73%
rename from Tests/AssemblerLib/LocalToGlobalIndexMap.cpp
rename to Tests/NumLib/LocalToGlobalIndexMap.cpp
index d9f6e7e2b77913e9891b2aea0a274b7979afad92..2105ab2404888dcf2436e25eab7f45b4fb5a702c 100644
--- a/Tests/AssemblerLib/LocalToGlobalIndexMap.cpp
+++ b/Tests/NumLib/LocalToGlobalIndexMap.cpp
@@ -12,7 +12,7 @@
 
 #include <gtest/gtest.h>
 
-#include "AssemblerLib/LocalToGlobalIndexMap.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
 #include "MeshLib/MeshGenerators/MeshGenerator.h"
 
 #include "MeshLib/MeshSearch/NodeSearch.h"
@@ -20,11 +20,11 @@
 #include "MeshLib/Mesh.h"
 
 
-class AssemblerLibLocalToGlobalIndexMapTest : public ::testing::Test
+class NumLibLocalToGlobalIndexMapTest : public ::testing::Test
 {
 
 public:
-    AssemblerLibLocalToGlobalIndexMapTest()
+    NumLibLocalToGlobalIndexMapTest()
     {
         mesh.reset(MeshLib::MeshGenerator::generateLineMesh(1.0, mesh_size));
         nodesSubset.reset(new MeshLib::MeshSubset(*mesh, &mesh->getNodes()));
@@ -44,21 +44,21 @@ protected:
     static std::size_t const comp1_id = 1;
     std::vector<std::unique_ptr<MeshLib::MeshSubsets>> components;
 
-    std::unique_ptr<AssemblerLib::LocalToGlobalIndexMap const> dof_map;
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap const> dof_map;
 };
 
 #ifndef USE_PETSC
-TEST_F(AssemblerLibLocalToGlobalIndexMapTest, NumberOfRowsByComponent)
+TEST_F(NumLibLocalToGlobalIndexMapTest, NumberOfRowsByComponent)
 #else
-TEST_F(AssemblerLibLocalToGlobalIndexMapTest, DISABLED_NumberOfRowsByComponent)
+TEST_F(NumLibLocalToGlobalIndexMapTest, DISABLED_NumberOfRowsByComponent)
 #endif
 {
     // need to store the size because the components will be moved into the
     // DOF-table.
     std::size_t components_size = components.size();
 
-    dof_map.reset(new AssemblerLib::LocalToGlobalIndexMap(std::move(components),
-        AssemblerLib::ComponentOrder::BY_COMPONENT));
+    dof_map.reset(new NumLib::LocalToGlobalIndexMap(std::move(components),
+        NumLib::ComponentOrder::BY_COMPONENT));
 
     // There must be as many rows as nodes in the input times the number of
     // components.
@@ -66,17 +66,17 @@ TEST_F(AssemblerLibLocalToGlobalIndexMapTest, DISABLED_NumberOfRowsByComponent)
 }
 
 #ifndef USE_PETSC
-TEST_F(AssemblerLibLocalToGlobalIndexMapTest, NumberOfRowsByLocation)
+TEST_F(NumLibLocalToGlobalIndexMapTest, NumberOfRowsByLocation)
 #else
-TEST_F(AssemblerLibLocalToGlobalIndexMapTest, DISABLED_NumberOfRowsByLocation)
+TEST_F(NumLibLocalToGlobalIndexMapTest, DISABLED_NumberOfRowsByLocation)
 #endif
 {
     // need to store the size because the components will be moved into the
     // DOF-table.
     std::size_t components_size = components.size();
 
-    dof_map.reset(new AssemblerLib::LocalToGlobalIndexMap(std::move(components),
-        AssemblerLib::ComponentOrder::BY_LOCATION));
+    dof_map.reset(new NumLib::LocalToGlobalIndexMap(std::move(components),
+        NumLib::ComponentOrder::BY_LOCATION));
 
     // There must be as many rows as nodes in the input times the number of
     // components.
@@ -84,13 +84,13 @@ TEST_F(AssemblerLibLocalToGlobalIndexMapTest, DISABLED_NumberOfRowsByLocation)
 }
 
 #ifndef USE_PETSC
-TEST_F(AssemblerLibLocalToGlobalIndexMapTest, SubsetByComponent)
+TEST_F(NumLibLocalToGlobalIndexMapTest, SubsetByComponent)
 #else
-TEST_F(AssemblerLibLocalToGlobalIndexMapTest, DISABLED_SubsetByComponent)
+TEST_F(NumLibLocalToGlobalIndexMapTest, DISABLED_SubsetByComponent)
 #endif
 {
-    dof_map.reset(new AssemblerLib::LocalToGlobalIndexMap(std::move(components),
-        AssemblerLib::ComponentOrder::BY_COMPONENT));
+    dof_map.reset(new NumLib::LocalToGlobalIndexMap(std::move(components),
+        NumLib::ComponentOrder::BY_COMPONENT));
 
     // Select some elements from the full mesh.
     std::array<std::size_t, 3> const ids = {{ 0, 5, 8 }};
@@ -106,7 +106,7 @@ TEST_F(AssemblerLibLocalToGlobalIndexMapTest, DISABLED_SubsetByComponent)
     auto selected_component = std::unique_ptr<MeshLib::MeshSubsets>{
         new MeshLib::MeshSubsets{selected_subset.get()}};
 
-    auto dof_map_subset = std::unique_ptr<AssemblerLib::LocalToGlobalIndexMap>{
+    auto dof_map_subset = std::unique_ptr<NumLib::LocalToGlobalIndexMap>{
         dof_map->deriveBoundaryConstrainedMap(0,  // variable id
                                               1,  // component id
                                               std::move(selected_component),
diff --git a/Tests/AssemblerLib/LocalToGlobalIndexMapMultiComponent.cpp b/Tests/NumLib/LocalToGlobalIndexMapMultiComponent.cpp
similarity index 91%
rename from Tests/AssemblerLib/LocalToGlobalIndexMapMultiComponent.cpp
rename to Tests/NumLib/LocalToGlobalIndexMapMultiComponent.cpp
index 0184d26fc2efa95fe60bc2818a190b91a1561be0..7be34c60b3b4430e63166f4ff1c1b62141bd8ca4 100644
--- a/Tests/AssemblerLib/LocalToGlobalIndexMapMultiComponent.cpp
+++ b/Tests/NumLib/LocalToGlobalIndexMapMultiComponent.cpp
@@ -10,7 +10,7 @@
 
 #include <memory>
 
-#include "AssemblerLib/VectorMatrixAssembler.h"
+#include "NumLib/Assembler/VectorMatrixAssembler.h"
 #include "MeshLib/Elements/Element.h"
 #include "MeshLib/Location.h"
 #include "MeshLib/Mesh.h"
@@ -23,15 +23,15 @@
 #include "MeshGeoToolsLib/MeshNodeSearcher.h"
 #include "MeshGeoToolsLib/BoundaryElementsSearcher.h"
 
-namespace AL = AssemblerLib;
+namespace AL = NumLib;
 namespace MeL = MeshLib;
 namespace MGTL = MeshGeoToolsLib;
 
-class AssemblerLibLocalToGlobalIndexMapMultiDOFTest : public ::testing::Test
+class NumLibLocalToGlobalIndexMapMultiDOFTest : public ::testing::Test
 {
 public:
     static const std::size_t mesh_subdivs = 4;
-    AssemblerLibLocalToGlobalIndexMapMultiDOFTest()
+    NumLibLocalToGlobalIndexMapMultiDOFTest()
     {
         mesh.reset(MeL::MeshGenerator::generateRegularQuadMesh(2.0, mesh_subdivs));
         mesh_items_all_nodes.reset(new MeL::MeshSubset(*mesh, &mesh->getNodes()));
@@ -71,7 +71,7 @@ public:
             mesh_items_all_nodes->getIntersectionByNodes(nodes));
     }
 
-    ~AssemblerLibLocalToGlobalIndexMapMultiDOFTest()
+    ~NumLibLocalToGlobalIndexMapMultiDOFTest()
     {
         for (auto e : boundary_elements)
             delete e;
@@ -143,7 +143,7 @@ struct ComputeGlobalIndexByLocation
 
 
 template <AL::ComponentOrder ComponentOrder>
-void AssemblerLibLocalToGlobalIndexMapMultiDOFTest::test(
+void NumLibLocalToGlobalIndexMapMultiDOFTest::test(
     const unsigned num_components,
     const unsigned selected_component,
     std::function<std::size_t(std::size_t, std::size_t)> const
@@ -218,9 +218,9 @@ void assert_equal(AL::LocalToGlobalIndexMap const& dof1, AL::LocalToGlobalIndexM
 }
 
 #ifndef USE_PETSC
-TEST_F(AssemblerLibLocalToGlobalIndexMapMultiDOFTest, Test1Comp)
+TEST_F(NumLibLocalToGlobalIndexMapMultiDOFTest, Test1Comp)
 #else
-TEST_F(AssemblerLibLocalToGlobalIndexMapMultiDOFTest, DISABLED_Test1Comp)
+TEST_F(NumLibLocalToGlobalIndexMapMultiDOFTest, DISABLED_Test1Comp)
 #endif
 {
     unsigned const num_components = 1;
@@ -240,9 +240,9 @@ TEST_F(AssemblerLibLocalToGlobalIndexMapMultiDOFTest, DISABLED_Test1Comp)
 }
 
 #ifndef USE_PETSC
-TEST_F(AssemblerLibLocalToGlobalIndexMapMultiDOFTest, TestMultiCompByComponent)
+TEST_F(NumLibLocalToGlobalIndexMapMultiDOFTest, TestMultiCompByComponent)
 #else
-TEST_F(AssemblerLibLocalToGlobalIndexMapMultiDOFTest, DISABLED_TestMultiCompByComponent)
+TEST_F(NumLibLocalToGlobalIndexMapMultiDOFTest, DISABLED_TestMultiCompByComponent)
 #endif
 {
     unsigned const num_components = 5;
@@ -253,9 +253,9 @@ TEST_F(AssemblerLibLocalToGlobalIndexMapMultiDOFTest, DISABLED_TestMultiCompByCo
 }
 
 #ifndef USE_PETSC
-TEST_F(AssemblerLibLocalToGlobalIndexMapMultiDOFTest, TestMultiCompByLocation)
+TEST_F(NumLibLocalToGlobalIndexMapMultiDOFTest, TestMultiCompByLocation)
 #else
-TEST_F(AssemblerLibLocalToGlobalIndexMapMultiDOFTest, DISABLED_TestMultiCompByLocation)
+TEST_F(NumLibLocalToGlobalIndexMapMultiDOFTest, DISABLED_TestMultiCompByLocation)
 #endif
 {
     unsigned const num_components = 5;
diff --git a/Tests/NumLib/ODEs.h b/Tests/NumLib/ODEs.h
index b693bf6efdbefa6aa5f65dae1bdd477fa83d8cb6..bb3c642e6f985e0e8b6485c5acd23bdfbd272a29 100644
--- a/Tests/NumLib/ODEs.h
+++ b/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
diff --git a/Tests/AssemblerLib/SteadyDiffusion2DExample1.h b/Tests/NumLib/SteadyDiffusion2DExample1.h
similarity index 98%
rename from Tests/AssemblerLib/SteadyDiffusion2DExample1.h
rename to Tests/NumLib/SteadyDiffusion2DExample1.h
index e5a4737c35e369a5b99e588cd33ecf5bbe54eafa..5e72dce9b35b50b7e7ca6ed638feb91127c9ba19 100644
--- a/Tests/AssemblerLib/SteadyDiffusion2DExample1.h
+++ b/Tests/NumLib/SteadyDiffusion2DExample1.h
@@ -49,7 +49,7 @@ template<typename IndexType>struct SteadyDiffusion2DExample1
             // already stored in the _localA matrix.
         }
 
-        void addToGlobal(AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
+        void addToGlobal(NumLib::LocalToGlobalIndexMap::RowColumnIndices const& indices,
                          GlobalMatrix& /*M*/, GlobalMatrix& K, GlobalVector& b) const
         {
             K.add(indices, *_localA);
diff --git a/Tests/NumLib/TestExtrapolation.cpp b/Tests/NumLib/TestExtrapolation.cpp
index c26a3c7ffce7265b2fd728fb9a8dd37771812462..c573d4e70c42a836f4e4676fdc253f96944e4bc5 100644
--- a/Tests/NumLib/TestExtrapolation.cpp
+++ b/Tests/NumLib/TestExtrapolation.cpp
@@ -10,12 +10,12 @@
 #include <random>
 #include <gtest/gtest.h>
 
-#include "AssemblerLib/VectorMatrixAssembler.h"
+#include "NumLib/DOF/MatrixProviderUser.h"
+#include "NumLib/DOF/MatrixVectorTraits.h"
+#include "NumLib/DOF/UnifiedMatrixSetters.h"
+#include "NumLib/Assembler/VectorMatrixAssembler.h"
 
 #include "MathLib/LinAlg/BLAS.h"
-#include "MathLib/LinAlg/MatrixProviderUser.h"
-#include "MathLib/LinAlg/MatrixVectorTraits.h"
-#include "MathLib/LinAlg/UnifiedMatrixSetters.h"
 
 #include "MeshLib/MeshGenerators/MeshGenerator.h"
 
@@ -24,7 +24,7 @@
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 
-#include "ProcessLib/NumericsConfig.h"
+#include "NumLib/NumericsConfig.h"
 #include "ProcessLib/Utils/LocalDataInitializer.h"
 #include "ProcessLib/Utils/CreateLocalAssemblers.h"
 #include "ProcessLib/Utils/InitShapeMatrices.h"
@@ -156,7 +156,7 @@ public:
     using GlobalVector = typename GlobalSetup::VectorType;
 
     using LocalAssembler = LocalAssemblerDataInterface<GlobalMatrix, GlobalVector>;
-    using GlobalAssembler = AssemblerLib::VectorMatrixAssembler<
+    using GlobalAssembler = NumLib::VectorMatrixAssembler<
         GlobalMatrix, GlobalVector, LocalAssembler,
         // The exact tag does not matter here.
         NumLib::ODESystemTag::FirstOrderImplicitQuasilinear>;
@@ -175,9 +175,9 @@ public:
         all_mesh_subsets.emplace_back(
                     new MeshLib::MeshSubsets{&_mesh_subset_all_nodes});
 
-        _dof_table.reset(new AssemblerLib::LocalToGlobalIndexMap(
+        _dof_table.reset(new NumLib::LocalToGlobalIndexMap(
               std::move(all_mesh_subsets),
-              AssemblerLib::ComponentOrder::BY_COMPONENT));
+              NumLib::ComponentOrder::BY_COMPONENT));
 
         // Passing _dof_table works, because this process has only one variable
         // and the variable has exactly one component.
@@ -196,7 +196,7 @@ public:
         {
             auto inner_cb = [&loc_asm, property](
                 std::vector<double> const& local_x,
-                AssemblerLib::LocalToGlobalIndexMap::RowColumnIndices const&
+                NumLib::LocalToGlobalIndexMap::RowColumnIndices const&
             ) {
                 loc_asm.interpolateNodalValuesToIntegrationPoints(local_x, property);
             };
@@ -253,7 +253,7 @@ private:
     unsigned const _integration_order;
 
     MeshLib::MeshSubset _mesh_subset_all_nodes;
-    std::unique_ptr<AssemblerLib::LocalToGlobalIndexMap> _dof_table;
+    std::unique_ptr<NumLib::LocalToGlobalIndexMap> _dof_table;
 
     std::unique_ptr<GlobalAssembler> _global_assembler;
     std::vector<std::unique_ptr<LocalAssembler>> _local_assemblers;
diff --git a/Tests/AssemblerLib/TestMeshComponentMap.cpp b/Tests/NumLib/TestMeshComponentMap.cpp
similarity index 84%
rename from Tests/AssemblerLib/TestMeshComponentMap.cpp
rename to Tests/NumLib/TestMeshComponentMap.cpp
index 326593eccea25ec6cff683e5f3b4be7b6ba9b3be..d4ae505f1111e4fe595878d93ee1b6b195847e8e 100644
--- a/Tests/AssemblerLib/TestMeshComponentMap.cpp
+++ b/Tests/NumLib/TestMeshComponentMap.cpp
@@ -14,20 +14,20 @@
 #include <memory>
 #include <vector>
 
-#include "AssemblerLib/MeshComponentMap.h"
+#include "NumLib/DOF/MeshComponentMap.h"
 
 #include "MeshLib/MeshGenerators/MeshGenerator.h"
 #include "MeshLib/MeshSubsets.h"
 
-class AssemblerLibMeshComponentMapTest : public ::testing::Test
+class NumLibMeshComponentMapTest : public ::testing::Test
 {
     public:
     typedef MeshLib::MeshItemType MeshItemType;
     typedef MeshLib::Location Location;
-    typedef AssemblerLib::MeshComponentMap MeshComponentMap;
+    typedef NumLib::MeshComponentMap MeshComponentMap;
 
     public:
-    AssemblerLibMeshComponentMapTest()
+    NumLibMeshComponentMapTest()
         : mesh(nullptr), nodesSubset(nullptr), cmap(nullptr)
     {
         mesh = MeshLib::MeshGenerator::generateLineMesh(1.0, mesh_size);
@@ -38,7 +38,7 @@ class AssemblerLibMeshComponentMapTest : public ::testing::Test
         components.emplace_back(new MeshLib::MeshSubsets{nodesSubset});
     }
 
-    ~AssemblerLibMeshComponentMapTest()
+    ~NumLibMeshComponentMapTest()
     {
         delete cmap;
         delete nodesSubset;
@@ -68,16 +68,16 @@ class AssemblerLibMeshComponentMapTest : public ::testing::Test
 };
 
 #ifndef USE_PETSC
-TEST_F(AssemblerLibMeshComponentMapTest, CheckOrderByComponent)
+TEST_F(NumLibMeshComponentMapTest, CheckOrderByComponent)
 #else
-TEST_F(AssemblerLibMeshComponentMapTest, DISABLED_CheckOrderByComponent)
+TEST_F(NumLibMeshComponentMapTest, DISABLED_CheckOrderByComponent)
 #endif
 {
     // - Entries in the vector are arranged in the order of a component type and then node ID
     // - For example, x=[(node 0, comp 0) (node 1, comp 0) ... (node n, comp0), (node 0, comp1) ... ]
 
     cmap = new MeshComponentMap(components,
-        AssemblerLib::ComponentOrder::BY_COMPONENT);
+        NumLib::ComponentOrder::BY_COMPONENT);
 
     ASSERT_EQ(2 * mesh->getNNodes(), cmap->dofSizeWithGhosts());
     for (std::size_t i = 0; i < mesh_size; i++)
@@ -96,16 +96,16 @@ TEST_F(AssemblerLibMeshComponentMapTest, DISABLED_CheckOrderByComponent)
 }
 
 #ifndef USE_PETSC
-TEST_F(AssemblerLibMeshComponentMapTest, CheckOrderByLocation)
+TEST_F(NumLibMeshComponentMapTest, CheckOrderByLocation)
 #else
-TEST_F(AssemblerLibMeshComponentMapTest, DISABLED_CheckOrderByLocation)
+TEST_F(NumLibMeshComponentMapTest, DISABLED_CheckOrderByLocation)
 #endif
 {
     // - Entries in the vector are arranged in the order of node ID and then a component type
     // - For example, x=[(node 0, comp 0) (node 0, comp 1) ... (node n, comp0), (node n, comp1) ]
 
     cmap = new MeshComponentMap(components,
-        AssemblerLib::ComponentOrder::BY_LOCATION);
+        NumLib::ComponentOrder::BY_LOCATION);
 
     ASSERT_EQ(2 * mesh->getNNodes(), cmap->dofSizeWithGhosts());
     for (std::size_t i = 0; i < mesh_size; i++)
@@ -124,13 +124,13 @@ TEST_F(AssemblerLibMeshComponentMapTest, DISABLED_CheckOrderByLocation)
 }
 
 #ifndef USE_PETSC
-TEST_F(AssemblerLibMeshComponentMapTest, OutOfRangeAccess)
+TEST_F(NumLibMeshComponentMapTest, OutOfRangeAccess)
 #else
-TEST_F(AssemblerLibMeshComponentMapTest, DISABLED_OutOfRangeAccess)
+TEST_F(NumLibMeshComponentMapTest, DISABLED_OutOfRangeAccess)
 #endif
 {
     cmap = new MeshComponentMap(components,
-        AssemblerLib::ComponentOrder::BY_COMPONENT);
+        NumLib::ComponentOrder::BY_COMPONENT);
 
     ASSERT_EQ(MeshComponentMap::nop, cmap->getGlobalIndex(
         Location(mesh->getID(), MeshItemType::Node, mesh_size + 1), comp0_id));
@@ -143,13 +143,13 @@ TEST_F(AssemblerLibMeshComponentMapTest, DISABLED_OutOfRangeAccess)
 }
 
 #ifndef USE_PETSC
-TEST_F(AssemblerLibMeshComponentMapTest, SubsetOfNodesByComponent)
+TEST_F(NumLibMeshComponentMapTest, SubsetOfNodesByComponent)
 #else
-TEST_F(AssemblerLibMeshComponentMapTest, DISABLED_SubsetOfNodesByComponent)
+TEST_F(NumLibMeshComponentMapTest, DISABLED_SubsetOfNodesByComponent)
 #endif
 {
     cmap = new MeshComponentMap(components,
-        AssemblerLib::ComponentOrder::BY_COMPONENT);
+        NumLib::ComponentOrder::BY_COMPONENT);
 
     // Select some nodes from the full mesh.
     std::array<std::size_t, 3> const ids = {{ 0, 5, 9 }};
@@ -180,13 +180,13 @@ TEST_F(AssemblerLibMeshComponentMapTest, DISABLED_SubsetOfNodesByComponent)
 }
 
 #ifndef USE_PETSC
-TEST_F(AssemblerLibMeshComponentMapTest, SubsetOfNodesByLocation)
+TEST_F(NumLibMeshComponentMapTest, SubsetOfNodesByLocation)
 #else
-TEST_F(AssemblerLibMeshComponentMapTest, DISABLED_SubsetOfNodesByLocation)
+TEST_F(NumLibMeshComponentMapTest, DISABLED_SubsetOfNodesByLocation)
 #endif
 {
     cmap = new MeshComponentMap(components,
-        AssemblerLib::ComponentOrder::BY_LOCATION);
+        NumLib::ComponentOrder::BY_LOCATION);
 
     // Select some nodes from the full mesh.
     std::array<std::size_t, 3> const ids = {{ 0, 5, 9 }};
diff --git a/Tests/NumLib/TestODEInt.cpp b/Tests/NumLib/TestODEInt.cpp
index 5541c4e201ff6e1037130ed689a5af8466b1c32d..eedde108ee7b8d2188b366cb77471e2b2d5da868 100644
--- a/Tests/NumLib/TestODEInt.cpp
+++ b/Tests/NumLib/TestODEInt.cpp
@@ -7,7 +7,7 @@
 
 #include "BaseLib/BuildInfo.h"
 #include "NumLib/ODESolver/TimeLoopSingleODE.h"
-#include "ProcessLib/NumericsConfig.h"
+#include "NumLib/NumericsConfig.h"
 #include "ODEs.h"
 
 
diff --git a/Tests/AssemblerLib/TestSerialExecutor.cpp b/Tests/NumLib/TestSerialExecutor.cpp
similarity index 88%
rename from Tests/AssemblerLib/TestSerialExecutor.cpp
rename to Tests/NumLib/TestSerialExecutor.cpp
index 0f47f417ce9cfa52144550ba6ff6d1fff5ac2e8c..0ae72d9ee14658377bda0ab4c95df029a78f4259 100644
--- a/Tests/AssemblerLib/TestSerialExecutor.cpp
+++ b/Tests/NumLib/TestSerialExecutor.cpp
@@ -13,10 +13,10 @@
 #include <vector>
 #include <functional>
 #include <numeric>
-#include "AssemblerLib/SerialExecutor.h"
+#include "NumLib/Assembler/SerialExecutor.h"
 
 template <typename ContainerElement_>
-class AssemblerLibSerialExecutor : public ::testing::Test
+class NumLibSerialExecutor : public ::testing::Test
 {
 public:
     using ContainerElement = ContainerElement_;
@@ -69,9 +69,9 @@ public:
 
 typedef ::testing::Types<int> TestCases;
 
-TYPED_TEST_CASE(AssemblerLibSerialExecutor, TestCases);
+TYPED_TEST_CASE(NumLibSerialExecutor, TestCases);
 
-TYPED_TEST(AssemblerLibSerialExecutor, ContainerArgument)
+TYPED_TEST(NumLibSerialExecutor, ContainerArgument)
 {
     using Elem         = typename TestFixture::ContainerElement;
     using Container    = typename TestFixture::Container;
@@ -84,14 +84,14 @@ TYPED_TEST(AssemblerLibSerialExecutor, ContainerArgument)
                     TestFixture::subtractFromReferenceStatic(value, index, ref_inner);
                 };
 
-            AssemblerLib::SerialExecutor::executeDereferenced(
+            NumLib::SerialExecutor::executeDereferenced(
                 cb_static, ctnr, ref);
         }
     );
 
     TestFixture::test(
         [this](PtrContainer const& ctnr, Container& ref) {
-            AssemblerLib::SerialExecutor::executeMemberDereferenced(
+            NumLib::SerialExecutor::executeMemberDereferenced(
                 *static_cast<TestFixture*>(this), &TestFixture::subtractFromReference,
                 ctnr, ref
             );
diff --git a/Tests/AssemblerLib/TestSerialLinearSolver.cpp b/Tests/NumLib/TestSerialLinearSolver.cpp
similarity index 94%
rename from Tests/AssemblerLib/TestSerialLinearSolver.cpp
rename to Tests/NumLib/TestSerialLinearSolver.cpp
index 402244719e7c896ad2e9e6c7789756d11f1820ef..88e1cbd2e4f2279325a3fe6e56bb5db867583b45 100644
--- a/Tests/AssemblerLib/TestSerialLinearSolver.cpp
+++ b/Tests/NumLib/TestSerialLinearSolver.cpp
@@ -16,7 +16,7 @@
 
 #include <gtest/gtest.h>
 
-#include "AssemblerLib/VectorMatrixAssembler.h"
+#include "NumLib/Assembler/VectorMatrixAssembler.h"
 
 #include "MathLib/LinAlg/ApplyKnownSolution.h"
 #include "MathLib/LinAlg/Solvers/GaussAlgorithm.h"
@@ -30,12 +30,12 @@
 #include "MeshLib/Mesh.h"
 #include "MeshLib/Node.h"
 
-#include "ProcessLib/NumericsConfig.h"
+#include "NumLib/NumericsConfig.h"
 
 #include "../TestTools.h"
 #include "SteadyDiffusion2DExample1.h"
 
-TEST(AssemblerLibSerialLinearSolver, Steady2DdiffusionQuadElem)
+TEST(NumLibSerialLinearSolver, Steady2DdiffusionQuadElem)
 {
     // example
     using Example = SteadyDiffusion2DExample1<GlobalIndexType>;
@@ -58,8 +58,8 @@ TEST(AssemblerLibSerialLinearSolver, Steady2DdiffusionQuadElem)
     // define a mesh item composition in a vector
     std::vector<std::unique_ptr<MeshLib::MeshSubsets>> vec_comp_dis;
     vec_comp_dis.emplace_back(new MeshLib::MeshSubsets{&mesh_items_all_nodes});
-    AssemblerLib::LocalToGlobalIndexMap local_to_global_index_map(
-        std::move(vec_comp_dis), AssemblerLib::ComponentOrder::BY_COMPONENT);
+    NumLib::LocalToGlobalIndexMap local_to_global_index_map(
+        std::move(vec_comp_dis), NumLib::ComponentOrder::BY_COMPONENT);
 
     //--------------------------------------------------------------------------
     // Construct a linear system
@@ -102,7 +102,7 @@ TEST(AssemblerLibSerialLinearSolver, Steady2DdiffusionQuadElem)
 
     // TODO in the future use simpler NumLib::ODESystemTag
     // Local and global assemblers.
-    typedef AssemblerLib::VectorMatrixAssembler<
+    typedef NumLib::VectorMatrixAssembler<
             GlobalMatrix, GlobalVector, LocalAssembler,
             NumLib::ODESystemTag::FirstOrderImplicitQuasilinear> GlobalAssembler;
 
diff --git a/Tests/AssemblerLib/TestVectorMatrixBuilder.cpp b/Tests/NumLib/TestVectorMatrixBuilder.cpp
similarity index 74%
rename from Tests/AssemblerLib/TestVectorMatrixBuilder.cpp
rename to Tests/NumLib/TestVectorMatrixBuilder.cpp
index a47171805fbfc2cc1cf83e3e1c2891717c83c513..dbc1d7811df8ca65e90d58c774836204e245c39f 100644
--- a/Tests/AssemblerLib/TestVectorMatrixBuilder.cpp
+++ b/Tests/NumLib/TestVectorMatrixBuilder.cpp
@@ -10,26 +10,26 @@
 #include <memory>
 #include <gtest/gtest.h>
 
-#include "AssemblerLib/MeshComponentMap.h"
+#include "NumLib/DOF/MeshComponentMap.h"
 
 #include "MeshLib/MeshGenerators/MeshGenerator.h"
 #include "MeshLib/MeshSubsets.h"
 
-#include "AssemblerLib/VectorMatrixBuilder.h"
+#include "NumLib/Assembler/VectorMatrixBuilder.h"
 
 template <typename Builder>
-class AssemblerLibVectorMatrixBuilder : public ::testing::Test
+class NumLibVectorMatrixBuilder : public ::testing::Test
 {
     public:
     typedef MeshLib::MeshItemType MeshItemType;
     typedef MeshLib::Location Location;
-    typedef AssemblerLib::MeshComponentMap MeshComponentMap;
+    typedef NumLib::MeshComponentMap MeshComponentMap;
 
     typedef typename Builder::VectorType VectorType;
     typedef typename Builder::MatrixType MatrixType;
 
     public:
-    AssemblerLibVectorMatrixBuilder()
+    NumLibVectorMatrixBuilder()
         : mesh(nullptr), nodesSubset(nullptr), cmap(nullptr)
     {
         mesh = MeshLib::MeshGenerator::generateLineMesh(1.0, mesh_size);
@@ -40,10 +40,10 @@ class AssemblerLibVectorMatrixBuilder : public ::testing::Test
         components.emplace_back(new MeshLib::MeshSubsets{nodesSubset});
 
         cmap = new MeshComponentMap(components,
-            AssemblerLib::ComponentOrder::BY_COMPONENT);
+            NumLib::ComponentOrder::BY_COMPONENT);
     }
 
-    ~AssemblerLibVectorMatrixBuilder()
+    ~NumLibVectorMatrixBuilder()
     {
         delete cmap;
         delete nodesSubset;
@@ -58,12 +58,12 @@ class AssemblerLibVectorMatrixBuilder : public ::testing::Test
     MeshComponentMap const* cmap;
 };
 
-TYPED_TEST_CASE_P(AssemblerLibVectorMatrixBuilder);
+TYPED_TEST_CASE_P(NumLibVectorMatrixBuilder);
 
 #ifndef USE_PETSC
-TYPED_TEST_P(AssemblerLibVectorMatrixBuilder, createVector)
+TYPED_TEST_P(NumLibVectorMatrixBuilder, createVector)
 #else
-TYPED_TEST_P(AssemblerLibVectorMatrixBuilder, DISABLED_createVector)
+TYPED_TEST_P(NumLibVectorMatrixBuilder, DISABLED_createVector)
 #endif
 {
     typedef typename TestFixture::VectorType V;
@@ -77,9 +77,9 @@ TYPED_TEST_P(AssemblerLibVectorMatrixBuilder, DISABLED_createVector)
 }
 
 #ifndef USE_PETSC
-TYPED_TEST_P(AssemblerLibVectorMatrixBuilder, createMatrix)
+TYPED_TEST_P(NumLibVectorMatrixBuilder, createMatrix)
 #else
-TYPED_TEST_P(AssemblerLibVectorMatrixBuilder, DISABLED_createMatrix)
+TYPED_TEST_P(NumLibVectorMatrixBuilder, DISABLED_createMatrix)
 #endif
 {
     typedef typename TestFixture::MatrixType M;
@@ -94,10 +94,10 @@ TYPED_TEST_P(AssemblerLibVectorMatrixBuilder, DISABLED_createMatrix)
 }
 
 #ifndef USE_PETSC
-REGISTER_TYPED_TEST_CASE_P(AssemblerLibVectorMatrixBuilder,
+REGISTER_TYPED_TEST_CASE_P(NumLibVectorMatrixBuilder,
     createVector, createMatrix);
 #else
-REGISTER_TYPED_TEST_CASE_P(AssemblerLibVectorMatrixBuilder,
+REGISTER_TYPED_TEST_CASE_P(NumLibVectorMatrixBuilder,
     DISABLED_createVector, DISABLED_createMatrix);
 #endif
 
@@ -118,17 +118,17 @@ REGISTER_TYPED_TEST_CASE_P(AssemblerLibVectorMatrixBuilder,
 
 typedef ::testing::Types
     <
-      AssemblerLib::VectorMatrixBuilder<
+      NumLib::VectorMatrixBuilder<
       MathLib::EigenMatrix, MathLib::EigenVector>
 #ifdef USE_LIS
-    , AssemblerLib::VectorMatrixBuilder<
+    , NumLib::VectorMatrixBuilder<
         MathLib::LisMatrix, MathLib::LisVector>
 #endif  // USE_LIS
 #ifdef USE_PETSC
-    , AssemblerLib::VectorMatrixBuilder<
+    , NumLib::VectorMatrixBuilder<
         MathLib::PETScMatrix, MathLib::PETScVector>
 #endif  // USE_PETSC
     > TestTypes;
 
-INSTANTIATE_TYPED_TEST_CASE_P(templated, AssemblerLibVectorMatrixBuilder,
+INSTANTIATE_TYPED_TEST_CASE_P(templated, NumLibVectorMatrixBuilder,
     TestTypes);
diff --git a/Tests/testrunner.cpp b/Tests/testrunner.cpp
index 0e9991d41d76b5cb4c74b88a19376a190f82374c..b14fe810ef19a7ea905c8860e3d731d66ecbc8d1 100644
--- a/Tests/testrunner.cpp
+++ b/Tests/testrunner.cpp
@@ -19,9 +19,9 @@
 
 #include "Applications/ApplicationsLib/LogogSetup.h"
 #include "Applications/ApplicationsLib/LinearSolverLibrarySetup.h"
+#include "NumLib/DOF/GlobalMatrixProviders.h"
 #include "BaseLib/BuildInfo.h"
 #include "BaseLib/TemplateLogogFormatterSuppressedGCC.h"
-#include "MathLib/LinAlg/GlobalMatrixProviders.h"
 
 #ifdef OGS_BUILD_GUI
 #include <QApplication>