diff --git a/NumLib/Fem/Integration/IntegrationGaussPrism.h b/NumLib/Fem/Integration/IntegrationGaussPrism.h
index d634e46a19bfb6cb5fc60ba4eec11fe5da06a02e..181932ac8aa424e5bd4ebd77625e2bc0ca46fcb9 100644
--- a/NumLib/Fem/Integration/IntegrationGaussPrism.h
+++ b/NumLib/Fem/Integration/IntegrationGaussPrism.h
@@ -29,24 +29,24 @@ public:
      *
      * @param order     integration order (default 2)
      */
-    explicit IntegrationGaussPrism(std::size_t order = 2)
+    explicit IntegrationGaussPrism(unsigned order = 2)
     : _order(2), _n_sampl_pt(0)
     {
         this->setIntegrationOrder(order);
     }
 
     /// Change the integration order.
-    void setIntegrationOrder(std::size_t /*order*/)
+    void setIntegrationOrder(unsigned /*order*/)
     {
         _order = 2; // fixed
         _n_sampl_pt = getNumberOfPoints(_order);
     }
 
     /// return current integration order.
-    std::size_t getIntegrationOrder() const {return _order;}
+    unsigned getIntegrationOrder() const {return _order;}
 
     /// return the number of sampling points
-    std::size_t getNumberOfPoints() const {return _n_sampl_pt;}
+    unsigned getNumberOfPoints() const {return _n_sampl_pt;}
 
     /**
      * get coordinates of a integration point
@@ -54,7 +54,7 @@ public:
      * @param igp      The integration point index
      * @return a weighted point
      */
-    WeightedPoint getWeightedPoint(std::size_t igp)
+    WeightedPoint getWeightedPoint(unsigned igp) const
     {
         return getWeightedPoint(getIntegrationOrder(), igp);
     }
@@ -67,7 +67,7 @@ public:
      * @return weight
      */
     static WeightedPoint
-    getWeightedPoint(std::size_t /*order*/, std::size_t igp)
+    getWeightedPoint(unsigned /*order*/, unsigned igp)
     {
         const unsigned gp_r = igp % 3;
         const unsigned gp_t = (unsigned)(igp/3);
@@ -81,7 +81,7 @@ public:
 
     template <typename Method>
     static WeightedPoint
-    getWeightedPoint(std::size_t igp)
+    getWeightedPoint(unsigned igp)
     {
         return WeightedPoint(Method::X[igp], Method::W[igp]);
     }
@@ -93,16 +93,16 @@ public:
      * @param order    the number of integration points
      * @return the number of points
      */
-    static std::size_t
-    getNumberOfPoints(std::size_t order)
+    static unsigned
+    getNumberOfPoints(unsigned order)
     {
         if (order==2) return 6;
         return 0;
     }
 
 private:
-    std::size_t _order;
-    std::size_t _n_sampl_pt;
+    unsigned _order;
+    unsigned _n_sampl_pt;
 };
 
 }
diff --git a/NumLib/Fem/Integration/IntegrationGaussPyramid.h b/NumLib/Fem/Integration/IntegrationGaussPyramid.h
index 0cc01a1c2383b93ef770d0f351ff910726ae7f8e..9fd9925875d9c651bf3bc49dbd3151b10849bf63 100644
--- a/NumLib/Fem/Integration/IntegrationGaussPyramid.h
+++ b/NumLib/Fem/Integration/IntegrationGaussPyramid.h
@@ -28,24 +28,24 @@ public:
      *
      * @param order     integration order (default 2)
      */
-    explicit IntegrationGaussPyramid(std::size_t order = 2)
+    explicit IntegrationGaussPyramid(unsigned order = 2)
     : _order(order), _n_sampl_pt(0)
     {
         this->setIntegrationOrder(order);
     }
 
     /// Change the integration order.
-    void setIntegrationOrder(std::size_t order)
+    void setIntegrationOrder(unsigned order)
     {
         _order = order;
         _n_sampl_pt = getNumberOfPoints(_order);
     }
 
     /// return current integration order.
-    std::size_t getIntegrationOrder() const {return _order;}
+    unsigned getIntegrationOrder() const {return _order;}
 
     /// return the number of sampling points
-    std::size_t getNumberOfPoints() const {return _n_sampl_pt;}
+    unsigned getNumberOfPoints() const {return _n_sampl_pt;}
 
     /**
      * get coordinates of a integration point
@@ -53,7 +53,7 @@ public:
      * @param igp      The integration point index
      * @return a weighted point
      */
-    WeightedPoint getWeightedPoint(std::size_t igp)
+    WeightedPoint getWeightedPoint(unsigned igp) const
     {
         return getWeightedPoint(getIntegrationOrder(), igp);
     }
@@ -66,7 +66,7 @@ public:
      * @return weight
      */
     static WeightedPoint
-    getWeightedPoint(std::size_t order, std::size_t igp)
+    getWeightedPoint(unsigned order, unsigned igp)
     {
         switch (order)
         {
@@ -79,7 +79,7 @@ public:
 
     template <typename Method>
     static WeightedPoint
-    getWeightedPoint(std::size_t igp)
+    getWeightedPoint(unsigned igp)
     {
         return WeightedPoint(Method::X[igp], Method::W[igp]);
     }
@@ -91,8 +91,8 @@ public:
      * @param order    the number of integration points
      * @return the number of points
      */
-    static std::size_t
-    getNumberOfPoints(std::size_t order)
+    static unsigned
+    getNumberOfPoints(unsigned order)
     {
         switch (order)
         {
@@ -104,8 +104,8 @@ public:
     }
 
 private:
-    std::size_t _order;
-    std::size_t _n_sampl_pt;
+    unsigned _order;
+    unsigned _n_sampl_pt;
 };
 
 }
diff --git a/NumLib/Fem/Integration/IntegrationGaussRegular-impl.h b/NumLib/Fem/Integration/IntegrationGaussRegular-impl.h
index b1200cd891cc41acd2fe0c8470cf31e4e700cf72..610f5f16cf9cc5e490be6a849dbaf55f9b503911 100644
--- a/NumLib/Fem/Integration/IntegrationGaussRegular-impl.h
+++ b/NumLib/Fem/Integration/IntegrationGaussRegular-impl.h
@@ -16,46 +16,46 @@ namespace NumLib
 {
 
 template <>
-inline std::array<std::size_t, 1>
-IntegrationGaussRegular<1>::getPositionIndices(std::size_t /*order*/, std::size_t igp)
+inline std::array<unsigned, 1>
+IntegrationGaussRegular<1>::getPositionIndices(unsigned /*order*/, unsigned igp)
 {
-    std::array<std::size_t, 1> result;
+    std::array<unsigned, 1> result;
     result[0] = igp;
     return result;
 }
 
 template <>
-inline std::array<std::size_t, 2>
-IntegrationGaussRegular<2>::getPositionIndices(std::size_t order, std::size_t igp)
+inline std::array<unsigned, 2>
+IntegrationGaussRegular<2>::getPositionIndices(unsigned order, unsigned igp)
 {
     assert(igp < order*order);
-    std::array<std::size_t, 2> result;
+    std::array<unsigned, 2> result;
     result[0] = igp / order;
     result[1] = igp % order;
     return result;
 }
 
 template <>
-inline std::array<std::size_t, 3>
-IntegrationGaussRegular<3>::getPositionIndices(std::size_t order, std::size_t igp)
+inline std::array<unsigned, 3>
+IntegrationGaussRegular<3>::getPositionIndices(unsigned order, unsigned igp)
 {
     assert(igp < order*order*order);
-    std::size_t const gp_r = igp / (order * order);
-    std::size_t const gp_s = igp % (order * order);
-    std::array<std::size_t, 3> result;
+    unsigned const gp_r = igp / (order * order);
+    unsigned const gp_s = igp % (order * order);
+    std::array<unsigned, 3> result;
     result[0] = gp_r;
     result[1] = gp_s / order;
     result[2] = gp_s % order;
     return result;
 }
 
-template <std::size_t N_DIM>
-inline
-MathLib::TemplateWeightedPoint<double,double,N_DIM>
-IntegrationGaussRegular<N_DIM>::getWeightedPoint(std::size_t order, std::size_t igp)
+template <unsigned N_DIM>
+inline MathLib::TemplateWeightedPoint<double, double, N_DIM>
+IntegrationGaussRegular<N_DIM>::getWeightedPoint(unsigned order,
+                                                 unsigned igp)
 {
     assert(igp < std::pow(order, N_DIM));
-    std::array<std::size_t, N_DIM> const pos = getPositionIndices(order, igp);
+    std::array<unsigned, N_DIM> const pos = getPositionIndices(order, igp);
 
     switch (order)
     {
@@ -68,11 +68,11 @@ IntegrationGaussRegular<N_DIM>::getWeightedPoint(std::size_t order, std::size_t
     return MathLib::TemplateWeightedPoint<double, double, N_DIM>(std::array<double, N_DIM>(), 0);
 }
 
-template <std::size_t N_DIM>
+template <unsigned N_DIM>
 template <typename Method>
-inline
-MathLib::TemplateWeightedPoint<double, double, N_DIM>
-IntegrationGaussRegular<N_DIM>::getWeightedPoint(std::array<std::size_t, N_DIM> const& pos)
+inline MathLib::TemplateWeightedPoint<double, double, N_DIM>
+IntegrationGaussRegular<N_DIM>::getWeightedPoint(
+    std::array<unsigned, N_DIM> const& pos)
 {
     std::array<double, N_DIM> coords;
     double weight = 1;
diff --git a/NumLib/Fem/Integration/IntegrationGaussRegular.h b/NumLib/Fem/Integration/IntegrationGaussRegular.h
index 6840d3edcf7ea287bfacffba6d13f108467b55b2..4d349297f86ce77133be9ff04f245b7add36e683 100644
--- a/NumLib/Fem/Integration/IntegrationGaussRegular.h
+++ b/NumLib/Fem/Integration/IntegrationGaussRegular.h
@@ -26,7 +26,7 @@ namespace NumLib
 /// Gauss quadrature rule for regular shape elements: line, quad and hex.
 ///
 /// \tparam N_DIM    Spatial dimension
-template <std::size_t N_DIM>
+template <unsigned N_DIM>
 class IntegrationGaussRegular
 {
     typedef typename MathLib::TemplateWeightedPoint<double, double, N_DIM>
@@ -36,31 +36,31 @@ public:
     /// order.
     ///
     /// @param order     integration order (default 2)
-    explicit IntegrationGaussRegular(std::size_t order = 2)
+    explicit IntegrationGaussRegular(unsigned order = 2)
     : _order(order), _n_sampl_pt(0)
     {
         this->setIntegrationOrder(order);
     }
 
     /// Change the integration order.
-    void setIntegrationOrder(std::size_t order)
+    void setIntegrationOrder(unsigned order)
     {
-        this->_n_sampl_pt = static_cast<std::size_t>(std::pow(order, N_DIM));
+        this->_n_sampl_pt = static_cast<unsigned>(std::pow(order, N_DIM));
         this->_order = order;
     }
 
     /// return current integration order.
-    std::size_t getIntegrationOrder() const {return _order;}
+    unsigned getIntegrationOrder() const {return _order;}
 
     /// return the number of sampling points
-    std::size_t getNumberOfPoints() const {return _n_sampl_pt;}
+    unsigned getNumberOfPoints() const {return _n_sampl_pt;}
 
     /// Get coordinates of the integration point.
     ///
     /// @param igp       The integration point index
     /// @return a weighted point
     WeightedPoint
-    getWeightedPoint(std::size_t igp) const
+    getWeightedPoint(unsigned igp) const
     {
         return getWeightedPoint(getIntegrationOrder(), igp);
     }
@@ -70,7 +70,7 @@ public:
     /// @param order     The number of integration points
     /// @param igp       The integration point index
     /// @return  a tuple of position indexes
-    static std::array<std::size_t, N_DIM> getPositionIndices(std::size_t order, std::size_t igp);
+    static std::array<unsigned, N_DIM> getPositionIndices(unsigned order, unsigned igp);
 
     /// Get coordinates of the integration point.
     ///
@@ -78,7 +78,7 @@ public:
     /// @param igp       The integration point index
     /// @return a weighted point
     static WeightedPoint
-    getWeightedPoint(std::size_t order, std::size_t igp);
+    getWeightedPoint(unsigned order, unsigned igp);
 
 private:
     /// Computes weighted point using given integration method.
@@ -88,11 +88,11 @@ private:
     template <typename Method>
     static
     WeightedPoint
-    getWeightedPoint(std::array<std::size_t, N_DIM> const& pos);
+    getWeightedPoint(std::array<unsigned, N_DIM> const& pos);
 
 private:
-    std::size_t _order;
-    std::size_t _n_sampl_pt;
+    unsigned _order;
+    unsigned _n_sampl_pt;
 };
 
 } // NumLib
diff --git a/NumLib/Fem/Integration/IntegrationGaussTet.h b/NumLib/Fem/Integration/IntegrationGaussTet.h
index 9a2a0fe59df717e197af0e99b2fd45e5583bf54d..17c97d9a49e4a335de3f67536ff42aef0b1b4039 100644
--- a/NumLib/Fem/Integration/IntegrationGaussTet.h
+++ b/NumLib/Fem/Integration/IntegrationGaussTet.h
@@ -28,24 +28,24 @@ public:
      *
      * @param order     integration order (default 2)
      */
-    explicit IntegrationGaussTet(std::size_t order = 2)
+    explicit IntegrationGaussTet(unsigned order = 2)
     : _order(order), _n_sampl_pt(0)
     {
         this->setIntegrationOrder(order);
     }
 
     /// Change the integration order.
-    void setIntegrationOrder(std::size_t order)
+    void setIntegrationOrder(unsigned order)
     {
         _order = order;
         _n_sampl_pt = getNumberOfPoints(order);
     }
 
     /// return current integration order.
-    std::size_t getIntegrationOrder() const {return _order;}
+    unsigned getIntegrationOrder() const {return _order;}
 
     /// return the number of sampling points
-    std::size_t getNumberOfPoints() const {return _n_sampl_pt;}
+    unsigned getNumberOfPoints() const {return _n_sampl_pt;}
 
     /**
      * get coordinates of a integration point
@@ -53,7 +53,7 @@ public:
      * @param igp      The integration point index
      * @return a weighted point
      */
-    WeightedPoint getWeightedPoint(std::size_t igp)
+    WeightedPoint getWeightedPoint(unsigned igp) const
     {
         return getWeightedPoint(getIntegrationOrder(), igp);
     }
@@ -66,7 +66,7 @@ public:
      * @return weight
      */
     static WeightedPoint
-    getWeightedPoint(std::size_t order, std::size_t igp)
+    getWeightedPoint(unsigned order, unsigned igp)
     {
         switch (order)
         {
@@ -79,7 +79,7 @@ public:
 
     template <typename Method>
     static WeightedPoint
-    getWeightedPoint(std::size_t igp)
+    getWeightedPoint(unsigned igp)
     {
         return WeightedPoint(Method::X[igp], Method::W[igp]);
     }
@@ -91,8 +91,8 @@ public:
      * @param order    the number of integration points
      * @return the number of points
      */
-    static std::size_t
-    getNumberOfPoints(std::size_t order)
+    static unsigned
+    getNumberOfPoints(unsigned order)
     {
         switch (order)
         {
@@ -104,8 +104,8 @@ public:
     }
 
 private:
-    std::size_t _order;
-    std::size_t _n_sampl_pt;
+    unsigned _order;
+    unsigned _n_sampl_pt;
 };
 
 }
diff --git a/NumLib/Fem/Integration/IntegrationGaussTri.h b/NumLib/Fem/Integration/IntegrationGaussTri.h
index 0d3ec3cc71827c6c03303365e597fd64cc793c57..cefd971ea512b5bd1b35aa42ee35373fffd2d423 100644
--- a/NumLib/Fem/Integration/IntegrationGaussTri.h
+++ b/NumLib/Fem/Integration/IntegrationGaussTri.h
@@ -44,24 +44,24 @@ public:
      *
      * @param order     integration order (default 2)
      */
-    explicit IntegrationGaussTri(std::size_t order = 2)
+    explicit IntegrationGaussTri(unsigned order = 2)
     : _order(order), _n_sampl_pt(0)
     {
         this->setIntegrationOrder(order);
     }
 
     /// Change the integration order.
-    void setIntegrationOrder(std::size_t order)
+    void setIntegrationOrder(unsigned order)
     {
         _order = order;
         _n_sampl_pt = getNumberOfPoints(order);
     }
 
     /// return current integration order.
-    std::size_t getIntegrationOrder() const {return _order;}
+    unsigned getIntegrationOrder() const {return _order;}
 
     /// return the number of sampling points
-    std::size_t getNumberOfPoints() const {return _n_sampl_pt;}
+    unsigned getNumberOfPoints() const {return _n_sampl_pt;}
 
     /**
      * get coordinates of a integration point
@@ -69,7 +69,7 @@ public:
      * @param igp      The integration point index
      * @return a weighted point
      */
-    WeightedPoint getWeightedPoint(std::size_t igp)
+    WeightedPoint getWeightedPoint(unsigned igp) const
     {
         return getWeightedPoint(getIntegrationOrder(), igp);
     }
@@ -82,7 +82,7 @@ public:
      * @return weight
      */
     static WeightedPoint
-    getWeightedPoint(std::size_t order, std::size_t igp)
+    getWeightedPoint(unsigned order, unsigned igp)
     {
         switch (order)
         {
@@ -95,7 +95,7 @@ public:
 
     template <typename Method>
     static WeightedPoint
-    getWeightedPoint(std::size_t igp)
+    getWeightedPoint(unsigned igp)
     {
         return WeightedPoint(Method::X[igp], 0.5 * Method::W[igp]);
     }
@@ -107,8 +107,8 @@ public:
      * @param order    the number of integration points
      * @return the number of points
      */
-    static std::size_t
-    getNumberOfPoints(std::size_t order)
+    static unsigned
+    getNumberOfPoints(unsigned order)
     {
         switch (order)
         {
@@ -120,8 +120,8 @@ public:
     }
 
 private:
-    std::size_t _order;
-    std::size_t _n_sampl_pt;
+    unsigned _order;
+    unsigned _n_sampl_pt;
 };
 
 }
diff --git a/NumLib/Fem/Integration/IntegrationPoint.h b/NumLib/Fem/Integration/IntegrationPoint.h
index ee2c3b5af79f9a317e5f1d7408ee4a358e79e0cf..7f8b8790350f47e21ea415abbabde7ef20edb051 100644
--- a/NumLib/Fem/Integration/IntegrationPoint.h
+++ b/NumLib/Fem/Integration/IntegrationPoint.h
@@ -22,23 +22,23 @@ class IntegrationPoint
 
 public:
     /// IntegrationPoint constructor for given order.
-    explicit IntegrationPoint(std::size_t /* order */)
+    explicit IntegrationPoint(unsigned /* order */)
     {
     }
 
     /// Change the integration order.
-    void setIntegrationOrder(std::size_t /* order */)
+    void setIntegrationOrder(unsigned /* order */)
     {
     }
 
     /// Return current integration order.
-    std::size_t getIntegrationOrder() const
+    unsigned getIntegrationOrder() const
     {
         return 0;
     }
 
     /// Return the number of sampling points.
-    std::size_t getNumberOfPoints() const
+    unsigned getNumberOfPoints() const
     {
         return 1;
     }
@@ -47,7 +47,7 @@ public:
     ///
     /// \param igp      The integration point index.
     /// \return a weighted point.
-    WeightedPoint getWeightedPoint(std::size_t igp)
+    WeightedPoint getWeightedPoint(unsigned igp) const
     {
         return getWeightedPoint(getIntegrationOrder(), igp);
     }
@@ -57,14 +57,14 @@ public:
     /// \param order    the number of integration points.
     /// \param igp      the sampling point id.
     /// \return weight
-    static WeightedPoint getWeightedPoint(std::size_t /* order */,
-                                          std::size_t /*igp*/)
+    static WeightedPoint getWeightedPoint(unsigned /* order */,
+                                          unsigned /*igp*/)
     {
         return WeightedPoint({{1}}, 1);
     }
 
     template <typename Method>
-    static WeightedPoint getWeightedPoint(std::size_t /*igp*/)
+    static WeightedPoint getWeightedPoint(unsigned /*igp*/)
     {
         return WeightedPoint({{1}}, 1);
     }
@@ -73,7 +73,7 @@ public:
     ///
     /// \param order    the number of integration points
     /// \return the number of points.
-    static std::size_t getNumberOfPoints(std::size_t /* order */)
+    static unsigned getNumberOfPoints(unsigned /* order */)
     {
         return 1;
     }
diff --git a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h
index 544e4236b2a8ce167da765870ed796014d3d2ed8..bac5f8e35a495dfcfe98319dbfa54f10bac69d96 100644
--- a/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/GenericNaturalBoundaryConditionLocalAssembler.h
@@ -39,17 +39,16 @@ protected:
 public:
     GenericNaturalBoundaryConditionLocalAssembler(
         MeshLib::Element const& e, unsigned const integration_order)
-        : _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
+        : _integration_method(integration_order),
+          _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
                                             IntegrationMethod, GlobalDim>(
-              e, integration_order)),
-          _integration_order(integration_order)
+              e, _integration_method))
     {
     }
 
 protected:
-    std::vector<typename ShapeMatricesType::ShapeMatrices> const
-        _shape_matrices;
-    unsigned const _integration_order;
+    IntegrationMethod const _integration_method;
+    std::vector<typename ShapeMatricesType::ShapeMatrices> const _shape_matrices;
 };
 
 }  // ProcessLib
diff --git a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
index 663e79d582b231b38f335ddf2600ee5fc1078b1f..fd9a5d43e9beff7eff9f6e0d5623d366a7e894b1 100644
--- a/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/NeumannBoundaryConditionLocalAssembler.h
@@ -46,18 +46,17 @@ public:
     {
         _local_rhs.setZero();
 
-        IntegrationMethod integration_method(Base::_integration_order);
-        std::size_t const n_integration_points =
-            integration_method.getNumberOfPoints();
+        unsigned const n_integration_points =
+            Base::_integration_method.getNumberOfPoints();
 
         SpatialPosition pos;
         pos.setElementID(id);
 
-        for (std::size_t ip(0); ip < n_integration_points; ip++)
+        for (unsigned ip = 0; ip < n_integration_points; ip++)
         {
             pos.setIntegrationPoint(ip);
             auto const& sm = Base::_shape_matrices[ip];
-            auto const& wp = integration_method.getWeightedPoint(ip);
+            auto const& wp = Base::_integration_method.getWeightedPoint(ip);
             _local_rhs.noalias() += sm.N * _neumann_bc_parameter(t, pos)[0] *
                                     sm.detJ * wp.getWeight();
         }
diff --git a/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h b/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
index cdeb1089c9f1e19660acb5bf699203089b44e8f2..a6322168e37758d3a058cdf6684af416ff7f6497 100644
--- a/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
+++ b/ProcessLib/BoundaryCondition/RobinBoundaryConditionLocalAssembler.h
@@ -51,18 +51,17 @@ public:
         _local_K.setZero();
         _local_rhs.setZero();
 
-        IntegrationMethod integration_method(Base::_integration_order);
-        std::size_t const n_integration_points =
-            integration_method.getNumberOfPoints();
+        unsigned const n_integration_points =
+            Base::_integration_method.getNumberOfPoints();
 
         SpatialPosition pos;
         pos.setElementID(id);
 
-        for (std::size_t ip = 0; ip < n_integration_points; ++ip)
+        for (unsigned ip = 0; ip < n_integration_points; ++ip)
         {
             pos.setIntegrationPoint(ip);
             auto const& sm = Base::_shape_matrices[ip];
-            auto const& wp = integration_method.getWeightedPoint(ip);
+            auto const& wp = Base::_integration_method.getWeightedPoint(ip);
 
             double const alpha = _data.alpha(t, pos)[0];
             double const u_0 = _data.u_0(t, pos)[0];
diff --git a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
index 870bfe7482a6e970a77e67f1088ccc758c33c1a2..d0a87e982d587b3907b649e97ac211140e02187c 100644
--- a/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
+++ b/ProcessLib/GroundwaterFlow/GroundwaterFlowFEM.h
@@ -67,14 +67,15 @@ public:
                        std::size_t const local_matrix_size,
                        unsigned const integration_order,
                        GroundwaterFlowProcessData const& process_data)
-        : _element(element)
-        , _shape_matrices(
-              initShapeMatrices<ShapeFunction, ShapeMatricesType, IntegrationMethod, GlobalDim>(
-                  element, integration_order))
-        , _process_data(process_data)
-        , _localA(local_matrix_size, local_matrix_size) // TODO narrowing conversion
-        , _localRhs(local_matrix_size)
-        , _integration_order(integration_order)
+        : _element(element),
+          _process_data(process_data),
+          _localA(local_matrix_size,
+                  local_matrix_size),  // TODO narrowing conversion
+          _localRhs(local_matrix_size),
+          _integration_method(integration_order),
+          _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
+                                            IntegrationMethod, GlobalDim>(
+              element, _integration_method))
     {
         // This assertion is valid only if all nodal d.o.f. use the same shape matrices.
         assert(local_matrix_size == ShapeFunction::NPOINTS * NUM_NODAL_DOF);
@@ -88,8 +89,8 @@ public:
         _localA.setZero();
         _localRhs.setZero();
 
-        IntegrationMethod integration_method(_integration_order);
-        unsigned const n_integration_points = integration_method.getNumberOfPoints();
+        unsigned const n_integration_points =
+            _integration_method.getNumberOfPoints();
 
         SpatialPosition pos;
         pos.setElementID(_element.getID());
@@ -98,7 +99,7 @@ public:
         {
             pos.setIntegrationPoint(ip);
             auto const& sm = _shape_matrices[ip];
-            auto const& wp = integration_method.getWeightedPoint(ip);
+            auto const& wp = _integration_method.getWeightedPoint(ip);
             auto const k = _process_data.hydraulic_conductivity(t, pos)[0];
 
             _localA.noalias() += sm.dNdx.transpose() * k * sm.dNdx *
@@ -151,13 +152,13 @@ public:
 
 private:
     MeshLib::Element const& _element;
-    std::vector<ShapeMatrices> _shape_matrices;
     GroundwaterFlowProcessData const& _process_data;
 
     NodalMatrixType _localA;
     NodalVectorType _localRhs;
 
-    unsigned const _integration_order;
+    IntegrationMethod const _integration_method;
+    std::vector<ShapeMatrices> _shape_matrices;
 
     std::vector<std::vector<double>> _darcy_velocities
         = std::vector<std::vector<double>>(
diff --git a/ProcessLib/TES/TESLocalAssembler-impl.h b/ProcessLib/TES/TESLocalAssembler-impl.h
index bb3fc62238571c38f95de1457f5c2fb39341ae9d..a99009406a9cc943684af4d6b458f6094a797e7f 100644
--- a/ProcessLib/TES/TESLocalAssembler-impl.h
+++ b/ProcessLib/TES/TESLocalAssembler-impl.h
@@ -106,15 +106,15 @@ namespace TES
 {
 template <typename ShapeFunction_, typename IntegrationMethod_,
           unsigned GlobalDim>
-TESLocalAssembler<
-    ShapeFunction_, IntegrationMethod_,
-    GlobalDim>::TESLocalAssembler(MeshLib::Element const& e,
-                                  std::size_t const /*local_matrix_size*/,
-                                  unsigned const integration_order,
-                                  AssemblyParams const& asm_params)
-    : _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
+TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::
+    TESLocalAssembler(MeshLib::Element const& e,
+                      std::size_t const /*local_matrix_size*/,
+                      unsigned const integration_order,
+                      AssemblyParams const& asm_params)
+    : _integration_method(integration_order),
+      _shape_matrices(initShapeMatrices<ShapeFunction, ShapeMatricesType,
                                         IntegrationMethod_, GlobalDim>(
-          e, integration_order)),
+          e, _integration_method)),
       _d(asm_params,
          // TODO narrowing conversion
          static_cast<const unsigned>(
@@ -125,8 +125,7 @@ TESLocalAssembler<
                ShapeFunction::NPOINTS * NODAL_DOF),
       _local_K(ShapeFunction::NPOINTS * NODAL_DOF,
                ShapeFunction::NPOINTS * NODAL_DOF),
-      _local_b(ShapeFunction::NPOINTS * NODAL_DOF),
-      _integration_order(integration_order)
+      _local_b(ShapeFunction::NPOINTS * NODAL_DOF)
 {
 }
 
@@ -142,15 +141,15 @@ void TESLocalAssembler<ShapeFunction_, IntegrationMethod_, GlobalDim>::
     _local_K.setZero();
     _local_b.setZero();
 
-    IntegrationMethod_ integration_method(_integration_order);
-    unsigned const n_integration_points = integration_method.getNumberOfPoints();
+    unsigned const n_integration_points =
+        _integration_method.getNumberOfPoints();
 
     _d.preEachAssemble();
 
-    for (std::size_t ip(0); ip < n_integration_points; ip++)
+    for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
         auto const& sm = _shape_matrices[ip];
-        auto const& wp = integration_method.getWeightedPoint(ip);
+        auto const& wp = _integration_method.getWeightedPoint(ip);
         auto const weight = wp.getWeight();
 
         _d.assembleIntegrationPoint(ip, local_x, sm.N, sm.dNdx, sm.J, sm.detJ,
diff --git a/ProcessLib/TES/TESLocalAssembler.h b/ProcessLib/TES/TESLocalAssembler.h
index 1826b2c6f4cc2fec5448b5559c80aa37d3c914b5..29fd69c792783a08bc76aa1ce55c94d377493c34 100644
--- a/ProcessLib/TES/TESLocalAssembler.h
+++ b/ProcessLib/TES/TESLocalAssembler.h
@@ -107,6 +107,8 @@ public:
     std::vector<double> const& getIntPtDarcyVelocityZ(
         std::vector<double>& /*cache*/) const override;
 private:
+    IntegrationMethod_ const _integration_method;
+
     std::vector<ShapeMatrices> _shape_matrices;
 
     using LAT = LocalAssemblerTraits<ShapeMatricesType, ShapeFunction::NPOINTS,
@@ -128,9 +130,6 @@ private:
     NodalMatrixType _local_M;
     NodalMatrixType _local_K;
     NodalVectorType _local_b;
-
-    // TODO Use the value from Process
-    unsigned const _integration_order;
 };
 
 }  // namespace TES
diff --git a/ProcessLib/Utils/InitShapeMatrices.h b/ProcessLib/Utils/InitShapeMatrices.h
index 08e1f06d5226af793ea2ab623a3de16f7c65d907..6fbf385775f96d6e54298ff401d3a926b6ec62c7 100644
--- a/ProcessLib/Utils/InitShapeMatrices.h
+++ b/ProcessLib/Utils/InitShapeMatrices.h
@@ -17,12 +17,10 @@
 
 namespace ProcessLib
 {
-
-
-template<typename ShapeFunction, typename ShapeMatricesType, typename IntegrationMethod,
-         unsigned GlobalDim>
-std::vector<typename ShapeMatricesType::ShapeMatrices>
-initShapeMatrices(MeshLib::Element const& e, unsigned integration_order)
+template <typename ShapeFunction, typename ShapeMatricesType,
+          typename IntegrationMethod, unsigned GlobalDim>
+std::vector<typename ShapeMatricesType::ShapeMatrices> initShapeMatrices(
+    MeshLib::Element const& e, IntegrationMethod const& integration_method)
 {
     std::vector<typename ShapeMatricesType::ShapeMatrices> shape_matrices;
 
@@ -31,11 +29,10 @@ initShapeMatrices(MeshLib::Element const& e, unsigned integration_order)
 
     FemType fe(*static_cast<const typename ShapeFunction::MeshElement*>(&e));
 
-    IntegrationMethod integration_method(integration_order);
-    std::size_t const n_integration_points = integration_method.getNumberOfPoints();
+    unsigned const n_integration_points = integration_method.getNumberOfPoints();
 
     shape_matrices.reserve(n_integration_points);
-    for (std::size_t 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);
         fe.computeShapeFunctions(
diff --git a/Tests/NumLib/TestExtrapolation.cpp b/Tests/NumLib/TestExtrapolation.cpp
index 2d3685d4908b21a4aa922ab770fcdab80a055c06..aee0e030dd1ce15a3efc358366950b17189babb7 100644
--- a/Tests/NumLib/TestExtrapolation.cpp
+++ b/Tests/NumLib/TestExtrapolation.cpp
@@ -77,8 +77,8 @@ public:
         : _shape_matrices(
               ProcessLib::initShapeMatrices<ShapeFunction, ShapeMatricesType,
                                             IntegrationMethod, GlobalDim>(
-                  e, integration_order))
-        , _int_pt_values(_shape_matrices.size())
+                  e, IntegrationMethod{integration_order})),
+          _int_pt_values(_shape_matrices.size())
     {
     }