diff --git a/Applications/FileIO/GocadIO/GocadSGridReader.cpp b/Applications/FileIO/GocadIO/GocadSGridReader.cpp
index d33c8ec4ec8199bedb74d2c9c2682231d2ecc288..af78fa944e81fad425fc1d0ffb773ea5033ebbfa 100644
--- a/Applications/FileIO/GocadIO/GocadSGridReader.cpp
+++ b/Applications/FileIO/GocadIO/GocadSGridReader.cpp
@@ -715,7 +715,7 @@ void GocadSGridReader::applySplitInformation(
 
 void GocadSGridReader::modifyElement(MeshLib::Element* hex,
                                      MeshLib::Node const* node2sub,
-                                     MeshLib::Node* substitute_node) const
+                                     MeshLib::Node* substitute_node)
 {
     // get the node pointers of the cell
     MeshLib::Node* const* hex_nodes(hex->getNodes());
diff --git a/Applications/FileIO/GocadIO/GocadSGridReader.h b/Applications/FileIO/GocadIO/GocadSGridReader.h
index 726599896b7f9e9cd85bfc36c4ac8577f3c23579..3e2715450c8bddac0759a3ac6f69c524bcbae3a8 100644
--- a/Applications/FileIO/GocadIO/GocadSGridReader.h
+++ b/Applications/FileIO/GocadIO/GocadSGridReader.h
@@ -74,8 +74,9 @@ private:
     void applySplitInformation(
         std::vector<MeshLib::Node*>& nodes,
         std::vector<MeshLib::Element*> const& elements) const;
-    void modifyElement(MeshLib::Element* hex, MeshLib::Node const* node2sub,
-                       MeshLib::Node* substitute_node) const;
+    static void modifyElement(MeshLib::Element* hex,
+                              MeshLib::Node const* node2sub,
+                              MeshLib::Node* substitute_node);
 
     void addFaceSetQuad(
         GocadNode* face_set_node, std::size_t face_set_number,
diff --git a/Applications/FileIO/SHPInterface.h b/Applications/FileIO/SHPInterface.h
index 885229be52ba532c49716541abf094e56932ae24..f9f63bbd69610ed4f8e50ad0b335654cbf45ecce 100644
--- a/Applications/FileIO/SHPInterface.h
+++ b/Applications/FileIO/SHPInterface.h
@@ -61,7 +61,8 @@ public:
     }
 
     /// Reads the header of the shape file.
-    bool readSHPInfo(const std::string &filename, int &shapeType, int &numberOfEntities);
+    static bool readSHPInfo(const std::string& filename, int& shapeType,
+                            int& numberOfEntities);
 
     /// Reads data from the shape file.
     void readSHPFile(const std::string& filename, OGSType choice,
diff --git a/Applications/FileIO/TetGenInterface.cpp b/Applications/FileIO/TetGenInterface.cpp
index 39ae415b1ae1d0cc6e3e20e80461e52dae7b684d..09ff26c2ad0ee8661660140ca8a96f31ff03fc17 100644
--- a/Applications/FileIO/TetGenInterface.cpp
+++ b/Applications/FileIO/TetGenInterface.cpp
@@ -293,11 +293,11 @@ bool TetGenInterface::readNodesFromStream (std::ifstream &ins,
     return false;
 }
 
-bool TetGenInterface::parseNodesFileHeader(std::string &line,
-                                           std::size_t &n_nodes,
-                                           std::size_t &dim,
-                                           std::size_t &n_attributes,
-                                           bool &boundary_markers) const
+bool TetGenInterface::parseNodesFileHeader(std::string& line,
+                                           std::size_t& n_nodes,
+                                           std::size_t& dim,
+                                           std::size_t& n_attributes,
+                                           bool& boundary_markers)
 {
     std::list<std::string> pnt_header = BaseLib::splitString(line, ' ');
     if (pnt_header.empty())
@@ -437,10 +437,10 @@ bool TetGenInterface::readElementsFromStream(
     return false;
 }
 
-bool TetGenInterface::parseElementsFileHeader(std::string &line,
+bool TetGenInterface::parseElementsFileHeader(std::string& line,
                                               std::size_t& n_tets,
                                               std::size_t& n_nodes_per_tet,
-                                              bool& region_attribute) const
+                                              bool& region_attribute)
 {
     std::size_t pos_beg;
     std::size_t pos_end;
@@ -574,10 +574,11 @@ bool TetGenInterface::parseElements(std::ifstream& ins,
     return true;
 }
 
-bool TetGenInterface::writeTetGenSmesh(const std::string &file_name,
-                                       const GeoLib::GEOObjects &geo_objects,
-                                       const std::string &geo_name,
-                                       const std::vector<GeoLib::Point> &attribute_points) const
+bool TetGenInterface::writeTetGenSmesh(
+    const std::string& file_name,
+    const GeoLib::GEOObjects& geo_objects,
+    const std::string& geo_name,
+    const std::vector<GeoLib::Point>& attribute_points)
 {
     std::vector<GeoLib::Point*> const*const points = geo_objects.getPointVec(geo_name);
     std::vector<GeoLib::Surface*> const*const surfaces = geo_objects.getSurfaceVec(geo_name);
@@ -793,7 +794,10 @@ void TetGenInterface::write3dElements(std::ofstream &out,
     out.seekp(after_elems_pos);
 }
 
-void TetGenInterface::writeElementToFacets(std::ofstream &out, const MeshLib::Element &element, unsigned &element_count, std::string const& matId) const
+void TetGenInterface::writeElementToFacets(std::ofstream& out,
+                                           const MeshLib::Element& element,
+                                           unsigned& element_count,
+                                           std::string const& matId)
 {
     element_count++;
     if (element.getGeomType() == MeshLib::MeshElemType::TRIANGLE)
diff --git a/Applications/FileIO/TetGenInterface.h b/Applications/FileIO/TetGenInterface.h
index cab781b13ce96a09ed015a94fb12acaccc4f5b8b..690fece898f123f4f24a9cb5cb06760ff9b86b9b 100644
--- a/Applications/FileIO/TetGenInterface.h
+++ b/Applications/FileIO/TetGenInterface.h
@@ -68,10 +68,11 @@ public:
      * @param attribute_points  attribute points containing material IDs (if the vector is empty no attributes are written).
      * @return returns true on success and false otherwise.
      */
-    bool writeTetGenSmesh(const std::string &file_name,
-                          const GeoLib::GEOObjects &geo_objects,
-                          const std::string &geo_name,
-                          const std::vector<GeoLib::Point> &attribute_points) const;
+    static bool writeTetGenSmesh(
+        const std::string& file_name,
+        const GeoLib::GEOObjects& geo_objects,
+        const std::string& geo_name,
+        const std::vector<GeoLib::Point>& attribute_points);
 
     /**
      * Writes the geometry of a given name to TetGen smesh-file.
@@ -121,11 +122,11 @@ private:
      * @param boundary_markers  have the nodes boundary information (output)
      * @return true, if the file header is read, false if the method detects an error
      */
-    bool parseNodesFileHeader(std::string &line,
-                              std::size_t &n_nodes,
-                              std::size_t &dim,
-                              std::size_t &n_attributes,
-                              bool &boundary_markers) const;
+    static bool parseNodesFileHeader(std::string& line,
+                                     std::size_t& n_nodes,
+                                     std::size_t& dim,
+                                     std::size_t& n_attributes,
+                                     bool& boundary_markers);
     /**
      * method parses the lines reading the nodes from TetGen nodes file
      * @param ins      the input stream (input)
@@ -162,10 +163,10 @@ private:
      * @param region_attribute  is on output true, if there
      * @return
      */
-    bool parseElementsFileHeader(std::string &line,
-                                 std::size_t &n_tets,
-                                 std::size_t &n_nodes_per_tet,
-                                 bool &region_attribute) const;
+    static bool parseElementsFileHeader(std::string& line,
+                                        std::size_t& n_tets,
+                                        std::size_t& n_nodes_per_tet,
+                                        bool& region_attribute);
     /**
      * Method parses the tetrahedras and put them in the element vector of the mesh class.
      * @param ins the input stream
@@ -206,10 +207,10 @@ private:
                          std::vector<MeshLib::Node> &attribute_points) const;
 
     /// Writes facet information from a 2D element to the stream and increments the total element count accordingly
-    void writeElementToFacets(std::ofstream &out,
-                              const MeshLib::Element &element,
-                              unsigned &element_count,
-                              std::string const& matId) const;
+    static void writeElementToFacets(std::ofstream& out,
+                                     const MeshLib::Element& element,
+                                     unsigned& element_count,
+                                     std::string const& matId);
 
     /// the value is true if the indexing is zero based, else false
     bool _zero_based_idx{false};
diff --git a/ChemistryLib/PhreeqcIOData/Output.h b/ChemistryLib/PhreeqcIOData/Output.h
index 50ec32e9bf44bb40a98d9251417e66602eccae77..9c7db8dd9afd2d8a080a13a20be80954c9f2d78f 100644
--- a/ChemistryLib/PhreeqcIOData/Output.h
+++ b/ChemistryLib/PhreeqcIOData/Output.h
@@ -29,14 +29,14 @@ public:
     {
     }
 
-    int getNumberOfItemsInDisplay()
+    static int getNumberOfItemsInDisplay()
     {
         return display_simulation_id + display_state + display_solution_id +
                display_distance + display_current_time + display_time_step +
                display_pH + display_pe;
     }
 
-    int getNumberOfDroppedItems()
+    static int getNumberOfDroppedItems()
     {
         return display_simulation_id + display_state + display_solution_id +
                display_distance + display_current_time + display_time_step;
diff --git a/GeoLib/QuadTree.h b/GeoLib/QuadTree.h
index 589e906e44017138f981ce8e96016c9ce690bbfa..b9cbdef9563dbbee13ab5b7a2892917f8755ca43 100644
--- a/GeoLib/QuadTree.h
+++ b/GeoLib/QuadTree.h
@@ -555,7 +555,7 @@ private:
         _is_leaf = false;
     }
 
-    bool needToRefine (QuadTree<POINT>* node)
+    static bool needToRefine (QuadTree<POINT>* node)
     {
         QuadTree<POINT>* north_neighbor (node->getNorthNeighbor ());
         if (north_neighbor != nullptr)
diff --git a/MaterialLib/Adsorption/Adsorption.cpp b/MaterialLib/Adsorption/Adsorption.cpp
index 05b92b2496a06b08c0ef2698ad5615462a5fa7ff..af782a26f2c5b3887c3724108f0020d5289b8a22 100644
--- a/MaterialLib/Adsorption/Adsorption.cpp
+++ b/MaterialLib/Adsorption/Adsorption.cpp
@@ -100,7 +100,8 @@ double AdsorptionReaction::getReactionRate(const double p_Ads, const double T_Ad
 }
 
 // Evaluate adsorbtion potential A
-double AdsorptionReaction::getPotential(const double p_Ads, double T_Ads, const double M_Ads) const
+double AdsorptionReaction::getPotential(const double p_Ads, double T_Ads,
+                                        const double M_Ads)
 {
     double A = MaterialLib::PhysicalConstant::IdealGasConstant * T_Ads * log(getEquilibriumVapourPressure(T_Ads)/p_Ads) / (M_Ads*1.e3); // in kJ/kg = J/g
     return A;
diff --git a/MaterialLib/Adsorption/Adsorption.h b/MaterialLib/Adsorption/Adsorption.h
index fdb0f04061a76f8475e2f0dd171309e075718463..3b1390fee32f16250ccbb664f34ecf187caf4ee3 100644
--- a/MaterialLib/Adsorption/Adsorption.h
+++ b/MaterialLib/Adsorption/Adsorption.h
@@ -45,7 +45,8 @@ protected:
     virtual double dCharacteristicCurve(const double A) const = 0;
 
 private:
-    double getPotential(const double p_Ads, const double T_Ads, const double M_Ads) const;
+    static double getPotential(const double p_Ads, const double T_Ads,
+                               const double M_Ads);
     double getEntropy(const double T_Ads, const double A) const;
 };
 
diff --git a/MaterialLib/Fluid/GibbsFreeEnergy/DimensionLessGibbsFreeEnergyRegion1.cpp b/MaterialLib/Fluid/GibbsFreeEnergy/DimensionLessGibbsFreeEnergyRegion1.cpp
index 353a49f6aeee2b01b8a67ecce7fd22f3557024de..b562a7b5af21ce3bbc2cdd4b3582b2a4b53f11d5 100644
--- a/MaterialLib/Fluid/GibbsFreeEnergy/DimensionLessGibbsFreeEnergyRegion1.cpp
+++ b/MaterialLib/Fluid/GibbsFreeEnergy/DimensionLessGibbsFreeEnergyRegion1.cpp
@@ -43,7 +43,7 @@ static const double ji[34] = {
     3,  17, -4, 0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41};
 
 double DimensionLessGibbsFreeEnergyRegion1::get_gamma(const double tau,
-                                                      const double pi) const
+                                                      const double pi)
 {
     double val = 0.;
     for (int i = 0; i < 34; i++)
@@ -55,7 +55,7 @@ double DimensionLessGibbsFreeEnergyRegion1::get_gamma(const double tau,
 }
 
 double DimensionLessGibbsFreeEnergyRegion1::get_dgamma_dtau(
-    const double tau, const double pi) const
+    const double tau, const double pi)
 {
     double val = 0.;
     for (int i = 0; i < 34; i++)
@@ -68,7 +68,7 @@ double DimensionLessGibbsFreeEnergyRegion1::get_dgamma_dtau(
 }
 
 double DimensionLessGibbsFreeEnergyRegion1::get_dgamma_dtau_dtau(
-    const double tau, const double pi) const
+    const double tau, const double pi)
 {
     double val = 0.;
     for (int i = 0; i < 34; i++)
@@ -81,7 +81,7 @@ double DimensionLessGibbsFreeEnergyRegion1::get_dgamma_dtau_dtau(
 }
 
 double DimensionLessGibbsFreeEnergyRegion1::get_dgamma_dpi(
-    const double tau, const double pi) const
+    const double tau, const double pi)
 {
     double val = 0.;
     for (int i = 0; i < 34; i++)
@@ -94,7 +94,7 @@ double DimensionLessGibbsFreeEnergyRegion1::get_dgamma_dpi(
 }
 
 double DimensionLessGibbsFreeEnergyRegion1::get_dgamma_dpi_dpi(
-    const double tau, const double pi) const
+    const double tau, const double pi)
 {
     double val = 0.;
     for (int i = 0; i < 34; i++)
@@ -107,7 +107,7 @@ double DimensionLessGibbsFreeEnergyRegion1::get_dgamma_dpi_dpi(
 }
 
 double DimensionLessGibbsFreeEnergyRegion1::get_dgamma_dtau_dpi(
-    const double tau, const double pi) const
+    const double tau, const double pi)
 {
     double val = 0.;
     for (int i = 0; i < 34; i++)
diff --git a/MaterialLib/Fluid/GibbsFreeEnergy/DimensionLessGibbsFreeEnergyRegion1.h b/MaterialLib/Fluid/GibbsFreeEnergy/DimensionLessGibbsFreeEnergyRegion1.h
index 79b122b816e236df5df8cac47d47cc17f353038c..8191e0fa4aa9cff0b9397c803fd1bfa1b1af52c4 100644
--- a/MaterialLib/Fluid/GibbsFreeEnergy/DimensionLessGibbsFreeEnergyRegion1.h
+++ b/MaterialLib/Fluid/GibbsFreeEnergy/DimensionLessGibbsFreeEnergyRegion1.h
@@ -41,7 +41,7 @@ public:
      * @param tau Dimension less pressure
      * @return    The value
      */
-    double get_gamma(const double tau, const double pi) const;
+    static double get_gamma(const double tau, const double pi);
 
     /**
      * Get the 1st order partial derivative of the dimension less Gibbs free
@@ -51,7 +51,7 @@ public:
      * @param tau Dimension less pressure
      * @return    The value
      */
-    double get_dgamma_dtau(const double tau, const double pi) const;
+    static double get_dgamma_dtau(const double tau, const double pi);
 
     /**
      * Get the 2nd order partial derivative of the dimension less Gibbs free
@@ -61,7 +61,7 @@ public:
      * @param tau Dimension less pressure
      * @return    The value
      */
-    double get_dgamma_dtau_dtau(const double tau, const double pi) const;
+    static double get_dgamma_dtau_dtau(const double tau, const double pi);
 
     /**
      * Get the 1st order partial derivative of the dimension less Gibbs free
@@ -71,7 +71,7 @@ public:
      * @param tau Dimension less pressure
      * @return    The value
      */
-    double get_dgamma_dpi(const double tau, const double pi) const;
+    static double get_dgamma_dpi(const double tau, const double pi);
 
     /**
      * Get the 2nd order partial derivative of the dimension less Gibbs free
@@ -81,7 +81,7 @@ public:
      * @param tau Dimension less pressure
      * @return    The value
      */
-    double get_dgamma_dpi_dpi(const double tau, const double pi) const;
+    static double get_dgamma_dpi_dpi(const double tau, const double pi);
 
     /**
      * Get the 2nd order partial derivative of the dimension less Gibbs free
@@ -91,7 +91,7 @@ public:
      * @param tau Dimension less pressure
      * @return    The value
      */
-    double get_dgamma_dtau_dpi(const double tau, const double pi) const;
+    static double get_dgamma_dtau_dpi(const double tau, const double pi);
 };
 
 }  // namespace Fluid
diff --git a/MaterialLib/Fluid/WaterVaporProperties/WaterVaporProperties.cpp b/MaterialLib/Fluid/WaterVaporProperties/WaterVaporProperties.cpp
index e3b65fd5ca4d8b68e3c8b2a9b64870d6cf423901..9740e08c03f9223d6abc3e1657672e405abbd461 100644
--- a/MaterialLib/Fluid/WaterVaporProperties/WaterVaporProperties.cpp
+++ b/MaterialLib/Fluid/WaterVaporProperties/WaterVaporProperties.cpp
@@ -76,7 +76,7 @@ double WaterVaporProperties::calculatedDensityNonwetdT(
 double WaterVaporProperties::getWaterVaporEnthalpySimple(const double temperature,
     const double heat_capacity_water_vapor,
     const double /*pressure*/,
-    const double /*latent_heat_evaporation*/) const
+    const double /*latent_heat_evaporation*/)
 {
     return heat_capacity_water_vapor * (temperature - CelsiusZeroInKelvin) +
         h_wg;
diff --git a/MaterialLib/Fluid/WaterVaporProperties/WaterVaporProperties.h b/MaterialLib/Fluid/WaterVaporProperties/WaterVaporProperties.h
index 198ec8a76832d018e78a97d2a59d65d437033c28..36cf485e0c4ea9b6513b79f5f5fee5bebb894a87 100644
--- a/MaterialLib/Fluid/WaterVaporProperties/WaterVaporProperties.h
+++ b/MaterialLib/Fluid/WaterVaporProperties/WaterVaporProperties.h
@@ -53,10 +53,10 @@ public:
                                  const double T,
                                  const double mass_density_water) const;
     /// Specific enthalpy of water vapor
-    double getWaterVaporEnthalpySimple(const double temperature,
+    static double getWaterVaporEnthalpySimple(const double temperature,
         const double heat_capacity_water_vapor,
         const double /*pressure*/,
-        const double /*latent_heat_evaporation*/) const;
+        const double /*latent_heat_evaporation*/);
 private:
     const double& _water_mol_mass = PhysicalConstant::MolarMass::Water;
     const double& _air_mol_mass = PhysicalConstant::MolarMass::Air;
diff --git a/MaterialLib/MPL/Medium.cpp b/MaterialLib/MPL/Medium.cpp
index cb57bb45f09e4a7fc8bb84f7ae573851b5b55d5c..94669996766ffb316f21f2cacdc79655fbc3a341 100644
--- a/MaterialLib/MPL/Medium.cpp
+++ b/MaterialLib/MPL/Medium.cpp
@@ -63,7 +63,7 @@ std::size_t Medium::numberOfPhases() const
     return phases_.size();
 }
 
-std::string Medium::description() const
+std::string Medium::description()
 {
     return "medium";
 }
diff --git a/MaterialLib/MPL/Medium.h b/MaterialLib/MPL/Medium.h
index f9c96f58d72d9db3436f12a6d5aa95403c3b0e0b..da4a99cc6cc378b1346efea7360de872126e2e00 100644
--- a/MaterialLib/MPL/Medium.h
+++ b/MaterialLib/MPL/Medium.h
@@ -48,7 +48,7 @@ public:
     std::size_t numberOfPhases() const;
 
     /// Short description of the medium.
-    std::string description() const;
+    static std::string description();
 
     template <typename T>
     T value(PropertyType const p) const
diff --git a/MathLib/LinAlg/Eigen/EigenMatrix.h b/MathLib/LinAlg/Eigen/EigenMatrix.h
index d766e6c22923e5d8ed56fd30ced936a890ba6507..9b85b65f693bceb5d98d41f4f5dad02af4445aee 100644
--- a/MathLib/LinAlg/Eigen/EigenMatrix.h
+++ b/MathLib/LinAlg/Eigen/EigenMatrix.h
@@ -57,7 +57,7 @@ public:
     IndexType getNumberOfColumns() const { return _mat.cols(); }
 
     /// return a start index of the active data range
-    IndexType getRangeBegin() const  { return 0; }
+    static constexpr IndexType getRangeBegin() { return 0; }
 
     /// return an end index of the active data range
     IndexType getRangeEnd() const  { return getNumberOfRows(); }
@@ -130,7 +130,7 @@ public:
     }
 
     /// return always true, i.e. the matrix is always ready for use
-    bool isAssembled() const { return true; }
+    static constexpr bool isAssembled() { return true; }
 
     /// printout this matrix for debugging
     void write(const std::string &filename) const
diff --git a/MathLib/LinAlg/Eigen/EigenVector.h b/MathLib/LinAlg/Eigen/EigenVector.h
index b6ddcb7a66e613ba78e8e19bfb4d2baf01a3dbc4..a300a4ce6c02ae8da1d5310a119f4c9c55499a5f 100644
--- a/MathLib/LinAlg/Eigen/EigenVector.h
+++ b/MathLib/LinAlg/Eigen/EigenVector.h
@@ -47,7 +47,7 @@ public:
     IndexType size() const { return static_cast<IndexType>(_vec.size()); }
 
     /// return a start index of the active data range
-    IndexType getRangeBegin() const { return 0;}
+    static constexpr IndexType getRangeBegin() { return 0; }
 
     /// return an end index of the active data range
     IndexType getRangeEnd() const { return size(); }
diff --git a/MeshGeoToolsLib/BoundaryElementsAlongPolyline.cpp b/MeshGeoToolsLib/BoundaryElementsAlongPolyline.cpp
index fe5c8b07fabf38c9939d2e61f38f4a7ba9321149..7e984d237be64faafedfc333a37ee93682ad0bff 100644
--- a/MeshGeoToolsLib/BoundaryElementsAlongPolyline.cpp
+++ b/MeshGeoToolsLib/BoundaryElementsAlongPolyline.cpp
@@ -102,7 +102,7 @@ BoundaryElementsAlongPolyline::~BoundaryElementsAlongPolyline()
 
 bool BoundaryElementsAlongPolyline::includesAllEdgeNodeIDs(
     const std::vector<std::size_t>& vec_node_ids, const MeshLib::Element& edge,
-    std::vector<std::size_t>& edge_node_distances) const
+    std::vector<std::size_t>& edge_node_distances)
 {
     unsigned j = 0;
     for (; j < edge.getNumberOfBaseNodes(); j++)
diff --git a/MeshGeoToolsLib/BoundaryElementsAlongPolyline.h b/MeshGeoToolsLib/BoundaryElementsAlongPolyline.h
index 3b6272331b05d9d2d8f467c0b12e5cbe195b2206..4adbc473ed4415765b888ef13589ee97c9ad5fe2 100644
--- a/MeshGeoToolsLib/BoundaryElementsAlongPolyline.h
+++ b/MeshGeoToolsLib/BoundaryElementsAlongPolyline.h
@@ -75,10 +75,10 @@ private:
      * the beginning of the given node ID vector
      * @return true if all element nodes are included in the vector
      */
-    bool includesAllEdgeNodeIDs(
+    static bool includesAllEdgeNodeIDs(
         const std::vector<std::size_t>& vec_node_ids,
         const MeshLib::Element& edge,
-        std::vector<std::size_t>& edge_node_distances) const;
+        std::vector<std::size_t>& edge_node_distances);
 
     /**
      * Modify node ordering of an edge so that its first node is closer to the
diff --git a/MeshLib/IO/Legacy/MeshIO.cpp b/MeshLib/IO/Legacy/MeshIO.cpp
index f587824c594b30f7ec9b28b96811a7ed3fb8282e..3b9a66614410e8c9ad887992c1fb59da39a94939 100644
--- a/MeshLib/IO/Legacy/MeshIO.cpp
+++ b/MeshLib/IO/Legacy/MeshIO.cpp
@@ -159,7 +159,7 @@ MeshLib::Mesh* MeshIO::loadMeshFromFile(const std::string& file_name)
     return nullptr;
 }
 
-std::size_t MeshIO::readMaterialID(std::istream & in) const
+std::size_t MeshIO::readMaterialID(std::istream & in)
 {
     unsigned index;
     unsigned material_id;
@@ -381,7 +381,7 @@ void MeshIO::writeElements(
     }
 }
 
-std::string MeshIO::ElemType2StringOutput(const MeshLib::MeshElemType t) const
+std::string MeshIO::ElemType2StringOutput(const MeshLib::MeshElemType t)
 {
     if (t == MeshLib::MeshElemType::LINE)
     {
diff --git a/MeshLib/IO/Legacy/MeshIO.h b/MeshLib/IO/Legacy/MeshIO.h
index adf90df9c36d523bca58dae95f0aefb4441f4c2d..802ad6d387afdf21ca378eda07559e4abf6a899e 100644
--- a/MeshLib/IO/Legacy/MeshIO.h
+++ b/MeshLib/IO/Legacy/MeshIO.h
@@ -55,10 +55,10 @@ private:
     void writeElements(std::vector<MeshLib::Element*> const& ele_vec,
                        MeshLib::PropertyVector<int> const* const material_ids,
                        std::ostream& out) const;
-    std::size_t readMaterialID(std::istream & in) const;
+    static std::size_t readMaterialID(std::istream & in);
     MeshLib::Element* readElement(
         std::istream& in, const std::vector<MeshLib::Node*>& nodes) const;
-    std::string ElemType2StringOutput(const MeshLib::MeshElemType t) const;
+    static std::string ElemType2StringOutput(const MeshLib::MeshElemType t);
 
     const MeshLib::Mesh* _mesh{nullptr};
 
diff --git a/MeshLib/MeshEditing/MeshRevision.cpp b/MeshLib/MeshEditing/MeshRevision.cpp
index 5f48663b0a34d115d5c5a67dc092e89b4093b0e7..1ad29941837f61eb209028add34b80b621562bcd 100644
--- a/MeshLib/MeshEditing/MeshRevision.cpp
+++ b/MeshLib/MeshEditing/MeshRevision.cpp
@@ -213,7 +213,8 @@ std::vector<MeshLib::Node*> MeshRevision::constructNewNodesArray(const std::vect
     return new_nodes;
 }
 
-unsigned MeshRevision::getNumberOfUniqueNodes(MeshLib::Element const*const element) const
+unsigned MeshRevision::getNumberOfUniqueNodes(
+    MeshLib::Element const* const element)
 {
     unsigned const nNodes(element->getNumberOfBaseNodes());
     unsigned count(nNodes);
@@ -654,7 +655,7 @@ unsigned MeshRevision::reducePrism(MeshLib::Element const*const org_elem,
                     // triangle edge collapsed
                     const unsigned i_offset = (i>2) ? i-3 : i+3;
                     const unsigned j_offset = (i>2) ? j-3 : j+3;
-                    const unsigned k = this->lutPrismThirdNode(i,j);
+                    const unsigned k = lutPrismThirdNode(i,j);
                     if (k == std::numeric_limits<unsigned>::max())
                     {
                         ERR("Unexpected error during prism reduction.");
@@ -806,8 +807,9 @@ MeshLib::Element* MeshRevision::constructFourNodeElement(
     return nullptr;
 }
 
-unsigned MeshRevision::findPyramidTopNode(MeshLib::Element const& element,
-    std::array<std::size_t,4> const& base_node_ids) const
+unsigned MeshRevision::findPyramidTopNode(
+    MeshLib::Element const& element,
+    std::array<std::size_t, 4> const& base_node_ids)
 {
     const std::size_t nNodes (element.getNumberOfBaseNodes());
     for (std::size_t i=0; i<nNodes; ++i)
@@ -828,13 +830,13 @@ unsigned MeshRevision::findPyramidTopNode(MeshLib::Element const& element,
     return std::numeric_limits<unsigned>::max(); // should never be reached if called correctly
 }
 
-unsigned MeshRevision::lutHexDiametralNode(unsigned id) const
+unsigned MeshRevision::lutHexDiametralNode(unsigned id)
 {
     return _hex_diametral_nodes[id];
 }
 
 std::array<unsigned, 4> MeshRevision::lutHexCuttingQuadNodes(unsigned id1,
-                                                             unsigned id2) const
+                                                             unsigned id2)
 {
     std::array<unsigned,4> nodes{};
     if      (id1==0 && id2==1) { nodes[0]=3; nodes[1]=2; nodes[2]=5; nodes[3]=4; }
@@ -868,23 +870,23 @@ std::array<unsigned, 4> MeshRevision::lutHexCuttingQuadNodes(unsigned id1,
 std::pair<unsigned, unsigned> MeshRevision::lutHexBackNodes(unsigned i,
                                                             unsigned j,
                                                             unsigned k,
-                                                            unsigned l) const
+                                                            unsigned l)
 {
     // collapsed edges are *not* connected
     std::pair<unsigned, unsigned> back(std::numeric_limits<unsigned>::max(), std::numeric_limits<unsigned>::max());
-    if      (this->lutHexDiametralNode(i) == k) { back.first = i; back.second = this->lutHexDiametralNode(l); }
-    else if (this->lutHexDiametralNode(i) == l) { back.first = i; back.second = this->lutHexDiametralNode(k); }
-    else if (this->lutHexDiametralNode(j) == k) { back.first = j; back.second = this->lutHexDiametralNode(l); }
-    else if (this->lutHexDiametralNode(j) == l) { back.first = j; back.second = this->lutHexDiametralNode(k); }
+    if      (lutHexDiametralNode(i) == k) { back.first = i; back.second = lutHexDiametralNode(l); }
+    else if (lutHexDiametralNode(i) == l) { back.first = i; back.second = lutHexDiametralNode(k); }
+    else if (lutHexDiametralNode(j) == k) { back.first = j; back.second = lutHexDiametralNode(l); }
+    else if (lutHexDiametralNode(j) == l) { back.first = j; back.second = lutHexDiametralNode(k); }
     // collapsed edges *are* connected
-    else if (i==k) { back.first = this->lutHexDiametralNode(l); back.second = j; }
-    else if (i==l) { back.first = this->lutHexDiametralNode(k); back.second = j; }
-    else if (j==k) { back.first = this->lutHexDiametralNode(l); back.second = i; }
-    else if (j==l) { back.first = this->lutHexDiametralNode(k); back.second = i; }
+    else if (i==k) { back.first = lutHexDiametralNode(l); back.second = j; }
+    else if (i==l) { back.first = lutHexDiametralNode(k); back.second = j; }
+    else if (j==k) { back.first = lutHexDiametralNode(l); back.second = i; }
+    else if (j==l) { back.first = lutHexDiametralNode(k); back.second = i; }
     return back;
 }
 
-unsigned MeshRevision::lutPrismThirdNode(unsigned id1, unsigned id2) const
+unsigned MeshRevision::lutPrismThirdNode(unsigned id1, unsigned id2)
 {
     if ((id1 == 0 && id2 == 1) || (id1 == 1 && id2 == 2))
     {
diff --git a/MeshLib/MeshEditing/MeshRevision.h b/MeshLib/MeshEditing/MeshRevision.h
index b23aba3ce7b990b5154c11bd71edbb3e1b916680..379c6bbafa4f3e1c75457ff3d914c01bd9356e56 100644
--- a/MeshLib/MeshEditing/MeshRevision.h
+++ b/MeshLib/MeshEditing/MeshRevision.h
@@ -69,7 +69,8 @@ private:
         const std::vector<std::size_t> &id_map) const;
 
     /// Calculates the number of unique nodes in an element (i.e. uncollapsed nodes)
-    unsigned getNumberOfUniqueNodes(MeshLib::Element const*const element) const;
+    static unsigned getNumberOfUniqueNodes(
+        MeshLib::Element const* const element);
 
     /// Resets the node IDs of the source mesh (needs to be called after everything is done).
     void resetNodeIDs();
@@ -154,22 +155,23 @@ private:
                          unsigned min_elem_dim) const;
 
     // In an element with 5 unique nodes, return the node that will be the top of the resulting pyramid
-    unsigned findPyramidTopNode(MeshLib::Element const& element,
-        std::array<std::size_t,4> const& base_node_ids) const;
+    static unsigned findPyramidTopNode(MeshLib::Element const& element,
+        std::array<std::size_t,4> const& base_node_ids);
 
     /// Lookup-table for returning the diametral node id of the given node id in a Hex
-    unsigned lutHexDiametralNode(unsigned id) const;
+    static unsigned lutHexDiametralNode(unsigned id);
 
     /// Lookup-table for returning four nodes connected to the two nodes (id1, id2) forming an edge in a Hex
-    std::array<unsigned, 4> lutHexCuttingQuadNodes(unsigned id1,
-                                                   unsigned id2) const;
+    static std::array<unsigned, 4> lutHexCuttingQuadNodes(unsigned id1,
+                                                          unsigned id2);
 
     /// When a hex is subdivided into two prisms, this returns the nodes of the hex edge that will serve as the back of one of the prisms.
-    std::pair<unsigned, unsigned> lutHexBackNodes(unsigned i, unsigned j,
-                                                  unsigned k, unsigned l) const;
+    static std::pair<unsigned, unsigned> lutHexBackNodes(unsigned i, unsigned j,
+                                                         unsigned k,
+                                                         unsigned l);
 
     /// Lookup-table for returning the third node of bottom or top triangle given the other two
-    unsigned lutPrismThirdNode(unsigned id1, unsigned id2) const;
+    static unsigned lutPrismThirdNode(unsigned id1, unsigned id2);
 
     /// The original mesh used for constructing the class
     Mesh& _mesh;
diff --git a/MeshLib/MeshGenerators/LayeredMeshGenerator.h b/MeshLib/MeshGenerators/LayeredMeshGenerator.h
index d83234e872c8545b948cb0b2a09ece160f1ef15b..289bb9513c079894e1f5095b525dae72fb23a08e 100644
--- a/MeshLib/MeshGenerators/LayeredMeshGenerator.h
+++ b/MeshLib/MeshGenerators/LayeredMeshGenerator.h
@@ -87,7 +87,8 @@ protected:
                                    std::size_t new_node_id) const;
 
     /// Calculates a data-dependent epsilon value
-    double calcEpsilon(GeoLib::Raster const& low, GeoLib::Raster const& high);
+    static double calcEpsilon(GeoLib::Raster const& low,
+                              GeoLib::Raster const& high);
 
     double _elevation_epsilon{0.0001};
     double _minimum_thickness{std::numeric_limits<float>::epsilon()};
diff --git a/NumLib/Fem/Integration/IntegrationPoint.h b/NumLib/Fem/Integration/IntegrationPoint.h
index 69d54fdf91c5adf7577758554e1479f83f9e38fa..a59a00366fb5dcffcfa5146862c10482645b9480 100644
--- a/NumLib/Fem/Integration/IntegrationPoint.h
+++ b/NumLib/Fem/Integration/IntegrationPoint.h
@@ -29,24 +29,16 @@ public:
     }
 
     /// Change the integration order.
-    void setIntegrationOrder(unsigned /* order */)
-    {
-    }
+    static void setIntegrationOrder(unsigned /* order */) {}
 
     /// Return current integration order.
-    unsigned getIntegrationOrder() const
-    {
-        return 0;
-    }
+    static constexpr unsigned getIntegrationOrder() { return 0; }
 
     /// Return the number of sampling points.
-    unsigned getNumberOfPoints() const
-    {
-        return 1;
-    }
+    static constexpr unsigned getNumberOfPoints() { return 1; }
 
     /// \copydoc IntegrationGaussLegendreRegular::getWeightedPoint(unsigned) const
-    WeightedPoint getWeightedPoint(unsigned igp) const
+    static WeightedPoint getWeightedPoint(unsigned igp)
     {
         return getWeightedPoint(getIntegrationOrder(), igp);
     }
@@ -69,7 +61,7 @@ public:
     ///
     /// \param order    the number of integration points
     /// \return the number of points.
-    static unsigned getNumberOfPoints(unsigned order)
+    static constexpr unsigned getNumberOfPoints(unsigned order)
     {
         (void)order;
         return 1;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
index b6e753491bbb619ca6442d75829fbd4a934ad0c7..f88b30fc2086035e07982f43f17b6b07e128d6f5 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
@@ -135,7 +135,7 @@ std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>
 BHECommonCoaxial::getBHEInflowDirichletBCNodesAndComponents(
     std::size_t const top_node_id,
     std::size_t const /*bottom_node_id*/,
-    int const in_component_id) const
+    int const in_component_id)
 {
     return {std::make_pair(top_node_id, in_component_id),
             std::make_pair(top_node_id, in_component_id + 1)};
@@ -145,7 +145,7 @@ std::optional<
     std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>>
 BHECommonCoaxial::getBHEBottomDirichletBCNodesAndComponents(
     std::size_t const bottom_node_id, int const in_component_id,
-    int const out_component_id) const
+    int const out_component_id)
 {
     return {{std::make_pair(bottom_node_id, in_component_id),
              std::make_pair(bottom_node_id, out_component_id)}};
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
index 161f28d967ffe8b4b0727af053a25b1d7faf05d9..bd67c455acade3d9c424c87d3c75a7643ac6732f 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
@@ -49,17 +49,17 @@ public:
     static constexpr std::pair<int, int> inflow_outflow_bc_component_ids[] = {
         {0, 1}};
 
-    std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>
+    static std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>
     getBHEInflowDirichletBCNodesAndComponents(
         std::size_t const top_node_id,
         std::size_t const /*bottom_node_id*/,
-        int const in_component_id) const;
+        int const in_component_id);
 
-    std::optional<
+    static std::optional<
         std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>>
     getBHEBottomDirichletBCNodesAndComponents(std::size_t const bottom_node_id,
                                               int const in_component_id,
-                                              int const out_component_id) const;
+                                              int const out_component_id);
 
     std::array<double, number_of_unknowns> pipeHeatConductions() const;
 
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp
index 1b3ca1d6996ba0f8bc12c6e89b44eba6cd3b131a..f8015e89efc4db19e8ab8082c721d0fc74b39cb3 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp
@@ -151,7 +151,7 @@ std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>
 BHE_1P::getBHEInflowDirichletBCNodesAndComponents(
     std::size_t const top_node_id,
     std::size_t const bottom_node_id,
-    int const in_component_id) const
+    int const in_component_id)
 {
     return {std::make_pair(top_node_id, in_component_id),
             std::make_pair(bottom_node_id, in_component_id)};
@@ -162,7 +162,7 @@ std::optional<
 BHE_1P::getBHEBottomDirichletBCNodesAndComponents(
     std::size_t const /*bottom_node_id*/,
     int const /*in_component_id*/,
-    int const /*out_component_id*/) const
+    int const /*out_component_id*/)
 {
     return {};
 }
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1P.h b/ProcessLib/HeatTransportBHE/BHE/BHE_1P.h
index 6f401c9f5612bd3de1b5803b9ad0bce1cf2149e7..c5bbf6dec36ed442d6ee6044a5b8c40aa38c8e28 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1P.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1P.h
@@ -62,12 +62,12 @@ public:
               typename RMatrixType,
               typename RPiSMatrixType,
               typename RSMatrixType>
-    void assembleRMatrices(
+    static void assembleRMatrices(
         int const idx_bhe_unknowns,
         Eigen::MatrixBase<SingleUnknownMatrixType> const& matBHE_loc_R,
         Eigen::MatrixBase<RMatrixType>& R_matrix,
         Eigen::MatrixBase<RPiSMatrixType>& R_pi_s_matrix,
-        Eigen::MatrixBase<RSMatrixType>& R_s_matrix) const
+        Eigen::MatrixBase<RSMatrixType>& R_s_matrix)
     {
         // Here we are looping over two resistance terms
         // First PHI_fg is the resistance between pipe and grout
@@ -112,17 +112,17 @@ public:
     static constexpr std::pair<int, int> inflow_outflow_bc_component_ids[] = {
         {0, 1}};
 
-    std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>
+    static std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>
     getBHEInflowDirichletBCNodesAndComponents(std::size_t const top_node_id,
                                               std::size_t const bottom_node_id,
-                                              int const in_component_id) const;
+                                              int const in_component_id);
 
-    std::optional<
+    static std::optional<
         std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>>
     getBHEBottomDirichletBCNodesAndComponents(
         std::size_t const /*bottom_node_id*/,
         int const /*in_component_id*/,
-        int const /*out_component_id*/) const;
+        int const /*out_component_id*/);
 
 public:
     std::array<double, number_of_unknowns> crossSectionAreas() const;
@@ -139,7 +139,7 @@ private:
     std::array<double, number_of_unknowns> calcThermalResistances(
         double const Nu);
 
-    double compute_R_gs(double const chi, double const R_g);
+    static double compute_R_gs(double const chi, double const R_g);
 
 private:
     /// PHI_fg, PHI_gs;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
index e84b5c9e11de975499dec4db53c375e2e63d9d6e..7083e7212c051fafeb3bdb59ec2508cfc89830b8 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
@@ -232,7 +232,7 @@ std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>
 BHE_1U::getBHEInflowDirichletBCNodesAndComponents(
     std::size_t const top_node_id,
     std::size_t const /*bottom_node_id*/,
-    int const in_component_id) const
+    int const in_component_id)
 {
     return {std::make_pair(top_node_id, in_component_id),
             std::make_pair(top_node_id, in_component_id + 1)};
@@ -243,7 +243,7 @@ std::optional<
 BHE_1U::getBHEBottomDirichletBCNodesAndComponents(
     std::size_t const bottom_node_id,
     int const in_component_id,
-    int const out_component_id) const
+    int const out_component_id)
 {
     return {{std::make_pair(bottom_node_id, in_component_id),
              std::make_pair(bottom_node_id, out_component_id)}};
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
index 43b3451cd3daa44c7681496469462bd13a6b4c5e..11cee1e40b6c95e0aaefb0acd6d58f704445fd3f 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.h
@@ -64,12 +64,12 @@ public:
               typename RMatrixType,
               typename RPiSMatrixType,
               typename RSMatrixType>
-    void assembleRMatrices(
+    static void assembleRMatrices(
         int const idx_bhe_unknowns,
         Eigen::MatrixBase<SingleUnknownMatrixType> const& matBHE_loc_R,
         Eigen::MatrixBase<RMatrixType>& R_matrix,
         Eigen::MatrixBase<RPiSMatrixType>& R_pi_s_matrix,
-        Eigen::MatrixBase<RSMatrixType>& R_s_matrix) const
+        Eigen::MatrixBase<RSMatrixType>& R_s_matrix)
     {
         switch (idx_bhe_unknowns)
         {
@@ -139,16 +139,16 @@ public:
     static constexpr std::pair<int, int> inflow_outflow_bc_component_ids[] = {
         {0, 1}};
 
-    std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>
+    static std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>
     getBHEInflowDirichletBCNodesAndComponents(
         std::size_t const top_node_id,
         std::size_t const /*bottom_node_id*/,
-        int const in_component_id) const;
-    std::optional<
+        int const in_component_id);
+    static std::optional<
         std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>>
     getBHEBottomDirichletBCNodesAndComponents(std::size_t const bottom_node_id,
                                               int const in_component_id,
-                                              int const out_component_id) const;
+                                              int const out_component_id);
 
 public:
     std::array<double, number_of_unknowns> crossSectionAreas() const;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
index 82ec322b0042f6018bd259664d38869002033b31..806778eb4de82b13cc0362b52649ccd945c01719 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
@@ -262,7 +262,7 @@ std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>
 BHE_2U::getBHEInflowDirichletBCNodesAndComponents(
     std::size_t const top_node_id,
     std::size_t const /*bottom_node_id*/,
-    int const in_component_id) const
+    int const in_component_id)
 {
     return {std::make_pair(top_node_id, in_component_id),
             std::make_pair(top_node_id, in_component_id + 2)};
@@ -273,7 +273,7 @@ std::optional<
 BHE_2U::getBHEBottomDirichletBCNodesAndComponents(
     std::size_t const bottom_node_id,
     int const in_component_id,
-    int const out_component_id) const
+    int const out_component_id)
 {
     return {{std::make_pair(bottom_node_id, in_component_id),
              std::make_pair(bottom_node_id, out_component_id)}};
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h
index 048447da113c371b95fefc54867230e90d46edab..1ae7032bfb86f1ebd2cd49f0e811e00a95bfe1b3 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.h
@@ -64,12 +64,12 @@ public:
               typename RMatrixType,
               typename RPiSMatrixType,
               typename RSMatrixType>
-    void assembleRMatrices(
+    static void assembleRMatrices(
         int const idx_bhe_unknowns,
         Eigen::MatrixBase<SingleUnknownMatrixType> const& matBHE_loc_R,
         Eigen::MatrixBase<RMatrixType>& R_matrix,
         Eigen::MatrixBase<RPiSMatrixType>& R_pi_s_matrix,
-        Eigen::MatrixBase<RSMatrixType>& R_s_matrix) const
+        Eigen::MatrixBase<RSMatrixType>& R_s_matrix)
     {
         switch (idx_bhe_unknowns)
         {
@@ -201,17 +201,17 @@ public:
     static constexpr std::pair<int, int> inflow_outflow_bc_component_ids[] = {
         {0, 2}, {1, 3}};
 
-    std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>
+    static std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>
     getBHEInflowDirichletBCNodesAndComponents(
         std::size_t const top_node_id,
         std::size_t const /*bottom_node_id*/,
-        int const in_component_id) const;
+        int const in_component_id);
 
-    std::optional<
+    static std::optional<
         std::array<std::pair<std::size_t /*node_id*/, int /*component*/>, 2>>
     getBHEBottomDirichletBCNodesAndComponents(std::size_t const bottom_node_id,
                                               int const in_component_id,
-                                              int const out_component_id) const;
+                                              int const out_component_id);
 
 public:
     std::array<double, number_of_unknowns> crossSectionAreas() const;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.h b/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.h
index b3f1bafe499a2b8537a8ec254e012ac80a526291..5a9c6383839316ca49648c101e7691016c4bb6f8 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_CXA.h
@@ -61,12 +61,12 @@ public:
     template <int NPoints, typename SingleUnknownMatrixType,
               typename RMatrixType, typename RPiSMatrixType,
               typename RSMatrixType>
-    void assembleRMatrices(
+    static void assembleRMatrices(
         int const idx_bhe_unknowns,
         Eigen::MatrixBase<SingleUnknownMatrixType> const& matBHE_loc_R,
         Eigen::MatrixBase<RMatrixType>& R_matrix,
         Eigen::MatrixBase<RPiSMatrixType>& R_pi_s_matrix,
-        Eigen::MatrixBase<RSMatrixType>& R_s_matrix) const
+        Eigen::MatrixBase<RSMatrixType>& R_s_matrix)
     {
         switch (idx_bhe_unknowns)
         {
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.h b/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.h
index 949dbed35761628d1e0500c56da9d7095b551e01..462b5574735107029e680672e971f98ea159999d 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_CXC.h
@@ -59,12 +59,12 @@ public:
     template <int NPoints, typename SingleUnknownMatrixType,
               typename RMatrixType, typename RPiSMatrixType,
               typename RSMatrixType>
-    void assembleRMatrices(
+    static void assembleRMatrices(
         int const idx_bhe_unknowns,
         Eigen::MatrixBase<SingleUnknownMatrixType> const& matBHE_loc_R,
         Eigen::MatrixBase<RMatrixType>& R_matrix,
         Eigen::MatrixBase<RPiSMatrixType>& R_pi_s_matrix,
-        Eigen::MatrixBase<RSMatrixType>& R_s_matrix) const
+        Eigen::MatrixBase<RSMatrixType>& R_s_matrix)
     {
         switch (idx_bhe_unknowns)
         {
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture.h
index 20e4e6d78b00f793cf07805c4a47a069285247a3..94ae08a4e618484b65de891431051250b1bd49f9 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture.h
@@ -219,15 +219,15 @@ public:
     }
 
 private:
-    std::vector<double> const& getIntPtSigma(
-        std::vector<double>& cache, std::size_t const /*component*/) const
+    static std::vector<double> const& getIntPtSigma(
+        std::vector<double>& cache, std::size_t const /*component*/)
     {
         cache.resize(0);
         return cache;
     }
 
-    std::vector<double> const& getIntPtEpsilon(
-        std::vector<double>& cache, std::size_t const /*component*/) const
+    static std::vector<double> const& getIntPtEpsilon(
+        std::vector<double>& cache, std::size_t const /*component*/)
     {
         cache.resize(0);
 
diff --git a/ProcessLib/Output/Output.cpp b/ProcessLib/Output/Output.cpp
index 8f798302066bd3d14e3d97a3721f75a92df22f07..54340fd6ffb793b2aecc344156c95fd9c5d4b6d0 100644
--- a/ProcessLib/Output/Output.cpp
+++ b/ProcessLib/Output/Output.cpp
@@ -252,7 +252,7 @@ void Output::outputMesh(OutputFile const& output_file,
                         MeshLib::IO::PVDFile* const pvd_file,
                         MeshLib::Mesh const& mesh,
                         int const timestep,
-                        double const t) const
+                        double const t)
 {
     DBUG("output to {:s}", output_file.path);
 
diff --git a/ProcessLib/Output/Output.h b/ProcessLib/Output/Output.h
index 7cadc3603966d13a0bc5cf4bb0c3af705aa2533c..3756d838d6d6687261be993100c3557073b046ad 100644
--- a/ProcessLib/Output/Output.h
+++ b/ProcessLib/Output/Output.h
@@ -89,11 +89,11 @@ public:
 private:
     struct OutputFile;
 
-    void outputMesh(OutputFile const& output_file,
+    static void outputMesh(OutputFile const& output_file,
                     MeshLib::IO::PVDFile* const pvd_file,
                     MeshLib::Mesh const& mesh,
                     int const timestep,
-                    double const t) const;
+                    double const t);
 #ifdef OGS_USE_XDMF
     void outputMeshXdmf(OutputFile const& output_file,
                         MeshLib::Mesh const& mesh,
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPMaterialProperties.cpp b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPMaterialProperties.cpp
index 0010524d4664d78c784c5be3f84ed0940f436cc1..184ad125ec44c1ff14ff46c2b854af8c2cc9e471 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPMaterialProperties.cpp
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPMaterialProperties.cpp
@@ -129,8 +129,7 @@ ThermalTwoPhaseFlowWithPPMaterialProperties::getThermalConductivityWetSolid(
 double
 ThermalTwoPhaseFlowWithPPMaterialProperties::calculateUnsatHeatConductivity(
     double const /*t*/, ParameterLib::SpatialPosition const& /*x*/,
-    double const Sw, double const lambda_pm_dry,
-    double const lambda_pm_wet) const
+    double const Sw, double const lambda_pm_dry, double const lambda_pm_wet)
 {
     double lambda_pm =
         lambda_pm_dry + std::sqrt(Sw) * (lambda_pm_wet - lambda_pm_dry);
@@ -204,7 +203,7 @@ double
 ThermalTwoPhaseFlowWithPPMaterialProperties::getLiquidWaterEnthalpySimple(
     const double temperature,
     const double heat_capacity_liquid_water,
-    const double /*pl*/) const
+    const double /*pl*/)
 {
     return heat_capacity_liquid_water * (temperature - CelsiusZeroInKelvin);
 }
diff --git a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPMaterialProperties.h b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPMaterialProperties.h
index 9d56a96f51a70516d2e364c72083e9d5eaf5d68e..fa5af55df181a8bc15c1ac852ca89dfb6ced5f5b 100644
--- a/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPMaterialProperties.h
+++ b/ProcessLib/ThermalTwoPhaseFlowWithPP/ThermalTwoPhaseFlowWithPPMaterialProperties.h
@@ -75,12 +75,12 @@ public:
     double getThermalConductivityDrySolid(const double p, const double T) const;
     double getThermalConductivityWetSolid(const double p, const double T) const;
     /// Calculates the unsaturated heat conductivity
-    double calculateUnsatHeatConductivity(
+    static double calculateUnsatHeatConductivity(
         double const t,
         ParameterLib::SpatialPosition const& x,
         double const Sw,
         double const lambda_pm_dry,
-        double const lambda_pm_wet) const;
+        double const lambda_pm_wet);
     /// water vapor saturation pressure
     double calculateSaturatedVaporPressure(const double T) const;
     /// partial water vapor pressure in nonwetting phase
@@ -111,9 +111,9 @@ public:
                                 const double heat_capacity_dry_air,
                                 const double /*pg*/) const;
     /// Specific enthalpy of liquid water
-    double getLiquidWaterEnthalpySimple(const double temperature,
+    static double getLiquidWaterEnthalpySimple(const double temperature,
                                         const double heat_capacity_liquid_water,
-                                        const double /*pl*/) const;
+                                        const double /*pl*/);
     const MaterialLib::TwoPhaseFlowWithPP::TwoPhaseFlowWithPPMaterialProperties&
     getTwoPhaseMaterialModel() const
     {
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.cpp b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.cpp
index 12d82d7de99419e84e271ae4afdb0c38321aac1e..126341f582ae88ea182694b9e125f27710005274 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.cpp
@@ -272,7 +272,7 @@ void TwoPhaseFlowWithPrhoMaterialProperties::calculateJacobian(
 * for calculating molar fraction of light component in the liquid phase
 */
 double TwoPhaseFlowWithPrhoMaterialProperties::calculateEquilibiumRhoWetLight(
-    double const pg, double const Sw, double const rho_wet_h2) const
+    double const pg, double const Sw, double const rho_wet_h2)
 {
     double const rho_equilibrium_wet_h2 = pg * HenryConstantH2 * H2;
     return std::min(1 - Sw, rho_equilibrium_wet_h2 - rho_wet_h2);
@@ -283,7 +283,7 @@ double TwoPhaseFlowWithPrhoMaterialProperties::calculateEquilibiumRhoWetLight(
 */
 double TwoPhaseFlowWithPrhoMaterialProperties::calculateSaturation(
     double /*PL*/, double X, double Sw, double rho_wet_h2, double rho_nonwet_h2,
-    double /*T*/) const
+    double /*T*/)
 {
     return X - (Sw * rho_wet_h2 + (1 - Sw) * rho_nonwet_h2);
 }
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.h b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.h
index bde5e2505851b917b9c9d6988151dd735c029801..87305bac2b44eb60ffd2c588d8b81c149544ad2b 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.h
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.h
@@ -152,14 +152,15 @@ private:
     /** Complementary condition 1
     * for calculating molar fraction of light component in the liquid phase
     */
-    double calculateEquilibiumRhoWetLight(double const pg, double const Sw,
-                                          double const rho_wet_h2) const;
+    static double calculateEquilibiumRhoWetLight(double const pg,
+                                                 double const Sw,
+                                                 double const rho_wet_h2);
     /** Complementary condition 2
     * for calculating the saturation
     */
-    double calculateSaturation(double /*PL*/, double X, double Sw,
-                               double rho_wet_h2, double rho_nonwet_h2,
-                               double /*T*/) const;
+    static double calculateSaturation(double /*PL*/, double X, double Sw,
+                                      double rho_wet_h2, double rho_nonwet_h2,
+                                      double /*T*/);
     /**
     * Calculate the derivatives using the analytical way
     */