diff --git a/GeoLib/AnalyticalGeometry.cpp b/GeoLib/AnalyticalGeometry.cpp
index ecb614bc36be7cba2fd540bba4cb5a539f256064..713b9c14ac8150847bea7ed6fbe7cd4ad56420a5 100644
--- a/GeoLib/AnalyticalGeometry.cpp
+++ b/GeoLib/AnalyticalGeometry.cpp
@@ -442,35 +442,28 @@ void computeAndInsertAllIntersectionPoints(GeoLib::PointVec& pnt_vec,
     }
 }
 
-GeoLib::Polygon rotatePolygonToXY(GeoLib::Polygon const& polygon_in,
-                                  Eigen::Vector3d& plane_normal)
+std::tuple<std::vector<GeoLib::Point*>, Eigen::Vector3d>
+rotatePolygonPointsToXY(GeoLib::Polygon const& polygon_in)
 {
     // 1 copy all points
-    auto* polygon_pnts(new std::vector<GeoLib::Point*>);
+    std::vector<GeoLib::Point*> polygon_points;
+    polygon_points.reserve(polygon_in.getNumberOfPoints());
     for (std::size_t k(0); k < polygon_in.getNumberOfPoints(); k++)
     {
-        polygon_pnts->push_back(new GeoLib::Point(*(polygon_in.getPoint(k))));
+        polygon_points.push_back(new GeoLib::Point(*(polygon_in.getPoint(k))));
     }
 
     // 2 rotate points
-    double d_polygon;
-    std::tie(plane_normal, d_polygon) = GeoLib::getNewellPlane(*polygon_pnts);
+    auto [plane_normal, d_polygon] = GeoLib::getNewellPlane(polygon_points);
     Eigen::Matrix3d const rot_mat =
         GeoLib::computeRotationMatrixToXY(plane_normal);
-    GeoLib::rotatePoints(rot_mat, *polygon_pnts);
+    GeoLib::rotatePoints(rot_mat, polygon_points);
 
     // 3 set z coord to zero
-    std::for_each(polygon_pnts->begin(), polygon_pnts->end(),
+    std::for_each(polygon_points.begin(), polygon_points.end(),
                   [](GeoLib::Point* p) { (*p)[2] = 0.0; });
 
-    // 4 create new polygon
-    GeoLib::Polyline rot_polyline(*polygon_pnts);
-    for (std::size_t k(0); k < polygon_in.getNumberOfPoints(); k++)
-    {
-        rot_polyline.addPoint(k);
-    }
-    rot_polyline.addPoint(0);
-    return GeoLib::Polygon(rot_polyline);
+    return {polygon_points, plane_normal};
 }
 
 std::vector<MathLib::Point3d> lineSegmentIntersect2d(
diff --git a/GeoLib/AnalyticalGeometry.h b/GeoLib/AnalyticalGeometry.h
index 982531a6d392ee21e80f2afe073bfab8fb460470..7820093df65de329caced26451dbacc87d5db058 100644
--- a/GeoLib/AnalyticalGeometry.h
+++ b/GeoLib/AnalyticalGeometry.h
@@ -202,22 +202,21 @@ std::unique_ptr<Point> triangleLineIntersection(
  * pnt_vec. For each intersection an id is returned.  This id is used to split the two
  * intersecting straight line segments in four straight line segments.
  */
-void computeAndInsertAllIntersectionPoints(PointVec &pnt_vec,
-    std::vector<Polyline*> & plys);
+void computeAndInsertAllIntersectionPoints(PointVec& pnt_vec,
+                                           std::vector<Polyline*>& plys);
 
 /**
- * Function rotates a polygon to the xy plane. For this reason, (1) the points of
- * the given polygon are copied, (2) a so called Newell plane is computed
- * (getNewellPlane()) and the points are rotated, (3) for security the
- * \f$z\f$ coordinates of the rotated points are set to zero and finally, (4) a
- * new polygon is constructed using the rotated points.
+ * Function rotates a polygon to the xy plane. For this reason, (1) the points
+ * of the given polygon are copied, (2) a so called Newell plane is computed
+ * (getNewellPlane()) and the points are rotated, (3) for accuracy reasons the
+ * \f$z\f$ coordinates of the rotated points are set to zero
  * \see getNewellPlane()
  * @param polygon_in a copy of the polygon_in polygon will be rotated
- * @param plane_normal the normal of the original Newell plane
- * @return a rotated polygon
+ * @return vector of rotated points and normal based on the original Newell
+ * plane
  */
-Polygon rotatePolygonToXY(Polygon const& polygon_in,
-                                  Eigen::Vector3d& plane_normal);
+std::tuple<std::vector<GeoLib::Point*>, Eigen::Vector3d>
+rotatePolygonPointsToXY(GeoLib::Polygon const& polygon_in);
 
 /// Sorts the vector of segments such that the \f$i\f$-th segment is connected
 /// with the \f$i+1\f$st segment, i.e. the end point of the \f$i\f$-th segment
diff --git a/MeshGeoToolsLib/MeshEditing/MarkNodesOutsideOfPolygon.h b/MeshGeoToolsLib/MeshEditing/MarkNodesOutsideOfPolygon.h
index 603ef6053deb9f9fc956ca299ff68112731b541c..462f2403202282178ba7666094eb7ef6accb6ccc 100644
--- a/MeshGeoToolsLib/MeshEditing/MarkNodesOutsideOfPolygon.h
+++ b/MeshGeoToolsLib/MeshEditing/MarkNodesOutsideOfPolygon.h
@@ -23,9 +23,9 @@ namespace MeshGeoToolsLib
 std::vector<bool> markNodesOutSideOfPolygon(
     std::vector<MeshLib::Node*> const& nodes, GeoLib::Polygon const& polygon)
 {
-    // *** rotate polygon to xy_plane
-    Eigen::Vector3d normal;
-    GeoLib::Polygon rot_polygon(GeoLib::rotatePolygonToXY(polygon, normal));
+    // *** rotate polygon points to xy-plane
+    auto [rotated_polygon_points, normal] =
+        GeoLib::rotatePolygonPointsToXY(polygon);
 
     // *** rotate mesh nodes to xy-plane
     // 1 copy all mesh nodes to GeoLib::Points
@@ -41,13 +41,24 @@ std::vector<bool> markNodesOutSideOfPolygon(
     std::for_each(rotated_nodes.begin(), rotated_nodes.end(),
                   [](GeoLib::Point* p) { (*p)[2] = 0.0; });
 
-    // *** mark rotated nodes
     std::vector<bool> outside(rotated_nodes.size(), true);
-    for (std::size_t k(0); k < rotated_nodes.size(); k++)
+    // *** mark rotated nodes inside rotated polygon
     {
-        if (rot_polygon.isPntInPolygon(*(rotated_nodes[k])))
+        // create new polygon using the rotated points
+        GeoLib::Polyline rotated_polyline(rotated_polygon_points);
+        for (std::size_t k(0); k < polygon.getNumberOfPoints(); k++)
         {
-            outside[k] = false;
+            rotated_polyline.addPoint(k);
+        }
+        rotated_polyline.addPoint(0);
+        GeoLib::Polygon const rotated_polygon(rotated_polyline);
+
+        for (std::size_t k(0); k < rotated_nodes.size(); k++)
+        {
+            if (rotated_polygon.isPntInPolygon(*(rotated_nodes[k])))
+            {
+                outside[k] = false;
+            }
         }
     }
 
@@ -56,9 +67,7 @@ std::vector<bool> markNodesOutSideOfPolygon(
         delete rotated_node;
     }
 
-    std::vector<GeoLib::Point*>& rot_polygon_pnts(
-        const_cast<std::vector<GeoLib::Point*>&>(rot_polygon.getPointsVec()));
-    for (auto& rot_polygon_pnt : rot_polygon_pnts)
+    for (auto& rot_polygon_pnt : rotated_polygon_points)
     {
         delete rot_polygon_pnt;
     }