diff --git a/ProcessLib/Utils/InitShapeMatrices.h b/NumLib/Fem/InitShapeMatrices.h
similarity index 61%
rename from ProcessLib/Utils/InitShapeMatrices.h
rename to NumLib/Fem/InitShapeMatrices.h
index fe31adc7afa0193964e155058fa1494aabf59aca..ab5f8feba2104ac73ce95a0ed223326554073e69 100644
--- a/ProcessLib/Utils/InitShapeMatrices.h
+++ b/NumLib/Fem/InitShapeMatrices.h
@@ -11,14 +11,14 @@
 
 #include <vector>
 
+#include "FiniteElement/TemplateIsoparametric.h"
 #include "MeshLib/Elements/Element.h"
-#include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
 
-
-namespace ProcessLib
+namespace NumLib
 {
-template <typename ShapeFunction, typename ShapeMatricesType,
-          typename IntegrationMethod, unsigned GlobalDim>
+template <typename ShapeFunction, typename ShapeMatricesType, int GlobalDim,
+          ShapeMatrixType SelectedShapeMatrixType = ShapeMatrixType::ALL,
+          typename IntegrationMethod>
 std::vector<typename ShapeMatricesType::ShapeMatrices,
             Eigen::aligned_allocator<typename ShapeMatricesType::ShapeMatrices>>
 initShapeMatrices(MeshLib::Element const& e, bool is_axially_symmetric,
@@ -30,18 +30,19 @@ initShapeMatrices(MeshLib::Element const& e, bool is_axially_symmetric,
         shape_matrices;
 
     auto const fe =
-        NumLib::createIsoparametricFiniteElement<ShapeFunction,
-                                                 ShapeMatricesType>(e);
+        createIsoparametricFiniteElement<ShapeFunction, ShapeMatricesType>(e);
 
-    unsigned const n_integration_points = integration_method.getNumberOfPoints();
+    unsigned const n_integration_points =
+        integration_method.getNumberOfPoints();
 
     shape_matrices.reserve(n_integration_points);
-    for (unsigned ip = 0; ip < n_integration_points; ++ip) {
+    for (unsigned ip = 0; ip < n_integration_points; ++ip)
+    {
         shape_matrices.emplace_back(ShapeFunction::DIM, GlobalDim,
-                                     ShapeFunction::NPOINTS);
+                                    ShapeFunction::NPOINTS);
         fe.computeShapeFunctions(
-                integration_method.getWeightedPoint(ip).getCoords(),
-                shape_matrices[ip], GlobalDim, is_axially_symmetric);
+            integration_method.getWeightedPoint(ip).getCoords(),
+            shape_matrices[ip], GlobalDim, is_axially_symmetric);
     }
 
     return shape_matrices;
@@ -53,8 +54,7 @@ double interpolateXCoordinate(
     typename ShapeMatricesType::ShapeMatrices::ShapeType const& N)
 {
     auto const fe =
-        NumLib::createIsoparametricFiniteElement<ShapeFunction,
-                                                 ShapeMatricesType>(e);
+        createIsoparametricFiniteElement<ShapeFunction, ShapeMatricesType>(e);
 
     return fe.interpolateZerothCoordinate(N);
 }
@@ -65,10 +65,9 @@ std::array<double, 3> interpolateCoordinates(
     typename ShapeMatricesType::ShapeMatrices::ShapeType const& N)
 {
     auto const fe =
-        NumLib::createIsoparametricFiniteElement<ShapeFunction,
-                                                 ShapeMatricesType>(e);
+        createIsoparametricFiniteElement<ShapeFunction, ShapeMatricesType>(e);
 
     return fe.interpolateCoordinates(N);
 }
 
-}  // namespace ProcessLib
+}  // namespace NumLib
diff --git a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
index fe2750f51085bf99ef0bf8c1d5ae21c8eb57884f..e456bab90a0344c92f21c89f538951dec354c162 100644
--- a/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/ConstraintDirichletBoundaryConditionLocalAssembler.h
@@ -10,17 +10,14 @@
 
 #pragma once
 
-#include "NumLib/Fem/ShapeMatrixPolicy.h"
-
+#include "MeshLib/Elements/Elements.h"
+#include "MeshLib/Elements/MapBulkElementPoint.h"
+#include "MeshLib/Elements/Utils.h"
 #include "NumLib/DOF/DOFTableUtil.h"
-
+#include "NumLib/Fem/InitShapeMatrices.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/Process.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
-#include "MeshLib/Elements/MapBulkElementPoint.h"
-#include "MeshLib/Elements/Elements.h"
-#include "MeshLib/Elements/Utils.h"
 
 namespace ProcessLib
 {
@@ -86,27 +83,19 @@ public:
           _surface_element_normal(MeshLib::calculateNormalizedSurfaceNormal(
               _surface_element, *(bulk_mesh.getElements()[_bulk_element_id])))
     {
-        auto const fe = NumLib::createIsoparametricFiniteElement<
-            ShapeFunction, ShapeMatricesType>(_surface_element);
+        auto const shape_matrices =
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      GlobalDim, NumLib::ShapeMatrixType::N_J>(
+                _surface_element, is_axially_symmetric, _integration_method);
+
+        auto const bulk_face_id = bulk_ids[_surface_element.getID()].second;
 
         auto const n_integration_points =
             _integration_method.getNumberOfPoints();
-
-        auto const bulk_face_id = bulk_ids[_surface_element.getID()].second;
-        std::vector<
-            typename ShapeMatricesType::ShapeMatrices,
-            Eigen::aligned_allocator<typename ShapeMatricesType::ShapeMatrices>>
-            shape_matrices;
-        shape_matrices.reserve(n_integration_points);
         _ip_data.reserve(n_integration_points);
+
         for (unsigned ip = 0; ip < n_integration_points; ++ip)
         {
-            shape_matrices.emplace_back(ShapeFunction::DIM, GlobalDim,
-                                        ShapeFunction::NPOINTS);
-            fe.template computeShapeFunctions<NumLib::ShapeMatrixType::N_J>(
-                _integration_method.getWeightedPoint(ip).getCoords(),
-                shape_matrices[ip], GlobalDim, is_axially_symmetric);
-
             auto const& wp = _integration_method.getWeightedPoint(ip);
             auto bulk_element_point = MeshLib::getBulkElementPoint(
                 bulk_mesh, _bulk_element_id, bulk_face_id, wp);
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h
index 40a29e985c21b8b1da067d225ce0de4f89862b33..328fe5cb3ec41c2551f6180e8cc10d1dbf0af619 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h
@@ -12,8 +12,8 @@
 
 #include "MeshLib/Elements/Element.h"
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
 
 namespace ProcessLib
 {
@@ -60,9 +60,9 @@ private:
             ns_and_weights;
         ns_and_weights.reserve(integration_method.getNumberOfPoints());
 
-        auto sms = initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                                     IntegrationMethod, GlobalDim>(
-            e, is_axially_symmetric, integration_method);
+        auto sms = NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                             GlobalDim>(e, is_axially_symmetric,
+                                                        integration_method);
         for (unsigned ip = 0; ip < sms.size(); ++ip)
         {
             auto& sm = sms[ip];
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
index 6b2db58606c8ba1575e4ae3e5e6c8ae4e2a256a7..9b47168948b7b56d768d2de2b5bcc03d72c9ed97 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
@@ -12,9 +12,9 @@
 
 #include "GenericNaturalBoundaryConditionLocalAssembler.h"
 #include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/Parameter.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
 
 namespace ProcessLib
 {
@@ -78,7 +78,8 @@ public:
             ParameterLib::SpatialPosition const position{
                 boost::none, Base::_element.getID(), ip,
                 MathLib::Point3d(
-                    interpolateCoordinates<ShapeFunction, ShapeMatricesType>(
+                    NumLib::interpolateCoordinates<ShapeFunction,
+                                                   ShapeMatricesType>(
                         Base::_element, N))};
 
             if (_data.integral_measure)
diff --git a/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h
index 8ed02b06f0136add80ca7754fc72319a1efd691f..d32103068b7f7a8bf82dc700b64cee0f52bfb104 100644
--- a/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/NormalTractionBoundaryConditionLocalAssembler.h
@@ -77,8 +77,8 @@ public:
         _ip_data.reserve(n_integration_points);
 
         auto const shape_matrices_u =
-            initShapeMatrices<ShapeFunctionDisplacement, ShapeMatricesType,
-                              IntegrationMethod, GlobalDim>(
+            NumLib::initShapeMatrices<ShapeFunctionDisplacement,
+                                      ShapeMatricesType, GlobalDim>(
                 e, is_axially_symmetric, _integration_method);
 
         GlobalDimVectorType element_normal(GlobalDim);
diff --git a/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
index d7b6c7af557d68cc2384c7be9aa42569b3be0a8a..380fb039dc9f36ea658c085b0610659e6f0cd4fe 100644
--- a/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
@@ -75,7 +75,8 @@ public:
             ParameterLib::SpatialPosition const position{
                 boost::none, Base::_element.getID(), ip,
                 MathLib::Point3d(
-                    interpolateCoordinates<ShapeFunction, ShapeMatricesType>(
+                    NumLib::interpolateCoordinates<ShapeFunction,
+                                                   ShapeMatricesType>(
                         Base::_element, N))};
 
             double integral_measure = 1.0;
diff --git a/ProcessLib/ComponentTransport/ComponentTransportFEM.h b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
index 2425bfefe4182adc1dfd7513ceed32606da5f894..2c3b5f5d5162e0dc9f9f3f8f4489271502dcdf8b 100644
--- a/ProcessLib/ComponentTransport/ComponentTransportFEM.h
+++ b/ProcessLib/ComponentTransport/ComponentTransportFEM.h
@@ -23,13 +23,13 @@
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "NumLib/Function/Interpolation.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
 #include "ProcessLib/ProcessVariable.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
 
 namespace ProcessLib
 {
@@ -146,9 +146,9 @@ public:
         _ip_data.reserve(n_integration_points);
 
         auto const shape_matrices =
-            initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                              IntegrationMethod, GlobalDim>(
-                element, is_axially_symmetric, _integration_method);
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      GlobalDim>(element, is_axially_symmetric,
+                                                 _integration_method);
 
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
diff --git a/ProcessLib/Deformation/MaterialForces.h b/ProcessLib/Deformation/MaterialForces.h
index a6fc394803d0f8114652223bb1123bd7d6c6a9ff..74f7c7e26c6be37e95c939bcb6d8799dce565d92 100644
--- a/ProcessLib/Deformation/MaterialForces.h
+++ b/ProcessLib/Deformation/MaterialForces.h
@@ -17,8 +17,8 @@
 #include "MathLib/LinAlg/LinAlg.h"
 #include "MathLib/LinAlg/MatrixVectorTraits.h"
 #include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "ProcessLib/Deformation/GMatrix.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
 
 namespace ProcessLib
 {
@@ -59,7 +59,7 @@ std::vector<double> const& getMaterialForces(
         auto const& psi = _ip_data[ip].free_energy_density;
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+            NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
                 element, _ip_data[ip].N);
 
         // For the 2D case the 33-component is needed (and the four entries
diff --git a/ProcessLib/HT/HTFEM.h b/ProcessLib/HT/HTFEM.h
index a48d136b84f29f5c1be19e6e6183717da19a3c5a..701d448e400a11492301f1a26807013def616d8b 100644
--- a/ProcessLib/HT/HTFEM.h
+++ b/ProcessLib/HT/HTFEM.h
@@ -13,19 +13,16 @@
 #include <Eigen/Dense>
 #include <vector>
 
+#include "HTLocalAssemblerInterface.h"
 #include "HTProcessData.h"
-
 #include "MaterialLib/MPL/Medium.h"
 #include "MaterialLib/MPL/Utils/FormEffectiveThermalConductivity.h"
-
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/Parameter.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
-#include "HTLocalAssemblerInterface.h"
 
 namespace ProcessLib
 {
@@ -69,9 +66,9 @@ public:
         _ip_data.reserve(n_integration_points);
 
         auto const shape_matrices =
-            initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                              IntegrationMethod, GlobalDim>(
-                element, is_axially_symmetric, _integration_method);
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      GlobalDim>(element, is_axially_symmetric,
+                                                 _integration_method);
 
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
diff --git a/ProcessLib/HT/MonolithicHTFEM.h b/ProcessLib/HT/MonolithicHTFEM.h
index 516bb58972f79aa5538c2f3aa75d1288f6d2c7c1..b36e3476950da60c14db43e67a26bb5d5f74b4d1 100644
--- a/ProcessLib/HT/MonolithicHTFEM.h
+++ b/ProcessLib/HT/MonolithicHTFEM.h
@@ -14,18 +14,16 @@
 #include <Eigen/Dense>
 #include <vector>
 
+#include "HTFEM.h"
 #include "HTProcessData.h"
 #include "MaterialLib/MPL/Medium.h"
 #include "MaterialLib/MPL/Utils/FormEigenTensor.h"
-
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/Parameter.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
-#include "HTFEM.h"
 
 namespace ProcessLib
 {
diff --git a/ProcessLib/HT/StaggeredHTFEM.h b/ProcessLib/HT/StaggeredHTFEM.h
index 5dca497bfd94808e6cdf018f1e2c8601e902710d..5230997fb59a2928aa295336ab349b10a90c7f5e 100644
--- a/ProcessLib/HT/StaggeredHTFEM.h
+++ b/ProcessLib/HT/StaggeredHTFEM.h
@@ -14,14 +14,13 @@
 #include <Eigen/Dense>
 #include <vector>
 
+#include "HTFEM.h"
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/Parameter.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
-#include "HTFEM.h"
 
 namespace ProcessLib
 {
diff --git a/ProcessLib/HeatConduction/HeatConductionFEM.h b/ProcessLib/HeatConduction/HeatConductionFEM.h
index dc676f5214950a3f872ed6c26d05066a9790d7dc..72740f7c502cfb3b2758efefea8b072c22a194e4 100644
--- a/ProcessLib/HeatConduction/HeatConductionFEM.h
+++ b/ProcessLib/HeatConduction/HeatConductionFEM.h
@@ -19,10 +19,10 @@
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
 
 namespace ProcessLib
 {
@@ -79,9 +79,10 @@ public:
         : _element(element),
           _process_data(process_data),
           _integration_method(integration_order),
-          _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                                            IntegrationMethod, GlobalDim>(
-              element, is_axially_symmetric, _integration_method)),
+          _shape_matrices(
+              NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                        GlobalDim>(
+                  element, is_axially_symmetric, _integration_method)),
           _heat_fluxes(
               GlobalDim,
               std::vector<double>(_integration_method.getNumberOfPoints()))
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
index 548e2efa719983d27f364c424de1fcb7c5a41f09..03dbc0e02f5362aea46c620310ed9c510fb382c0 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerBHE-impl.h
@@ -12,7 +12,7 @@
 
 #include "HeatTransportBHELocalAssemblerBHE.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 
 namespace ProcessLib
 {
@@ -40,9 +40,9 @@ HeatTransportBHELocalAssemblerBHE<ShapeFunction, IntegrationMethod, BHEType>::
     _secondary_data.N.resize(n_integration_points);
 
     auto const shape_matrices =
-        initShapeMatrices<ShapeFunction, ShapeMatricesType, IntegrationMethod,
-                          3 /* GlobalDim */>(e, is_axially_symmetric,
-                                             _integration_method);
+        NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                  3 /* GlobalDim */>(e, is_axially_symmetric,
+                                                     _integration_method);
 
     // ip data initialization
     for (unsigned ip = 0; ip < n_integration_points; ip++)
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil-impl.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil-impl.h
index 2a796f104d9421fe8590f91c27b3a179e8a8a43d..e6c3bee521a09f265ad1e7914976809002dc917b 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil-impl.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil-impl.h
@@ -13,15 +13,14 @@
 #include <valarray>
 #include <vector>
 
+#include "HeatTransportBHEProcessAssemblerInterface.h"
 #include "MaterialLib/MPL/Medium.h"
 #include "MaterialLib/MPL/Utils/FormEigenTensor.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "NumLib/Function/Interpolation.h"
 #include "ProcessLib/HeatTransportBHE/HeatTransportBHEProcessData.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
-#include "HeatTransportBHEProcessAssemblerInterface.h"
 #include "SecondaryData.h"
 
 namespace ProcessLib
@@ -45,9 +44,10 @@ HeatTransportBHELocalAssemblerSoil<ShapeFunction, IntegrationMethod>::
     _ip_data.reserve(n_integration_points);
     _secondary_data.N.resize(n_integration_points);
 
-    _shape_matrices = initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                                        IntegrationMethod, 3 /* GlobalDim */>(
-        e, is_axially_symmetric, _integration_method);
+    _shape_matrices =
+        NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                  3 /* GlobalDim */>(e, is_axially_symmetric,
+                                                     _integration_method);
 
     ParameterLib::SpatialPosition x_position;
     x_position.setElementID(_element_id);
diff --git a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h
index 536d1148555402705252a8f94c0cfebae586e9a8..1e1aa0a469728d31896d228f34a5a5334e6e080e 100644
--- a/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h
+++ b/ProcessLib/HeatTransportBHE/LocalAssemblers/HeatTransportBHELocalAssemblerSoil.h
@@ -12,15 +12,13 @@
 
 #include <vector>
 
+#include "HeatTransportBHEProcessAssemblerInterface.h"
+#include "IntegrationPointDataSoil.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ProcessLib/Deformation/BMatrixPolicy.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
 #include "ProcessLib/HeatTransportBHE/HeatTransportBHEProcessData.h"
-
-#include "HeatTransportBHEProcessAssemblerInterface.h"
-#include "IntegrationPointDataSoil.h"
 #include "SecondaryData.h"
 
 namespace ProcessLib
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index 90401c561904bdd83dbac9ab064492b5e890e602..d4668b975f6d2ad422f05ff64f8fe99f1b2feb06 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -50,14 +50,14 @@ HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
     _secondary_data.N_u.resize(n_integration_points);
 
     auto const shape_matrices_u =
-        initShapeMatrices<ShapeFunctionDisplacement,
-                          ShapeMatricesTypeDisplacement, IntegrationMethod,
-                          DisplacementDim>(e, is_axially_symmetric,
-                                           _integration_method);
+        NumLib::initShapeMatrices<ShapeFunctionDisplacement,
+                                  ShapeMatricesTypeDisplacement,
+                                  DisplacementDim>(e, is_axially_symmetric,
+                                                   _integration_method);
 
     auto const shape_matrices_p =
-        initShapeMatrices<ShapeFunctionPressure, ShapeMatricesTypePressure,
-                          IntegrationMethod, DisplacementDim>(
+        NumLib::initShapeMatrices<ShapeFunctionPressure,
+                                  ShapeMatricesTypePressure, DisplacementDim>(
             e, is_axially_symmetric, _integration_method);
 
     auto const& solid_material =
@@ -213,9 +213,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
         auto const& dNdx_p = _ip_data[ip].dNdx_p;
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunctionDisplacement,
-                                   ShapeMatricesTypeDisplacement>(_element,
-                                                                  N_u);
+            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+                                           ShapeMatricesTypeDisplacement>(
+                _element, N_u);
         auto const B =
             LinearBMatrix::computeBMatrix<DisplacementDim,
                                           ShapeFunctionDisplacement::NPOINTS,
@@ -541,9 +541,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
 
         auto& eps = _ip_data[ip].eps;
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunctionDisplacement,
-                                   ShapeMatricesTypeDisplacement>(_element,
-                                                                  N_u);
+            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+                                           ShapeMatricesTypeDisplacement>(
+                _element, N_u);
         auto const B =
             LinearBMatrix::computeBMatrix<DisplacementDim,
                                           ShapeFunctionDisplacement::NPOINTS,
@@ -614,9 +614,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
         auto const& N_p = _ip_data[ip].N_p;
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunctionDisplacement,
-                                   ShapeMatricesTypeDisplacement>(_element,
-                                                                  N_u);
+            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+                                           ShapeMatricesTypeDisplacement>(
+                _element, N_u);
         auto const B =
             LinearBMatrix::computeBMatrix<DisplacementDim,
                                           ShapeFunctionDisplacement::NPOINTS,
@@ -725,9 +725,9 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
         auto const& dNdx_u = _ip_data[ip].dNdx_u;
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunctionDisplacement,
-                                   ShapeMatricesTypeDisplacement>(_element,
-                                                                  N_u);
+            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+                                           ShapeMatricesTypeDisplacement>(
+                _element, N_u);
         auto const B =
             LinearBMatrix::computeBMatrix<DisplacementDim,
                                           ShapeFunctionDisplacement::NPOINTS,
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
index e939b2be387e68fcc948465d78788231a657d3cd..bc35308e3d1e4a422f910aeb399046e275aeb695 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM.h
@@ -13,20 +13,19 @@
 #include <memory>
 #include <vector>
 
+#include "HydroMechanicsProcessData.h"
+#include "LocalAssemblerInterface.h"
 #include "MaterialLib/PhysicalConstant.h"
 #include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
 #include "MathLib/KelvinVector.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/Deformation/BMatrixPolicy.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
-#include "HydroMechanicsProcessData.h"
-#include "LocalAssemblerInterface.h"
 
 namespace ProcessLib
 {
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
index ad016b1f60bf71674f436fa82cba1dd4946c55bb..4b789bf9cedea88791d87036ef6d644ac32b41ac 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerFracture-impl.h
@@ -11,11 +11,8 @@
 #pragma once
 
 #include "HydroMechanicsLocalAssemblerFracture.h"
-
 #include "MaterialLib/FractureModels/FractureIdentity2.h"
-
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "ProcessLib/LIE/Common/LevelSetFunction.h"
 
 namespace ProcessLib
@@ -62,15 +59,14 @@ HydroMechanicsLocalAssemblerFracture<ShapeFunctionDisplacement,
     _ip_data.reserve(n_integration_points);
 
     auto const shape_matrices_u =
-        initShapeMatrices<ShapeFunctionDisplacement,
-                          ShapeMatricesTypeDisplacement, IntegrationMethod,
-                          GlobalDim>(e, is_axially_symmetric,
-                                     integration_method);
+        NumLib::initShapeMatrices<ShapeFunctionDisplacement,
+                                  ShapeMatricesTypeDisplacement, GlobalDim>(
+            e, is_axially_symmetric, integration_method);
 
     auto const shape_matrices_p =
-        initShapeMatrices<ShapeFunctionPressure, ShapeMatricesTypePressure,
-                          IntegrationMethod, GlobalDim>(e, is_axially_symmetric,
-                                                        integration_method);
+        NumLib::initShapeMatrices<ShapeFunctionPressure,
+                                  ShapeMatricesTypePressure, GlobalDim>(
+            e, is_axially_symmetric, integration_method);
 
     auto const& frac_prop = *_process_data.fracture_property;
 
diff --git a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
index 6a9c5211dd69d95286a3f356bec5abf7db7daa65..3166b50a934fd4252afa4fd63bc12e155487faa3 100644
--- a/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
+++ b/ProcessLib/LIE/HydroMechanics/LocalAssembler/HydroMechanicsLocalAssemblerMatrix-impl.h
@@ -11,14 +11,13 @@
 #pragma once
 
 #include "HydroMechanicsLocalAssemblerMatrix.h"
-
 #include "MaterialLib/PhysicalConstant.h"
 #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
 #include "MathLib/KelvinVector.h"
 #include "MeshLib/ElementStatus.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Function/Interpolation.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
 
 namespace ProcessLib
 {
@@ -53,15 +52,14 @@ HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
     _ip_data.reserve(n_integration_points);
 
     auto const shape_matrices_u =
-        initShapeMatrices<ShapeFunctionDisplacement,
-                          ShapeMatricesTypeDisplacement, IntegrationMethod,
-                          GlobalDim>(e, is_axially_symmetric,
-                                     integration_method);
+        NumLib::initShapeMatrices<ShapeFunctionDisplacement,
+                                  ShapeMatricesTypeDisplacement, GlobalDim>(
+            e, is_axially_symmetric, integration_method);
 
     auto const shape_matrices_p =
-        initShapeMatrices<ShapeFunctionPressure, ShapeMatricesTypePressure,
-                          IntegrationMethod, GlobalDim>(e, is_axially_symmetric,
-                                                        integration_method);
+        NumLib::initShapeMatrices<ShapeFunctionPressure,
+                                  ShapeMatricesTypePressure, GlobalDim>(
+            e, is_axially_symmetric, integration_method);
 
     auto& solid_material = MaterialLib::Solids::selectSolidConstitutiveRelation(
         _process_data.solid_materials, _process_data.material_ids, e.getID());
@@ -203,9 +201,9 @@ void HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
         auto const& H_u = ip_data.H_u;
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunctionDisplacement,
-                                   ShapeMatricesTypeDisplacement>(_element,
-                                                                  N_u);
+            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+                                           ShapeMatricesTypeDisplacement>(
+                _element, N_u);
         auto const B =
             LinearBMatrix::computeBMatrix<GlobalDim,
                                           ShapeFunctionDisplacement::NPOINTS,
@@ -339,9 +337,9 @@ void HydroMechanicsLocalAssemblerMatrix<ShapeFunctionDisplacement,
         auto const& dNdx_u = ip_data.dNdx_u;
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunctionDisplacement,
-                                   ShapeMatricesTypeDisplacement>(_element,
-                                                                  N_u);
+            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+                                           ShapeMatricesTypeDisplacement>(
+                _element, N_u);
         auto const B =
             LinearBMatrix::computeBMatrix<GlobalDim,
                                           ShapeFunctionDisplacement::NPOINTS,
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h
index e1c7c7490d6a0e54129d3f8215567f6777b43999..90070e4d6266c712379858623f637545fbf0d561 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h
@@ -10,15 +10,12 @@
 
 #pragma once
 
-#include "SmallDeformationLocalAssemblerFracture.h"
-
 #include <Eigen/Eigen>
 
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
-
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "ProcessLib/LIE/Common/LevelSetFunction.h"
+#include "SmallDeformationLocalAssemblerFracture.h"
 
 namespace ProcessLib
 {
@@ -43,9 +40,10 @@ SmallDeformationLocalAssemblerFracture<ShapeFunction, IntegrationMethod,
           dofIndex_to_localIndex),
       _process_data(process_data),
       _integration_method(integration_order),
-      _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                                        IntegrationMethod, DisplacementDim>(
-          e, is_axially_symmetric, _integration_method)),
+      _shape_matrices(
+          NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                    DisplacementDim>(e, is_axially_symmetric,
+                                                     _integration_method)),
       _element(e)
 {
     assert(_element.getDimension() == DisplacementDim - 1);
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h
index 85039a99fc5ba0202cfd1ad898e7fc230a79154e..ac39ac1631547221bf7550f44a4ea9c773d5cd39 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix-impl.h
@@ -10,27 +10,22 @@
 
 #pragma once
 
-#include "SmallDeformationLocalAssemblerMatrix.h"
-
+#include <Eigen/Eigen>
 #include <valarray>
 #include <vector>
 
-#include <Eigen/Eigen>
-
+#include "IntegrationPointDataMatrix.h"
 #include "MaterialLib/PhysicalConstant.h"
 #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
-
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ProcessLib/Deformation/BMatrixPolicy.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
 #include "ProcessLib/LIE/SmallDeformation/SmallDeformationProcessData.h"
-
-#include "IntegrationPointDataMatrix.h"
 #include "SecondaryData.h"
 #include "SmallDeformationLocalAssemblerInterface.h"
+#include "SmallDeformationLocalAssemblerMatrix.h"
 
 namespace ProcessLib
 {
@@ -60,9 +55,9 @@ SmallDeformationLocalAssemblerMatrix<ShapeFunction, IntegrationMethod,
     _secondary_data.N.resize(n_integration_points);
 
     auto const shape_matrices =
-        initShapeMatrices<ShapeFunction, ShapeMatricesType, IntegrationMethod,
-                          DisplacementDim>(e, is_axially_symmetric,
-                                           _integration_method);
+        NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                  DisplacementDim>(e, is_axially_symmetric,
+                                                   _integration_method);
 
     auto& solid_material = MaterialLib::Solids::selectSolidConstitutiveRelation(
         _process_data.solid_materials, _process_data.material_ids, e.getID());
@@ -132,8 +127,8 @@ void SmallDeformationLocalAssemblerMatrix<ShapeFunction, IntegrationMethod,
         auto const& N = _ip_data[ip].N;
         auto const& dNdx = _ip_data[ip].dNdx;
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(_element,
-                                                                     N);
+            NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+                _element, N);
         auto const B =
             LinearBMatrix::computeBMatrix<DisplacementDim,
                                           ShapeFunction::NPOINTS,
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix.h
index 8378ee6610c1600aea7a6fb91c678f0eb01e7196..654497f1fbcb93e27b2757cbb21d4d805d357df8 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrix.h
@@ -12,14 +12,12 @@
 
 #include <vector>
 
+#include "IntegrationPointDataMatrix.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ProcessLib/Deformation/BMatrixPolicy.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
 #include "ProcessLib/LIE/SmallDeformation/SmallDeformationProcessData.h"
-
-#include "IntegrationPointDataMatrix.h"
 #include "SecondaryData.h"
 #include "SmallDeformationLocalAssemblerInterface.h"
 
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
index 06b49f64877ea6b45e268af4c8a4b8895d575f0f..b423acabe3b83eabbc6067b350d6af14e0995bd6 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture-impl.h
@@ -10,34 +10,26 @@
 
 #pragma once
 
-#include "SmallDeformationLocalAssemblerMatrixNearFracture.h"
-
+#include <Eigen/Eigen>
 #include <valarray>
 #include <vector>
 
-#include <Eigen/Eigen>
-
+#include "IntegrationPointDataMatrix.h"
 #include "MaterialLib/PhysicalConstant.h"
-
 #include "MathLib/KelvinVector.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "MathLib/Point3d.h"
-
 #include "MeshLib/Node.h"
-
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
-
 #include "ProcessLib/Deformation/BMatrixPolicy.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
 #include "ProcessLib/LIE/Common/LevelSetFunction.h"
 #include "ProcessLib/LIE/Common/Utils.h"
-
-#include "IntegrationPointDataMatrix.h"
 #include "SecondaryData.h"
 #include "SmallDeformationLocalAssemblerInterface.h"
+#include "SmallDeformationLocalAssemblerMatrixNearFracture.h"
 
 namespace ProcessLib
 {
@@ -67,14 +59,12 @@ SmallDeformationLocalAssemblerMatrixNearFracture<ShapeFunction,
       _element(e),
       _is_axially_symmetric(is_axially_symmetric)
 {
-    std::vector<
-        ShapeMatrices,
-        Eigen::aligned_allocator<typename ShapeMatricesType::ShapeMatrices>>
-        shape_matrices = initShapeMatrices<ShapeFunction,
-                                           ShapeMatricesType,
-                                           IntegrationMethod,
-                                           DisplacementDim>(
-            e, is_axially_symmetric, _integration_method);
+    std::vector<ShapeMatrices, Eigen::aligned_allocator<
+                                   typename ShapeMatricesType::ShapeMatrices>>
+        shape_matrices =
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      DisplacementDim>(e, is_axially_symmetric,
+                                                       _integration_method);
 
     unsigned const n_integration_points =
         _integration_method.getNumberOfPoints();
@@ -234,8 +224,8 @@ void SmallDeformationLocalAssemblerMatrixNearFracture<
         }
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(_element,
-                                                                     N);
+            NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+                _element, N);
         auto const B =
             LinearBMatrix::computeBMatrix<DisplacementDim,
                                           ShapeFunction::NPOINTS,
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture.h
index 8e6a9e21691f3c8e488894fe11d23dd5f81e6cef..b596f06541165b3b7db207ce2c97df8c8dc6f50b 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerMatrixNearFracture.h
@@ -13,17 +13,14 @@
 #include <unordered_map>
 #include <vector>
 
+#include "IntegrationPointDataMatrix.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
-
 #include "ProcessLib/Deformation/BMatrixPolicy.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
 #include "ProcessLib/LIE/Common/FractureProperty.h"
 #include "ProcessLib/LIE/Common/JunctionProperty.h"
 #include "ProcessLib/LIE/SmallDeformation/SmallDeformationProcessData.h"
-
-#include "IntegrationPointDataMatrix.h"
 #include "SecondaryData.h"
 #include "SmallDeformationLocalAssemblerInterface.h"
 
diff --git a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
index 674e2cd2466348d828ee58b49b3d8929f9f49d0c..135236cb6664e028b3993ad421cf9b0acab29e93 100644
--- a/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
+++ b/ProcessLib/LiquidFlow/LiquidFlowLocalAssembler.h
@@ -14,20 +14,18 @@
 
 #include <vector>
 
+#include "LiquidFlowData.h"
 #include "MaterialLib/MPL/Medium.h"
 #include "MaterialLib/MPL/Utils/FormEffectiveThermalConductivity.h"
-
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
-#include "LiquidFlowData.h"
 
 namespace ProcessLib
 {
@@ -100,9 +98,9 @@ public:
         _ip_data.reserve(n_integration_points);
 
         auto const& shape_matrices =
-            initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                              IntegrationMethod, GlobalDim>(
-                element, is_axially_symmetric, _integration_method);
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      GlobalDim>(element, is_axially_symmetric,
+                                                 _integration_method);
 
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
index f549b4e9360cd65af269f578d429647f9513cfc2..d577aae2adb406f3e88e95e31bfc9d475425c2e2 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
@@ -87,8 +87,8 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
         auto const& dNdx = _ip_data[ip].dNdx;
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(_element,
-                                                                     N);
+            NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+                _element, N);
 
         auto const& B =
             LinearBMatrix::computeBMatrix<DisplacementDim,
@@ -185,8 +185,8 @@ void PhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
             double const d_ip = N.dot(d);
             double const degradation = d_ip * d_ip * (1 - k) + k;
             auto const x_coord =
-                interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
-                    _element, N);
+                NumLib::interpolateXCoordinate<ShapeFunction,
+                                               ShapeMatricesType>(_element, N);
             auto const& B = LinearBMatrix::computeBMatrix<
                 DisplacementDim, ShapeFunction::NPOINTS,
                 typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM.h b/ProcessLib/PhaseField/PhaseFieldFEM.h
index 58d02d3d817c3b7f5c3cc87bf5a53985e9fce23b..194aa829e7597895db01c9652922fc71498ba5c8 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM.h
@@ -19,12 +19,12 @@
 #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/SpatialPosition.h"
 #include "PhaseFieldProcessData.h"
 #include "ProcessLib/Deformation/BMatrixPolicy.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
 #include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
 
 namespace ProcessLib
@@ -169,9 +169,9 @@ public:
         }
 
         auto const shape_matrices =
-            initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                              IntegrationMethod, DisplacementDim>(
-                e, is_axially_symmetric, _integration_method);
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      DisplacementDim>(e, is_axially_symmetric,
+                                                       _integration_method);
 
         ParameterLib::SpatialPosition x_position;
         x_position.setElementID(_element.getID());
diff --git a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM-impl.h b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM-impl.h
index a1ab6fee86e870150066e5a0feb4bc7f0f5c270f..5ec8390e80865f51ffe6608ad89347744aa82124 100644
--- a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM-impl.h
+++ b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM-impl.h
@@ -37,9 +37,8 @@ LocalAssemblerData<ShapeFunction, IntegrationMethod, GlobalDim>::
     _ip_data.reserve(n_integration_points);
 
     auto const shape_matrices =
-        initShapeMatrices<ShapeFunction, ShapeMatricesType, IntegrationMethod,
-                          GlobalDim>(element, is_axially_symmetric,
-                                     _integration_method);
+        NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType, GlobalDim>(
+            element, is_axially_symmetric, _integration_method);
 
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
diff --git a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM.h b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM.h
index 4f76c6c267c8c3137eb5bde971238695c5408005..bdb72dabff579263987e048a5b571032307fb4ef 100644
--- a/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM.h
+++ b/ProcessLib/RichardsComponentTransport/RichardsComponentTransportFEM.h
@@ -18,11 +18,11 @@
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "NumLib/Function/Interpolation.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
 #include "RichardsComponentTransportProcessData.h"
 
 namespace ProcessLib
diff --git a/ProcessLib/RichardsFlow/RichardsFlowFEM.h b/ProcessLib/RichardsFlow/RichardsFlowFEM.h
index dd83b3c417b4b2aa82b97bc1de91d1a5d995c2d2..38e15f90655af33fb5971ba231085c63fa1133ba 100644
--- a/ProcessLib/RichardsFlow/RichardsFlowFEM.h
+++ b/ProcessLib/RichardsFlow/RichardsFlowFEM.h
@@ -21,12 +21,12 @@
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "NumLib/Function/Interpolation.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
 #include "RichardsFlowProcessData.h"
 
 namespace ProcessLib
@@ -115,9 +115,9 @@ public:
         _ip_data.reserve(n_integration_points);
 
         auto const shape_matrices =
-            initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                              IntegrationMethod, GlobalDim>(
-                element, is_axially_symmetric, _integration_method);
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      GlobalDim>(element, is_axially_symmetric,
+                                                 _integration_method);
 
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index da82b8ef3ff21bc01fd06e344ec95db2d1d1a715..3580a33a8c4fa1d3b41890a47e055c8aa62b4645 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -49,14 +49,14 @@ RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
     _secondary_data.N_u.resize(n_integration_points);
 
     auto const shape_matrices_u =
-        initShapeMatrices<ShapeFunctionDisplacement,
-                          ShapeMatricesTypeDisplacement, IntegrationMethod,
-                          DisplacementDim>(e, is_axially_symmetric,
-                                           _integration_method);
+        NumLib::initShapeMatrices<ShapeFunctionDisplacement,
+                                  ShapeMatricesTypeDisplacement,
+                                  DisplacementDim>(e, is_axially_symmetric,
+                                                   _integration_method);
 
     auto const shape_matrices_p =
-        initShapeMatrices<ShapeFunctionPressure, ShapeMatricesTypePressure,
-                          IntegrationMethod, DisplacementDim>(
+        NumLib::initShapeMatrices<ShapeFunctionPressure,
+                                  ShapeMatricesTypePressure, DisplacementDim>(
             e, is_axially_symmetric, _integration_method);
 
     auto const& solid_material =
@@ -303,9 +303,9 @@ void RichardsMechanicsLocalAssembler<
         auto const& dNdx_p = _ip_data[ip].dNdx_p;
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunctionDisplacement,
-                                   ShapeMatricesTypeDisplacement>(_element,
-                                                                  N_u);
+            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+                                           ShapeMatricesTypeDisplacement>(
+                _element, N_u);
         auto const B =
             LinearBMatrix::computeBMatrix<DisplacementDim,
                                           ShapeFunctionDisplacement::NPOINTS,
@@ -642,9 +642,9 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         auto const& dNdx_p = _ip_data[ip].dNdx_p;
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunctionDisplacement,
-                                   ShapeMatricesTypeDisplacement>(_element,
-                                                                  N_u);
+            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+                                           ShapeMatricesTypeDisplacement>(
+                _element, N_u);
         auto const B =
             LinearBMatrix::computeBMatrix<DisplacementDim,
                                           ShapeFunctionDisplacement::NPOINTS,
@@ -1313,9 +1313,9 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         variables[static_cast<int>(MPL::Variable::temperature)] = temperature;
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunctionDisplacement,
-                                   ShapeMatricesTypeDisplacement>(_element,
-                                                                  N_u);
+            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+                                           ShapeMatricesTypeDisplacement>(
+                _element, N_u);
         auto const B =
             LinearBMatrix::computeBMatrix<DisplacementDim,
                                           ShapeFunctionDisplacement::NPOINTS,
@@ -1387,9 +1387,9 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
         auto const& dNdx_u = _ip_data[ip].dNdx_u;
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunctionDisplacement,
-                                   ShapeMatricesTypeDisplacement>(_element,
-                                                                  N_u);
+            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+                                           ShapeMatricesTypeDisplacement>(
+                _element, N_u);
         auto const B =
             LinearBMatrix::computeBMatrix<DisplacementDim,
                                           ShapeFunctionDisplacement::NPOINTS,
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
index c161c0a7812d57455274993644c8916cab08bdb6..93aa12a8b5ca181c9a686b9870dd0b043b6c0556 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM.h
@@ -13,19 +13,18 @@
 #include <memory>
 #include <vector>
 
+#include "IntegrationPointData.h"
+#include "LocalAssemblerInterface.h"
 #include "MaterialLib/SolidModels/LinearElasticIsotropic.h"
 #include "MathLib/KelvinVector.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/Deformation/BMatrixPolicy.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
-#include "IntegrationPointData.h"
-#include "LocalAssemblerInterface.h"
 #include "RichardsMechanicsProcessData.h"
 
 namespace ProcessLib
@@ -129,10 +128,10 @@ public:
             {
                 ParameterLib::SpatialPosition const x_position{
                     boost::none, _element.getID(), ip,
-                    MathLib::Point3d(
-                        interpolateCoordinates<ShapeFunctionDisplacement,
-                                               ShapeMatricesTypeDisplacement>(
-                            _element, ip_data.N_u))};
+                    MathLib::Point3d(NumLib::interpolateCoordinates<
+                                     ShapeFunctionDisplacement,
+                                     ShapeMatricesTypeDisplacement>(
+                        _element, ip_data.N_u))};
 
                 ip_data.sigma_eff =
                     MathLib::KelvinVector::symmetricTensorToKelvinVector<
diff --git a/ProcessLib/SmallDeformation/SmallDeformationFEM.h b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
index fcc28dc600aa5b98476013a4a7c528081d51794d..00f24fd70dc19445a33c834b0960af52efc551fb 100644
--- a/ProcessLib/SmallDeformation/SmallDeformationFEM.h
+++ b/ProcessLib/SmallDeformation/SmallDeformationFEM.h
@@ -19,6 +19,7 @@
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/Deformation/BMatrixPolicy.h"
@@ -26,7 +27,6 @@
 #include "ProcessLib/Deformation/LinearBMatrix.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
 #include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
 #include "SmallDeformationProcessData.h"
 
@@ -125,9 +125,9 @@ public:
         _secondary_data.N.resize(n_integration_points);
 
         auto const shape_matrices =
-            initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                              IntegrationMethod, DisplacementDim>(
-                e, is_axially_symmetric, _integration_method);
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      DisplacementDim>(e, is_axially_symmetric,
+                                                       _integration_method);
 
         auto& solid_material =
             MaterialLib::Solids::selectSolidConstitutiveRelation(
@@ -206,9 +206,10 @@ public:
             {
                 ParameterLib::SpatialPosition const x_position{
                     boost::none, _element.getID(), ip,
-                    MathLib::Point3d(interpolateCoordinates<ShapeFunction,
-                                                            ShapeMatricesType>(
-                        _element, ip_data.N))};
+                    MathLib::Point3d(
+                        NumLib::interpolateCoordinates<ShapeFunction,
+                                                       ShapeMatricesType>(
+                            _element, ip_data.N))};
 
                 ip_data.sigma =
                     MathLib::KelvinVector::symmetricTensorToKelvinVector<
@@ -279,8 +280,8 @@ public:
             }
 
             auto const x_coord =
-                interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
-                    _element, N);
+                NumLib::interpolateXCoordinate<ShapeFunction,
+                                               ShapeMatricesType>(_element, N);
             auto const B = LinearBMatrix::computeBMatrix<
                 DisplacementDim, ShapeFunction::NPOINTS,
                 typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
diff --git a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalFEM.h b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalFEM.h
index 525053841db49cfcd7afd8ec2b4c33f7aaa6778a..a228f3314fa1d213837f9fedeb5c44443762ddb3 100644
--- a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalFEM.h
+++ b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalFEM.h
@@ -23,6 +23,7 @@
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "MeshLib/findElementsWithinRadius.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "NumLib/Function/Interpolation.h"
 #include "ParameterLib/Parameter.h"
@@ -30,7 +31,6 @@
 #include "ProcessLib/Deformation/Divergence.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
 #include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
 #include "SmallDeformationNonlocalProcessData.h"
 
@@ -91,9 +91,9 @@ public:
         _secondary_data.N.resize(n_integration_points);
 
         auto const shape_matrices =
-            initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                              IntegrationMethod, DisplacementDim>(
-                e, is_axially_symmetric, _integration_method);
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      DisplacementDim>(e, is_axially_symmetric,
+                                                       _integration_method);
 
         auto& solid_material =
             MaterialLib::Solids::selectSolidConstitutiveRelation(
@@ -337,8 +337,8 @@ public:
             auto const& dNdx = _ip_data[ip].dNdx;
 
             auto const x_coord =
-                interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
-                    _element, N);
+                NumLib::interpolateXCoordinate<ShapeFunction,
+                                               ShapeMatricesType>(_element, N);
             auto const B = LinearBMatrix::computeBMatrix<
                 DisplacementDim, ShapeFunction::NPOINTS,
                 typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
@@ -456,8 +456,8 @@ public:
             auto const& dNdx = _ip_data[ip].dNdx;
 
             auto const x_coord =
-                interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
-                    _element, N);
+                NumLib::interpolateXCoordinate<ShapeFunction,
+                                               ShapeMatricesType>(_element, N);
             auto const B = LinearBMatrix::computeBMatrix<
                 DisplacementDim, ShapeFunction::NPOINTS,
                 typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
@@ -595,8 +595,8 @@ public:
             auto const& dNdx = _ip_data[ip].dNdx;
 
             auto const x_coord =
-                interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
-                    _element, N);
+                NumLib::interpolateXCoordinate<ShapeFunction,
+                                               ShapeMatricesType>(_element, N);
             auto const B = LinearBMatrix::computeBMatrix<
                 DisplacementDim, ShapeFunction::NPOINTS,
                 typename BMatricesType::BMatrixType>(dNdx, N, x_coord,
diff --git a/ProcessLib/SourceTerms/LineSourceTermFEM.h b/ProcessLib/SourceTerms/LineSourceTermFEM.h
index f811a406f7abdaa082764b115160e23d0cb1e770..24fedf884034afb484ebd13f94f56d872c17279d 100644
--- a/ProcessLib/SourceTerms/LineSourceTermFEM.h
+++ b/ProcessLib/SourceTerms/LineSourceTermFEM.h
@@ -15,9 +15,9 @@
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
 #include "SourceTermIntegrationPointData.h"
 
 namespace ProcessLib
@@ -63,9 +63,9 @@ public:
             _integration_method.getNumberOfPoints();
 
         auto const shape_matrices =
-            initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                              IntegrationMethod, GlobalDim>(
-                _element, is_axially_symmetric, _integration_method);
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      GlobalDim>(_element, is_axially_symmetric,
+                                                 _integration_method);
 
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
@@ -94,8 +94,9 @@ public:
             ParameterLib::SpatialPosition const pos{
                 boost::none, _element.getID(), ip,
                 MathLib::Point3d(
-                    interpolateCoordinates<ShapeFunction, ShapeMatricesType>(
-                        _element, N))};
+                    NumLib::interpolateCoordinates<ShapeFunction,
+                                                   ShapeMatricesType>(_element,
+                                                                      N))};
             auto const st_val = _parameter(t, pos)[0];
 
             _local_rhs.noalias() += st_val * w * N;
diff --git a/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h b/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h
index afa1c555a3fb2651c9b1a3f5140a67dd214257bc..379ebd37d508423f5bf52ad73aa8d9799f261d0b 100644
--- a/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h
+++ b/ProcessLib/SourceTerms/Python/PythonSourceTermLocalAssembler.h
@@ -10,13 +10,12 @@
 
 #pragma once
 
-#include "PythonSourceTerm.h"
-
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
+#include "PythonSourceTerm.h"
 #include "PythonSourceTermLocalAssemblerInterface.h"
 #include "PythonSourceTermPythonSideInterface.h"
 
@@ -72,8 +71,8 @@ public:
             _integration_method.getNumberOfPoints();
 
         auto const shape_matrices =
-            initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                              IntegrationMethod, GlobalDim>(
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      GlobalDim>(
                 _element, _is_axially_symmetric, _integration_method);
 
         for (unsigned ip = 0; ip < n_integration_points; ip++)
diff --git a/ProcessLib/SourceTerms/VolumetricSourceTermFEM.h b/ProcessLib/SourceTerms/VolumetricSourceTermFEM.h
index fbeb33d42c8035fc8b80e7f15ef19fb5d33693af..e6e9030bc0b09a4cb6dd5b8e69207b7d26601f80 100644
--- a/ProcessLib/SourceTerms/VolumetricSourceTermFEM.h
+++ b/ProcessLib/SourceTerms/VolumetricSourceTermFEM.h
@@ -15,9 +15,9 @@
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
 #include "SourceTermIntegrationPointData.h"
 
 namespace ProcessLib
@@ -63,9 +63,9 @@ public:
             _integration_method.getNumberOfPoints();
 
         auto const shape_matrices =
-            initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                              IntegrationMethod, GlobalDim>(
-                _element, is_axially_symmetric, _integration_method);
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      GlobalDim>(_element, is_axially_symmetric,
+                                                 _integration_method);
 
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
@@ -94,8 +94,9 @@ public:
             ParameterLib::SpatialPosition const pos{
                 boost::none, _element.getID(), ip,
                 MathLib::Point3d(
-                    interpolateCoordinates<ShapeFunction, ShapeMatricesType>(
-                        _element, N))};
+                    NumLib::interpolateCoordinates<ShapeFunction,
+                                                   ShapeMatricesType>(_element,
+                                                                      N))};
             auto const st_val = _volumetric_source_term(t, pos)[0];
 
             _local_rhs.noalias() += st_val * w * N;
diff --git a/ProcessLib/SteadyStateDiffusion/SteadyStateDiffusionFEM.h b/ProcessLib/SteadyStateDiffusion/SteadyStateDiffusionFEM.h
index 834b1cf49a5a3094891830a01cce02f14d657c6a..2bcbabbc3b88d13422764267b7118cba6d1b5e2f 100644
--- a/ProcessLib/SteadyStateDiffusion/SteadyStateDiffusionFEM.h
+++ b/ProcessLib/SteadyStateDiffusion/SteadyStateDiffusionFEM.h
@@ -12,7 +12,6 @@
 
 #include <vector>
 
-#include "SteadyStateDiffusionData.h"
 #include "MaterialLib/MPL/Medium.h"
 #include "MaterialLib/MPL/Utils/FormEigenTensor.h"
 #include "MaterialLib/MPL/VariableType.h"
@@ -20,12 +19,13 @@
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "NumLib/Function/Interpolation.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
+#include "SteadyStateDiffusionData.h"
 
 namespace ProcessLib
 {
@@ -70,9 +70,10 @@ public:
         : _element(element),
           _process_data(process_data),
           _integration_method(integration_order),
-          _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                                            IntegrationMethod, GlobalDim>(
-              element, is_axially_symmetric, _integration_method))
+          _shape_matrices(
+              NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                        GlobalDim>(
+                  element, is_axially_symmetric, _integration_method))
     {
     }
 
diff --git a/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h b/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h
index 01b81bc20ee1398e059ecfeb30690af6791b5930..49c33741bc2e04a13c8d57f1150e74467ad2f523 100644
--- a/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h
+++ b/ProcessLib/SurfaceFlux/SurfaceFluxLocalAssembler.h
@@ -10,17 +10,14 @@
 
 #pragma once
 
-#include "NumLib/Fem/ShapeMatrixPolicy.h"
-
+#include "MeshLib/Elements/Elements.h"
+#include "MeshLib/Elements/FaceRule.h"
+#include "MeshLib/Elements/MapBulkElementPoint.h"
 #include "NumLib/DOF/DOFTableUtil.h"
-
+#include "NumLib/Fem/InitShapeMatrices.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/Process.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
-#include "MeshLib/Elements/MapBulkElementPoint.h"
-#include "MeshLib/Elements/FaceRule.h"
-#include "MeshLib/Elements/Elements.h"
 
 namespace ProcessLib
 {
@@ -73,25 +70,15 @@ public:
           _bulk_element_id(bulk_element_ids[surface_element.getID()]),
           _bulk_face_id(bulk_face_ids[surface_element.getID()])
     {
-        auto const fe = NumLib::createIsoparametricFiniteElement<
-            ShapeFunction, ShapeMatricesType>(_surface_element);
-
-        std::size_t const n_integration_points =
+        auto const n_integration_points =
             _integration_method.getNumberOfPoints();
 
-        std::vector<
-            typename ShapeMatricesType::ShapeMatrices,
-            Eigen::aligned_allocator<typename ShapeMatricesType::ShapeMatrices>>
-            shape_matrices;
-        shape_matrices.reserve(n_integration_points);
-        _detJ_times_integralMeasure.reserve(n_integration_points);
+        auto const shape_matrices =
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      GlobalDim, NumLib::ShapeMatrixType::N_J>(
+                _surface_element, is_axially_symmetric, _integration_method);
         for (std::size_t ip = 0; ip < n_integration_points; ++ip)
         {
-            shape_matrices.emplace_back(ShapeFunction::DIM, GlobalDim,
-                                        ShapeFunction::NPOINTS);
-            fe.template computeShapeFunctions<NumLib::ShapeMatrixType::N_J>(
-                _integration_method.getWeightedPoint(ip).getCoords(),
-                shape_matrices[ip], GlobalDim, is_axially_symmetric);
             _detJ_times_integralMeasure.push_back(
                 shape_matrices[ip].detJ * shape_matrices[ip].integralMeasure);
         }
diff --git a/ProcessLib/TES/TESLocalAssembler-impl.h b/ProcessLib/TES/TESLocalAssembler-impl.h
index 2e53aa07e2053ddafc2fb78b0ab34abf942fc1ac..14d8e773aa602a09e467fa1fac36052582c82625 100644
--- a/ProcessLib/TES/TESLocalAssembler-impl.h
+++ b/ProcessLib/TES/TESLocalAssembler-impl.h
@@ -13,10 +13,9 @@
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "NumLib/Function/Interpolation.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
 #include "TESLocalAssembler.h"
 #include "TESReactionAdaptor.h"
 
@@ -126,8 +125,8 @@ TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::
                       AssemblyParams const& asm_params)
     : _element(e),
       _integration_method(integration_order),
-      _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                                        IntegrationMethod_, GlobalDim>(
+      _shape_matrices(NumLib::initShapeMatrices<ShapeFunction,
+                                                ShapeMatricesType, GlobalDim>(
           e, is_axially_symmetric, _integration_method)),
       _d(asm_params, _integration_method.getNumberOfPoints(), GlobalDim)
 {
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler.h b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler.h
index 298252db51776618c849737951edb6ea35442387..739b6ef5af60cfea923e11ca95ac300e9851c3bb 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler.h
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPLocalAssembler.h
@@ -16,12 +16,11 @@
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
 #include "ThermalTwoPhaseFlowWithPPProcessData.h"
 
 namespace ProcessLib
@@ -114,9 +113,9 @@ public:
             _integration_method.getNumberOfPoints();
         _ip_data.reserve(n_integration_points);
         auto const shape_matrices =
-            initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                              IntegrationMethod, GlobalDim>(
-                element, is_axially_symmetric, _integration_method);
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      GlobalDim>(element, is_axially_symmetric,
+                                                 _integration_method);
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
             auto const& sm = shape_matrices[ip];
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
index 13c00ac9a70536643eddee2f961a4fea30022135..aadd7e619c9a4ebc223ad1e6ff9f9e9c1b6cb424 100644
--- a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
@@ -51,14 +51,14 @@ ThermoHydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
     _secondary_data.N_u.resize(n_integration_points);
 
     auto const shape_matrices_u =
-        initShapeMatrices<ShapeFunctionDisplacement,
-                          ShapeMatricesTypeDisplacement, IntegrationMethod,
-                          DisplacementDim>(e, is_axially_symmetric,
-                                           _integration_method);
+        NumLib::initShapeMatrices<ShapeFunctionDisplacement,
+                                  ShapeMatricesTypeDisplacement,
+                                  DisplacementDim>(e, is_axially_symmetric,
+                                                   _integration_method);
 
     auto const shape_matrices_p =
-        initShapeMatrices<ShapeFunctionPressure, ShapeMatricesTypePressure,
-                          IntegrationMethod, DisplacementDim>(
+        NumLib::initShapeMatrices<ShapeFunctionPressure,
+                                  ShapeMatricesTypePressure, DisplacementDim>(
             e, is_axially_symmetric, _integration_method);
 
     auto const& solid_material =
@@ -212,9 +212,9 @@ void ThermoHydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
         auto const T_int_pt = N_T.dot(T);
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunctionDisplacement,
-                                   ShapeMatricesTypeDisplacement>(_element,
-                                                                  N_u);
+            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+                                           ShapeMatricesTypeDisplacement>(
+                _element, N_u);
         auto const B =
             LinearBMatrix::computeBMatrix<DisplacementDim,
                                           ShapeFunctionDisplacement::NPOINTS,
@@ -571,9 +571,9 @@ void ThermoHydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
         auto const& dNdx_u = _ip_data[ip].dNdx_u;
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunctionDisplacement,
-                                   ShapeMatricesTypeDisplacement>(_element,
-                                                                  N_u);
+            NumLib::interpolateXCoordinate<ShapeFunctionDisplacement,
+                                           ShapeMatricesTypeDisplacement>(
+                _element, N_u);
         auto const B =
             LinearBMatrix::computeBMatrix<DisplacementDim,
                                           ShapeFunctionDisplacement::NPOINTS,
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM.h
index 5bddbba86d87fc4cc4a710322bb9ebc5afd2dc15..834c83f569b186725d35823e8dcfb1f8a32f9591 100644
--- a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM.h
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM.h
@@ -20,12 +20,12 @@
 #include "MathLib/KelvinVector.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/Deformation/BMatrixPolicy.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
 #include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
 #include "ThermoHydroMechanicsProcessData.h"
 
diff --git a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h
index e9a8afbd94591b543fada938a5dafedeae11e924..33e83927989449292f989ededf63fb818648a834 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h
+++ b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM-impl.h
@@ -94,8 +94,8 @@ void ThermoMechanicalPhaseFieldLocalAssembler<ShapeFunction, IntegrationMethod,
         auto const& dNdx = _ip_data[ip].dNdx;
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(_element,
-                                                                     N);
+            NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+                _element, N);
         auto const& B =
             LinearBMatrix::computeBMatrix<DisplacementDim,
                                           ShapeFunction::NPOINTS,
diff --git a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
index 0df4e3d795e795605ab4a295c4231e536396f6c8..47e0b3fc1b3342543ab59f4d03f6238f181741b7 100644
--- a/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
+++ b/ProcessLib/ThermoMechanicalPhaseField/ThermoMechanicalPhaseFieldFEM.h
@@ -19,11 +19,11 @@
 #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/SpatialPosition.h"
 #include "ProcessLib/Deformation/BMatrixPolicy.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
 #include "ProcessLib/Utils/SetOrGetIntegrationPointData.h"
 #include "ThermoMechanicalPhaseFieldProcessData.h"
 
@@ -184,11 +184,10 @@ public:
                 "support is implemented.");
         }
 
-
         auto const shape_matrices =
-            initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                              IntegrationMethod, DisplacementDim>(
-                e, is_axially_symmetric, _integration_method);
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      DisplacementDim>(e, is_axially_symmetric,
+                                                       _integration_method);
 
         ParameterLib::SpatialPosition x_position;
         x_position.setElementID(_element.getID());
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h
index 2ea0cd0fe075aaad6916fb120200ffe15c3aa964..ccee60157b59dd89de0e76899e9c37b124ac03ac 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM-impl.h
@@ -41,9 +41,9 @@ ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
     _secondary_data.N.resize(n_integration_points);
 
     auto const shape_matrices =
-        initShapeMatrices<ShapeFunction, ShapeMatricesType, IntegrationMethod,
-                          DisplacementDim>(e, is_axially_symmetric,
-                                           _integration_method);
+        NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                  DisplacementDim>(e, is_axially_symmetric,
+                                                   _integration_method);
 
     auto& solid_material = MaterialLib::Solids::selectSolidConstitutiveRelation(
         _process_data.solid_materials, _process_data.material_ids, e.getID());
@@ -174,8 +174,8 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
         auto const& dNdx = _ip_data[ip].dNdx;
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(_element,
-                                                                     N);
+            NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+                _element, N);
         auto const& B =
             LinearBMatrix::computeBMatrix<DisplacementDim,
                                           ShapeFunction::NPOINTS,
@@ -376,8 +376,8 @@ void ThermoMechanicsLocalAssembler<ShapeFunction, IntegrationMethod,
         auto const& dNdx = _ip_data[ip].dNdx;
 
         auto const x_coord =
-            interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(_element,
-                                                                     N);
+            NumLib::interpolateXCoordinate<ShapeFunction, ShapeMatricesType>(
+                _element, N);
         auto const& B =
             LinearBMatrix::computeBMatrix<DisplacementDim,
                                           ShapeFunction::NPOINTS,
diff --git a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
index afeba9d82e63518fc873c45cccd0c55cccf8f60c..ce76296bf3a9da62edc81a565cccdf9e48b2c961 100644
--- a/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
+++ b/ProcessLib/ThermoMechanics/ThermoMechanicsFEM.h
@@ -13,18 +13,17 @@
 #include <memory>
 #include <vector>
 
+#include "LocalAssemblerInterface.h"
 #include "MaterialLib/SolidModels/SelectSolidConstitutiveRelation.h"
 #include "MathLib/KelvinVector.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/Deformation/BMatrixPolicy.h"
 #include "ProcessLib/Deformation/LinearBMatrix.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
-#include "LocalAssemblerInterface.h"
 #include "ThermoMechanicsProcessData.h"
 
 namespace ProcessLib
@@ -166,9 +165,10 @@ public:
             {
                 ParameterLib::SpatialPosition const x_position{
                     boost::none, _element.getID(), ip,
-                    MathLib::Point3d(interpolateCoordinates<ShapeFunction,
-                                                            ShapeMatricesType>(
-                        _element, ip_data.N))};
+                    MathLib::Point3d(
+                        NumLib::interpolateCoordinates<ShapeFunction,
+                                                       ShapeMatricesType>(
+                            _element, ip_data.N))};
 
                 ip_data.sigma =
                     MathLib::KelvinVector::symmetricTensorToKelvinVector<
diff --git a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h
index f695f6e5bb3e5e45bc3b93e97ebd57991bf1ab26..4aca9768855db218668566d1e2e9befe787c9743 100644
--- a/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h
+++ b/ProcessLib/TwoPhaseFlowWithPP/TwoPhaseFlowWithPPLocalAssembler.h
@@ -15,12 +15,11 @@
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
 #include "TwoPhaseFlowWithPPProcessData.h"
 
 namespace ProcessLib
@@ -109,9 +108,9 @@ public:
             _integration_method.getNumberOfPoints();
         _ip_data.reserve(n_integration_points);
         auto const shape_matrices =
-            initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                              IntegrationMethod, GlobalDim>(
-                element, is_axially_symmetric, _integration_method);
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      GlobalDim>(element, is_axially_symmetric,
+                                                 _integration_method);
         for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
             auto const& sm = shape_matrices[ip];
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler.h b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler.h
index 023a5d9a420d14a55dabe6849798df7212771371..10ed053e4c3977d3f32e0be48947678590edcfe6 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler.h
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoLocalAssembler.h
@@ -11,16 +11,16 @@
 #pragma once
 
 #include <vector>
+
 #include "MaterialLib/PhysicalConstant.h"
 #include "MathLib/LinAlg/Eigen/EigenMapTools.h"
 #include "NumLib/Extrapolation/ExtrapolatableElement.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "ParameterLib/Parameter.h"
 #include "ProcessLib/LocalAssemblerInterface.h"
 #include "ProcessLib/LocalAssemblerTraits.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
 #include "TwoPhaseFlowWithPrhoProcessData.h"
 
 namespace ProcessLib
@@ -103,9 +103,10 @@ public:
         TwoPhaseFlowWithPrhoProcessData const& process_data)
         : _element(element),
           _integration_method(integration_order),
-          _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                                            IntegrationMethod, GlobalDim>(
-              element, is_axially_symmetric, _integration_method)),
+          _shape_matrices(
+              NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                        GlobalDim>(
+                  element, is_axially_symmetric, _integration_method)),
           _process_data(process_data),
           _saturation(
               std::vector<double>(_integration_method.getNumberOfPoints())),
diff --git a/Tests/MathLib/TestGaussLegendreIntegration.cpp b/Tests/MathLib/TestGaussLegendreIntegration.cpp
index 42c9f01311c09dea5de6b3d4b7924a6f146d35f1..7a474fe0810bba82590c007d33503537ac99ac7f 100644
--- a/Tests/MathLib/TestGaussLegendreIntegration.cpp
+++ b/Tests/MathLib/TestGaussLegendreIntegration.cpp
@@ -13,19 +13,14 @@
 #include <limits>
 
 #include "InfoLib/TestInfo.h"
-
 #include "MathLib/Integration/GaussLegendreTet.h"
-
 #include "MeshLib/Elements/Element.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
 #include "MeshLib/MeshSubset.h"
-
 #include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
-
 #include "ProcessLib/Utils/CreateLocalAssemblers.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
-
 #include "Tests/VectorUtils.h"
 
 namespace GaussLegendreTest
@@ -89,8 +84,8 @@ public:
         double integral = 0;
 
         auto const sms =
-            ProcessLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                                          IntegrationMethod, GlobalDim>(
+            NumLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                      GlobalDim>(
                 _e, false /*is_axially_symmetric*/,
                 IntegrationMethod{integration_order});
         IntegrationMethod integration_method{integration_order};
diff --git a/Tests/NumLib/TestExtrapolation.cpp b/Tests/NumLib/TestExtrapolation.cpp
index e32f6867d4d8fe4c398ab8ea8c96b8393c404cd4..4a383fd72a446f249b48fcc393db8088cb65a2ad 100644
--- a/Tests/NumLib/TestExtrapolation.cpp
+++ b/Tests/NumLib/TestExtrapolation.cpp
@@ -10,24 +10,20 @@
 #include <gtest/gtest.h>
 
 #include "MathLib/LinAlg/LinAlg.h"
-
 #include "MeshLib/IO/writeMeshToFile.h"
 #include "MeshLib/MeshGenerators/MeshGenerator.h"
-
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "NumLib/DOF/MatrixProviderUser.h"
 #include "NumLib/Extrapolation/ExtrapolatableElementCollection.h"
 #include "NumLib/Extrapolation/Extrapolator.h"
 #include "NumLib/Extrapolation/LocalLinearLeastSquaresExtrapolator.h"
 #include "NumLib/Fem/FiniteElement/TemplateIsoparametric.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
 #include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "NumLib/Function/Interpolation.h"
 #include "NumLib/NumericsConfig.h"
-
 #include "ProcessLib/Utils/CreateLocalAssemblers.h"
-#include "ProcessLib/Utils/InitShapeMatrices.h"
 #include "ProcessLib/Utils/LocalDataInitializer.h"
-
 #include "Tests/VectorUtils.h"
 
 namespace ExtrapolationTest
@@ -86,11 +82,9 @@ public:
                        std::size_t const /*local_matrix_size*/,
                        bool is_axially_symmetric,
                        unsigned const integration_order)
-        : _shape_matrices(
-              ProcessLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
-                                            IntegrationMethod, GlobalDim>(
-                  e, is_axially_symmetric,
-                  IntegrationMethod{integration_order})),
+        : _shape_matrices(NumLib::initShapeMatrices<
+                          ShapeFunction, ShapeMatricesType, GlobalDim>(
+              e, is_axially_symmetric, IntegrationMethod{integration_order})),
           _int_pt_values(_shape_matrices.size())
     {
     }
diff --git a/Tests/NumLib/TestFe.cpp b/Tests/NumLib/TestFe.cpp
index 2bb4c92950d6f7a1d8d385634c94d478432b0d24..395f0b369508a3ae99e9bf43797a18baa08baf3b 100644
--- a/Tests/NumLib/TestFe.cpp
+++ b/Tests/NumLib/TestFe.cpp
@@ -9,29 +9,26 @@
 
 #include <gtest/gtest.h>
 
-#include <vector>
-#include <cmath>
-
 #include <Eigen/Eigen>
+#include <cmath>
+#include <vector>
 
-#include "MeshLib/ElementCoordinatesMappingLocal.h"
-#include "MeshLib/CoordinateSystem.h"
-
-#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
-#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
-#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
-#include "NumLib/Fem/ShapeMatrixPolicy.h"
-
+#include "FeTestData/TestFeHEX8.h"
 #include "FeTestData/TestFeLINE2.h"
 #include "FeTestData/TestFeLINE2Y.h"
 #include "FeTestData/TestFeLINE3.h"
-#include "FeTestData/TestFeTRI3.h"
-#include "FeTestData/TestFeQUAD4.h"
-#include "FeTestData/TestFeHEX8.h"
-#include "FeTestData/TestFeTET4.h"
 #include "FeTestData/TestFePRISM6.h"
 #include "FeTestData/TestFePYRA5.h"
-
+#include "FeTestData/TestFeQUAD4.h"
+#include "FeTestData/TestFeTET4.h"
+#include "FeTestData/TestFeTRI3.h"
+#include "MeshLib/CoordinateSystem.h"
+#include "MeshLib/ElementCoordinatesMappingLocal.h"
+#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
+#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
+#include "NumLib/Fem/InitShapeMatrices.h"
+#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
+#include "NumLib/Fem/ShapeMatrixPolicy.h"
 #include "Tests/TestTools.h"
 
 using namespace NumLib;
@@ -49,6 +46,7 @@ struct TestCase
     using ShapeMatrixTypes = ShapeMatrixPolicy_<typename TestFeType::ShapeFunction, GlobalDim>;
     template <typename X>
     using ShapeMatrixPolicy = ShapeMatrixPolicy_<X, GlobalDim>;
+    using ShapeFunction = typename TestFeType::ShapeFunction;
 };
 
 using TestTypes =
@@ -78,6 +76,7 @@ class NumLibFemIsoTest : public ::testing::Test, public T::TestFeType
 {
  public:
      using ShapeMatrixTypes = typename T::ShapeMatrixTypes;
+     using ShapeFunction = typename T::ShapeFunction;
      using TestFeType = typename T::TestFeType;
      // Matrix types
      using NodalMatrix = typename ShapeMatrixTypes::NodalMatrixType;
@@ -191,21 +190,21 @@ TYPED_TEST_CASE(NumLibFemIsoTest, TestTypes);
 TYPED_TEST(NumLibFemIsoTest, CheckMassMatrix)
 {
     // Refer to typedefs in the fixture
-    using FeType = typename TestFixture::FeType;
     using NodalMatrix = typename TestFixture::NodalMatrix;
-    using ShapeMatricesType = typename TestFixture::ShapeMatricesType;
 
-    // create a finite element object
-    FeType fe(*this->mesh_element);
+    auto const shape_matrices =
+        NumLib::initShapeMatrices<typename TestFixture::ShapeFunction,
+                                  typename TestFixture::ShapeMatrixTypes,
+                                  TestFixture::global_dim,
+                                  ShapeMatrixType::N_J>(
+            *this->mesh_element, false /*is_axially_symmetric*/,
+            this->integration_method);
 
     // evaluate a mass matrix M = int{ N^T D N }dA_e
     NodalMatrix M = NodalMatrix::Zero(this->e_nnodes, this->e_nnodes);
-    ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
     for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
-        shape.setZero();
+        auto const& shape = shape_matrices[i];
         auto wp = this->integration_method.getWeightedPoint(i);
-        fe.template computeShapeFunctions<ShapeMatrixType::N_J>(
-            wp.getCoords(), shape, this->global_dim, false);
         M.noalias() += shape.N.transpose() * shape.N * shape.detJ * wp.getWeight();
     }
     //std::cout << "M=\n" << M;
@@ -215,21 +214,21 @@ TYPED_TEST(NumLibFemIsoTest, CheckMassMatrix)
 TYPED_TEST(NumLibFemIsoTest, CheckLaplaceMatrix)
 {
     // Refer to typedefs in the fixture
-    using FeType = typename TestFixture::FeType;
     using NodalMatrix = typename TestFixture::NodalMatrix;
-    using ShapeMatricesType = typename TestFixture::ShapeMatricesType;
 
-    // create a finite element object
-    FeType fe(*this->mesh_element);
+    auto const shape_matrices =
+        NumLib::initShapeMatrices<typename TestFixture::ShapeFunction,
+                                  typename TestFixture::ShapeMatrixTypes,
+                                  TestFixture::global_dim,
+                                  ShapeMatrixType::DNDX>(
+            *this->mesh_element, false /*is_axially_symmetric*/,
+            this->integration_method);
 
     // evaluate a Laplace matrix K = int{ dNdx^T D dNdx }dA_e
     NodalMatrix K = NodalMatrix::Zero(this->e_nnodes, this->e_nnodes);
-    ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
     for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
-        shape.setZero();
+        auto const& shape = shape_matrices[i];
         auto wp = this->integration_method.getWeightedPoint(i);
-        fe.template computeShapeFunctions<ShapeMatrixType::DNDX>(
-            wp.getCoords(), shape, this->global_dim, false);
         K.noalias() += shape.dNdx.transpose() * this->globalD * shape.dNdx * shape.detJ * wp.getWeight();
     }
     //std::cout << "K=\n" << K << std::endl;
@@ -239,21 +238,21 @@ TYPED_TEST(NumLibFemIsoTest, CheckLaplaceMatrix)
 TYPED_TEST(NumLibFemIsoTest, CheckMassLaplaceMatrices)
 {
     // Refer to typedefs in the fixture
-    using FeType = typename TestFixture::FeType;
     using NodalMatrix = typename TestFixture::NodalMatrix;
-    using ShapeMatricesType = typename TestFixture::ShapeMatricesType;
 
-    // create a finite element object
-    FeType fe(*this->mesh_element);
+    auto const shape_matrices =
+        NumLib::initShapeMatrices<typename TestFixture::ShapeFunction,
+                                  typename TestFixture::ShapeMatrixTypes,
+                                  TestFixture::global_dim>(
+            *this->mesh_element, false /*is_axially_symmetric*/,
+            this->integration_method);
 
     // evaluate both mass and laplace matrices at once
     NodalMatrix M = NodalMatrix::Zero(this->e_nnodes, this->e_nnodes);
     NodalMatrix K = NodalMatrix::Zero(this->e_nnodes, this->e_nnodes);
-    ShapeMatricesType shape(this->dim, this->global_dim, this->e_nnodes);
     for (std::size_t i=0; i < this->integration_method.getNumberOfPoints(); i++) {
-        shape.setZero();
+        auto const& shape = shape_matrices[i];
         auto wp = this->integration_method.getWeightedPoint(i);
-        fe.computeShapeFunctions(wp.getCoords(), shape, this->global_dim, false);
         M.noalias() += shape.N.transpose() * shape.N * shape.detJ * wp.getWeight();
         K.noalias() += shape.dNdx.transpose() * (this->globalD * shape.dNdx) * shape.detJ * wp.getWeight();
     }