diff --git a/Applications/DataExplorer/VtkVis/VisualizationWidget.cpp b/Applications/DataExplorer/VtkVis/VisualizationWidget.cpp
index f9906b998713c7669120123f2add200b96e1438b..e944415da9a7d9fe634dc5c5a7deda9b1dedbb92 100644
--- a/Applications/DataExplorer/VtkVis/VisualizationWidget.cpp
+++ b/Applications/DataExplorer/VtkVis/VisualizationWidget.cpp
@@ -121,7 +121,7 @@ void VisualizationWidget::showAll(int x, int y, int z)
     vtkCamera* cam = _vtkRender->GetActiveCamera();
     double* fp = cam->GetFocalPoint();
     double* p = cam->GetPosition();
-    double dist = sqrt(vtkMath::Distance2BetweenPoints(p, fp));
+    double dist = std::sqrt(vtkMath::Distance2BetweenPoints(p, fp));
     cam->SetPosition(fp[0]+(x*dist), fp[1]+(y*dist), fp[2]+(z*dist));
 
     if (x != 0 || y != 0)
diff --git a/Applications/DataExplorer/VtkVis/VtkTextureOnSurfaceFilter.cpp b/Applications/DataExplorer/VtkVis/VtkTextureOnSurfaceFilter.cpp
index d11ac19727c2f2bb05ec5da4882df2866aa6c7cf..58a4c51d44ad1f4ebde681938e9be5bba01028ee 100644
--- a/Applications/DataExplorer/VtkVis/VtkTextureOnSurfaceFilter.cpp
+++ b/Applications/DataExplorer/VtkVis/VtkTextureOnSurfaceFilter.cpp
@@ -92,9 +92,9 @@ int VtkTextureOnSurfaceFilter::RequestData( vtkInformation* request,
             points->GetPoint(i, coords);
             GeoLib::Point* pnt2 = new GeoLib::Point(coords);
             if (i<173)
-                dist += sqrt(MathLib::sqrDist(pnt, pnt2));
+                dist += std::sqrt(MathLib::sqrDist(pnt, pnt2));
             else
-                dist -= sqrt(MathLib::sqrDist(pnt, pnt2));
+                dist -= std::sqrt(MathLib::sqrDist(pnt, pnt2));
         }
         points->GetPoint(i, coords);
         double x = MathLib::normalize(0, 8404, dist);
diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index 4feae3e931a4c0ac6c573a80664369236f14f5b4..678f2db61cb63667fc46058f5f9b4bef2887276b 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -249,8 +249,9 @@ bool lineSegmentIntersect(GeoLib::LineSegment const& s0,
     GeoLib::Point const p0(a[0]+rhs[0]*v[0], a[1]+rhs[0]*v[1], a[2]+rhs[0]*v[2]);
     GeoLib::Point const p1(c[0]+rhs[1]*w[0], c[1]+rhs[1]*w[1], c[2]+rhs[1]*w[2]);
 
-    double const min_dist(sqrt(MathLib::sqrDist(p0, p1)));
-    double const min_seg_len(std::min(sqrt(sqr_len_v), sqrt(sqr_len_w)));
+    double const min_dist(std::sqrt(MathLib::sqrDist(p0, p1)));
+    double const min_seg_len(
+        std::min(std::sqrt(sqr_len_v), std::sqrt(sqr_len_w)));
     if (min_dist < min_seg_len * 1e-6) {
         s[0] = 0.5 * (p0[0] + p1[0]);
         s[1] = 0.5 * (p0[1] + p1[1]);
diff --git a/GeoLib/GEOObjects.h b/GeoLib/GEOObjects.h
index a7dc195abb1173b039a64fefc96732dd77f800c1..21edcbcc23bf2f3b0d3b6ae12c508a1ce8dd15ba 100644
--- a/GeoLib/GEOObjects.h
+++ b/GeoLib/GEOObjects.h
@@ -90,11 +90,12 @@ public:
      * @param pnt_id_name_map names corresponding to the points
      * @param eps relative tolerance value for testing of point uniqueness
      */
-    void addPointVec(std::unique_ptr<std::vector<Point*>> points,
-                     std::string& name,
-                     std::unique_ptr<std::map<std::string, std::size_t>>
-                         pnt_id_name_map = nullptr,
-                     double eps = sqrt(std::numeric_limits<double>::epsilon()));
+    void addPointVec(
+        std::unique_ptr<std::vector<Point*>> points,
+        std::string& name,
+        std::unique_ptr<std::map<std::string, std::size_t>> pnt_id_name_map =
+            nullptr,
+        double eps = std::sqrt(std::numeric_limits<double>::epsilon()));
 
     /**
      * Returns the point vector with the given name.
diff --git a/GeoLib/Grid.h b/GeoLib/Grid.h
index 8e6c075ab67285ba97c159b22d2d5604413713a0..2bb39f87e4170726a44e47e59faba04a9c433446 100644
--- a/GeoLib/Grid.h
+++ b/GeoLib/Grid.h
@@ -431,7 +431,7 @@ POINT* Grid<POINT>::getNearestPoint(P const& pnt) const
     std::array<double,6> dists(getPointCellBorderDistances(pnt, coords));
 
     if (calcNearestPointInGridCell(pnt, coords, sqr_min_dist, nearest_pnt)) {
-        double min_dist(sqrt(sqr_min_dist));
+        double min_dist(std::sqrt(sqr_min_dist));
         if (dists[0] >= min_dist && dists[1] >= min_dist
             && dists[2] >= min_dist && dists[3] >= min_dist
             && dists[4] >= min_dist && dists[5] >= min_dist) {
@@ -478,7 +478,7 @@ POINT* Grid<POINT>::getNearestPoint(P const& pnt) const
         } // end while
     } // end else
 
-    double len(sqrt(MathLib::sqrDist(pnt, *nearest_pnt)));
+    double len(std::sqrt(MathLib::sqrDist(pnt, *nearest_pnt)));
     // search all other grid cells within the cube with the edge nodes
     std::vector<std::vector<POINT*> const*> vecs_of_pnts(
         getPntVecsOfGridCellsIntersectingCube(pnt, len));
diff --git a/GeoLib/IO/XmlIO/Qt/XmlStnInterface.cpp b/GeoLib/IO/XmlIO/Qt/XmlStnInterface.cpp
index 465c4cf1fa0e0700a8bc4b8c365aa4527948f279..df88156fb89f290b38c4995a6d32d0480f1a81c9 100644
--- a/GeoLib/IO/XmlIO/Qt/XmlStnInterface.cpp
+++ b/GeoLib/IO/XmlIO/Qt/XmlStnInterface.cpp
@@ -208,7 +208,9 @@ void XmlStnInterface::readStratigraphy( const QDomNode &stratRoot,
             /* add other horizon features here */
 
             double depth (horizon.attribute("z").toDouble());
-            if (fabs(depth - depth_check) > std::numeric_limits<double>::epsilon()) // skip soil-layer if its thickness is zero
+            if (std::abs(depth - depth_check) >
+                std::numeric_limits<double>::
+                    epsilon())  // skip soil-layer if its thickness is zero
             {
                 borehole->addSoilLayer(horizon.attribute("x").toDouble(),
                                        horizon.attribute("y").toDouble(),
@@ -326,7 +328,7 @@ void XmlStnInterface::writeBoreholeData(QDomDocument &doc,
     boreholeTag.appendChild(stationDepthTag);
     QDomText stationDepthText = doc.createTextNode(QString::number(borehole->getDepth(), 'f'));
     stationDepthTag.appendChild(stationDepthText);
-    if (fabs(borehole->getDate()) > 0)
+    if (std::abs(borehole->getDate()) > 0)
     {
         QDomElement stationDateTag = doc.createElement("bdate");
         boreholeTag.appendChild(stationDateTag);
diff --git a/GeoLib/Raster.cpp b/GeoLib/Raster.cpp
index 7336958d5def40d875c52278044664ae7256d3fd..48623dd40387a6842e21d5b7fd178e290a7af35d 100644
--- a/GeoLib/Raster.cpp
+++ b/GeoLib/Raster.cpp
@@ -88,8 +88,8 @@ double Raster::interpolateValueAtPoint(MathLib::Point3d const& pnt) const
     double const yIdx (std::floor(yPos));    //  so not to over- or underflow.
 
     // weights for bilinear interpolation
-    double const xShift = std::fabs((xPos - xIdx) - 0.5);
-    double const yShift = std::fabs((yPos - yIdx) - 0.5);
+    double const xShift = std::abs((xPos - xIdx) - 0.5);
+    double const yShift = std::abs((yPos - yIdx) - 0.5);
     std::array<double,4> weight = {{ (1-xShift)*(1-yShift), xShift*(1-yShift), xShift*yShift, (1-xShift)*yShift }};
 
     // neighbors to include in interpolation
@@ -119,7 +119,8 @@ double Raster::interpolateValueAtPoint(MathLib::Point3d const& pnt) const
         }
 
         // remove no data values
-        if (std::fabs(pix_val[j] - _header.no_data) < std::numeric_limits<double>::epsilon())
+        if (std::abs(pix_val[j] - _header.no_data) <
+            std::numeric_limits<double>::epsilon())
         {
             weight[j] = 0;
             no_data_count++;
diff --git a/MathLib/Curve/PiecewiseLinearMonotonicCurve.cpp b/MathLib/Curve/PiecewiseLinearMonotonicCurve.cpp
index 6ddc1026b29ca4ee23a45dde458bdbf92a1b348d..fdd73775ea0ba013e6a0c9a6389e7992be3cd737 100644
--- a/MathLib/Curve/PiecewiseLinearMonotonicCurve.cpp
+++ b/MathLib/Curve/PiecewiseLinearMonotonicCurve.cpp
@@ -34,7 +34,7 @@ bool PiecewiseLinearMonotonicCurve::isStrongMonotonic() const
 {
     const double gradient0 = getDerivative(_supp_pnts[0]);
 
-    if (std::fabs(gradient0) < std::numeric_limits<double>::min())
+    if (std::abs(gradient0) < std::numeric_limits<double>::min())
     {
         return false;
     }
diff --git a/MathLib/Nonlinear/Root1D.h b/MathLib/Nonlinear/Root1D.h
index edc37751936911a8a877ac39b45fd2495bce48aa..a2e9f8f273b2a3e9492aac0bbd51f6b4e05def69 100644
--- a/MathLib/Nonlinear/Root1D.h
+++ b/MathLib/Nonlinear/Root1D.h
@@ -112,7 +112,7 @@ public:
     }
 
     //! Returns the size of the current search interval.
-    double getRange() const { return std::fabs(_a - _b); }
+    double getRange() const { return std::abs(_a - _b); }
 
 private:
     Function const& _f;
diff --git a/MathLib/TemplatePoint.h b/MathLib/TemplatePoint.h
index a1831f748699eaea8a325faf3d20773140b9e8a9..a106db21677e48dbb8b3994f9e5b6275be820b81 100644
--- a/MathLib/TemplatePoint.h
+++ b/MathLib/TemplatePoint.h
@@ -155,8 +155,8 @@ bool lessEq(TemplatePoint<T, DIM> const& a, TemplatePoint<T, DIM> const& b,
 {
     auto coordinateIsLargerEps = [&eps](T const u, T const v) -> bool
     {
-        return std::fabs(u - v) > eps * std::min(std::fabs(v), std::fabs(u)) &&
-               std::fabs(u - v) > eps;
+        return std::abs(u - v) > eps * std::min(std::abs(v), std::abs(u)) &&
+               std::abs(u - v) > eps;
     };
 
     for (std::size_t i = 0; i < DIM; ++i)
diff --git a/MeshLib/Elements/LineRule2.cpp b/MeshLib/Elements/LineRule2.cpp
index 373a8ef4b8b4146b72a8ad6b94469281638ef281..fbc833b43d5663e7b905f473b7670c8cab7eaff7 100644
--- a/MeshLib/Elements/LineRule2.cpp
+++ b/MeshLib/Elements/LineRule2.cpp
@@ -21,7 +21,7 @@ const unsigned LineRule2::edge_nodes[1][2] =
 
 double LineRule2::computeVolume(Node const* const* _nodes)
 {
-    return sqrt(MathLib::sqrDist(*_nodes[0], *_nodes[1]));
+    return std::sqrt(MathLib::sqrDist(*_nodes[0], *_nodes[1]));
 }
 
 bool LineRule2::isPntInElement(Node const* const* nodes,
diff --git a/MeshLib/MeshGenerators/MeshLayerMapper.cpp b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
index b333efd6723574b21c951a5814e51376c00b6003..9f9a1e895709b996f043fbe5d902cd8b3188a0ee 100644
--- a/MeshLib/MeshGenerators/MeshLayerMapper.cpp
+++ b/MeshLib/MeshGenerators/MeshLayerMapper.cpp
@@ -301,7 +301,7 @@ bool MeshLayerMapper::layerMapping(MeshLib::Mesh const& mesh,
 
         double elevation (raster.getValueAtPoint(*nodes[i]));
         constexpr double eps = std::numeric_limits<double>::epsilon();
-        if (std::fabs(elevation - header.no_data) < eps)
+        if (std::abs(elevation - header.no_data) < eps)
         {
             if (ignore_nodata)
             {
diff --git a/MeshLib/MeshQuality/ElementSizeMetric.cpp b/MeshLib/MeshQuality/ElementSizeMetric.cpp
index 245b695b2ba89c9154440a06f22697d2c316fdb9..47eaecae48d018abc9e8544d9cbb00e41afd64e9 100644
--- a/MeshLib/MeshQuality/ElementSizeMetric.cpp
+++ b/MeshLib/MeshQuality/ElementSizeMetric.cpp
@@ -59,7 +59,7 @@ std::size_t ElementSizeMetric::calc1dQuality()
         double area(std::numeric_limits<double>::max());
         _element_quality_metric[k] = elements[k]->getContent();
         if (_element_quality_metric[k] <
-            sqrt(fabs(std::numeric_limits<double>::epsilon())))
+            std::sqrt(std::abs(std::numeric_limits<double>::epsilon())))
         {
             error_count++;
         }
@@ -93,7 +93,7 @@ std::size_t ElementSizeMetric::calc2dQuality()
             continue;
         }
         double const area = elem.getContent();
-        if (area < sqrt(fabs(std::numeric_limits<double>::epsilon())))
+        if (area < std::sqrt(std::abs(std::numeric_limits<double>::epsilon())))
         {
             error_count++;
         }
@@ -128,7 +128,7 @@ std::size_t ElementSizeMetric::calc3dQuality()
         }
 
         double const volume (elem.getContent());
-        if (volume < sqrt(fabs(std::numeric_limits<double>::epsilon())))
+        if (volume < sqrt(std::abs(std::numeric_limits<double>::epsilon())))
         {
             error_count++;
         }
diff --git a/MeshLib/MeshQuality/RadiusEdgeRatioMetric.cpp b/MeshLib/MeshQuality/RadiusEdgeRatioMetric.cpp
index 3c7b733a5b6805b99cb12070b8d27e156f791708..e9c7524d14d14107dc5e8cda88238fcb600e7ad2 100644
--- a/MeshLib/MeshQuality/RadiusEdgeRatioMetric.cpp
+++ b/MeshLib/MeshQuality/RadiusEdgeRatioMetric.cpp
@@ -36,7 +36,7 @@ void RadiusEdgeRatioMetric::calculateQuality ()
         std::copy_n(elem.getNodes(), n_nodes, pnts.begin());
         GeoLib::MinimalBoundingSphere const s(pnts);
         auto const& [min, max] = computeSqrEdgeLengthRange(elem);
-        _element_quality_metric[k] = sqrt(min)/(2*s.getRadius());
+        _element_quality_metric[k] = std::sqrt(min) / (2 * s.getRadius());
     }
 }
 
diff --git a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp
index ab69b5d7f39f89d83ce04a8e4e49bf6e1075aac6..86a318c2c700405c1c784b819a68dc5d36323888 100644
--- a/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp
+++ b/NumLib/TimeStepping/Algorithms/EvolutionaryPIDcontroller.cpp
@@ -116,7 +116,7 @@ double EvolutionaryPIDcontroller::limitStepSize(
         // If the last time step was rejected and the new predicted time step
         // size is identical to that of the previous rejected step, the new
         // step size is then reduced by half.
-        if (std::fabs(limited_h - _ts_current.dt()) <
+        if (std::abs(limited_h - _ts_current.dt()) <
             std::numeric_limits<double>::min())
         {
             limited_h = std::max(_h_min, 0.5 * limited_h);
diff --git a/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp b/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp
index bd4d547a02c9462cf0f01ca2a9bf2fbb9fc6c9ea..c6923916a1c720af9bf98c70f0a8904cf4238f12 100644
--- a/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp
+++ b/NumLib/TimeStepping/Algorithms/IterationNumberBasedTimeStepping.cpp
@@ -69,7 +69,7 @@ std::tuple<bool, double> IterationNumberBasedTimeStepping::next(
         double dt = getNextTimeStepSize();
         // In case it is the first time be rejected, re-computed dt again with
         // current dt
-        if (std::fabs(dt - _ts_current.dt()) <
+        if (std::abs(dt - _ts_current.dt()) <
             std::numeric_limits<double>::epsilon())
         {
             // time step was rejected, keep dt for the next dt computation.
diff --git a/ProcessLib/Output/Output.cpp b/ProcessLib/Output/Output.cpp
index c2b4e1f9af02d11fbba479e32be5876ce2b200a9..4cd99e4386ecbebfd9954d7a6ebc84de6ccbd5b3 100644
--- a/ProcessLib/Output/Output.cpp
+++ b/ProcessLib/Output/Output.cpp
@@ -64,7 +64,7 @@ bool Output::shallDoOutput(int timestep, double const t)
     auto const fixed_output_time = std::lower_bound(
         cbegin(_fixed_output_times), cend(_fixed_output_times), t);
     if ((fixed_output_time != cend(_fixed_output_times)) &&
-        (std::fabs(*fixed_output_time - t) <
+        (std::abs(*fixed_output_time - t) <
          std::numeric_limits<double>::epsilon()))
     {
         return true;
diff --git a/ProcessLib/PhaseField/PhaseFieldProcess.cpp b/ProcessLib/PhaseField/PhaseFieldProcess.cpp
index 7bb76091ac23e106d45b3864a3a3153774f7fcec..86083c05caa11cbb4ffedddc2fc1a50b1bba8d14 100644
--- a/ProcessLib/PhaseField/PhaseFieldProcess.cpp
+++ b/ProcessLib/PhaseField/PhaseFieldProcess.cpp
@@ -319,7 +319,7 @@ void PhaseFieldProcess<DisplacementDim>::postNonLinearSolverConcreteProcess(
             _process_data.pressure =
                 _process_data.injected_volume / _process_data.crack_volume;
             _process_data.pressure_error =
-                std::fabs(_process_data.pressure_old - _process_data.pressure) /
+                std::abs(_process_data.pressure_old - _process_data.pressure) /
                 _process_data.pressure;
             INFO("Internal pressure: {:g} and Pressure error: {:.4e}",
                  _process_data.pressure, _process_data.pressure_error);
diff --git a/ProcessLib/TimeLoop.cpp b/ProcessLib/TimeLoop.cpp
index a8ad37f24efa725020fa0c29f2ea837f4af90381..37603ef24985dda5a787e0d3b49a660215960d19 100644
--- a/ProcessLib/TimeLoop.cpp
+++ b/ProcessLib/TimeLoop.cpp
@@ -458,7 +458,7 @@ double TimeLoop::computeTimeStepping(const double prev_dt, double& t,
     dt = NumLib::possiblyClampDtToNextFixedTime(t, dt,
                                                 _output->getFixedOutputTimes());
     // Check whether the time stepping is stabilized
-    if (std::fabs(dt - prev_dt) < std::numeric_limits<double>::epsilon())
+    if (std::abs(dt - prev_dt) < std::numeric_limits<double>::epsilon())
     {
         if (_last_step_rejected)
         {
@@ -656,7 +656,7 @@ bool TimeLoop::loop()
         }
 
         dt = computeTimeStepping(prev_dt, t, accepted_steps, rejected_steps);
-        if (std::fabs(t - _end_time) < std::numeric_limits<double>::epsilon() ||
+        if (std::abs(t - _end_time) < std::numeric_limits<double>::epsilon() ||
             t + dt > _end_time)
         {
             break;