diff --git a/Applications/DataHolderLib/Project.cpp b/Applications/DataHolderLib/Project.cpp
index d5d7da788c6..477a68a6ed5 100644
--- a/Applications/DataHolderLib/Project.cpp
+++ b/Applications/DataHolderLib/Project.cpp
@@ -77,7 +77,7 @@ bool Project::getUniqueName(std::string& name) const
             cpName = cpName + "-" + std::to_string(count);
         }
 
-        for (auto& mesh : _mesh_vec)
+        for (auto const& mesh : _mesh_vec)
         {
             if (cpName == mesh->getName())
             {
diff --git a/Applications/FileIO/GMSInterface.cpp b/Applications/FileIO/GMSInterface.cpp
index 408b9df1b10..d3aedea1e57 100644
--- a/Applications/FileIO/GMSInterface.cpp
+++ b/Applications/FileIO/GMSInterface.cpp
@@ -147,10 +147,10 @@ void GMSInterface::writeBoreholesToGMS(
 
     for (auto station_as_point : *stations)
     {
-        auto* station = static_cast<GeoLib::StationBorehole*>(station_as_point);
-        std::vector<GeoLib::Point*> profile = station->getProfile();
-        std::vector<std::string> soilNames = station->getSoilNames();
-        // std::size_t idx = 0;
+        auto const* station =
+            static_cast<GeoLib::StationBorehole*>(station_as_point);
+        std::vector<GeoLib::Point*> const& profile = station->getProfile();
+        std::vector<std::string> const& soilNames = station->getSoilNames();
         std::string current_soil_name;
 
         std::size_t nLayers = profile.size();
diff --git a/Applications/FileIO/Gmsh/GMSHPolygonTree.cpp b/Applications/FileIO/Gmsh/GMSHPolygonTree.cpp
index fd683000c69..f1e17af9d6f 100644
--- a/Applications/FileIO/Gmsh/GMSHPolygonTree.cpp
+++ b/Applications/FileIO/Gmsh/GMSHPolygonTree.cpp
@@ -461,7 +461,7 @@ void GMSHPolygonTree::writeAdditionalPointData(std::size_t& pnt_id_offset,
                                                std::size_t sfc_number,
                                                std::ostream& out) const
 {
-    if (auto* adaptive_mesh_density =
+    if (auto const* adaptive_mesh_density =
             dynamic_cast<GMSHAdaptiveMeshDensity*>(&_mesh_density_strategy))
     {
         std::vector<GeoLib::Point*> steiner_pnts;
diff --git a/Applications/FileIO/GocadIO/GocadAsciiReader.cpp b/Applications/FileIO/GocadIO/GocadAsciiReader.cpp
index 3d53896a17e..32e21db227d 100644
--- a/Applications/FileIO/GocadIO/GocadAsciiReader.cpp
+++ b/Applications/FileIO/GocadIO/GocadAsciiReader.cpp
@@ -355,7 +355,7 @@ bool parseNodes(std::ifstream& in,
 
 /// Parses the segments of the current line
 bool parseLineSegments(std::ifstream& in,
-                       std::vector<MeshLib::Node*>& nodes,
+                       std::vector<MeshLib::Node*> const& nodes,
                        std::vector<MeshLib::Element*>& elems,
                        std::map<std::size_t, std::size_t> const& node_id_map,
                        MeshLib::Properties& mesh_prop)
diff --git a/Applications/FileIO/GocadIO/GocadSGridReader.cpp b/Applications/FileIO/GocadIO/GocadSGridReader.cpp
index 3a3523b3105..aec0bbc81bb 100644
--- a/Applications/FileIO/GocadIO/GocadSGridReader.cpp
+++ b/Applications/FileIO/GocadIO/GocadSGridReader.cpp
@@ -711,9 +711,9 @@ void GocadSGridReader::applySplitInformation(
     }
 }
 
-void GocadSGridReader::modifyElement(MeshLib::Element* hex,
-                                     MeshLib::Node const* node2sub,
-                                     MeshLib::Node* substitute_node)
+void GocadSGridReader::modifyElement(MeshLib::Element const* const hex,
+                                     MeshLib::Node const* const node2sub,
+                                     MeshLib::Node* const 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 610c90dc08f..25ed1a197f0 100644
--- a/Applications/FileIO/GocadIO/GocadSGridReader.h
+++ b/Applications/FileIO/GocadIO/GocadSGridReader.h
@@ -74,7 +74,7 @@ class GocadSGridReader final
     void applySplitInformation(
         std::vector<MeshLib::Node*>& nodes,
         std::vector<MeshLib::Element*> const& elements) const;
-    static void modifyElement(MeshLib::Element* hex,
+    static void modifyElement(MeshLib::Element const* hex,
                               MeshLib::Node const* node2sub,
                               MeshLib::Node* substitute_node);
 
diff --git a/Applications/FileIO/Legacy/createSurface.cpp b/Applications/FileIO/Legacy/createSurface.cpp
index c21bd3a2097..05de5ce25a4 100644
--- a/Applications/FileIO/Legacy/createSurface.cpp
+++ b/Applications/FileIO/Legacy/createSurface.cpp
@@ -95,7 +95,7 @@ bool createSurface(GeoLib::Polyline const& ply,
         WARN("Call to '{:s}' returned non-zero value {:d}.", gmsh_command,
              gmsh_return_value);
     }
-    auto surface_mesh = MeshLib::IO::readMeshFromFile(msh_file.string());
+    auto const* surface_mesh = MeshLib::IO::readMeshFromFile(msh_file.string());
     if (!surface_mesh)
     {
         WARN("The surface mesh could not be created.");
diff --git a/Applications/Utils/MeshEdit/NodeReordering.cpp b/Applications/Utils/MeshEdit/NodeReordering.cpp
index 636872eafb6..98d660ec22e 100644
--- a/Applications/Utils/MeshEdit/NodeReordering.cpp
+++ b/Applications/Utils/MeshEdit/NodeReordering.cpp
@@ -140,11 +140,9 @@ void reorderNonlinearNodes(MeshLib::Mesh& mesh)
     }
 
     BaseLib::makeVectorUnique(base_nodes,
-                              [](MeshLib::Node* a, MeshLib::Node* b)
-                              { return a->getID() < b->getID(); });
+                              MeshLib::idsComparator<MeshLib::Node*>);
     BaseLib::makeVectorUnique(nonlinear_nodes,
-                              [](MeshLib::Node* a, MeshLib::Node* b)
-                              { return a->getID() < b->getID(); });
+                              MeshLib::idsComparator<MeshLib::Node*>);
 
     std::vector<MeshLib::Node*>& allnodes =
         const_cast<std::vector<MeshLib::Node*>&>(mesh.getNodes());
diff --git a/ChemistryLib/CreateChemicalSolverInterface.cpp b/ChemistryLib/CreateChemicalSolverInterface.cpp
index 4cd871b5ba6..9da6a484ab5 100644
--- a/ChemistryLib/CreateChemicalSolverInterface.cpp
+++ b/ChemistryLib/CreateChemicalSolverInterface.cpp
@@ -158,7 +158,7 @@ createChemicalSolverInterface<ChemicalSolver::PhreeqcKernel>(
     auto const ls_name =
         //! \ogs_file_param{prj__chemical_system__linear_solver}
         config.getConfigParameter<std::string>("linear_solver");
-    auto& linear_solver = BaseLib::getOrError(
+    auto const& linear_solver = BaseLib::getOrError(
         linear_solvers, ls_name,
         "A linear solver with the given name does not exist.");
 
diff --git a/ChemistryLib/PhreeqcIO.cpp b/ChemistryLib/PhreeqcIO.cpp
index 05768035a42..115649588e0 100644
--- a/ChemistryLib/PhreeqcIO.cpp
+++ b/ChemistryLib/PhreeqcIO.cpp
@@ -636,7 +636,7 @@ std::ostream& operator<<(std::ostream& os, PhreeqcIO const& phreeqc_io)
     return os;
 }
 
-void PhreeqcIO::callPhreeqc()
+void PhreeqcIO::callPhreeqc() const
 {
     INFO("Phreeqc: Executing chemical calculation.");
     if (RunFile(phreeqc_instance_id, _phreeqc_input_file.c_str()) != IPQ_OK)
@@ -826,11 +826,12 @@ std::istream& operator>>(std::istream& in, PhreeqcIO& phreeqc_io)
                     auto const& secondary_variables =
                         user_punch->secondary_variables;
                     // Update values of secondary variables
-                    auto& secondary_variable = BaseLib::findElementOrError(
-                        secondary_variables.begin(), secondary_variables.end(),
-                        compare_by_name,
-                        "Could not find secondary variable '" + item_name +
-                            "'.");
+                    auto const& secondary_variable =
+                        BaseLib::findElementOrError(
+                            secondary_variables.begin(),
+                            secondary_variables.end(), compare_by_name,
+                            "Could not find secondary variable '" + item_name +
+                                "'.");
                     (*secondary_variable.value)[chemical_system_id] =
                         accepted_items[item_id];
                     break;
@@ -896,13 +897,14 @@ void PhreeqcIO::computeSecondaryVariable(
     std::size_t const ele_id,
     std::vector<GlobalIndexType> const& chemical_system_indices)
 {
-    for (auto& kinetic_reactant : _chemical_system->kinetic_reactants)
+    for (auto const& kinetic_reactant : _chemical_system->kinetic_reactants)
     {
         (*kinetic_reactant.mesh_prop_molality)[ele_id] =
             averageReactantMolality(kinetic_reactant, chemical_system_indices);
     }
 
-    for (auto& equilibrium_reactant : _chemical_system->equilibrium_reactants)
+    for (auto const& equilibrium_reactant :
+         _chemical_system->equilibrium_reactants)
     {
         (*equilibrium_reactant.mesh_prop_molality)[ele_id] =
             averageReactantMolality(equilibrium_reactant,
diff --git a/ChemistryLib/PhreeqcIO.h b/ChemistryLib/PhreeqcIO.h
index c97befcbb43..a33e18acc47 100644
--- a/ChemistryLib/PhreeqcIO.h
+++ b/ChemistryLib/PhreeqcIO.h
@@ -97,7 +97,7 @@ class PhreeqcIO final : public ChemicalSolverInterface
 private:
     void writeInputsToFile(double const dt);
 
-    void callPhreeqc();
+    void callPhreeqc() const;
 
     void readOutputsFromFile();
 
diff --git a/ChemistryLib/PhreeqcKernel.cpp b/ChemistryLib/PhreeqcKernel.cpp
index 5a475800ff3..614b7a912a9 100644
--- a/ChemistryLib/PhreeqcKernel.cpp
+++ b/ChemistryLib/PhreeqcKernel.cpp
@@ -192,6 +192,7 @@ static bool isHydrogen(std::string_view const element)
 void PhreeqcKernel::setAqueousSolutions(
     std::vector<GlobalVector*> const& process_solutions)
 {
+    assert(!process_solutions.empty());
     // Loop over chemical systems
     std::size_t const num_chemical_systems = process_solutions[0]->size();
     for (std::size_t chemical_system_id = 0;
diff --git a/ChemistryLib/PhreeqcKernelData/CreateKineticReactant.cpp b/ChemistryLib/PhreeqcKernelData/CreateKineticReactant.cpp
index 543d1599484..8a1ee4c640b 100644
--- a/ChemistryLib/PhreeqcKernelData/CreateKineticReactant.cpp
+++ b/ChemistryLib/PhreeqcKernelData/CreateKineticReactant.cpp
@@ -43,7 +43,7 @@ std::unique_ptr<Kinetics> createKineticReactants(
             //! \ogs_file_param{prj__chemical_system__kinetic_reactants__kinetic_reactant__initial_amount}
             reactant_config.getConfigParameter<double>("initial_amount");
 
-        auto amount = MeshLib::getOrCreateMeshProperty<double>(
+        auto* const amount = MeshLib::getOrCreateMeshProperty<double>(
             const_cast<MeshLib::Mesh&>(mesh),
             name,
             MeshLib::MeshItemType::Node,
diff --git a/GeoLib/GEOObjects.cpp b/GeoLib/GEOObjects.cpp
index bad7ec9a252..ad16d68e75d 100644
--- a/GeoLib/GEOObjects.cpp
+++ b/GeoLib/GEOObjects.cpp
@@ -222,9 +222,9 @@ const PolylineVec* GEOObjects::getPolylineVecObj(const std::string& name) const
 bool GEOObjects::removePolylineVec(std::string const& name)
 {
     _callbacks->removePolylineVec(name);
-    auto it =
-        std::find_if(_ply_vecs.begin(), _ply_vecs.end(),
-                     [&name](PolylineVec* v) { return v->getName() == name; });
+    auto it = std::find_if(_ply_vecs.begin(), _ply_vecs.end(),
+                           [&name](PolylineVec const* const v)
+                           { return v->getName() == name; });
     if (it != _ply_vecs.end())
     {
         delete *it;
@@ -342,7 +342,7 @@ void GEOObjects::getStationVectorNames(std::vector<std::string>& names) const
 std::vector<std::string> GEOObjects::getGeometryNames() const
 {
     std::vector<std::string> names;
-    for (auto point : _pnt_vecs)
+    for (auto const point : _pnt_vecs)
     {
         if (point->getType() == PointVec::PointType::POINT)
         {
@@ -385,10 +385,7 @@ int GEOObjects::mergeGeometries(std::vector<std::string> const& geo_names,
 
     std::vector<std::size_t> pnt_offsets(n_geo_names, 0);
 
-    if (!mergePoints(geo_names, merged_geo_name, pnt_offsets))
-    {
-        return 1;
-    }
+    mergePoints(geo_names, merged_geo_name, pnt_offsets);
 
     mergePolylines(geo_names, merged_geo_name, pnt_offsets);
 
@@ -397,7 +394,7 @@ int GEOObjects::mergeGeometries(std::vector<std::string> const& geo_names,
     return 0;
 }
 
-bool GEOObjects::mergePoints(std::vector<std::string> const& geo_names,
+void GEOObjects::mergePoints(std::vector<std::string> const& geo_names,
                              std::string& merged_geo_name,
                              std::vector<std::size_t>& pnt_offsets)
 {
@@ -442,7 +439,6 @@ bool GEOObjects::mergePoints(std::vector<std::string> const& geo_names,
 
     addPointVec(std::move(merged_points), merged_geo_name,
                 std::move(merged_pnt_names), 1e-6);
-    return true;
 }
 
 void GEOObjects::mergePolylines(std::vector<std::string> const& geo_names,
diff --git a/GeoLib/GEOObjects.h b/GeoLib/GEOObjects.h
index 39eb3ea762a..a132ef82062 100644
--- a/GeoLib/GEOObjects.h
+++ b/GeoLib/GEOObjects.h
@@ -323,9 +323,8 @@ class GEOObjects final
      * @param merged_geo_name The (new) name of the geometry resulting from
      * merging.
      * @param pnt_offsets offsets in the merged vector storing the points
-     * @return true, if merging the points succeeded, else false
      */
-    bool mergePoints(std::vector<std::string> const& geo_names,
+    void mergePoints(std::vector<std::string> const& geo_names,
                      std::string& merged_geo_name,
                      std::vector<std::size_t>& pnt_offsets);
 
diff --git a/GeoLib/Grid.h b/GeoLib/Grid.h
index 1660a9a5b8d..cf6a4b7f5d0 100644
--- a/GeoLib/Grid.h
+++ b/GeoLib/Grid.h
@@ -204,10 +204,9 @@ Grid<POINT>::Grid(InputIterator first, InputIterator last,
         {max[0] - min[0], max[1] - min[1], max[2] - min[2]}};
 
     // enlarge delta
-    for (auto& d : delta)
-    {
-        d = std::nextafter(d, std::numeric_limits<double>::max());
-    }
+    constexpr double direction = std::numeric_limits<double>::max();
+    std::transform(begin(delta), end(delta), begin(delta),
+                   [](double const d) { return std::nextafter(d, direction); });
 
     assert(n_pnts > 0);
     initNumberOfSteps(max_num_per_grid_cell, static_cast<std::size_t>(n_pnts),
diff --git a/GeoLib/Polygon.cpp b/GeoLib/Polygon.cpp
index 583a8eb20b9..d7569bd63c6 100644
--- a/GeoLib/Polygon.cpp
+++ b/GeoLib/Polygon.cpp
@@ -380,7 +380,7 @@ bool Polygon::getNextIntersectionPointPolygonLine(
     }
     else
     {
-        for (auto polygon : _simple_polygon_list)
+        for (auto const* polygon : _simple_polygon_list)
         {
             for (auto seg_it(polygon->begin()); seg_it != polygon->end();
                  ++seg_it)
diff --git a/GeoLib/Polyline.cpp b/GeoLib/Polyline.cpp
index c41414bc201..a4ebcbc51ab 100644
--- a/GeoLib/Polyline.cpp
+++ b/GeoLib/Polyline.cpp
@@ -170,13 +170,6 @@ LineSegment Polyline::getSegment(std::size_t i) const
                        _ply_pnts[_ply_pnt_ids[i + 1]], false);
 }
 
-LineSegment Polyline::getSegment(std::size_t i)
-{
-    assert(i < getNumberOfSegments());
-    return LineSegment(_ply_pnts[_ply_pnt_ids[i]],
-                       _ply_pnts[_ply_pnt_ids[i + 1]], false);
-}
-
 void Polyline::setPointID(std::size_t idx, std::size_t id)
 {
     assert(idx < _ply_pnt_ids.size());
diff --git a/GeoLib/Polyline.h b/GeoLib/Polyline.h
index 60790b717f7..5aaf9a4ef09 100644
--- a/GeoLib/Polyline.h
+++ b/GeoLib/Polyline.h
@@ -222,7 +222,6 @@ class Polyline : public GeoObject
     std::vector<std::size_t> _ply_pnt_ids;
 
     LineSegment getSegment(std::size_t i) const;
-    LineSegment getSegment(std::size_t i);
 };
 
 bool containsEdge (const Polyline& ply, std::size_t id0, std::size_t id1);
diff --git a/GeoLib/QuadTree.h b/GeoLib/QuadTree.h
index 520d837c048..b1504c02e04 100644
--- a/GeoLib/QuadTree.h
+++ b/GeoLib/QuadTree.h
@@ -115,15 +115,11 @@ template <typename POINT> class QuadTree
             return false;
         }
 
-        if (!_is_leaf) {
-            for (auto& child : _children)
-            {
-                if (child->addPoint(pnt))
-                {
-                    return true;
-                }
-            }
-            return false;
+        if (!_is_leaf)
+        {
+            return std::any_of(begin(_children), end(_children),
+                               [&pnt](auto* child)
+                               { return child->addPoint(pnt); });
         }
 
         // check if point is already in quadtree
@@ -540,7 +536,7 @@ template <typename POINT> class QuadTree
         _is_leaf = false;
     }
 
-    static bool needToRefine (QuadTree<POINT>* node)
+    static bool needToRefine(QuadTree<POINT> const* const node)
     {
         QuadTree<POINT>* north_neighbor (node->getNorthNeighbor ());
         if (north_neighbor != nullptr)
diff --git a/GeoLib/Station.cpp b/GeoLib/Station.cpp
index 374357aeb9c..8bdc51596bc 100644
--- a/GeoLib/Station.cpp
+++ b/GeoLib/Station.cpp
@@ -22,12 +22,13 @@
 
 namespace GeoLib
 {
-Station::Station(double x, double y, double z, std::string name)
+Station::Station(double const x, double const y, double const z,
+                 std::string name)
     : Point(x, y, z), _name(std::move(name))
 {
 }
 
-Station::Station(Point* coords, std::string name)
+Station::Station(Point const* const coords, std::string name)
     : Point(*coords), _name(std::move(name))
 {
 }
diff --git a/GeoLib/Station.h b/GeoLib/Station.h
index b0a3a076138..a4b51060bf4 100644
--- a/GeoLib/Station.h
+++ b/GeoLib/Station.h
@@ -48,7 +48,7 @@ class Station : public Point
     explicit Station(double x = 0.0, double y = 0.0, double z = 0.0,
                      std::string name = "");
 
-    explicit Station(Point* coords, std::string name = "");
+    explicit Station(Point const* coords, std::string name = "");
 
     /**
      * Constructor copies the source object
diff --git a/MaterialLib/MPL/CheckMaterialSpatialDistributionMap.h b/MaterialLib/MPL/CheckMaterialSpatialDistributionMap.h
index f3d78a977fc..7d4ec3b2a94 100644
--- a/MaterialLib/MPL/CheckMaterialSpatialDistributionMap.h
+++ b/MaterialLib/MPL/CheckMaterialSpatialDistributionMap.h
@@ -28,10 +28,8 @@ void checkMaterialSpatialDistributionMap(
     ContainerLiquid const& required_properties_liquid_phase,
     ContainerGas const& required_properties_gas_phase)
 {
-    for (auto const& element : mesh.getElements())
+    for (auto const element_id : mesh.getElements() | MeshLib::views::ids)
     {
-        auto const element_id = element->getID();
-
         auto const& medium = *media_map.getMedium(element_id);
         if (!required_properties_medium.empty())
         {
diff --git a/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp b/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp
index 99a45ba589b..7d256940f31 100644
--- a/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp
+++ b/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp
@@ -120,10 +120,9 @@ PermeabilityMohrCoulombFailureIndexModel<DisplacementDim>::value(
     if (f >= 1.0)
     {
         const double exp_value = std::exp(b_ * f);
-        for (auto& k_i : k_data)
-        {
-            k_i = std::min(k_i + kr_ * exp_value, k_max_);
-        }
+        std::transform(begin(k_data), end(k_data), begin(k_data),
+                       [&](double const k_i)
+                       { return std::min(k_i + kr_ * exp_value, k_max_); });
     }
 
     // Local coordinate transformation is only applied for the case that the
diff --git a/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp b/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp
index dfd4a138f26..666db44232a 100644
--- a/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp
+++ b/MeshGeoToolsLib/AppendLinesAlongPolyline.cpp
@@ -45,7 +45,7 @@ std::unique_ptr<MeshLib::Mesh> appendLinesAlongPolylines(
     // for each polyline
     for (std::size_t k(0); k < n_ply; k++)
     {
-        auto* ply = ply_vec.getVector()[k];
+        auto const* const ply = ply_vec.getVector()[k];
 
         // search nodes on the polyline
         MeshGeoToolsLib::MeshNodesAlongPolyline mshNodesAlongPoly(
diff --git a/MeshGeoToolsLib/MeshNodesAlongSurface.cpp b/MeshGeoToolsLib/MeshNodesAlongSurface.cpp
index d9a2223845d..f4c9fed2d00 100644
--- a/MeshGeoToolsLib/MeshNodesAlongSurface.cpp
+++ b/MeshGeoToolsLib/MeshNodesAlongSurface.cpp
@@ -35,7 +35,7 @@ MeshNodesAlongSurface::MeshNodesAlongSurface(MeshLib::Mesh const& mesh,
     // loop over all nodes
     for (std::size_t i = 0; i < n_nodes; i++)
     {
-        auto* node = mesh_nodes[i];
+        auto const* const node = mesh_nodes[i];
         if (!sfc.isPntInBoundingVolume(*node, epsilon_radius))
         {
             continue;
diff --git a/MeshLib/Elements/TemplateElement-impl.h b/MeshLib/Elements/TemplateElement-impl.h
index c0aa0f8fd06..3f45956b32f 100644
--- a/MeshLib/Elements/TemplateElement-impl.h
+++ b/MeshLib/Elements/TemplateElement-impl.h
@@ -19,8 +19,8 @@ TemplateElement<ELEMENT_RULE>::TemplateElement(Node* nodes[n_all_nodes],
 {
     std::copy_n(nodes, n_all_nodes, std::begin(_nodes));
     delete[] nodes;
-    this->_neighbors = new Element*[getNumberOfNeighbors()];
-    std::fill(this->_neighbors, this->_neighbors + getNumberOfNeighbors(),
+    this->_neighbors = new Element*[ELEMENT_RULE::n_neighbors];
+    std::fill(this->_neighbors, this->_neighbors + ELEMENT_RULE::n_neighbors,
               nullptr);
 
     this->space_dimension_ = ELEMENT_RULE::dimension;
@@ -31,8 +31,9 @@ TemplateElement<ELEMENT_RULE>::TemplateElement(
     std::array<Node*, n_all_nodes> const& nodes, std::size_t id)
     : Element(id), _nodes{nodes}
 {
-    this->_neighbors = new Element*[getNumberOfNeighbors()];
-    std::fill(this->_neighbors, this->_neighbors + getNumberOfNeighbors(), nullptr);
+    this->_neighbors = new Element*[ELEMENT_RULE::n_neighbors];
+    std::fill(this->_neighbors, this->_neighbors + ELEMENT_RULE::n_neighbors,
+              nullptr);
 
     this->space_dimension_ = ELEMENT_RULE::dimension;
 }
@@ -42,8 +43,8 @@ TemplateElement<ELEMENT_RULE>::TemplateElement(
     TemplateElement<ELEMENT_RULE> const& e)
     : Element(e.getID()), _nodes{e._nodes}
 {
-    this->_neighbors = new Element*[getNumberOfNeighbors()];
-    for (unsigned i = 0; i < getNumberOfNeighbors(); i++)
+    this->_neighbors = new Element*[ELEMENT_RULE::n_neighbors];
+    for (unsigned i = 0; i < ELEMENT_RULE::n_neighbors; i++)
     {
         this->_neighbors[i] = e._neighbors[i];
     }
diff --git a/MeshLib/Elements/Utils.h b/MeshLib/Elements/Utils.h
index 68804dd7004..309c2de0846 100644
--- a/MeshLib/Elements/Utils.h
+++ b/MeshLib/Elements/Utils.h
@@ -34,9 +34,7 @@ inline std::vector<Node*> getBaseNodes(std::vector<Element*> const& elements)
                   std::back_inserter(base_nodes));
     }
 
-    BaseLib::makeVectorUnique(base_nodes, [](Node const* a, Node* b) {
-        return a->getID() < b->getID();
-    });
+    BaseLib::makeVectorUnique(base_nodes, MeshLib::idsComparator<Node*>);
 
     return base_nodes;
 }
diff --git a/MeshLib/Mesh.cpp b/MeshLib/Mesh.cpp
index 3e644ed03ac..ac05a6129fa 100644
--- a/MeshLib/Mesh.cpp
+++ b/MeshLib/Mesh.cpp
@@ -186,7 +186,7 @@ std::pair<double, double> minMaxEdgeLength(
     using limits = std::numeric_limits<double>;
     auto const bounds = ranges::accumulate(
         elements, std::pair{limits::infinity(), -limits::infinity()}, min_max,
-        [](Element* const e) { return computeSqrEdgeLengthRange(*e); });
+        [](Element const* const e) { return computeSqrEdgeLengthRange(*e); });
 
     return {std::sqrt(bounds.first), std::sqrt(bounds.second)};
 }
@@ -325,7 +325,7 @@ std::vector<std::vector<Node*>> calculateNodesConnectedByElements(
         // Make nodes unique and sorted by their ids.
         // This relies on the node's id being equivalent to it's address.
         std::sort(adjacent_nodes.begin(), adjacent_nodes.end(),
-                  [](Node* a, Node* b) { return a->getID() < b->getID(); });
+                  idsComparator<Node*>);
         auto const last =
             std::unique(adjacent_nodes.begin(), adjacent_nodes.end());
         adjacent_nodes.erase(last, adjacent_nodes.end());
diff --git a/MeshLib/Mesh.h b/MeshLib/Mesh.h
index 871a036f060..807561b4085 100644
--- a/MeshLib/Mesh.h
+++ b/MeshLib/Mesh.h
@@ -195,6 +195,21 @@ bool isBaseNode(Node const& node,
 std::pair<double, double> minMaxEdgeLength(
     std::vector<Element*> const& elements);
 
+/// Lexicographic comparison of ids of two objects of type T. T can be a pointer
+/// or a value type.
+template <typename T>
+bool idsComparator(T const a, T const b)
+{
+    if constexpr (std::is_pointer_v<T>)
+    {
+        return a->getID() < b->getID();
+    }
+    else
+    {
+        return a.getID() < b.getID();
+    }
+}
+
 /// MeshLib specific, lazy, non-owning, non-mutating, composable range views.
 namespace views
 {
diff --git a/MeshLib/MeshSearch/ElementSearch.cpp b/MeshLib/MeshSearch/ElementSearch.cpp
index 7df4438e3c0..4e85a2fc404 100644
--- a/MeshLib/MeshSearch/ElementSearch.cpp
+++ b/MeshLib/MeshSearch/ElementSearch.cpp
@@ -10,6 +10,11 @@
 
 #include "ElementSearch.h"
 
+#include <range/v3/algorithm/any_of.hpp>
+#include <range/v3/range/conversion.hpp>
+#include <range/v3/view/filter.hpp>
+#include <range/v3/view/take.hpp>
+
 #include "BaseLib/Algorithm.h"
 #include "MeshLib/Elements/Element.h"
 #include "MeshLib/Node.h"
@@ -21,22 +26,13 @@ ElementSearch::ElementSearch(const MeshLib::Mesh& mesh) : _mesh(mesh) {}
 template <typename Container, typename Predicate>
 std::vector<std::size_t> filter(Container const& container, Predicate const& p)
 {
-    std::vector<std::size_t> matchedIDs;
-    std::size_t i = 0;
-    for (auto value : container)
-    {
-        if (p(value))
-        {
-            matchedIDs.push_back(i);
-        }
-        i++;
-    }
-    return matchedIDs;
+    return ranges::views::filter(container, p) | views::ids |
+           ranges::to<std::vector>;
 }
 
 std::size_t ElementSearch::searchByElementType(MeshElemType eleType)
 {
-    auto matchedIDs = filter(_mesh.getElements(), [&](MeshLib::Element* e)
+    auto matchedIDs = filter(_mesh.getElements(), [&](MeshLib::Element const* e)
                              { return e->getGeomType() == eleType; });
 
     this->updateUnion(matchedIDs);
@@ -45,8 +41,9 @@ std::size_t ElementSearch::searchByElementType(MeshElemType eleType)
 
 std::size_t ElementSearch::searchByContent(double eps)
 {
-    auto matchedIDs = filter(_mesh.getElements(), [&eps](MeshLib::Element* e)
-                             { return e->getContent() < eps; });
+    auto matchedIDs =
+        filter(_mesh.getElements(), [&eps](MeshLib::Element const* e)
+               { return e->getContent() < eps; });
 
     this->updateUnion(matchedIDs);
     return matchedIDs.size();
@@ -54,20 +51,15 @@ std::size_t ElementSearch::searchByContent(double eps)
 
 std::size_t ElementSearch::searchByBoundingBox(GeoLib::AABB const& aabb)
 {
-    auto matchedIDs =
-        filter(_mesh.getElements(),
-               [&aabb](MeshLib::Element* e)
-               {
-                   std::size_t const nElemNodes(e->getNumberOfBaseNodes());
-                   for (std::size_t n = 0; n < nElemNodes; ++n)
-                   {
-                       if (aabb.containsPoint(*e->getNode(n), 0))
-                       {
-                           return true;  // any node of element is in aabb.
-                       }
-                   }
-                   return false;  // no nodes of element are in aabb.
-               });
+    auto matchedIDs = filter(
+        _mesh.getElements(),
+        [&aabb](MeshLib::Element const* e)
+        {
+            // any node of element is in aabb.
+            return ranges::any_of(
+                e->nodes() | ranges::views::take(e->getNumberOfBaseNodes()),
+                [&aabb](auto const* n) { return aabb.containsPoint(*n, 0); });
+        });
 
     this->updateUnion(matchedIDs);
     return matchedIDs.size();
diff --git a/MeshLib/Utils/DuplicateMeshComponents.cpp b/MeshLib/Utils/DuplicateMeshComponents.cpp
index 0f17843ed00..2d58d886f80 100644
--- a/MeshLib/Utils/DuplicateMeshComponents.cpp
+++ b/MeshLib/Utils/DuplicateMeshComponents.cpp
@@ -120,7 +120,7 @@ std::vector<Element*> cloneElements(std::vector<Element*> const& elements)
     cloned_elements.reserve(elements.size());
     std::transform(begin(elements), end(elements),
                    std::back_inserter(cloned_elements),
-                   [](Element* const e) { return e->clone(); });
+                   [](Element const* const e) { return e->clone(); });
     return cloned_elements;
 }
 
diff --git a/MeshLib/Utils/GetElementRotationMatrices.cpp b/MeshLib/Utils/GetElementRotationMatrices.cpp
index 0c07fa3fc33..a117a731843 100644
--- a/MeshLib/Utils/GetElementRotationMatrices.cpp
+++ b/MeshLib/Utils/GetElementRotationMatrices.cpp
@@ -24,7 +24,7 @@ std::vector<Eigen::MatrixXd> getElementRotationMatrices(
 {
     std::vector<Eigen::MatrixXd> element_rotation_matrices;
     element_rotation_matrices.reserve(elements.size());
-    for (auto const& element : elements)
+    for (auto const* const element : elements)
     {
         int const element_dimension = static_cast<int>(element->getDimension());
         if (element_dimension == space_dimension)
diff --git a/MeshToolsLib/MeshGenerators/VoxelGridFromMesh.cpp b/MeshToolsLib/MeshGenerators/VoxelGridFromMesh.cpp
index d6dcfe2c4c7..5cfdec7d22d 100644
--- a/MeshToolsLib/MeshGenerators/VoxelGridFromMesh.cpp
+++ b/MeshToolsLib/MeshGenerators/VoxelGridFromMesh.cpp
@@ -146,7 +146,7 @@ bool checkDyncast(MeshLib::Mesh& mesh,
 }
 
 void mapMeshArraysOntoGrid(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
-                           std::unique_ptr<MeshLib::Mesh>& grid)
+                           std::unique_ptr<MeshLib::Mesh> const& grid)
 {
     assert(mesh != nullptr);
     assert(grid != nullptr);
@@ -166,4 +166,4 @@ void mapMeshArraysOntoGrid(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
         }
     }
 }
-}  // namespace MeshToolsLib::MeshGenerator::VoxelGridFromMesh
\ No newline at end of file
+}  // namespace MeshToolsLib::MeshGenerator::VoxelGridFromMesh
diff --git a/MeshToolsLib/MeshGenerators/VoxelGridFromMesh.h b/MeshToolsLib/MeshGenerators/VoxelGridFromMesh.h
index ea43a90c1c1..55a61f8aa6b 100644
--- a/MeshToolsLib/MeshGenerators/VoxelGridFromMesh.h
+++ b/MeshToolsLib/MeshGenerators/VoxelGridFromMesh.h
@@ -48,5 +48,5 @@ bool removeUnusedGridCells(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
                            std::unique_ptr<MeshLib::Mesh>& grid);
 // map the cell data of mesh to voxelgrid
 void mapMeshArraysOntoGrid(vtkSmartPointer<vtkUnstructuredGrid> const& mesh,
-                           std::unique_ptr<MeshLib::Mesh>& grid);
-}  // namespace MeshToolsLib::MeshGenerator::VoxelGridFromMesh
\ No newline at end of file
+                           std::unique_ptr<MeshLib::Mesh> const& grid);
+}  // namespace MeshToolsLib::MeshGenerator::VoxelGridFromMesh
diff --git a/MeshToolsLib/MeshQuality/MeshValidation.cpp b/MeshToolsLib/MeshQuality/MeshValidation.cpp
index 7381bacc4f9..3db81708983 100644
--- a/MeshToolsLib/MeshQuality/MeshValidation.cpp
+++ b/MeshToolsLib/MeshQuality/MeshValidation.cpp
@@ -171,7 +171,7 @@ MeshValidation::ElementErrorCodeOutput(
 {
     const auto nErrorFlags(
         static_cast<std::size_t>(ElementErrorFlag::MaxValue));
-    ElementErrorFlag flags[nErrorFlags] = {
+    const ElementErrorFlag flags[nErrorFlags] = {
         ElementErrorFlag::ZeroVolume, ElementErrorFlag::NonCoplanar,
         ElementErrorFlag::NonConvex, ElementErrorFlag::NodeOrder};
     const std::size_t nElements(error_codes.size());
diff --git a/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
index c030be70292..7235edd5f51 100644
--- a/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
+++ b/NumLib/TimeStepping/Algorithms/FixedTimeStepping.cpp
@@ -171,7 +171,7 @@ FixedTimeStepping::FixedTimeStepping(
         auto const delta_t = repeat_dt_pairs.back().second;
         auto const repeat =
             static_cast<std::size_t>(std::ceil((_t_end - t_curr) / delta_t));
-        t_curr = addTimeIncrement(_dt_vector, repeat, delta_t, t_curr);
+        addTimeIncrement(_dt_vector, repeat, delta_t, t_curr);
     }
 
     incorporateFixedTimesForOutput(_t_initial, _t_end, _dt_vector,
diff --git a/ProcessLib/BoundaryConditionAndSourceTerm/ConstraintDirichletBoundaryCondition.cpp b/ProcessLib/BoundaryConditionAndSourceTerm/ConstraintDirichletBoundaryCondition.cpp
index b9593bd78b8..acf8af3f23b 100644
--- a/ProcessLib/BoundaryConditionAndSourceTerm/ConstraintDirichletBoundaryCondition.cpp
+++ b/ProcessLib/BoundaryConditionAndSourceTerm/ConstraintDirichletBoundaryCondition.cpp
@@ -114,15 +114,13 @@ void ConstraintDirichletBoundaryCondition::preTimestep(
     DBUG(
         "ConstraintDirichletBoundaryCondition::preTimestep: computing flux "
         "constraints");
-    for (auto const* boundary_element : _bc_mesh.getElements())
+    for (auto const id : _bc_mesh.getElements() | MeshLib::views::ids)
     {
-        _flux_values[boundary_element->getID()] =
-            _local_assemblers[boundary_element->getID()]->integrate(
-                x, t,
-                [this](std::size_t const element_id,
-                       MathLib::Point3d const& pnt, double const t,
-                       std::vector<GlobalVector*> const& x)
-                { return _getFlux(element_id, pnt, t, x); });
+        _flux_values[id] = _local_assemblers[id]->integrate(
+            x, t,
+            [this](std::size_t const element_id, MathLib::Point3d const& pnt,
+                   double const t, std::vector<GlobalVector*> const& x)
+            { return _getFlux(element_id, pnt, t, x); });
     }
 }
 
diff --git a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
index 5446d3f9112..022adcac2c3 100644
--- a/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
+++ b/ProcessLib/ComponentTransport/CreateComponentTransportProcess.cpp
@@ -47,10 +47,8 @@ void checkMPLProperties(
         MaterialPropertyLib::PropertyType::decay_rate,
         MaterialPropertyLib::PropertyType::pore_diffusion};
 
-    for (auto const& element : mesh.getElements())
+    for (auto const& element_id : mesh.getElements() | MeshLib::views::ids)
     {
-        auto const element_id = element->getID();
-
         auto const& medium = *media_map.getMedium(element_id);
         checkRequiredProperties(medium, required_properties_medium);
 
diff --git a/ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h b/ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h
index 2ecd9bad534..2a8be7faba1 100644
--- a/ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h
+++ b/ProcessLib/Deformation/SolidMaterialInternalToSecondaryVariables.h
@@ -137,7 +137,7 @@ template <typename InternalVariable>
 class CollectIntegrationPointDataForExtrapolation final
 {
 public:
-    CollectIntegrationPointDataForExtrapolation(
+    explicit CollectIntegrationPointDataForExtrapolation(
         InternalVariablesCollection<InternalVariable>&&
             internal_variables_collection)
         : internal_variables_collection_(
@@ -228,7 +228,7 @@ template <typename InternalVariable>
 class CollectIntegrationPointDataForIpWriter final
 {
 public:
-    CollectIntegrationPointDataForIpWriter(
+    explicit CollectIntegrationPointDataForIpWriter(
         InternalVariablesCollection<InternalVariable>&&
             internal_variables_collection)
         : internal_variables_collection_(internal_variables_collection)
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
index 98b24d5c432..7ab35cdf540 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.cpp
@@ -104,7 +104,7 @@ BHECommonCoaxial::pipeAdvectionVectors(
 
 std::array<double, BHECommonCoaxial::number_of_unknowns>
 BHECommonCoaxial::calcThermalResistances(double const Nu_inner_pipe,
-                                         double const Nu_annulus_pipe)
+                                         double const Nu_annulus_pipe) const
 {
     // thermal resistances due to advective flow of refrigerant in the pipes
     auto const R_advective = calculateAdvectiveThermalResistance(
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
index ea4fbf1dde0..e1876563a1c 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BHECommonCoaxial.h
@@ -43,7 +43,7 @@ class BHECommonCoaxial : public BHECommon
     double updateFlowRateAndTemperature(double T_out, double current_time);
 
     std::array<double, number_of_unknowns> calcThermalResistances(
-        double const Nu_inner_pipe, double const Nu_annulus_pipe);
+        double const Nu_inner_pipe, double const Nu_annulus_pipe) const;
 
     std::array<double, number_of_unknowns> pipeHeatCapacities() const;
 
diff --git a/ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp b/ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp
index 2b460264534..1ca62da3133 100644
--- a/ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/MeshUtils.cpp
@@ -24,9 +24,9 @@ std::vector<MeshLib::Element*> extractOneDimensionalElements(
 {
     std::vector<MeshLib::Element*> one_dimensional_elements;
 
-    copy_if(begin(elements), end(elements),
-            back_inserter(one_dimensional_elements),
-            [](MeshLib::Element* e) { return e->getDimension() == 1; });
+    copy_if(
+        begin(elements), end(elements), back_inserter(one_dimensional_elements),
+        [](MeshLib::Element const* const e) { return e->getDimension() == 1; });
 
     return one_dimensional_elements;
 }
@@ -80,7 +80,7 @@ BHEMeshData getBHEDataInMesh(MeshLib::Mesh const& mesh)
         std::vector<MeshLib::Element*>& vec_elements = bhe_elements[bhe_id];
         copy_if(begin(all_bhe_elements), end(all_bhe_elements),
                 back_inserter(vec_elements),
-                [&](MeshLib::Element* e)
+                [&](MeshLib::Element const* const e)
                 { return material_ids[e->getID()] == bhe_mat_id; });
         DBUG("-> found {:d} elements on the BHE_{:d}", vec_elements.size(),
              bhe_id);
@@ -100,8 +100,8 @@ BHEMeshData getBHEDataInMesh(MeshLib::Mesh const& mesh)
             }
         }
         BaseLib::makeVectorUnique(vec_nodes,
-                                  [](MeshLib::Node* node1, MeshLib::Node* node2)
-                                  { return node1->getID() < node2->getID(); });
+                                  MeshLib::idsComparator<MeshLib::Node*>);
+
         DBUG("-> found {:d} nodes on the BHE_{:d}", vec_nodes.size(), bhe_id);
     }
 
diff --git a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
index d69ac4a5f40..a40c37308ae 100644
--- a/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
+++ b/ProcessLib/HydroMechanics/CreateHydroMechanicsProcess.cpp
@@ -214,9 +214,8 @@ std::unique_ptr<Process> createHydroMechanicsProcess(
                                                 MaterialPropertyLib::density};
     std::array const requiredSolidProperties = {MaterialPropertyLib::density};
 
-    for (auto const& element : mesh.getElements())
+    for (auto const& element_id : mesh.getElements() | MeshLib::views::ids)
     {
-        auto const element_id = element->getID();
         media_map->checkElementHasMedium(element_id);
         auto const& medium = *media_map->getMedium(element_id);
         checkRequiredProperties(medium, requiredMediumProperties);
diff --git a/ProcessLib/LIE/Common/MeshUtils.cpp b/ProcessLib/LIE/Common/MeshUtils.cpp
index eac158b064a..1144ec87088 100644
--- a/ProcessLib/LIE/Common/MeshUtils.cpp
+++ b/ProcessLib/LIE/Common/MeshUtils.cpp
@@ -264,7 +264,7 @@ void getFractureMatrixDataInMesh(
             vec_fracture_elements[frac_id];
         std::copy_if(all_fracture_elements.begin(), all_fracture_elements.end(),
                      std::back_inserter(vec_elements),
-                     [&](MeshLib::Element* e)
+                     [&](MeshLib::Element const* const e)
                      { return (*material_ids)[e->getID()] == frac_mat_id; });
         DBUG("-> found {:d} elements on the fracture {:d}", vec_elements.size(),
              frac_id);
@@ -275,7 +275,7 @@ void getFractureMatrixDataInMesh(
     for (unsigned frac_id = 0; frac_id < vec_fracture_mat_IDs.size(); frac_id++)
     {
         std::vector<MeshLib::Node*>& vec_nodes = vec_fracture_nodes[frac_id];
-        for (MeshLib::Element* e : vec_fracture_elements[frac_id])
+        for (MeshLib::Element const* const e : vec_fracture_elements[frac_id])
         {
             for (unsigned i = 0; i < e->getNumberOfNodes(); i++)
             {
@@ -287,8 +287,7 @@ void getFractureMatrixDataInMesh(
             }
         }
         BaseLib::makeVectorUnique(vec_nodes,
-                                  [](MeshLib::Node* node1, MeshLib::Node* node2)
-                                  { return node1->getID() < node2->getID(); });
+                                  MeshLib::idsComparator<MeshLib::Node*>);
         DBUG("-> found {:d} nodes on the fracture {:d}", vec_nodes.size(),
              frac_id);
     }
@@ -307,7 +306,7 @@ void getFractureMatrixDataInMesh(
         auto const& fracture_elements = vec_fracture_elements[fid];
         std::vector<MeshLib::Element*> vec_ele;
         // first, collect matrix elements
-        for (MeshLib::Element* e : fracture_elements)
+        for (MeshLib::Element const* const e : fracture_elements)
         {
             // it is sufficient to iterate over base nodes, because they are
             // already connected to all neighbours
@@ -334,8 +333,7 @@ void getFractureMatrixDataInMesh(
             }
         }
         BaseLib::makeVectorUnique(vec_ele,
-                                  [](MeshLib::Element* e1, MeshLib::Element* e2)
-                                  { return e1->getID() < e2->getID(); });
+                                  MeshLib::idsComparator<MeshLib::Element*>);
 
         // second, append fracture elements
         std::copy(fracture_elements.begin(), fracture_elements.end(),
diff --git a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h
index a1e2da74c28..b214c225e2d 100644
--- a/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h
+++ b/ProcessLib/LIE/SmallDeformation/LocalAssembler/SmallDeformationLocalAssemblerFracture-impl.h
@@ -281,7 +281,7 @@ void SmallDeformationLocalAssemblerFracture<ShapeFunction, DisplacementDim>::
         auto& C = ip_data.C;
         auto& state = *ip_data.material_state_variables;
         auto& b_m = ip_data.aperture;
-        auto& N = _secondary_data.N[ip];
+        auto const& N = _secondary_data.N[ip];
 
         auto const ip_physical_coords(computePhysicalCoordinates(_element, N));
         std::vector<double> const levelsets(duGlobalEnrichments(
@@ -328,7 +328,7 @@ void SmallDeformationLocalAssemblerFracture<ShapeFunction, DisplacementDim>::
 
     for (unsigned ip = 0; ip < n_integration_points; ip++)
     {
-        auto& ip_data = _ip_data[ip];
+        auto const& ip_data = _ip_data[ip];
         ele_b += ip_data.aperture;
         ele_w += ip_data.w;
         ele_sigma += ip_data.sigma;
diff --git a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
index a7f644f8b99..b2f079c5fbe 100644
--- a/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
+++ b/ProcessLib/LIE/SmallDeformation/SmallDeformationProcess.cpp
@@ -97,10 +97,10 @@ SmallDeformationProcess<DisplacementDim>::SmallDeformationProcess(
     }
 
     // set branches
-    for (auto& vec_branch_nodeID_matID : vec_branch_nodeID_matIDs)
+    for (auto const& [vec_branch_nodeID, matID] : vec_branch_nodeID_matIDs)
     {
-        auto master_matId = vec_branch_nodeID_matID.second[0];
-        auto slave_matId = vec_branch_nodeID_matID.second[1];
+        auto master_matId = matID[0];
+        auto slave_matId = matID[1];
         auto& master_frac =
             _process_data.fracture_properties
                 [_process_data._map_materialID_to_fractureID[master_matId]];
@@ -108,13 +108,11 @@ SmallDeformationProcess<DisplacementDim>::SmallDeformationProcess(
             _process_data.fracture_properties
                 [_process_data._map_materialID_to_fractureID[slave_matId]];
 
-        master_frac.branches_master.push_back(
-            createBranchProperty(*mesh.getNode(vec_branch_nodeID_matID.first),
-                                 master_frac, slave_frac));
+        master_frac.branches_master.push_back(createBranchProperty(
+            *mesh.getNode(vec_branch_nodeID), master_frac, slave_frac));
 
-        slave_frac.branches_slave.push_back(
-            createBranchProperty(*mesh.getNode(vec_branch_nodeID_matID.first),
-                                 master_frac, slave_frac));
+        slave_frac.branches_slave.push_back(createBranchProperty(
+            *mesh.getNode(vec_branch_nodeID), master_frac, slave_frac));
     }
 
     // set junctions
@@ -146,10 +144,10 @@ SmallDeformationProcess<DisplacementDim>::SmallDeformationProcess(
     for (unsigned i = 0; i < vec_junction_nodeID_matIDs.size(); i++)
     {
         auto node = mesh.getNode(vec_junction_nodeID_matIDs[i].first);
-        for (auto e : mesh.getElementsConnectedToNode(*node))
+        for (auto id :
+             mesh.getElementsConnectedToNode(*node) | MeshLib::views::ids)
         {
-            _process_data._vec_ele_connected_junctionIDs[e->getID()].push_back(
-                i);
+            _process_data._vec_ele_connected_junctionIDs[id].push_back(i);
         }
     }
 
@@ -205,7 +203,7 @@ void SmallDeformationProcess<DisplacementDim>::constructDofTable()
     std::vector<MeshLib::MeshSubset> all_mesh_subsets;
     std::generate_n(std::back_inserter(all_mesh_subsets), DisplacementDim,
                     [&]() { return *_mesh_subset_matrix_nodes; });
-    for (auto& ms : _mesh_subset_fracture_nodes)
+    for (auto const& ms : _mesh_subset_fracture_nodes)
     {
         std::generate_n(std::back_inserter(all_mesh_subsets),
                         DisplacementDim,
diff --git a/ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.cpp b/ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.cpp
index 71ab59ce439..a94e9842b0c 100644
--- a/ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.cpp
+++ b/ProcessLib/RichardsComponentTransport/CreateRichardsComponentTransportProcess.cpp
@@ -46,10 +46,8 @@ void checkMPLProperties(
         MaterialPropertyLib::PropertyType::decay_rate,
         MaterialPropertyLib::PropertyType::pore_diffusion};
 
-    for (auto const& element : mesh.getElements())
+    for (auto const element_id : mesh.getElements() | MeshLib::views::ids)
     {
-        auto const element_id = element->getID();
-
         auto const& medium = *media_map.getMedium(element_id);
         checkRequiredProperties(medium, required_properties_medium);
 
diff --git a/ProcessLib/StokesFlow/CreateStokesFlowProcess.cpp b/ProcessLib/StokesFlow/CreateStokesFlowProcess.cpp
index 91cea3678a8..955011e536e 100644
--- a/ProcessLib/StokesFlow/CreateStokesFlowProcess.cpp
+++ b/ProcessLib/StokesFlow/CreateStokesFlowProcess.cpp
@@ -33,9 +33,8 @@ void checkMPLProperties(
     std::array const required_properties_medium = {
         MaterialPropertyLib::PropertyType::permeability};
 
-    for (auto const& element : mesh.getElements())
+    for (auto const element_id : mesh.getElements() | MeshLib::views::ids)
     {
-        auto const element_id = element->getID();
         auto const& medium = *media_map.getMedium(element_id);
 
         if (use_stokes_brinkman_form)
diff --git a/ProcessLib/TES/TESLocalAssemblerInner-impl.h b/ProcessLib/TES/TESLocalAssemblerInner-impl.h
index 9465f1787a3..3faa6e50452 100644
--- a/ProcessLib/TES/TESLocalAssemblerInner-impl.h
+++ b/ProcessLib/TES/TESLocalAssemblerInner-impl.h
@@ -34,7 +34,7 @@ TESLocalAssemblerInner<Traits>::TESLocalAssemblerInner(
 
 template <typename Traits>
 Eigen::Matrix3d TESLocalAssemblerInner<Traits>::getMassCoeffMatrix(
-    const unsigned int_pt)
+    const unsigned int_pt) const
 {
     // TODO: Dalton's law property
     const double dxn_dxm = Adsorption::AdsorptionReaction::dMolarFraction(
@@ -102,7 +102,7 @@ TESLocalAssemblerInner<Traits>::getLaplaceCoeffMatrix(const unsigned /*int_pt*/,
 
 template <typename Traits>
 Eigen::Matrix3d TESLocalAssemblerInner<Traits>::getAdvectionCoeffMatrix(
-    const unsigned /*int_pt*/)
+    const unsigned /*int_pt*/) const
 {
     const double A_pp = 0.0;
     const double A_pT = 0.0;
@@ -126,7 +126,7 @@ Eigen::Matrix3d TESLocalAssemblerInner<Traits>::getAdvectionCoeffMatrix(
 
 template <typename Traits>
 Eigen::Matrix3d TESLocalAssemblerInner<Traits>::getContentCoeffMatrix(
-    const unsigned /*int_pt*/)
+    const unsigned /*int_pt*/) const
 {
     const double C_pp = 0.0;
     const double C_pT = 0.0;
diff --git a/ProcessLib/TES/TESLocalAssemblerInner.h b/ProcessLib/TES/TESLocalAssemblerInner.h
index c3cb1f65038..e60b01fe361 100644
--- a/ProcessLib/TES/TESLocalAssemblerInner.h
+++ b/ProcessLib/TES/TESLocalAssemblerInner.h
@@ -51,11 +51,11 @@ class TESLocalAssemblerInner
     TESFEMReactionAdaptor& getReactionAdaptor() { return *_d.reaction_adaptor; }
     TESLocalAssemblerData const& getData() const { return _d; }
 private:
-    Eigen::Matrix3d getMassCoeffMatrix(const unsigned int_pt);
+    Eigen::Matrix3d getMassCoeffMatrix(const unsigned int_pt) const;
     typename Traits::LaplaceMatrix getLaplaceCoeffMatrix(const unsigned int_pt,
                                                          const unsigned dim);
-    Eigen::Matrix3d getAdvectionCoeffMatrix(const unsigned int_pt);
-    Eigen::Matrix3d getContentCoeffMatrix(const unsigned int_pt);
+    Eigen::Matrix3d getAdvectionCoeffMatrix(const unsigned int_pt) const;
+    Eigen::Matrix3d getContentCoeffMatrix(const unsigned int_pt) const;
     Eigen::Vector3d getRHSCoeffVector(const unsigned int_pt);
 
     void preEachAssembleIntegrationPoint(
diff --git a/ProcessLib/TES/TESOGS5MaterialModels.h b/ProcessLib/TES/TESOGS5MaterialModels.h
index 158b63f0174..7fe7ba41e5c 100644
--- a/ProcessLib/TES/TESOGS5MaterialModels.h
+++ b/ProcessLib/TES/TESOGS5MaterialModels.h
@@ -237,7 +237,7 @@ struct FluidHeatConductivityN2
                             loop1_term<2>(T) + loop1_term<3>(T) +
                             loop1_term<4>(T) + loop1_term<5>(T) +
                             loop1_term<6>(T);
-        const double temp(std::exp((f[8] / T)) - 1);
+        const double temp = std::expm1(f[8] / T);
         const double c_v0 =
             R *
             (sum1 + ((f[7] * (f[8] / T) * (f[8] / T) * (std::exp((f[8] / T)))) /
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.cpp b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.cpp
index 8266768f021..f4db9e39c0a 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.cpp
@@ -241,8 +241,8 @@ void TwoPhaseFlowWithPrhoMaterialProperties::calculateResidual(
     const double rho_h2_nonwet = pg * H2 / IdealGasConstant / T;
 
     // calculating residual
-    res(0) = calculateEquilibiumRhoWetLight(pg, Sw, rho_h2_wet);
-    res(1) = calculateSaturation(pl, X, Sw, rho_h2_wet, rho_h2_nonwet, T);
+    res[0] = calculateEquilibiumRhoWetLight(pg, Sw, rho_h2_wet);
+    res[1] = calculateSaturation(pl, X, Sw, rho_h2_wet, rho_h2_nonwet, T);
 }
 
 void TwoPhaseFlowWithPrhoMaterialProperties::calculateJacobian(
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.h b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.h
index a5294d8998a..6934b11891b 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.h
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.h
@@ -89,19 +89,18 @@ class TwoPhaseFlowWithPrhoMaterialProperties
     double getGasDensity(const double p, const double T) const;
     double getGasViscosity(const double p, const double T) const;
     double getLiquidViscosity(const double p, const double T) const;
-    bool computeConstitutiveRelation(
-        double const t,
-        ParameterLib::SpatialPosition const& x_position,
-        const int material_id,
-        double const pg,
-        double const X,
-        double const T,
-        double& Sw,
-        double& X_m,
-        double& dsw_dpg,
-        double& dsw_dX,
-        double& dxm_dpg,
-        double& dxm_dX);
+    bool computeConstitutiveRelation(double const t,
+                                     ParameterLib::SpatialPosition const& x,
+                                     const int material_id,
+                                     double const pg,
+                                     double const X,
+                                     double const T,
+                                     double& Sw,
+                                     double& X_m,
+                                     double& dsw_dpg,
+                                     double& dsw_dX,
+                                     double& dxm_dpg,
+                                     double& dxm_dX);
 
 protected:
     std::unique_ptr<MaterialLib::Fluid::FluidProperty> _liquid_density;
diff --git a/Tests/MeshLib/TestPrismQualityCriteria.cpp b/Tests/MeshLib/TestPrismQualityCriteria.cpp
index daf9cbc5fe3..09eaaea94fe 100644
--- a/Tests/MeshLib/TestPrismQualityCriteria.cpp
+++ b/Tests/MeshLib/TestPrismQualityCriteria.cpp
@@ -58,7 +58,7 @@ TEST_F(PrismElementQuality, ElementSize)
     for (auto const element_quality : element_quality_vector)
     {
         ASSERT_NEAR(expected_value, element_quality,
-                    std::numeric_limits<double>::epsilon());
+                    2 * std::numeric_limits<double>::epsilon());
     }
 }