diff --git a/Applications/DataExplorer/DataView/GEOModels.cpp b/Applications/DataExplorer/DataView/GEOModels.cpp
index 5b1e901ec7dbafb22295fb5b0930c4a10528bb23..379bae091ab233edad9f96c4e6899a51d9ed85c7 100644
--- a/Applications/DataExplorer/DataView/GEOModels.cpp
+++ b/Applications/DataExplorer/DataView/GEOModels.cpp
@@ -17,6 +17,7 @@
 
 #include <logog/include/logog.hpp>
 
+#include "Applications/FileIO/Legacy/createSurface.h"
 #include "GeoTreeModel.h"
 #include "StationTreeModel.h"
 
@@ -148,6 +149,13 @@ void GEOModels::removeSurfaceVec(std::string const& name)
     _geoModel->removeGeoList(name, GeoLib::GEOTYPE::SURFACE);
 }
 
+void GEOModels::renameGeometry(std::string const& old_name,
+                               std::string const& new_name)
+{
+    _geoModel->renameGeometry(old_name, new_name);
+    updateGeometry(new_name);
+}
+
 void GEOModels::connectPolylineSegments(
     const std::string& geoName,
     std::vector<std::size_t> const& indexlist,
@@ -175,22 +183,34 @@ void GEOModels::connectPolylineSegments(
             // insert result in a new vector of polylines (because the GEOObjects::appendPolylines()-method requires a vector)
             std::vector<GeoLib::Polyline*> connected_ply;
 
+            connected_ply.push_back(new_line);
+            _geo_objects.appendPolylineVec(connected_ply, geoName);
+
             if (closePly)
             {
                 new_line->closePolyline();
 
                 if (triangulatePly)
                 {
-                    std::vector<GeoLib::Surface*> new_sfc;
-                    new_sfc.push_back(GeoLib::Surface::createSurface(*new_line));
-                    _geo_objects.appendSurfaceVec(new_sfc, geoName);
+                    INFO(
+                        "Creating a surface by triangulation of the polyline "
+                        "...");
+                    if (FileIO::createSurface(*new_line, _geo_objects, geoName))
+                    {
+                        INFO("\t done");
+                    }
+                    else
+                    {
+                        WARN(
+                            "\t Creating a surface by triangulation of the "
+                            "polyline failed.");
+                    }
+                    plyVec = _geo_objects.getPolylineVecObj(geoName);
                 }
             }
 
-            connected_ply.push_back(new_line);
             if (!ply_name.empty())
                 plyVec->setNameOfElementByID(polylines->size(), ply_name);
-            _geo_objects.appendPolylineVec(connected_ply, geoName);
         }
         else
             OGSError::box("Error connecting polyines.");
diff --git a/Applications/DataExplorer/DataView/GEOModels.h b/Applications/DataExplorer/DataView/GEOModels.h
index 54981a89ab2fb13a42cdb65e806d3dfffce13b55..4e21183c05174ec72fa27729568c6303ebe04365 100644
--- a/Applications/DataExplorer/DataView/GEOModels.h
+++ b/Applications/DataExplorer/DataView/GEOModels.h
@@ -76,6 +76,9 @@ public slots:
     void appendSurfaceVec(std::string const& name);
     void removeSurfaceVec(std::string const& name);
 
+    void renameGeometry(std::string const& old_name,
+                        std::string const& new_name);
+
     /// Adds the name 'new_name' for the geo-object specified by the parameters
     void addNameForElement(std::string const& geometry_name,
                            GeoLib::GEOTYPE const object_type,
@@ -166,6 +169,12 @@ public:
         _geo_models.removeSurfaceVec(name);
     };
 
+    void renameGeometry(std::string const& old_name,
+                        std::string const& new_name) override
+    {
+        _geo_models.renameGeometry(old_name, new_name);
+    }
+
 private:
     GEOModels& _geo_models;
 };
diff --git a/Applications/DataExplorer/DataView/GeoTreeModel.cpp b/Applications/DataExplorer/DataView/GeoTreeModel.cpp
index 70a685fec67a5962b74c9ecce401e87f653353fb..6f61745a3d65c6114f349d338a8ef213a7bd1968 100644
--- a/Applications/DataExplorer/DataView/GeoTreeModel.cpp
+++ b/Applications/DataExplorer/DataView/GeoTreeModel.cpp
@@ -274,10 +274,31 @@ void GeoTreeModel::addChildren(GeoObjectListItem* sfcList,
     INFO("%d surfaces added.", end_index - start_index);
 }
 
+void GeoTreeModel::renameGeometry(std::string const& old_name,
+                                  std::string const& new_name)
+{
+    for (auto tree_item_entry : _lists)
+    {
+        if (old_name == tree_item_entry->data(0).toString().toStdString())
+        {
+            QVariant new_entry(QString::fromStdString(new_name));
+            tree_item_entry->setData(0, new_entry);
+            break;
+        }
+    }
+    for (auto tree_item_entry : _lists)
+    {
+        if (new_name == tree_item_entry->data(0).toString().toStdString())
+        {
+            INFO("Found tree_item_entry with name '%s'.", new_name.c_str());
+        }
+    }
+}
+
 /**
  * Removes the TreeItem with the given name including all its children
  */
-void GeoTreeModel::removeGeoList(const std::string &name, GeoLib::GEOTYPE type)
+void GeoTreeModel::removeGeoList(const std::string& name, GeoLib::GEOTYPE type)
 {
     for (std::size_t i = 0; i < _lists.size(); i++)
         if (name == _lists[i]->data(0).toString().toStdString())
diff --git a/Applications/DataExplorer/DataView/GeoTreeModel.h b/Applications/DataExplorer/DataView/GeoTreeModel.h
index 7f8f45210f6c13675d41141f38d5619becc97e2c..47062aa7c12d96147d18169c292aaddbae78c6b8 100644
--- a/Applications/DataExplorer/DataView/GeoTreeModel.h
+++ b/Applications/DataExplorer/DataView/GeoTreeModel.h
@@ -67,6 +67,9 @@ public:
      */
     void removeGeoList(const std::string &name, GeoLib::GEOTYPE type);
 
+    void renameGeometry(std::string const& old_name,
+                        std::string const& new_name);
+
     void setNameForItem(const std::string &name, GeoLib::GEOTYPE type, std::size_t id, std::string item_name);
 
     /*
diff --git a/Applications/DataExplorer/mainwindow.cpp b/Applications/DataExplorer/mainwindow.cpp
index 329b90cc730fda6c5f96ecc0c5ddc039425124ba..7425d5d6ef3306bb626cf102ce6d323a86f7ff5f 100644
--- a/Applications/DataExplorer/mainwindow.cpp
+++ b/Applications/DataExplorer/mainwindow.cpp
@@ -35,6 +35,7 @@
 #include "Applications/FileIO/GMSInterface.h"
 #include "Applications/FileIO/Gmsh/GMSHInterface.h"
 #include "Applications/FileIO/Gmsh/GmshReader.h"
+#include "Applications/FileIO/Legacy/OGSIOVer4.h"
 #include "Applications/FileIO/PetrelInterface.h"
 #include "Applications/FileIO/TetGenInterface.h"
 #include "Applications/FileIO/XmlIO/Qt/XmlGspInterface.h"
@@ -43,7 +44,6 @@
 #include "BaseLib/FileTools.h"
 #include "BaseLib/Histogram.h"
 #include "GeoLib/DuplicateGeometry.h"
-#include "GeoLib/IO/Legacy/OGSIOVer4.h"
 #include "GeoLib/IO/XmlIO/Qt/XmlGmlInterface.h"
 #include "GeoLib/IO/XmlIO/Qt/XmlStnInterface.h"
 #include "GeoLib/Raster.h"
@@ -450,9 +450,9 @@ void MainWindow::loadFile(ImportFileType::type t, const QString &fileName)
         {
             std::string unique_name;
             std::vector<std::string> errors;
-            if (!GeoLib::IO::Legacy::readGLIFileV4(fileName.toStdString(),
-                                                   _project.getGEOObjects(),
-                                                   unique_name, errors))
+            if (!FileIO::Legacy::readGLIFileV4(fileName.toStdString(),
+                                               _project.getGEOObjects(),
+                                               unique_name, errors))
             {
                 for (auto& error : errors)
                     OGSError::box(QString::fromStdString(error));
@@ -773,8 +773,8 @@ void MainWindow::writeGeometryToFile(QString gliName, QString fileName)
     QFileInfo fi(fileName);
     if (fi.suffix().toLower() == "gli")
     {
-        GeoLib::IO::Legacy::writeAllDataToGLIFileV4(fileName.toStdString(),
-                                                    _project.getGEOObjects());
+        FileIO::Legacy::writeAllDataToGLIFileV4(fileName.toStdString(),
+                                                _project.getGEOObjects());
         return;
     }
 #endif
diff --git a/Applications/FileIO/CMakeLists.txt b/Applications/FileIO/CMakeLists.txt
index b1c52ec754c10108ed8357ab08421c324df5e7bb..00b7fadad86ba5a7327503dc416e1891d61a30b3 100644
--- a/Applications/FileIO/CMakeLists.txt
+++ b/Applications/FileIO/CMakeLists.txt
@@ -1,5 +1,6 @@
 GET_SOURCE_FILES(SOURCES)
 APPEND_SOURCE_FILES(SOURCES Gmsh)
+APPEND_SOURCE_FILES(SOURCES Legacy)
 
 if(NOT Shapelib_FOUND)
     list(REMOVE_ITEM SOURCES SHPInterface.h SHPInterface.cpp)
diff --git a/GeoLib/IO/Legacy/OGSIOVer4.cpp b/Applications/FileIO/Legacy/OGSIOVer4.cpp
similarity index 91%
rename from GeoLib/IO/Legacy/OGSIOVer4.cpp
rename to Applications/FileIO/Legacy/OGSIOVer4.cpp
index 715efec64244bf3eb4fed11625bbaee01f001b64..a668c281a4c4445cdae209da1d024b3ee52d30cc 100644
--- a/GeoLib/IO/Legacy/OGSIOVer4.cpp
+++ b/Applications/FileIO/Legacy/OGSIOVer4.cpp
@@ -20,11 +20,12 @@
 
 #include <logog/include/logog.hpp>
 
+#include "Applications/FileIO/Legacy/createSurface.h"
+
 #include "BaseLib/FileTools.h"
 #include "BaseLib/StringTools.h"
 
 #include "GeoLib/AnalyticalGeometry.h"
-#include "GeoLib/EarClippingTriangulation.h"
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/Point.h"
 #include "GeoLib/PointVec.h"
@@ -36,9 +37,7 @@
 
 #include "GeoLib/IO/TINInterface.h"
 
-namespace GeoLib
-{
-namespace IO
+namespace FileIO
 {
 namespace Legacy {
 
@@ -268,21 +267,21 @@ std::string readPolylines(std::istream &in, std::vector<GeoLib::Polyline*>* ply_
    01/2010 TF signatur modification, reimplementation
 **************************************************************************/
 /** read a single Surface */
-std::string readSurface(std::istream &in,
-                        std::vector<GeoLib::Polygon*> &polygon_vec,
-                        std::vector<GeoLib::Surface*> &sfc_vec,
-                        std::map<std::string,std::size_t>& sfc_names,
-                        const std::vector<GeoLib::Polyline*> &ply_vec,
+std::string readSurface(std::istream& in,
+                        std::vector<GeoLib::Polygon*>& polygon_vec,
+                        std::vector<GeoLib::Surface*>& sfc_vec,
+                        std::map<std::string, std::size_t>& sfc_names,
+                        const std::vector<GeoLib::Polyline*>& ply_vec,
                         const std::map<std::string, std::size_t>& ply_vec_names,
-                        GeoLib::PointVec &pnt_vec,
-                        std::string const& path, std::vector<std::string>& errors)
+                        GeoLib::PointVec& pnt_vec, std::string const& path,
+                        std::vector<std::string>& errors)
 {
     std::string line;
     GeoLib::Surface* sfc(nullptr);
 
     int type (-1);
     std::string name;
-    std::size_t ply_id (0); // std::numeric_limits<std::size_t>::max());
+    std::size_t ply_id (std::numeric_limits<std::size_t>::max());
 
     do {
         in >> line;
@@ -385,13 +384,14 @@ std::string readSurface(std::istream &in,
    05/2004 CC Modification
    01/2010 TF changed signature of function, big modifications
 **************************************************************************/
-std::string readSurfaces(std::istream &in,
-                         std::vector<GeoLib::Surface*> &sfc_vec,
-                         std::map<std::string, std::size_t>& sfc_names,
-                         const std::vector<GeoLib::Polyline*> &ply_vec,
-                         const std::map<std::string,std::size_t>& ply_vec_names,
-                         GeoLib::PointVec & pnt_vec,
-                         const std::string &path, std::vector<std::string>& errors)
+std::string readSurfaces(
+    std::istream& in, std::vector<GeoLib::Surface*>& sfc_vec,
+    std::map<std::string, std::size_t>& sfc_names,
+    const std::vector<GeoLib::Polyline*>& ply_vec,
+    const std::map<std::string, std::size_t>& ply_vec_names,
+    GeoLib::PointVec& pnt_vec, const std::string& path,
+    std::vector<std::string>& errors, GeoLib::GEOObjects& geo,
+    std::string const& unique_name)
 {
     if (!in.good())
     {
@@ -405,25 +405,23 @@ std::string readSurfaces(std::istream &in,
     while (!in.eof() && !in.fail() && tag.find("#SURFACE") != std::string::npos)
     {
         std::size_t n_polygons (polygon_vec.size());
-        tag = readSurface(in,
-                          polygon_vec,
-                          sfc_vec,
-                          sfc_names,
-                          ply_vec,
-                          ply_vec_names,
-                          pnt_vec,
-                          path,
-                          errors);
+        tag = readSurface(in, polygon_vec, sfc_vec, sfc_names, ply_vec,
+                          ply_vec_names, pnt_vec, path, errors);
         if (n_polygons < polygon_vec.size())
         {
-            // subdivide polygon in simple polygons
-            GeoLib::Surface* sfc(GeoLib::Surface::createSurface(
-                                         *(dynamic_cast<GeoLib::Polyline*> (polygon_vec
-                                                                            [
-                                                                                    polygon_vec
-                                                                                    .
-                                                                                    size() - 1]))));
-            sfc_vec.push_back(sfc);
+            INFO("Creating a surface by triangulation of the polyline ...");
+            if (FileIO::createSurface(
+                    *(dynamic_cast<GeoLib::Polyline*>(polygon_vec.back())), geo,
+                    unique_name))
+            {
+                INFO("\t done");
+            }
+            else
+            {
+                WARN(
+                    "\t Creating a surface by triangulation of the polyline "
+                    "failed.");
+            }
         }
     }
     for (auto & k : polygon_vec)
@@ -487,19 +485,36 @@ bool readGLIFileV4(const std::string& fname,
     else
         INFO("GeoLib::readGLIFile(): tag #POLYLINE not found.");
 
+    if (!ply_vec->empty())
+        geo.addPolylineVec(std::move(ply_vec), unique_name,
+                           std::move(ply_names));
+
+    // Since ply_names is a unique_ptr and is given to the GEOObject instance
+    // geo it is not usable anymore. For this reason a copy is necessary.
+    std::map<std::string, std::size_t> ply_names_copy;
+    if (geo.getPolylineVecObj(unique_name))
+    {
+        ply_names_copy = std::map<std::string, std::size_t>{
+            geo.getPolylineVecObj(unique_name)->getNameIDMapBegin(),
+            geo.getPolylineVecObj(unique_name)->getNameIDMapEnd()};
+    }
+
     auto sfc_vec = std::make_unique<std::vector<GeoLib::Surface*>>();
     auto sfc_names = std::make_unique<std::map<std::string, std::size_t>>();
     if (tag.find("#SURFACE") != std::string::npos && in)
     {
         INFO("GeoLib::readGLIFile(): read surfaces from stream.");
+
         tag = readSurfaces(in,
                            *sfc_vec,
                            *sfc_names,
-                           *ply_vec,
-                           *ply_names,
+                           *geo.getPolylineVec(unique_name),
+                           ply_names_copy,
                            point_vec,
                            path,
-                           errors);
+                           errors,
+                           geo,
+                           unique_name);
         INFO("GeoLib::readGLIFile(): \tok, %d surfaces read.", sfc_vec->size());
     }
     else
@@ -507,11 +522,6 @@ bool readGLIFileV4(const std::string& fname,
 
     in.close();
 
-    if (!ply_vec->empty())
-        geo.addPolylineVec(
-            std::move(ply_vec), unique_name,
-            std::move(ply_names));  // KR: insert into GEOObjects if not empty
-
     if (!sfc_vec->empty())
         geo.addSurfaceVec(
             std::move(sfc_vec), unique_name,
@@ -688,5 +698,4 @@ void writeAllDataToGLIFileV4 (const std::string& fname, const GeoLib::GEOObjects
 }
 
 }
-} // end namespace IO
-} // end namespace GeoLib
+} // end namespace FileIO
diff --git a/GeoLib/IO/Legacy/OGSIOVer4.h b/Applications/FileIO/Legacy/OGSIOVer4.h
similarity index 93%
rename from GeoLib/IO/Legacy/OGSIOVer4.h
rename to Applications/FileIO/Legacy/OGSIOVer4.h
index 8e1f9ab992ba7434d9ab9284ec70bce149152df8..d79adbe5e16fae3f2818f0004cef35841e3f998c 100644
--- a/GeoLib/IO/Legacy/OGSIOVer4.h
+++ b/Applications/FileIO/Legacy/OGSIOVer4.h
@@ -20,11 +20,12 @@
 namespace GeoLib
 {
 class GEOObjects;
+}
 
-namespace IO
+namespace FileIO
+{
+namespace Legacy
 {
-namespace Legacy {
-
 /** Interface for handling geometry from OGS-5 and below (*.gli files) */
 
 /** Reads geometric objects from file in gli format */
@@ -42,5 +43,4 @@ void writeGLIFileV4 (const std::string& fname,
 void writeAllDataToGLIFileV4 (const std::string& fname, const GeoLib::GEOObjects& geo);
 
 }
-} // end namespace IO
-} // end namespace GeoLib
+} // end namespace FileIO
diff --git a/Applications/FileIO/Legacy/createSurface.cpp b/Applications/FileIO/Legacy/createSurface.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bb4a0808a7ad05190cbd6726fc8d2b84e3ad30cd
--- /dev/null
+++ b/Applications/FileIO/Legacy/createSurface.cpp
@@ -0,0 +1,121 @@
+/**
+ *
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include <cstdio>
+#include <list>
+#include <memory>
+
+#include <logog/include/logog.hpp>
+
+#include "createSurface.h"
+
+#include "Applications/FileIO/Gmsh/GMSHInterface.h"
+#include "Applications/FileIO/Gmsh/GmshReader.h"
+
+#include "GeoLib/GEOObjects.h"
+#include "GeoLib/Point.h"
+#include "GeoLib/Polygon.h"
+#include "GeoLib/Polyline.h"
+#include "GeoLib/Surface.h"
+
+#include "MeshLib/convertMeshToGeo.h"
+#include "MeshLib/Mesh.h"
+
+namespace FileIO
+{
+bool createSurface(GeoLib::Polyline const& ply,
+                   GeoLib::GEOObjects& geometries,
+                   std::string const& geometry_name)
+{
+    if (!ply.isClosed())
+    {
+        WARN("Error in createSurface() - Polyline is not closed.");
+        return false;
+    }
+
+    if (ply.getNumberOfPoints() <= 2)
+    {
+        WARN(
+            "Error in createSurface() - Polyline consists of less "
+            "than three points and therefore cannot be triangulated.");
+        return false;
+    }
+
+    // create new GEOObjects and insert a copy of the polyline
+    auto polyline_points = std::make_unique<std::vector<GeoLib::Point*>>();
+    GeoLib::GEOObjects geo;
+    auto ply_points = ply.getPointsVec();
+    for (auto p : ply_points)
+        polyline_points->push_back(new GeoLib::Point(*p));
+    std::string ply_name = "temporary_polyline_name";
+    geo.addPointVec(std::move(polyline_points), ply_name);
+    auto polyline =
+        std::make_unique<GeoLib::Polyline>(*geo.getPointVec(ply_name));
+    for (std::size_t k(0); k < ply.getNumberOfPoints(); ++k)
+        polyline->addPoint(ply.getPointID(k));
+    auto polylines = std::make_unique<std::vector<GeoLib::Polyline*>>();
+    polylines->push_back(polyline.release());
+    geo.addPolylineVec(std::move(polylines), ply_name);
+
+    // use GMSHInterface to create a mesh from the closed polyline
+    std::vector<std::string> geo_names;
+    geo.getGeometryNames(geo_names);
+    FileIO::GMSH::GMSHInterface gmsh_io(
+        geo, false, FileIO::GMSH::MeshDensityAlgorithm::FixedMeshDensity, 0.0,
+        0.0, 0.0, geo_names, false, false);
+    gmsh_io.setPrecision(std::numeric_limits<double>::digits10);
+
+    char file_base_name_c[L_tmpnam];
+    if (! std::tmpnam(file_base_name_c))
+    {
+       OGS_FATAL("Could not create unique file name.");
+    }
+    std::string const file_base_name(file_base_name_c);
+    gmsh_io.writeToFile(file_base_name + ".geo");
+    std::string gmsh_command =
+        "gmsh -2 -algo meshadapt \"" + file_base_name + ".geo\"";
+    gmsh_command += " -o \"" + file_base_name + ".msh\"";
+    std::system(gmsh_command.c_str());
+    auto surface_mesh =
+        FileIO::GMSH::readGMSHMesh("\"" + file_base_name + ".msh\"");
+    if (!surface_mesh)
+    {
+        WARN("The surface mesh could not be created.");
+        return false;
+    }
+    if (std::remove((file_base_name + ".geo").c_str()) !=0)
+        WARN("Could not remove temporary file '%s'.",
+            (file_base_name + ".geo").c_str());
+    if (std::remove((file_base_name + ".msh").c_str()) !=0)
+        WARN("Could not remove temporary file '%s'.",
+            (file_base_name + ".msh").c_str());
+
+    // convert the surface mesh into a geometric surface
+    if (!MeshLib::convertMeshToGeo(*surface_mesh, geometries,
+                                   std::numeric_limits<double>::epsilon()))
+    {
+        WARN("The surface mesh could not be converted to a geometry.");
+        return false;
+    }
+    std::string merged_geometry_name("geometry_with_surfaces");
+    geometries.mergeGeometries({geometry_name, surface_mesh->getName()},
+                               merged_geometry_name);
+    geometries.removeSurfaceVec(geometry_name);
+    geometries.removePolylineVec(geometry_name);
+    geometries.removePointVec(geometry_name);
+    geometries.removeSurfaceVec(surface_mesh->getName());
+    geometries.removePolylineVec(surface_mesh->getName());
+    geometries.removePointVec(surface_mesh->getName());
+    geometries.renameGeometry(merged_geometry_name, geometry_name);
+
+    return true;
+}
+
+} // end namespace
diff --git a/Applications/FileIO/Legacy/createSurface.h b/Applications/FileIO/Legacy/createSurface.h
new file mode 100644
index 0000000000000000000000000000000000000000..a51ba6042755305ca77152d880fcac753764bf0f
--- /dev/null
+++ b/Applications/FileIO/Legacy/createSurface.h
@@ -0,0 +1,30 @@
+/**
+ *
+ * \copyright
+ * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+namespace GeoLib
+{
+class Polyline;
+class GEOObjects;
+}
+
+namespace FileIO
+{
+/// Creates a plane surface from the given polyline. The polyline have to be
+/// closed, i.e. the first and the last point have to be the identical. The
+/// triangulation of the polyline is done by the finite element meshing tool
+/// Gmsh. Finally, the resulting mesh is converted into a GeoLib::Surface which
+/// is inserted into the \c GeoLib::GEOObjects instance \c geometries using the
+/// name \c geometry_name.
+bool createSurface(GeoLib::Polyline const& polyline,
+                   GeoLib::GEOObjects& geometries,
+                   std::string const& geometry_name);
+}
diff --git a/Applications/FileIO/SHPInterface.cpp b/Applications/FileIO/SHPInterface.cpp
index 81424262032901a4651da544ad82838e1eb371ec..250c474e7f68cfddb030c845104688273e77d9dc 100644
--- a/Applications/FileIO/SHPInterface.cpp
+++ b/Applications/FileIO/SHPInterface.cpp
@@ -19,6 +19,8 @@
 
 #include <logog/include/logog.hpp>
 
+#include "Applications/FileIO/Legacy/createSurface.h"
+
 #include "GeoLib/AnalyticalGeometry.h"
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/Polyline.h"
@@ -176,17 +178,20 @@ void SHPInterface::readPolygons(const SHPHandle &hSHP, int numberOfElements, con
     auto const polylines = _geoObjects.getPolylineVec(listName);
     auto sfc_vec = std::make_unique<std::vector<GeoLib::Surface*>>();
 
-    for (auto const* polyline : *polylines) {
-        GeoLib::Surface* sfc(GeoLib::Surface::createSurface(*polyline));
-        if (sfc)
-            sfc_vec->push_back(sfc);
-        else {
-            WARN("SHPInterface::readPolygons(): Could not triangulate polygon.")
+    for (auto const* polyline : *polylines)
+    {
+        INFO("Creating a surface by triangulation of the polyline ...");
+        if (FileIO::createSurface(*polyline, _geoObjects, listName))
+        {
+            INFO("\t done");
+        }
+        else
+        {
+            WARN(
+                "\t Creating a surface by triangulation of the polyline "
+                "failed.");
         }
     }
-
-    if (!sfc_vec->empty())
-        _geoObjects.addSurfaceVec(std::move(sfc_vec), listName);
 }
 
 bool SHPInterface::write2dMeshToSHP(const std::string &file_name, const MeshLib::Mesh &mesh)
diff --git a/GeoLib/IO/readGeometryFromFile.cpp b/Applications/FileIO/readGeometryFromFile.cpp
similarity index 87%
rename from GeoLib/IO/readGeometryFromFile.cpp
rename to Applications/FileIO/readGeometryFromFile.cpp
index 20d50b3a04138c6cb3276f5785b2e9b3bd5acd7d..3a6480a254ae2c9884f7308ff0b54aed74b8150c 100644
--- a/GeoLib/IO/readGeometryFromFile.cpp
+++ b/Applications/FileIO/readGeometryFromFile.cpp
@@ -16,13 +16,11 @@
 #include "BaseLib/FileTools.h"
 
 #include "GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.h"
-#include "GeoLib/IO/Legacy/OGSIOVer4.h"
+#include "Legacy/OGSIOVer4.h"
 
 #include "GeoLib/GEOObjects.h"
 
-namespace GeoLib
-{
-namespace IO
+namespace FileIO
 {
 void
 readGeometryFromFile(std::string const& fname, GeoLib::GEOObjects & geo_objs)
@@ -34,7 +32,7 @@ readGeometryFromFile(std::string const& fname, GeoLib::GEOObjects & geo_objs)
     } else {
         std::vector<std::string> errors;
         std::string geo_name; // geo_name is output of the reading function
-        GeoLib::IO::Legacy::readGLIFileV4(fname, geo_objs, geo_name, errors);
+        FileIO::Legacy::readGLIFileV4(fname, geo_objs, geo_name, errors);
     }
 
     std::vector<std::string> geo_names;
@@ -45,4 +43,3 @@ readGeometryFromFile(std::string const& fname, GeoLib::GEOObjects & geo_objs)
     }
 }
 }
-}
diff --git a/GeoLib/IO/readGeometryFromFile.h b/Applications/FileIO/readGeometryFromFile.h
similarity index 92%
rename from GeoLib/IO/readGeometryFromFile.h
rename to Applications/FileIO/readGeometryFromFile.h
index 15c40fdccd49d873ea7da1370c156374b24e6807..36257a999a9550f10f1999b46dfb8050058430cc 100644
--- a/GeoLib/IO/readGeometryFromFile.h
+++ b/Applications/FileIO/readGeometryFromFile.h
@@ -17,11 +17,8 @@ namespace GeoLib
     class GEOObjects;
 }
 
-namespace GeoLib
-{
-namespace IO
+namespace FileIO
 {
     void
     readGeometryFromFile(std::string const& fname, GeoLib::GEOObjects & geo_objs);
 }
-}
diff --git a/GeoLib/IO/writeGeometryToFile.cpp b/Applications/FileIO/writeGeometryToFile.cpp
similarity index 87%
rename from GeoLib/IO/writeGeometryToFile.cpp
rename to Applications/FileIO/writeGeometryToFile.cpp
index 831342ff6d053bd12cdce97621cf4fe4d048fef9..48f696ae885c09d7a9a04225c4aace23f71902c1 100644
--- a/GeoLib/IO/writeGeometryToFile.cpp
+++ b/Applications/FileIO/writeGeometryToFile.cpp
@@ -13,13 +13,11 @@
 #include "BaseLib/FileTools.h"
 
 #include "GeoLib/IO/XmlIO/Boost/BoostXmlGmlInterface.h"
-#include "GeoLib/IO/Legacy/OGSIOVer4.h"
+#include "Legacy/OGSIOVer4.h"
 
 #include "GeoLib/GEOObjects.h"
 
-namespace GeoLib
-{
-namespace IO
+namespace FileIO
 {
 void writeGeometryToFile(std::string const& geo_name,
     GeoLib::GEOObjects& geo_objs, std::string const& fname)
@@ -30,11 +28,10 @@ void writeGeometryToFile(std::string const& geo_name,
         xml.setNameForExport(geo_name);
         xml.writeToFile(fname);
     } else if (extension == "gli" || extension == "GLI") {
-        GeoLib::IO::Legacy::writeGLIFileV4(fname, geo_name, geo_objs);
+        FileIO::Legacy::writeGLIFileV4(fname, geo_name, geo_objs);
     } else {
         ERR("Writing of geometry failed, since it was not possible to determine"
             " the required format from file extension.");
     }
 }
 }
-}
diff --git a/GeoLib/IO/writeGeometryToFile.h b/Applications/FileIO/writeGeometryToFile.h
similarity index 97%
rename from GeoLib/IO/writeGeometryToFile.h
rename to Applications/FileIO/writeGeometryToFile.h
index a669cee45bce21b9fda8196e15d470a2a81e6a4d..4c1e39f4d2da7db8ece4c1124a7e88639c845cd3 100644
--- a/GeoLib/IO/writeGeometryToFile.h
+++ b/Applications/FileIO/writeGeometryToFile.h
@@ -15,8 +15,9 @@
 namespace GeoLib
 {
 class GEOObjects;
+}
 
-namespace IO
+namespace FileIO
 {
 
 /// Write geometry given by the \c geo_objs object and specified by the name
@@ -26,4 +27,3 @@ namespace IO
 void writeGeometryToFile(std::string const& geo_name,
     GeoLib::GEOObjects& geo_objs, std::string const& fname);
 }
-}
diff --git a/Applications/Utils/FileConverter/CMakeLists.txt b/Applications/Utils/FileConverter/CMakeLists.txt
index 1a486cca875886c968546f80b3b83e94159dea94..f466c7239a3e9c0a971084700ef8ef2cc2029ce6 100644
--- a/Applications/Utils/FileConverter/CMakeLists.txt
+++ b/Applications/Utils/FileConverter/CMakeLists.txt
@@ -13,7 +13,7 @@ endif ()
 
 add_executable(convertGEO convertGEO.cpp)
 set_target_properties(convertGEO PROPERTIES FOLDER Utilities)
-target_link_libraries(convertGEO GeoLib)
+target_link_libraries(convertGEO GeoLib ApplicationsFileIO)
 
 add_executable(generateMatPropsFromMatID generateMatPropsFromMatID.cpp )
 target_link_libraries(generateMatPropsFromMatID MeshLib)
diff --git a/Applications/Utils/FileConverter/convertGEO.cpp b/Applications/Utils/FileConverter/convertGEO.cpp
index 34fda9bb085882ce67dd69c318557dd02d02744c..43e61d9d1b57740b37b2b75e3eab6b31cdc67274 100644
--- a/Applications/Utils/FileConverter/convertGEO.cpp
+++ b/Applications/Utils/FileConverter/convertGEO.cpp
@@ -15,9 +15,8 @@
 #include "BaseLib/BuildInfo.h"
 
 #include "Applications/ApplicationsLib/LogogSetup.h"
-
-#include "GeoLib/IO/readGeometryFromFile.h"
-#include "GeoLib/IO/writeGeometryToFile.h"
+#include "Applications/FileIO/readGeometryFromFile.h"
+#include "Applications/FileIO/writeGeometryToFile.h"
 #include "GeoLib/GEOObjects.h"
 
 
@@ -39,12 +38,12 @@ int main (int argc, char* argv[])
     cmd.parse(argc, argv);
 
     GeoLib::GEOObjects geoObjects;
-    GeoLib::IO::readGeometryFromFile(argInputFileName.getValue(), geoObjects);
+    FileIO::readGeometryFromFile(argInputFileName.getValue(), geoObjects);
     std::vector<std::string> geo_names;
     geoObjects.getGeometryNames(geo_names);
     assert(geo_names.size() == 1);
 
-    GeoLib::IO::writeGeometryToFile(geo_names[0], geoObjects, argOutputFileName.getValue());
+    FileIO::writeGeometryToFile(geo_names[0], geoObjects, argOutputFileName.getValue());
 
     return EXIT_SUCCESS;
 }
diff --git a/Applications/Utils/GeoTools/CMakeLists.txt b/Applications/Utils/GeoTools/CMakeLists.txt
index 82cbb066974e94c626a3ced9d527b6cfc4e43502..8072f7d96cfe148e4c6e2ee00aaabb64db7031fa 100644
--- a/Applications/Utils/GeoTools/CMakeLists.txt
+++ b/Applications/Utils/GeoTools/CMakeLists.txt
@@ -1,6 +1,6 @@
 if(Qt5XmlPatterns_FOUND)
     add_executable(TriangulatePolyline TriangulatePolyline.cpp)
-    target_link_libraries(TriangulatePolyline GeoLib)
+    target_link_libraries(TriangulatePolyline GeoLib ApplicationsFileIO)
     set_target_properties(TriangulatePolyline PROPERTIES FOLDER Utilities)
     install(TARGETS
         TriangulatePolyline
diff --git a/Applications/Utils/GeoTools/TriangulatePolyline.cpp b/Applications/Utils/GeoTools/TriangulatePolyline.cpp
index acd65f60e9b6f7b116aabc08c1a710643f14479e..3b3674eed523328d2b6d94c584bb979b5e719ac5 100644
--- a/Applications/Utils/GeoTools/TriangulatePolyline.cpp
+++ b/Applications/Utils/GeoTools/TriangulatePolyline.cpp
@@ -16,6 +16,7 @@
 #include <tclap/CmdLine.h>
 
 #include "Applications/ApplicationsLib/LogogSetup.h"
+#include "Applications/FileIO/Legacy/createSurface.h"
 
 #include "BaseLib/BuildInfo.h"
 
@@ -95,19 +96,18 @@ int main(int argc, char *argv[])
             return EXIT_FAILURE;
     }
 
-    // create surface
-    INFO ("Triangulating surface...");
-    auto new_sfc = std::make_unique<std::vector<GeoLib::Surface*>>();
-    new_sfc->push_back(GeoLib::Surface::createSurface(*line));
-
-    GeoLib::SurfaceVec* sfc_vec (geo_objects.getSurfaceVecObj(geo_names[0]));
-    if (sfc_vec == nullptr)
+    INFO ("Creating a surface by triangulation of the polyline ...");
+    if (FileIO::createSurface(*line, geo_objects, geo_names[0]))
     {
-        geo_objects.addSurfaceVec(std::move(new_sfc), geo_names[0]);
-        sfc_vec = geo_objects.getSurfaceVecObj(geo_names[0]);
+        INFO("\t done");
     }
     else
-        geo_objects.appendSurfaceVec(*new_sfc, geo_names[0]);
+    {
+        WARN(
+            "\t Creating a surface by triangulation of the polyline "
+            "failed.");
+    }
+    GeoLib::SurfaceVec* sfc_vec(geo_objects.getSurfaceVecObj(geo_names[0]));
     std::size_t const sfc_id = geo_objects.getSurfaceVec(geo_names[0])->size() - 1;
     std::string const surface_name (polyline_name + "_surface");
     for (std::size_t i=1;;++i)
diff --git a/Applications/Utils/MeshEdit/CMakeLists.txt b/Applications/Utils/MeshEdit/CMakeLists.txt
index 8e9a765551c0c6e9f28a685813e6094dda34a91c..ebda53b715190e720bfa74a91cbfae54b9553463 100644
--- a/Applications/Utils/MeshEdit/CMakeLists.txt
+++ b/Applications/Utils/MeshEdit/CMakeLists.txt
@@ -26,7 +26,11 @@ ADD_VTK_DEPENDENCY(MoveMesh)
 set_target_properties(MoveMesh PROPERTIES FOLDER Utilities)
 
 add_executable(appendLinesAlongPolyline appendLinesAlongPolyline.cpp)
-target_link_libraries(appendLinesAlongPolyline MeshGeoToolsLib MeshLib)
+target_link_libraries(appendLinesAlongPolyline
+    MeshGeoToolsLib
+    MeshLib
+    ApplicationsFileIO
+)
 ADD_VTK_DEPENDENCY(appendLinesAlongPolyline)
 set_target_properties(appendLinesAlongPolyline PROPERTIES FOLDER Utilities)
 
@@ -47,7 +51,10 @@ set_target_properties(reviseMesh PROPERTIES FOLDER Utilities)
 
 add_executable(ResetPropertiesInPolygonalRegion
     ResetPropertiesInPolygonalRegion.cpp)
-target_link_libraries(ResetPropertiesInPolygonalRegion MeshLib)
+target_link_libraries(ResetPropertiesInPolygonalRegion
+    MeshLib
+    ApplicationsFileIO
+)
 set_target_properties(ResetPropertiesInPolygonalRegion
     PROPERTIES FOLDER Utilities)
 
@@ -64,7 +71,11 @@ set_target_properties(createLayeredMeshFromRasters PROPERTIES FOLDER Utilities)
 
 add_executable(CreateBoundaryConditionsAlongPolylines
     CreateBoundaryConditionsAlongPolylines.cpp )
-target_link_libraries(CreateBoundaryConditionsAlongPolylines MeshGeoToolsLib MeshLib)
+target_link_libraries(CreateBoundaryConditionsAlongPolylines
+    MeshGeoToolsLib
+    MeshLib
+    ApplicationsFileIO
+)
 ADD_VTK_DEPENDENCY(CreateBoundaryConditionsAlongPolylines)
 set_target_properties(CreateBoundaryConditionsAlongPolylines
     PROPERTIES FOLDER Utilities)
diff --git a/Applications/Utils/MeshEdit/CreateBoundaryConditionsAlongPolylines.cpp b/Applications/Utils/MeshEdit/CreateBoundaryConditionsAlongPolylines.cpp
index 7dfe2bdead3a37f4cdb141d0c96bf5e575fd9ef5..5f587fafea40708b6c26d77f554c9e793ae7b84f 100644
--- a/Applications/Utils/MeshEdit/CreateBoundaryConditionsAlongPolylines.cpp
+++ b/Applications/Utils/MeshEdit/CreateBoundaryConditionsAlongPolylines.cpp
@@ -17,6 +17,8 @@
 #include <tclap/CmdLine.h>
 
 #include "Applications/ApplicationsLib/LogogSetup.h"
+#include "Applications/FileIO/readGeometryFromFile.h"
+#include "Applications/FileIO/writeGeometryToFile.h"
 
 #include "BaseLib/FileTools.h"
 
@@ -25,8 +27,6 @@
 
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/Point.h"
-#include "GeoLib/IO/readGeometryFromFile.h"
-#include "GeoLib/IO/writeGeometryToFile.h"
 
 #include "MeshLib/Mesh.h"
 #include "MeshLib/Node.h"
@@ -90,9 +90,9 @@ void writeBCsAndGeometry(GeoLib::GEOObjects& geometry_sets,
 {
     if (write_gml) {
         INFO("write points to \"%s.gml\".", geo_name.c_str());
-        GeoLib::IO::writeGeometryToFile(geo_name, geometry_sets, out_fname+".gml");
+        FileIO::writeGeometryToFile(geo_name, geometry_sets, out_fname+".gml");
     }
-    GeoLib::IO::writeGeometryToFile(geo_name, geometry_sets, out_fname+".gli");
+    FileIO::writeGeometryToFile(geo_name, geometry_sets, out_fname+".gli");
 
     bool liquid_flow(false);
     if (bc_type == "LIQUID_FLOW")
@@ -177,7 +177,7 @@ int main (int argc, char* argv[])
 
     // *** read geometry
     GeoLib::GEOObjects geometries;
-    GeoLib::IO::readGeometryFromFile(geometry_fname.getValue(), geometries);
+    FileIO::readGeometryFromFile(geometry_fname.getValue(), geometries);
 
     std::string geo_name;
     {
diff --git a/Applications/Utils/MeshEdit/ResetPropertiesInPolygonalRegion.cpp b/Applications/Utils/MeshEdit/ResetPropertiesInPolygonalRegion.cpp
index 68d9b2e949f38de46cdb3ca522cdcbf93846419d..f3e2a1cb2b1fb4458fbf4d40ef90387ebc4bfb2e 100644
--- a/Applications/Utils/MeshEdit/ResetPropertiesInPolygonalRegion.cpp
+++ b/Applications/Utils/MeshEdit/ResetPropertiesInPolygonalRegion.cpp
@@ -16,6 +16,7 @@
 #include <tclap/CmdLine.h>
 
 #include "Applications/ApplicationsLib/LogogSetup.h"
+#include "Applications/FileIO/readGeometryFromFile.h"
 
 #include "MeshLib/IO/readMeshFromFile.h"
 #include "MeshLib/IO/writeMeshToFile.h"
@@ -23,7 +24,6 @@
 #include "GeoLib/AnalyticalGeometry.h"
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/Polygon.h"
-#include "GeoLib/IO/readGeometryFromFile.h"
 
 #include "MathLib/Vector3.h"
 #include "MathLib/LinAlg/Dense/DenseMatrix.h"
@@ -153,7 +153,7 @@ int main (int argc, char* argv[])
 
     // *** read geometry
     GeoLib::GEOObjects geometries;
-    GeoLib::IO::readGeometryFromFile(geometry_fname.getValue(), geometries);
+    FileIO::readGeometryFromFile(geometry_fname.getValue(), geometries);
 
     std::string geo_name;
     {
diff --git a/Applications/Utils/MeshEdit/appendLinesAlongPolyline.cpp b/Applications/Utils/MeshEdit/appendLinesAlongPolyline.cpp
index f5fe60c6dd59f7561ea5d40e14811c8368a8dcf3..ca4b1cec87900d02d7f3ca1752d1fca0f06f34f0 100644
--- a/Applications/Utils/MeshEdit/appendLinesAlongPolyline.cpp
+++ b/Applications/Utils/MeshEdit/appendLinesAlongPolyline.cpp
@@ -9,12 +9,12 @@
 #include <tclap/CmdLine.h>
 
 #include "Applications/ApplicationsLib/LogogSetup.h"
+#include "Applications/FileIO/readGeometryFromFile.h"
 
 #include "BaseLib/FileTools.h"
 
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/PolylineVec.h"
-#include "GeoLib/IO/readGeometryFromFile.h"
 
 #include "MeshGeoToolsLib/AppendLinesAlongPolyline.h"
 
@@ -44,7 +44,7 @@ int main (int argc, char* argv[])
 
     // read GEO objects
     GeoLib::GEOObjects geo_objs;
-    GeoLib::IO::readGeometryFromFile(geoFileArg.getValue(), geo_objs);
+    FileIO::readGeometryFromFile(geoFileArg.getValue(), geo_objs);
 
     std::vector<std::string> geo_names;
     geo_objs.getGeometryNames (geo_names);
diff --git a/Applications/Utils/MeshGeoTools/CMakeLists.txt b/Applications/Utils/MeshGeoTools/CMakeLists.txt
index 56cb517ec382bfe6b54667ef996c9e73a1f32f05..fde60d00cc200e27bbdf3335801f9d4f317f5f98 100644
--- a/Applications/Utils/MeshGeoTools/CMakeLists.txt
+++ b/Applications/Utils/MeshGeoTools/CMakeLists.txt
@@ -4,6 +4,7 @@ TARGET_LINK_LIBRARIES( ComputeSurfaceNodeIDsInPolygonalRegion
     MeshLib
     GeoLib
     MathLib
+    ApplicationsFileIO
     ${OGS_VTK_REQUIRED_LIBS}
 )
 
diff --git a/Applications/Utils/MeshGeoTools/ComputeSurfaceNodeIDsInPolygonalRegion.cpp b/Applications/Utils/MeshGeoTools/ComputeSurfaceNodeIDsInPolygonalRegion.cpp
index 27d81440963b7e38dcf13c4232ef65d63dff8aa2..f790b3c85370032d5082f1a2cd9f597cd3776cb3 100644
--- a/Applications/Utils/MeshGeoTools/ComputeSurfaceNodeIDsInPolygonalRegion.cpp
+++ b/Applications/Utils/MeshGeoTools/ComputeSurfaceNodeIDsInPolygonalRegion.cpp
@@ -17,15 +17,12 @@
 #include <tclap/CmdLine.h>
 
 #include "Applications/ApplicationsLib/LogogSetup.h"
-
+#include "Applications/FileIO/readGeometryFromFile.h"
 #include "BaseLib/Error.h"
 #include "BaseLib/FileTools.h"
-
-#include "MeshLib/IO/readMeshFromFile.h"
-
 #include "GeoLib/GEOObjects.h"
 #include "GeoLib/Polygon.h"
-#include "GeoLib/IO/readGeometryFromFile.h"
+#include "MeshLib/IO/readMeshFromFile.h"
 
 #include "MathLib/Vector3.h"
 
@@ -89,7 +86,7 @@ int main (int argc, char* argv[])
     INFO("Mesh read: %u nodes, %u elements.", mesh->getNumberOfNodes(), mesh->getNumberOfElements());
 
     GeoLib::GEOObjects geo_objs;
-    GeoLib::IO::readGeometryFromFile(geo_in.getValue(), geo_objs);
+    FileIO::readGeometryFromFile(geo_in.getValue(), geo_objs);
     std::vector<std::string> geo_names;
     geo_objs.getGeometryNames(geo_names);
     INFO("Geometry \"%s\" read: %u points, %u polylines.",
diff --git a/Applications/Utils/OGSFileConverter/CMakeLists.txt b/Applications/Utils/OGSFileConverter/CMakeLists.txt
index a3810f531f4b383ad5c5f1d7274cfff052d156e7..3cf643d8d9e9f9a81ba1dcb3e8e3d59217913345 100644
--- a/Applications/Utils/OGSFileConverter/CMakeLists.txt
+++ b/Applications/Utils/OGSFileConverter/CMakeLists.txt
@@ -21,7 +21,7 @@ add_library(OGSFileConverterLib
 )
 target_link_libraries(OGSFileConverterLib
     PUBLIC QtBase MathLib
-    INTERFACE MeshLib
+    INTERFACE MeshLib ApplicationsFileIO
 )
 ADD_VTK_DEPENDENCY(OGSFileConverterLib)
 
@@ -29,6 +29,7 @@ add_executable(OGSFileConverter main.cpp)
 
 target_link_libraries(OGSFileConverter
     OGSFileConverterLib
+    ApplicationsFileIO
 )
 ADD_VTK_DEPENDENCY(OGSFileConverter)
 
diff --git a/Applications/Utils/OGSFileConverter/OGSFileConverter.cpp b/Applications/Utils/OGSFileConverter/OGSFileConverter.cpp
index fee9e4b17cce57ecbcac00bf6350c92828abb4c3..60139a9ec74f5260d50adbe77969cfb823ed1569 100644
--- a/Applications/Utils/OGSFileConverter/OGSFileConverter.cpp
+++ b/Applications/Utils/OGSFileConverter/OGSFileConverter.cpp
@@ -19,11 +19,11 @@
 #include <QFileInfo>
 
 #include "Applications/DataExplorer/Base/OGSError.h"
+#include "Applications/FileIO/Legacy/OGSIOVer4.h"
 #include "BaseLib/FileTools.h"
 #include "BaseLib/StringTools.h"
 #include "FileListDialog.h"
 #include "GeoLib/GEOObjects.h"
-#include "GeoLib/IO/Legacy/OGSIOVer4.h"
 #include "GeoLib/IO/XmlIO/Qt/XmlGmlInterface.h"
 #include "MeshLib/IO/Legacy/MeshIO.h"
 #include "MeshLib/IO/VtkIO/VtuInterface.h"
@@ -60,7 +60,7 @@ void OGSFileConverter::convertGML2GLI(const QStringList &input, const QString &o
         }
         std::vector<std::string> geo_names;
         geo_objects.getGeometryNames(geo_names);
-        GeoLib::IO::Legacy::writeGLIFileV4(output_str, geo_names[0], geo_objects);
+        FileIO::Legacy::writeGLIFileV4(output_str, geo_names[0], geo_objects);
         geo_objects.removeSurfaceVec(geo_names[0]);
         geo_objects.removePolylineVec(geo_names[0]);
         geo_objects.removePointVec(geo_names[0]);
@@ -87,7 +87,7 @@ void OGSFileConverter::convertGLI2GML(const QStringList &input, const QString &o
         std::string unique_name;
         std::vector<std::string> errors;
 
-        GeoLib::IO::Legacy::readGLIFileV4(input_string.toStdString(),
+        FileIO::Legacy::readGLIFileV4(input_string.toStdString(),
                                           geo_objects, unique_name, errors);
         if (errors.empty() ||
             (errors.size() == 1 &&
diff --git a/GeoLib/CMakeLists.txt b/GeoLib/CMakeLists.txt
index f7b51618a7a2056ed0d03bb82baef72418e71c56..957d6c9043fc2f6999db05989831fa2bdac38ccf 100644
--- a/GeoLib/CMakeLists.txt
+++ b/GeoLib/CMakeLists.txt
@@ -2,7 +2,6 @@
 GET_SOURCE_FILES(SOURCES)
 
 APPEND_SOURCE_FILES(SOURCES IO)
-APPEND_SOURCE_FILES(SOURCES IO/Legacy)
 
 APPEND_SOURCE_FILES(SOURCES IO/XmlIO/Rapid)
 APPEND_SOURCE_FILES(SOURCES IO/XmlIO/Boost)
diff --git a/GeoLib/EarClippingTriangulation.cpp b/GeoLib/EarClippingTriangulation.cpp
deleted file mode 100644
index ff94ca1fb9f982622539d5d1cbc2e8fee3c4baa2..0000000000000000000000000000000000000000
--- a/GeoLib/EarClippingTriangulation.cpp
+++ /dev/null
@@ -1,298 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date   2011-02-23
- * \brief  Implementation of the EarClippingTriangulation class.
- *
- * \copyright
- * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#include "EarClippingTriangulation.h"
-
-#include "BaseLib/uniqueInsert.h"
-
-#include "MathLib/GeometricBasics.h"
-
-#include "Polygon.h"
-#include "Triangle.h"
-#include "Point.h"
-
-namespace GeoLib
-{
-EarClippingTriangulation::EarClippingTriangulation(const GeoLib::Polygon* polygon,
-        std::list<GeoLib::Triangle> &triangles, bool rot)
-{
-    copyPolygonPoints (polygon);
-
-    if (rot) {
-        rotatePointsToXY (_pnts);
-        ensureCWOrientation ();
-    }
-
-    initVertexList ();
-    initLists ();
-    clipEars ();
-
-    std::vector<GeoLib::Point*> const& ref_pnts_vec (polygon->getPointsVec());
-    std::list<GeoLib::Triangle>::const_iterator it (_triangles.begin());
-    if (_original_orient == GeoLib::CW) {
-        while (it != _triangles.end()) {
-            const std::size_t i0 (polygon->getPointID ((*it)[0]));
-            const std::size_t i1 (polygon->getPointID ((*it)[1]));
-            const std::size_t i2 (polygon->getPointID ((*it)[2]));
-            triangles.emplace_back(ref_pnts_vec, i0, i1, i2);
-            ++it;
-        }
-    } else {
-        std::size_t n_pnts (_pnts.size() - 1);
-        while (it != _triangles.end()) {
-            const std::size_t i0 (polygon->getPointID (n_pnts - (*it)[0]));
-            const std::size_t i1 (polygon->getPointID (n_pnts - (*it)[1]));
-            const std::size_t i2 (polygon->getPointID (n_pnts - (*it)[2]));
-            triangles.emplace_back(ref_pnts_vec, i0, i1, i2);
-            ++it;
-        }
-    }
-}
-
-EarClippingTriangulation::~EarClippingTriangulation()
-{
-    const std::size_t n_pnts (_pnts.size());
-    for (std::size_t k(0); k<n_pnts; k++) {
-        delete _pnts[k];
-    }
-}
-
-void EarClippingTriangulation::copyPolygonPoints (const GeoLib::Polygon* polygon)
-{
-    // copy points - last point is identical to the first
-    std::size_t n_pnts (polygon->getNumberOfPoints() - 1);
-    for (std::size_t k(0); k < n_pnts; k++) {
-        _pnts.push_back (new GeoLib::Point (*(polygon->getPoint(k))));
-    }
-}
-
-void EarClippingTriangulation::ensureCWOrientation ()
-{
-    std::size_t n_pnts (_pnts.size());
-    // get the left most upper point
-    std::size_t min_x_max_y_idx (0); // for orientation check
-    for (std::size_t k(0); k<n_pnts; k++) {
-        if ((*(_pnts[k]))[0] <= (*(_pnts[min_x_max_y_idx]))[0]) {
-            if ((*(_pnts[k]))[0] < (*(_pnts[min_x_max_y_idx]))[0]) {
-                min_x_max_y_idx = k;
-            } else {
-                if ((*(_pnts[k]))[1] > (*(_pnts[min_x_max_y_idx]))[1]) {
-                    min_x_max_y_idx = k;
-                }
-            }
-        }
-    }
-    // determine orientation
-    if (0 < min_x_max_y_idx && min_x_max_y_idx < n_pnts-1) {
-        _original_orient = GeoLib::getOrientation (
-            _pnts[min_x_max_y_idx-1], _pnts[min_x_max_y_idx], _pnts[min_x_max_y_idx+1]);
-    } else {
-        if (0 == min_x_max_y_idx) {
-            _original_orient = GeoLib::getOrientation (_pnts[n_pnts-1], _pnts[0], _pnts[1]);
-        } else {
-            _original_orient = GeoLib::getOrientation (_pnts[n_pnts-2], _pnts[n_pnts-1], _pnts[0]);
-        }
-    }
-    if (_original_orient == GeoLib::CCW) {
-        // switch orientation
-        for (std::size_t k(0); k<n_pnts/2; k++) {
-            std::swap (_pnts[k], _pnts[n_pnts-1-k]);
-        }
-    }
-}
-
-bool EarClippingTriangulation::isEar(std::size_t v0, std::size_t v1, std::size_t v2) const
-{
-    for (unsigned long v : _vertex_list)
-    {
-        if (v != v0 && v != v1 && v != v2)
-        {
-            if (MathLib::isPointInTriangle(*_pnts[v], *_pnts[v0], *_pnts[v1],
-                                           *_pnts[v2]))
-            {
-                return false;
-            }
-        }
-    }
-    return true;
-}
-
-void EarClippingTriangulation::initVertexList ()
-{
-    std::size_t n_pnts (_pnts.size());
-    for (std::size_t k(0); k<n_pnts; k++)
-        _vertex_list.push_back (k);
-}
-
-void EarClippingTriangulation::initLists ()
-{
-    // go through points checking ccw, cw or collinear order and identifying ears
-    std::list<std::size_t>::iterator it (_vertex_list.begin()), prev(_vertex_list.end()), next;
-    --prev;
-    next = it;
-    ++next;
-    GeoLib::Orientation orientation;
-    bool first_run(true); // saves special handling of the last case with identical code
-    while (_vertex_list.size() >= 3 && first_run) {
-        if (next == _vertex_list.end()) {
-            first_run = false;
-            next = _vertex_list.begin();
-        }
-        orientation  = getOrientation (_pnts[*prev], _pnts[*it], _pnts[*next]);
-        if (orientation == GeoLib::COLLINEAR) {
-            WARN("EarClippingTriangulation::initLists(): collinear points (%f, %f, %f), (%f, %f, %f), (%f, %f, %f)",
-                    (*_pnts[*prev])[0], (*_pnts[*prev])[1], (*_pnts[*prev])[2],
-                    (*_pnts[*it])[0], (*_pnts[*it])[1], (*_pnts[*it])[2],
-                    (*_pnts[*next])[0], (*_pnts[*next])[1], (*_pnts[*next])[2]);
-            it = _vertex_list.erase (it);
-            ++next;
-        } else {
-            if (orientation == GeoLib::CW) {
-                _convex_vertex_list.push_back (*it);
-                if (isEar (*prev, *it, *next))
-                    _ear_list.push_back (*it);
-            }
-            prev = it;
-            it = next;
-            ++next;
-        }
-    }
-}
-
-void EarClippingTriangulation::clipEars()
-{
-    std::list<std::size_t>::iterator it, prev, next;
-    // *** clip an ear
-    while (_vertex_list.size() > 3) {
-        // pop ear from list
-        std::size_t ear = _ear_list.front();
-        _ear_list.pop_front();
-        // remove ear tip from _convex_vertex_list
-        _convex_vertex_list.remove(ear);
-
-        // remove ear from vertex_list, apply changes to _ear_list, _convex_vertex_list
-        bool nfound(true);
-        it = _vertex_list.begin();
-        prev = _vertex_list.end();
-        --prev;
-        while (nfound && it != _vertex_list.end()) {
-            if (*it == ear) {
-                nfound = false;
-                it = _vertex_list.erase(it); // remove ear tip
-                next = it;
-                if (next == _vertex_list.end()) {
-                    next = _vertex_list.begin();
-                    prev = _vertex_list.end();
-                    --prev;
-                }
-                // add triangle
-                _triangles.emplace_back(_pnts, *prev, *next, ear);
-
-                // check the orientation of prevprev, prev, next
-                std::list<std::size_t>::iterator prevprev;
-                if (prev == _vertex_list.begin()) {
-                    prevprev = _vertex_list.end();
-                } else {
-                    prevprev = prev;
-                }
-                --prevprev;
-
-                // apply changes to _convex_vertex_list and _ear_list looking "backward"
-                GeoLib::Orientation orientation = GeoLib::getOrientation(_pnts[*prevprev], _pnts[*prev],
-                        _pnts[*next]);
-                if (orientation == GeoLib::CW) {
-                    BaseLib::uniquePushBack(_convex_vertex_list, *prev);
-                    // prev is convex
-                    if (isEar(*prevprev, *prev, *next)) {
-                        // prev is an ear tip
-                        BaseLib::uniquePushBack(_ear_list, *prev);
-                    } else {
-                        // if necessary remove prev
-                        _ear_list.remove(*prev);
-                    }
-                } else {
-                    // prev is not convex => reflex or collinear
-                    _convex_vertex_list.remove(*prev);
-                    _ear_list.remove(*prev);
-                    if (orientation == GeoLib::COLLINEAR) {
-                        prev = _vertex_list.erase(prev);
-                        if (prev == _vertex_list.begin()) {
-                            prev = _vertex_list.end();
-                            --prev;
-                        } else {
-                            --prev;
-                        }
-                    }
-                }
-
-                // check the orientation of prev, next, nextnext
-                std::list<std::size_t>::iterator nextnext,
-                        help_it(_vertex_list.end());
-                --help_it;
-                if (next == help_it) {
-                    nextnext = _vertex_list.begin();
-                } else {
-                    nextnext = next;
-                    ++nextnext;
-                }
-
-                // apply changes to _convex_vertex_list and _ear_list looking "forward"
-                orientation = getOrientation(_pnts[*prev], _pnts[*next],
-                        _pnts[*nextnext]);
-                if (orientation == GeoLib::CW) {
-                    BaseLib::uniquePushBack(_convex_vertex_list, *next);
-                    // next is convex
-                    if (isEar(*prev, *next, *nextnext)) {
-                        // next is an ear tip
-                        BaseLib::uniquePushBack(_ear_list, *next);
-                    } else {
-                        // if necessary remove *next
-                        _ear_list.remove(*next);
-                    }
-                } else {
-                    // next is not convex => reflex or collinear
-                    _convex_vertex_list.remove(*next);
-                    _ear_list.remove(*next);
-                    if (orientation == GeoLib::COLLINEAR) {
-                        next = _vertex_list.erase(next);
-                        if (next == _vertex_list.end())
-                            next = _vertex_list.begin();
-                    }
-                }
-            } else {
-                prev = it;
-                ++it;
-            }
-        }
-
-    }
-
-    // add last triangle
-    next = _vertex_list.begin();
-    prev = next;
-    ++next;
-    if (next == _vertex_list.end())
-        return;
-    it = next;
-    ++next;
-    if (next == _vertex_list.end())
-        return;
-
-    if (getOrientation(_pnts[*prev], _pnts[*it], _pnts[*next]) == GeoLib::CCW)
-        _triangles.emplace_back(_pnts, *prev, *it, *next);
-    else
-        _triangles.emplace_back(_pnts, *prev, *next, *it);
-}
-
-} // end namespace GeoLib
diff --git a/GeoLib/EarClippingTriangulation.h b/GeoLib/EarClippingTriangulation.h
deleted file mode 100644
index 839887766bd449edab857e77f1792e1ba286b216..0000000000000000000000000000000000000000
--- a/GeoLib/EarClippingTriangulation.h
+++ /dev/null
@@ -1,63 +0,0 @@
-/**
- * \file
- * \author Thomas Fischer
- * \date   2011-02-23
- * \brief  Definition of the EarClippingTriangulation class.
- *
- * \copyright
- * Copyright (c) 2012-2017, OpenGeoSys Community (http://www.opengeosys.org)
- *            Distributed under a Modified BSD License.
- *              See accompanying file LICENSE.txt or
- *              http://www.opengeosys.org/project/license
- *
- */
-
-#pragma once
-
-#include <list>
-#include <vector>
-
-#include "AnalyticalGeometry.h"
-
-namespace GeoLib
-{
-
-class Polygon;
-class Triangle;
-
-class EarClippingTriangulation
-{
-public:
-    EarClippingTriangulation(const GeoLib::Polygon* ply,
-                             std::list<GeoLib::Triangle> &triangles,
-                             bool rot = true);
-    virtual ~EarClippingTriangulation();
-private:
-    /**
-     * copies the points of the polygon to the vector _pnts
-     */
-    inline void copyPolygonPoints (const GeoLib::Polygon* polygon);
-    inline void ensureCWOrientation ();
-
-    inline bool isEar(std::size_t v0, std::size_t v1, std::size_t v2) const;
-
-    inline void initVertexList ();
-    inline void initLists ();
-    inline void clipEars ();
-
-    /**
-     * a copy of the polygon points
-     */
-    std::vector<GeoLib::Point*> _pnts;
-    std::list<std::size_t> _vertex_list;
-    std::list<std::size_t> _convex_vertex_list;
-    std::list<std::size_t> _ear_list;
-
-    /**
-     * triangles of the triangulation (maybe in the wrong orientation)
-     */
-    std::list<GeoLib::Triangle> _triangles;
-
-    GeoLib::Orientation _original_orient;
-};
-} // end namespace GeoLib
diff --git a/GeoLib/GEOObjects.cpp b/GeoLib/GEOObjects.cpp
index abc1caff65160c6f92984192d13092a2c5fc7ff1..80f663aaf2eb6aefee53ba60f0a6c257bbd26cae 100644
--- a/GeoLib/GEOObjects.cpp
+++ b/GeoLib/GEOObjects.cpp
@@ -511,6 +511,36 @@ void GEOObjects::mergeSurfaces(std::vector<std::string> const & geo_names,
     }
 }
 
+void GEOObjects::renameGeometry(std::string const& old_name,
+                                std::string const& new_name)
+{
+    _callbacks->renameGeometry(old_name, new_name);
+    for (auto* pnt_vec : _pnt_vecs)
+    {
+        if (pnt_vec->getName() == old_name)
+        {
+            pnt_vec->setName(new_name);
+            break;
+        }
+    }
+    for (auto* ply_vec : _ply_vecs)
+    {
+        if (ply_vec->getName() == old_name)
+        {
+            ply_vec->setName(new_name);
+            break;
+        }
+    }
+    for (auto* sfc_vec : _sfc_vecs)
+    {
+        if (sfc_vec->getName() == old_name)
+        {
+            sfc_vec->setName(new_name);
+            break;
+        }
+    }
+}
+
 const GeoLib::GeoObject* GEOObjects::getGeoObject(const std::string &geo_name,
                                                             GeoLib::GEOTYPE type,
                                                             const std::string &geo_obj_name) const
diff --git a/GeoLib/GEOObjects.h b/GeoLib/GEOObjects.h
index 3311627bd5f5f85d7df60260f1817891d02cfcfc..b9f2088fd3f6b6ba3423e7b0fe071040bab29667 100644
--- a/GeoLib/GEOObjects.h
+++ b/GeoLib/GEOObjects.h
@@ -74,6 +74,7 @@ public:
         virtual void addSurfaceVec(std::string const&){}
         virtual void appendSurfaceVec(std::string const&){}
         virtual void removeSurfaceVec(std::string const&){}
+        virtual void renameGeometry(std::string const&, std::string const&) {}
         virtual ~Callbacks() = default;
     };
 
@@ -192,8 +193,8 @@ public:
      * \return true if the surfaces are appended, false if the SurfaceVec with the
      * corresponding name does not exist and the surfaces are added.
      * */
-    bool appendSurfaceVec(const std::vector<Surface*> &surfaces,
-                                  const std::string &name);
+    bool appendSurfaceVec(const std::vector<Surface*>& surfaces,
+                          const std::string& name);
 
     /// Returns the surface vector with the given name as a const.
     const std::vector<Surface*>* getSurfaceVec(const std::string &name) const;
@@ -241,6 +242,13 @@ public:
      */
     int mergeGeometries(std::vector<std::string> const & names, std::string &merged_geo_name);
 
+    /// Renames an existing geometry, i.e. renames the internal PointVec,
+    /// PolylineVec and the SurfaceVec objects from \c old_name to \c new_name.
+    /// If no such PointVec, PolylineVec and SurfaceVec objects exist nothing
+    /// will happen.
+    void renameGeometry(std::string const& old_name,
+                        std::string const& new_name);
+
     /// Returns the geo object for a geometric item of the given name and type for the associated geometry.
     const GeoLib::GeoObject* getGeoObject(const std::string &geo_name,
                                           GeoLib::GEOTYPE type,
diff --git a/GeoLib/Surface.cpp b/GeoLib/Surface.cpp
index 2790f123fc91b140038ff2f8e5d032ac2c9db54f..cfe03c2be6541d2cc9865391d1868a9a82739972 100644
--- a/GeoLib/Surface.cpp
+++ b/GeoLib/Surface.cpp
@@ -22,7 +22,6 @@
 
 #include "Triangle.h"
 #include "Polyline.h"
-#include "EarClippingTriangulation.h"
 
 namespace GeoLib
 {
@@ -84,56 +83,6 @@ void Surface::addTriangle(std::size_t pnt_a,
     }
 }
 
-Surface* Surface::createSurface(const Polyline& ply)
-{
-    if (!ply.isClosed())
-    {
-        WARN("Error in Surface::createSurface() - Polyline is not closed.");
-        return nullptr;
-    }
-
-    if (ply.getNumberOfPoints() <= 2)
-    {
-        WARN(
-            "Error in Surface::createSurface() - Polyline consists of less "
-            "than three points and therefore cannot be triangulated.");
-        return nullptr;
-    }
-
-    // create empty surface
-    auto* sfc(new Surface(ply.getPointsVec()));
-
-    auto* polygon(new Polygon(ply));
-    polygon->computeListOfSimplePolygons();
-
-    // create surfaces from simple polygons
-    const std::list<GeoLib::Polygon*>& list_of_simple_polygons(
-        polygon->getListOfSimplePolygons());
-    for (auto simple_polygon : list_of_simple_polygons)
-    {
-        std::list<GeoLib::Triangle> triangles;
-        GeoLib::EarClippingTriangulation(simple_polygon, triangles);
-
-        // add Triangles to Surface
-        std::list<GeoLib::Triangle>::const_iterator it(triangles.begin());
-        while (it != triangles.end())
-        {
-            sfc->addTriangle((*it)[0], (*it)[1], (*it)[2]);
-            ++it;
-        }
-    }
-    delete polygon;
-    if (sfc->getNumberOfTriangles() == 0)
-    {
-        WARN(
-            "Surface::createSurface(): Triangulation does not contain any "
-            "triangle.");
-        delete sfc;
-        return nullptr;
-    }
-    return sfc;
-}
-
 std::size_t Surface::getNumberOfTriangles() const
 {
     return _sfc_triangles.size();
diff --git a/GeoLib/Surface.h b/GeoLib/Surface.h
index 54f3a455cde56eb3bb4270d6f0c5659bf4ec4218..97b6f6ec5fbbf0d9c01f5be6efb6d41f71feb081 100644
--- a/GeoLib/Surface.h
+++ b/GeoLib/Surface.h
@@ -49,9 +49,6 @@ public:
      * */
     void addTriangle(std::size_t pnt_a, std::size_t pnt_b, std::size_t pnt_c);
 
-    /// Triangulates a new surface based on closed polyline.
-    static Surface* createSurface(const Polyline& ply);
-
     /**
      * returns the number of triangles describing the Surface
      * */
diff --git a/Tests/GeoLib/IO/TestGLIReader.cpp b/Tests/GeoLib/IO/TestGLIReader.cpp
index 74e4069786702eed4db2fb5360caab09f111abfd..43d5a4f71418e2c82c6013428831d4b0072e6f97 100644
--- a/Tests/GeoLib/IO/TestGLIReader.cpp
+++ b/Tests/GeoLib/IO/TestGLIReader.cpp
@@ -16,7 +16,7 @@
 #include "gtest/gtest.h"
 
 #include "BaseLib/BuildInfo.h"
-#include "GeoLib/IO/Legacy/OGSIOVer4.h"
+#include "Applications/FileIO/Legacy/OGSIOVer4.h"
 #include "GeoLib/GEOObjects.h"
 
 class OGSIOVer4InterfaceTest : public ::testing::Test
@@ -55,7 +55,8 @@ TEST_F(OGSIOVer4InterfaceTest, SimpleTIN)
     GeoLib::GEOObjects geometries;
     std::vector<std::string> errors;
     std::string geometry_name("TestGeometry");
-    GeoLib::IO::Legacy::readGLIFileV4(_gli_fname, geometries, geometry_name, errors);
+    FileIO::Legacy::readGLIFileV4(_gli_fname, geometries, geometry_name,
+                                  errors);
 
     std::vector<GeoLib::Surface*> const*
         sfcs(geometries.getSurfaceVec(geometry_name));
@@ -78,7 +79,8 @@ TEST_F(OGSIOVer4InterfaceTest, StillCorrectTINWihtAdditionalValueAtEndOfLine)
     GeoLib::GEOObjects geometries;
     std::vector<std::string> errors;
     std::string geometry_name("TestGeometry");
-    GeoLib::IO::Legacy::readGLIFileV4(_gli_fname, geometries, geometry_name, errors);
+    FileIO::Legacy::readGLIFileV4(_gli_fname, geometries, geometry_name,
+                                  errors);
 
     std::vector<GeoLib::Surface*> const*
         sfcs(geometries.getSurfaceVec(geometry_name));
@@ -100,7 +102,8 @@ TEST_F(OGSIOVer4InterfaceTest, InvalidTIN_ZeroAreaTri)
     GeoLib::GEOObjects geometries;
     std::vector<std::string> errors;
     std::string geometry_name("TestGeometry");
-    GeoLib::IO::Legacy::readGLIFileV4(_gli_fname, geometries, geometry_name, errors);
+    FileIO::Legacy::readGLIFileV4(_gli_fname, geometries, geometry_name,
+                                  errors);
 
     std::vector<GeoLib::Surface*> const*
         sfcs(geometries.getSurfaceVec(geometry_name));
@@ -123,7 +126,8 @@ TEST_F(OGSIOVer4InterfaceTest, InvalidTIN_LineDoesNotStartWithID)
     GeoLib::GEOObjects geometries;
     std::vector<std::string> errors;
     std::string geometry_name("TestGeometry");
-    GeoLib::IO::Legacy::readGLIFileV4(_gli_fname, geometries, geometry_name, errors);
+    FileIO::Legacy::readGLIFileV4(_gli_fname, geometries, geometry_name,
+                                  errors);
 
     std::vector<GeoLib::Surface*> const*
         sfcs(geometries.getSurfaceVec(geometry_name));
@@ -145,7 +149,8 @@ TEST_F(OGSIOVer4InterfaceTest, InvalidTIN_PointIsMissing)
     GeoLib::GEOObjects geometries;
     std::vector<std::string> errors;
     std::string geometry_name("TestGeometry");
-    GeoLib::IO::Legacy::readGLIFileV4(_gli_fname, geometries, geometry_name, errors);
+    FileIO::Legacy::readGLIFileV4(_gli_fname, geometries, geometry_name,
+                                  errors);
 
     std::vector<GeoLib::Surface*> const*
         sfcs(geometries.getSurfaceVec(geometry_name));
@@ -166,7 +171,8 @@ TEST_F(OGSIOVer4InterfaceTest, InvalidTIN_CoordOfPointIsMissing)
     GeoLib::GEOObjects geometries;
     std::vector<std::string> errors;
     std::string geometry_name("TestGeometry");
-    GeoLib::IO::Legacy::readGLIFileV4(_gli_fname, geometries, geometry_name, errors);
+    FileIO::Legacy::readGLIFileV4(_gli_fname, geometries, geometry_name,
+                                  errors);
 
     std::vector<GeoLib::Surface*> const*
         sfcs(geometries.getSurfaceVec(geometry_name));
@@ -187,7 +193,8 @@ TEST_F(OGSIOVer4InterfaceTest, SimpleTIN_AdditionalEmptyLinesAtEnd)
     GeoLib::GEOObjects geometries;
     std::vector<std::string> errors;
     std::string geometry_name("TestGeometry");
-    GeoLib::IO::Legacy::readGLIFileV4(_gli_fname, geometries, geometry_name, errors);
+    FileIO::Legacy::readGLIFileV4(_gli_fname, geometries, geometry_name,
+                                  errors);
 
     std::vector<GeoLib::Surface*> const*
         sfcs(geometries.getSurfaceVec(geometry_name));
diff --git a/Tests/MeshLib/TestBoundaryElementSearch.cpp b/Tests/MeshLib/TestBoundaryElementSearch.cpp
index eb48b5f6039eba566a48cc6af44453516c086442..5f706051bb18deaa5445760e041994427094edf6 100644
--- a/Tests/MeshLib/TestBoundaryElementSearch.cpp
+++ b/Tests/MeshLib/TestBoundaryElementSearch.cpp
@@ -131,32 +131,27 @@ TEST_F(MeshLibBoundaryElementSearchInSimpleHexMesh, SurfaceSearch)
     pnts.push_back(new GeoLib::Point(_geometric_size, 0.0, _geometric_size));
     pnts.push_back(new GeoLib::Point(0.0, 0.0, _geometric_size));
 
-    GeoLib::Polyline ply_bottom(pnts);
-    ply_bottom.addPoint(0);
-    ply_bottom.addPoint(1);
-    ply_bottom.addPoint(2);
-    ply_bottom.addPoint(3);
-    ply_bottom.addPoint(0);
-    std::unique_ptr<GeoLib::Surface> sfc_bottom(GeoLib::Surface::createSurface(ply_bottom));
-
-    GeoLib::Polyline ply_front(pnts);
-    ply_front.addPoint(0);
-    ply_front.addPoint(1);
-    ply_front.addPoint(4);
-    ply_front.addPoint(5);
-    ply_front.addPoint(0);
-    std::unique_ptr<GeoLib::Surface> sfc_front(GeoLib::Surface::createSurface(ply_front));
+    GeoLib::Surface sfc_bottom(pnts);
+    sfc_bottom.addTriangle(0, 1, 2);
+    sfc_bottom.addTriangle(0, 2, 3);
+
+    GeoLib::Surface sfc_front(pnts);
+    sfc_front.addTriangle(0, 1, 4);
+    sfc_front.addTriangle(0, 4, 5);
 
     // perform search on the bottom surface
     MeshGeoToolsLib::MeshNodeSearcher mesh_node_searcher(
         *_hex_mesh,
         std::make_unique<MeshGeoToolsLib::HeuristicSearchLength>(*_hex_mesh),
         MeshGeoToolsLib::SearchAllNodes::Yes);
-    MeshGeoToolsLib::BoundaryElementsSearcher boundary_element_searcher(*_hex_mesh, mesh_node_searcher);
-    std::vector<MeshLib::Element*> const& found_faces_sfc_b(boundary_element_searcher.getBoundaryElements(*sfc_bottom));
+    MeshGeoToolsLib::BoundaryElementsSearcher boundary_element_searcher(
+        *_hex_mesh, mesh_node_searcher);
+    std::vector<MeshLib::Element*> const& found_faces_sfc_b(
+        boundary_element_searcher.getBoundaryElements(sfc_bottom));
     ASSERT_EQ(n_eles_2d, found_faces_sfc_b.size());
-    double sum_area_b = std::accumulate(found_faces_sfc_b.begin(), found_faces_sfc_b.end(), 0.0,
-                [](double v, MeshLib::Element*e){return v+e->getContent();});
+    double sum_area_b = std::accumulate(
+        found_faces_sfc_b.begin(), found_faces_sfc_b.end(), 0.0,
+        [](double v, MeshLib::Element* e) { return v + e->getContent(); });
     ASSERT_EQ(_geometric_size*_geometric_size, sum_area_b);
     auto connected_nodes_b = MeshLib::getUniqueNodes(found_faces_sfc_b);
     ASSERT_EQ(n_nodes_2d, connected_nodes_b.size());
@@ -164,10 +159,12 @@ TEST_F(MeshLibBoundaryElementSearchInSimpleHexMesh, SurfaceSearch)
         ASSERT_EQ(0.0, (*node)[2]); // check z coordinates
 
     // perform search on the front surface
-    std::vector<MeshLib::Element*> const& found_faces_sfc_f(boundary_element_searcher.getBoundaryElements(*sfc_front));
+    std::vector<MeshLib::Element*> const& found_faces_sfc_f(
+        boundary_element_searcher.getBoundaryElements(sfc_front));
     ASSERT_EQ(n_eles_2d, found_faces_sfc_f.size());
-    double sum_area_f = std::accumulate(found_faces_sfc_f.begin(), found_faces_sfc_f.end(), 0.0,
-                [](double v, MeshLib::Element*e){return v+e->getContent();});
+    double sum_area_f = std::accumulate(
+        found_faces_sfc_f.begin(), found_faces_sfc_f.end(), 0.0,
+        [](double v, MeshLib::Element* e) { return v + e->getContent(); });
     ASSERT_EQ(_geometric_size*_geometric_size, sum_area_f);
     auto connected_nodes_f = MeshLib::getUniqueNodes(found_faces_sfc_f);
     ASSERT_EQ(n_nodes_2d, connected_nodes_f.size());
diff --git a/Tests/MeshLib/TestMeshNodeSearch.cpp b/Tests/MeshLib/TestMeshNodeSearch.cpp
index 2beaddfdbcdd4016bddae0e41e024978c1a1c8c5..ccfcaed583d949ebad3341715919f959c7346bb7 100644
--- a/Tests/MeshLib/TestMeshNodeSearch.cpp
+++ b/Tests/MeshLib/TestMeshNodeSearch.cpp
@@ -227,30 +227,24 @@ TEST_F(MeshLibMeshNodeSearchInSimpleQuadMesh, SurfaceSearch)
         std::move(search_length), MeshGeoToolsLib::SearchAllNodes::Yes);
 
     // entire domain
-    GeoLib::Polyline ply0(pnts);
-    ply0.addPoint(0);
-    ply0.addPoint(1);
-    ply0.addPoint(2);
-    ply0.addPoint(3);
-    ply0.addPoint(0);
-    std::unique_ptr<GeoLib::Surface> sfc0(GeoLib::Surface::createSurface(ply0));
+    GeoLib::Surface sfc0(pnts);
+    sfc0.addTriangle(0, 1, 2);
+    sfc0.addTriangle(0, 2, 3);
 
-    std::vector<std::size_t> const& found_ids_sfc0(mesh_node_searcher.getMeshNodeIDsAlongSurface(*sfc0));
+    std::vector<std::size_t> const& found_ids_sfc0(
+        mesh_node_searcher.getMeshNodeIDsAlongSurface(sfc0));
 
     ASSERT_EQ(_quad_mesh->getNumberOfNodes(), found_ids_sfc0.size());
     for (std::size_t k(0); k<found_ids_sfc0.size(); k++)
         ASSERT_EQ(k, found_ids_sfc0[k]);
 
     // bottom half domain
-    GeoLib::Polyline ply1(pnts);
-    ply1.addPoint(0);
-    ply1.addPoint(1);
-    ply1.addPoint(4);
-    ply1.addPoint(5);
-    ply1.addPoint(0);
-    std::unique_ptr<GeoLib::Surface> sfc1(GeoLib::Surface::createSurface(ply1));
+    GeoLib::Surface sfc1(pnts);
+    sfc1.addTriangle(0, 1, 4);
+    sfc1.addTriangle(0, 4, 5);
 
-    std::vector<std::size_t> const& found_ids_sfc1(mesh_node_searcher.getMeshNodeIDsAlongSurface(*sfc1));
+    std::vector<std::size_t> const& found_ids_sfc1(
+        mesh_node_searcher.getMeshNodeIDsAlongSurface(sfc1));
 
     ASSERT_EQ(_quad_mesh->getNumberOfNodes()/2, found_ids_sfc1.size());
     for (std::size_t k(0); k<found_ids_sfc1.size(); k++)
@@ -282,44 +276,35 @@ TEST_F(MeshLibMeshNodeSearchInSimpleHexMesh, SurfaceSearch)
     const std::size_t n_nodes_3d = n_nodes_1d * n_nodes_1d * n_nodes_1d;
 
     // bottom surface
-    GeoLib::Polyline ply_bottom(pnts);
-    ply_bottom.addPoint(0);
-    ply_bottom.addPoint(1);
-    ply_bottom.addPoint(2);
-    ply_bottom.addPoint(3);
-    ply_bottom.addPoint(0);
-    std::unique_ptr<GeoLib::Surface> sfc_bottom(GeoLib::Surface::createSurface(ply_bottom));
-
-    std::vector<std::size_t> const& found_ids_sfc_b(mesh_node_searcher.getMeshNodeIDsAlongSurface(*sfc_bottom));
+    GeoLib::Surface sfc_bottom(pnts);
+    sfc_bottom.addTriangle(0, 1, 2);
+    sfc_bottom.addTriangle(0, 2, 3);
+
+    std::vector<std::size_t> const& found_ids_sfc_b(
+        mesh_node_searcher.getMeshNodeIDsAlongSurface(sfc_bottom));
     ASSERT_EQ(n_nodes_2d, found_ids_sfc_b.size());
     for (std::size_t k(0); k<found_ids_sfc_b.size(); k++)
         ASSERT_EQ(k, found_ids_sfc_b[k]);
 
     // top surface
-    GeoLib::Polyline ply_top(pnts);
-    ply_top.addPoint(4);
-    ply_top.addPoint(5);
-    ply_top.addPoint(6);
-    ply_top.addPoint(7);
-    ply_top.addPoint(4);
-    std::unique_ptr<GeoLib::Surface> sfc_top(GeoLib::Surface::createSurface(ply_top));
-
-    std::vector<std::size_t> const& found_ids_sfc_t(mesh_node_searcher.getMeshNodeIDsAlongSurface(*sfc_top));
+    GeoLib::Surface sfc_top(pnts);
+    sfc_top.addTriangle(4, 5, 6);
+    sfc_top.addTriangle(4, 6, 7);
+
+    std::vector<std::size_t> const& found_ids_sfc_t(
+        mesh_node_searcher.getMeshNodeIDsAlongSurface(sfc_top));
     ASSERT_EQ(n_nodes_2d, found_ids_sfc_t.size());
     const std::size_t offset_t = n_nodes_3d - n_nodes_2d;
     for (std::size_t k(0); k<found_ids_sfc_t.size(); k++)
         ASSERT_EQ(offset_t + k, found_ids_sfc_t[k]);
 
     // front
-    GeoLib::Polyline ply_front(pnts);
-    ply_front.addPoint(0);
-    ply_front.addPoint(1);
-    ply_front.addPoint(5);
-    ply_front.addPoint(4);
-    ply_front.addPoint(0);
-    std::unique_ptr<GeoLib::Surface> sfc_front(GeoLib::Surface::createSurface(ply_front));
-
-    std::vector<std::size_t> const& found_ids_sfc_f(mesh_node_searcher.getMeshNodeIDsAlongSurface(*sfc_front));
+    GeoLib::Surface sfc_front(pnts);
+    sfc_front.addTriangle(0, 1, 5);
+    sfc_front.addTriangle(0, 5, 4);
+
+    std::vector<std::size_t> const& found_ids_sfc_f(
+        mesh_node_searcher.getMeshNodeIDsAlongSurface(sfc_front));
     ASSERT_EQ(n_nodes_2d, found_ids_sfc_f.size());
     std::size_t cnt=0;
     for (std::size_t k(0); k<n_nodes_1d; k++) {
@@ -328,14 +313,12 @@ TEST_F(MeshLibMeshNodeSearchInSimpleHexMesh, SurfaceSearch)
     }
 
     // back
-    GeoLib::Polyline ply_back(pnts);
-    ply_back.addPoint(3);
-    ply_back.addPoint(2);
-    ply_back.addPoint(6);
-    ply_back.addPoint(7);
-    ply_back.addPoint(3);
-    std::unique_ptr<GeoLib::Surface> sfc_back(GeoLib::Surface::createSurface(ply_back));
-    std::vector<std::size_t> const& found_ids_sfc_back(mesh_node_searcher.getMeshNodeIDsAlongSurface(*sfc_back));
+    GeoLib::Surface sfc_back(pnts);
+    sfc_back.addTriangle(3, 2, 6);
+    sfc_back.addTriangle(3, 6, 7);
+
+    std::vector<std::size_t> const& found_ids_sfc_back(
+        mesh_node_searcher.getMeshNodeIDsAlongSurface(sfc_back));
     ASSERT_EQ(n_nodes_2d, found_ids_sfc_back.size());
     cnt = 0;
     const std::size_t y_offset = n_nodes_1d*(n_nodes_1d-1);
@@ -346,14 +329,12 @@ TEST_F(MeshLibMeshNodeSearchInSimpleHexMesh, SurfaceSearch)
     }
 
     // left
-    GeoLib::Polyline ply_left(pnts);
-    ply_left.addPoint(0);
-    ply_left.addPoint(3);
-    ply_left.addPoint(7);
-    ply_left.addPoint(4);
-    ply_left.addPoint(0);
-    std::unique_ptr<GeoLib::Surface> sfc_left(GeoLib::Surface::createSurface(ply_left));
-    std::vector<std::size_t> const& found_ids_sfc_left(mesh_node_searcher.getMeshNodeIDsAlongSurface(*sfc_left));
+    GeoLib::Surface sfc_left(pnts);
+    sfc_left.addTriangle(0, 3, 7);
+    sfc_left.addTriangle(0, 7, 4);
+
+    std::vector<std::size_t> const& found_ids_sfc_left(
+        mesh_node_searcher.getMeshNodeIDsAlongSurface(sfc_left));
     ASSERT_EQ(n_nodes_2d, found_ids_sfc_left.size());
     cnt = 0;
     for (std::size_t k(0); k<n_nodes_1d; k++) {
@@ -363,14 +344,12 @@ TEST_F(MeshLibMeshNodeSearchInSimpleHexMesh, SurfaceSearch)
     }
 
     // right
-    GeoLib::Polyline ply_right(pnts);
-    ply_right.addPoint(1);
-    ply_right.addPoint(2);
-    ply_right.addPoint(6);
-    ply_right.addPoint(5);
-    ply_right.addPoint(1);
-    std::unique_ptr<GeoLib::Surface> sfc_right(GeoLib::Surface::createSurface(ply_right));
-    std::vector<std::size_t> const& found_ids_sfc_right(mesh_node_searcher.getMeshNodeIDsAlongSurface(*sfc_right));
+    GeoLib::Surface sfc_right(pnts);
+    sfc_right.addTriangle(1, 2, 6);
+    sfc_right.addTriangle(1, 6, 5);
+
+    std::vector<std::size_t> const& found_ids_sfc_right(
+        mesh_node_searcher.getMeshNodeIDsAlongSurface(sfc_right));
     ASSERT_EQ(n_nodes_2d, found_ids_sfc_right.size());
     cnt = 0;
     for (std::size_t k(0); k<n_nodes_1d; k++) {
diff --git a/Tests/NumLib/TestDistribution.cpp b/Tests/NumLib/TestDistribution.cpp
index f8456f8410c3a96d56db6583867aca81c792b917..7a5c6481b7ae5796fdb68edda8224ec0f7d1eda2 100644
--- a/Tests/NumLib/TestDistribution.cpp
+++ b/Tests/NumLib/TestDistribution.cpp
@@ -62,7 +62,9 @@ public:
         plys->push_back(ply1);
 
         auto sfcs = std::make_unique<std::vector<GeoLib::Surface*>>();
-        _sfc1 = GeoLib::Surface::createSurface(*ply1);
+        _sfc1 = new GeoLib::Surface(*pnts);
+        _sfc1->addTriangle(0, 1, 2);
+        _sfc1->addTriangle(0, 2, 3);
         sfcs->push_back(_sfc1);
 
         _geo_objs.addPointVec(std::move(pnts), _project_name);
@@ -120,7 +122,9 @@ public:
         plys->push_back(ply1);
 
         auto sfcs = std::make_unique<std::vector<GeoLib::Surface*>>();
-        _sfc1 = GeoLib::Surface::createSurface(*ply1);
+        _sfc1 = new GeoLib::Surface(*pnts);
+        _sfc1->addTriangle(0, 3, 7);
+        _sfc1->addTriangle(0, 7, 4);
         sfcs->push_back(_sfc1);
 
         _geo_objs.addPointVec(std::move(pnts) ,_project_name);
diff --git a/Tests/NumLib/TestSpatialFunction.cpp b/Tests/NumLib/TestSpatialFunction.cpp
index b11a5bb977c73b48d94b426536765027e2a8779a..f73e9c233e1535472ebe44ba69890cc62023f1db 100644
--- a/Tests/NumLib/TestSpatialFunction.cpp
+++ b/Tests/NumLib/TestSpatialFunction.cpp
@@ -82,18 +82,16 @@ TEST(NumLib, SpatialFunctionInterpolationSurface)
     GeoLib::Point pt3(10.0, 10.0, 0.0);
     GeoLib::Point pt4(0.0, 10.0, 0.0);
     std::vector<GeoLib::Point*> pnts = {&pt1, &pt2, &pt3, &pt4};
-    GeoLib::Polyline ply0(pnts);
-    ply0.addPoint(0);
-    ply0.addPoint(1);
-    ply0.addPoint(2);
-    ply0.addPoint(3);
-    ply0.addPoint(0);
-    std::unique_ptr<GeoLib::Surface>  sfc1(GeoLib::Surface::createSurface(ply0));
+    GeoLib::Surface sfc1(pnts);
+    sfc1.addTriangle(0, 1, 2);
+    sfc1.addTriangle(0, 2, 3);
 
     // define a function
     const std::vector<std::size_t> vec_point_ids = {{0, 1, 2, 3}};
     const std::vector<double> vec_point_values = {{0., 100., 100., 0.}};
-    NumLib::LinearInterpolationOnSurface interpolate(*sfc1, vec_point_ids, vec_point_values, std::numeric_limits<double>::max());
+    NumLib::LinearInterpolationOnSurface interpolate(
+        sfc1, vec_point_ids, vec_point_values,
+        std::numeric_limits<double>::max());
 
     // normal
     for (unsigned k=0; k<10; k++) {