diff --git a/DESCRIPTION b/DESCRIPTION
index b845145259d9f6466b71c8106ab07871ef2e9a66..9cc1085ddba8483cb7b8ad6fde2b48d6b9e148f1 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -22,5 +22,6 @@ Imports:
     checkmate,
     tibble,
     R6,
-    stringr
+    stringr,
+    readr
 RoxygenNote: 7.1.1
diff --git a/NAMESPACE b/NAMESPACE
index bd01f7e413e772bc1f02e647e05350c48f2a6822..ad776bcc042535c0ece05efa7ec5ad15a0160e50 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,7 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
-export(data_to_xml)
-export(export_gml_to_file)
+export(export_xml_to_file)
+export(gml_data_to_xml)
 export(points_to_xml)
 export(polylines_to_xml)
 export(surfaces_to_xml)
diff --git a/R/create_gml.R b/R/create_gml.R
index c6222ad36931053df4c5db029fbbf334346ccb8f..4fa6a13b19941381f266440fb97fcb84cbc48c6c 100644
--- a/R/create_gml.R
+++ b/R/create_gml.R
@@ -1,13 +1,13 @@
 #This script contains various functions to turn data for a .gml file into the correct XML format
 
 #' Export function
-#' @param gml_data The .gml data (already in XML friendly format)
+#' @param xml_data The .gml data (already in XML friendly format)
 #' @param file_name The name of the .gml file to be written
 # @examples
-# export_gml_to_file(my_gml)
+# export_xml_to_file(my_gml)
 #' @export
-export_gml_to_file <- function(gml_data, file_name) {
-  doc <- xml2::as_xml_document(gml_data)
+export_xml_to_file <- function(xml_data, file_name) {
+  doc <- xml2::as_xml_document(xml_data)
   xml2::write_xml(doc, file_name, options = "format", encoding="ISO-8859-1")
   invisible()
 }
@@ -21,7 +21,7 @@ export_gml_to_file <- function(gml_data, file_name) {
 #' @return A XML document ready for export to a file
 # @examples (WIP)
 #' @export
-data_to_xml <- function(geo_name, points_tibble, polylines_list, surfaces_list) {
+gml_data_to_xml <- function(geo_name, points_tibble, polylines_list, surfaces_list) {
 
   #data_node <- xml2::read_xml('<OpenGeoSysGLI xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:ogs="http://www.opengeosys.org"/>')
 
@@ -42,14 +42,15 @@ data_to_xml <- function(geo_name, points_tibble, polylines_list, surfaces_list)
 #' @param point_tibble The specified tibble
 #' @return An XML node containing the points
 #' @examples
-#' my_tibble <- tibble::tibble(x = c(0, 0), y = c(0, 0), z = c(0, 1),is_origin = c(TRUE, FALSE))
+#' my_tibble <- tibble::tibble(x = c(0, 0), y = c(0, 0), z = c(0, 1), name = c("origin", ""))
 #' point_node <- points_to_xml(my_tibble)
 #' @export
 points_to_xml <- function(point_tibble) {
   points_node <- list(points = list())
+  has_names <- (length(point_tibble) == 4)
 
   for(i in 1:length(point_tibble[[1]])){
-    if(!point_tibble[[4]][[i]]) {
+    if(!has_names || point_tibble[[4]][[i]] == "") {
       points_node[[1]] <- c(points_node[[1]], list(point = structure(list(),
                                                                      id = (i-1),
                                                                      x = point_tibble[[1]][[i]],
@@ -61,7 +62,7 @@ points_to_xml <- function(point_tibble) {
                                                                      x = point_tibble[[1]][[i]],
                                                                      y = point_tibble[[2]][[i]],
                                                                      z = point_tibble[[3]][[i]],
-                                                                     name = "origin")))
+                                                                     name = point_tibble[[4]][[i]])))
     }
   }
 
@@ -79,8 +80,14 @@ polylines_to_xml <- function(polyline_list) {
   polylines_node <- list(polylines = list())
 
   for(i in 1:length(polyline_list)){
-    polylines_node[[1]] <- c(polylines_node[[1]], list(polyline = structure(list(pnt = list(polyline_list[[i]][[2]][[1]]),
-                                                                                 pnt = list(polyline_list[[i]][[2]][[2]])),
+
+    pnt_list <- list()
+
+    for(j in 1:length(polyline_list[[i]][[2]])) {
+      pnt_list <- c(pnt_list, list(pnt = list(polyline_list[[i]][[2]][[j]])))
+    }
+
+    polylines_node[[1]] <- c(polylines_node[[1]], list(polyline = structure(pnt_list,
                                                                             id = (i-1),
                                                                             name = polyline_list[[i]][[1]])))
   }
diff --git a/R/create_gml_R6.R b/R/create_gml_R6.R
index 1e0c7cc77604d07e007b6019976ffc1cdcacef58..e08168aacb8ece4df8ed93a9228a074eaf3c42be 100644
--- a/R/create_gml_R6.R
+++ b/R/create_gml_R6.R
@@ -223,7 +223,6 @@ GMLPoints <- R6::R6Class("GMLPoints",
                                  private$.ids <- self$validate_ids(gml_points)
                                  self$validate_dim(gml_points)
                                  self$validate_coordinates(gml_points)
-                                 self$validate_origin(gml_points)
 
                                  private$.gml_points <- gml_points
                              },
@@ -259,22 +258,6 @@ GMLPoints <- R6::R6Class("GMLPoints",
                                  }
                              },
 
-                             validate_origin = function(gml_points){
-                                 found_origin <- FALSE
-                                 for(i in 1:length(gml_points)){
-                                     if (gml_points[[i]]$is_origin == TRUE){
-                                         if(found_origin){
-                                             stop("More than one origin point detected", call. = FALSE)
-                                         }else{
-                                             found_origin <- TRUE
-                                         }
-                                     }
-                                 }
-                                 if(!found_origin){
-                                     stop("No origin point detected", call. = FALSE)
-                                 }
-                             },
-
                              as_list = function(...) {
                                  points_list <- list()
 
diff --git a/R/create_vtu.R b/R/create_vtu.R
new file mode 100644
index 0000000000000000000000000000000000000000..1befbdeb701ebbd27c5953dbda01fe064d069a31
--- /dev/null
+++ b/R/create_vtu.R
@@ -0,0 +1,16 @@
+#This script contains various functions to turn data for a .vtu file into the correct XML format
+
+vtu_data_to_xml <- function() {
+
+
+}
+
+
+
+vtu_piece_to_xml <- function(n_points, n_cells, point_data, cell_data, points, cells) {
+    piece_node <- list(piece = structure(list(),
+                                         NumberOfPoints = 0,
+                                         NumberOfCells = 0))
+    return(xml2::as_xml_document(piece_node))
+}
+
diff --git a/R/gml_prefabs.R b/R/gml_prefabs.R
new file mode 100644
index 0000000000000000000000000000000000000000..296f1004dc683dc49f1cb4b6affb815e744ead36
--- /dev/null
+++ b/R/gml_prefabs.R
@@ -0,0 +1,214 @@
+#This script contains some useful S3 classes for premade geometry
+
+#'Constructor based on HydroMechanics/IdealGas/flow_free_expansion/cube_1x1x1.gml
+new_r2ogs6_gml_cube_1x1x1 <- function() {
+    structure(list(
+        gml_geometry_name = "cube_1x1x1_geometry",
+
+        #Cube points (will become second XML child under root)
+        gml_points = tibble::tibble(x = c(0, 0, 0, 0, 1, 1, 1, 1),
+                                    y = c(0, 0, 1, 1, 0, 0, 1, 1),
+                                    z = c(0, 1, 1, 0, 0, 1, 1, 0),
+                                    name = c("origin", rep("", 7))),
+
+        #Cube polylines (will become third XML child under root)
+        gml_polylines = list(list(name = "front_left", c(0, 1)),
+                             list(name = "front_right", c(4, 5)),
+                             list(name = "front_bottom", c(0, 4)),
+                             list(name = "front_top", c(1, 5)),
+                             list(name = "bottom_left", c(0, 3)),
+                             list(name = "bottom_right", c(4, 7)),
+                             list(name = "top_left", c(1, 2)),
+                             list(name = "top_right", c(5, 6)),
+                             list(name = "back_left", c(2, 3)),
+                             list(name = "back_right", c(6, 7)),
+                             list(name = "back_bottom", c(3, 7)),
+                             list(name = "back_top", c(2, 6))),
+
+        #Cube surfaces (will become fourth XML child under root)
+        gml_surfaces = list(list(name = "left", c(0, 1, 2), c(0, 3, 2)),
+                            list(name = "right", c(4, 6, 5), c(4, 6, 7)),
+                            list(name = "top", c(1, 2, 5), c(5, 2, 6)),
+                            list(name = "bottom", c(0, 3, 4), c(4, 3, 7)),
+                            list(name = "front", c(0, 1, 4), c(4, 1, 5)),
+                            list(name = "back", c(2, 3, 6), c(6, 3, 7)))
+
+    ), class = "r2ogs6_gml_prefab")
+}
+
+
+#Another cube but without polylines
+
+
+#'Constructor based on HydroMechanics/IdealGas/flow_no_strain/square_1x1.gml
+new_r2ogs6_gml_square_1x1 <- function() {
+    structure(list(
+        gml_geometry_name = "square_1x1_geometry",
+        gml_points = tibble::tibble(x = c(0, 0, 1, 1),
+                                    y = c(0, 1, 0, 1),
+                                    z = c(0, 0, 0, 0),
+                                    name = c("origin", rep("", 3))),
+        gml_polylines = list(list(name = "left", c(0, 1)),
+                             list(name = "right", c(2, 3)),
+                             list(name = "bottom", c(0, 2)),
+                             list(name = "top", c(1, 3)))
+    ), class = "r2ogs6_gml_prefab")
+}
+
+
+#'Constructor based on HydroMechanics/IdealGas/flow_pressure_boundary/quad_1x10.gml
+new_r2ogs6_gml_quad_1x10 <- function() {
+    structure(list(
+        gml_geometry_name = "quad_1x10_geometry",
+        gml_points = tibble::tibble(x = c(0, 0, 10, 10),
+                                    y = c(0, 1, 0, 1),
+                                    z = c(0, 0, 0, 0),
+                                    name = c("origin", rep("", 3))),
+
+        gml_polylines = list(list(name = "left", c(0, 1)),
+                             list(name = "right", c(2, 3)),
+                             list(name = "bottom", c(0, 2)),
+                             list(name = "top", c(1, 3)))
+    ), class = "r2ogs6_gml_prefab")
+}
+
+
+#'Constructor based on HydroMechanics/StaggeredScheme/InjectionProduction1D/bar.gml
+new_r2ogs6_gml_bar <- function() {
+    structure(list(
+        gml_geometry_name = "bar",
+        gml_points = tibble::tibble(x = c(0, 10, 10, 0),
+                                    y = c(0, 0, 150, 150),
+                                    z = c(0, 0, 0, 0),
+                                    name = c(rep("", 4))),
+
+        gml_polylines = list(list(name = "left", c(0, 3)),
+                             list(name = "right", c(1, 2)),
+                             list(name = "bottom", c(0, 1)),
+                             list(name = "top", c(2, 3)))
+    ), class = "r2ogs6_gml_prefab")
+}
+
+
+#'Constructor based on HydroMechanics/Verification/hm1_1Dbeam.gml
+new_r2ogs6_gml_hm1_1Dbeam <- function() {
+    structure(list(
+        gml_geometry_name = "geometry",
+        gml_points = tibble::tibble(x = c(1, 1, 1, 1, 0, 0, 0, 0),
+                                    y = c(0.10000000000000001, 0, 0, 0.10000000000000001,
+                                          0, 0, 0.10000000000000001, 0.10000000000000001),
+                                    z = c(0, 0, 0.10000000000000001, 0.10000000000000001,
+                                          0, 0.10000000000000001, 0, 0.10000000000000001),
+                                    name = c("POINT0", "POINT1", "POINT2", "POINT3",
+                                             "POINT4", "POINT7", "POINT8", "POINT22")),
+
+        gml_polylines = list(list(name = "POLYLINE1", c(0, 1, 2, 3, 0)),
+                             list(name = "POLYLINE2", c(4, 1, 2, 5, 4)),
+                             list(name = "POLYLINE3", c(6, 0, 3, 7, 6)),
+                             list(name = "POLYLINE4", c(4, 1, 0, 6, 4)),
+                             list(name = "POLYLINE5", c(5, 2, 3, 7, 5)),
+                             list(name = "POLYLINE6", c(4, 6, 7, 5, 4)),
+                             list(name = "POLYLINE7", c(1, 0, 3, 2, 1))),
+
+        gml_surfaces = list(list(name = "SURFACE1", c(0, 2, 1), c(2, 0, 3)),
+                            list(name = "SURFACE2", c(4, 2, 1), c(2, 4, 5)),
+                            list(name = "SURFACE3", c(6, 3, 0), c(3, 6, 7)),
+                            list(name = "SURFACE4", c(4, 0, 1), c(0, 4, 6)),
+                            list(name = "SURFACE5", c(5, 3, 2), c(3, 5, 7)),
+                            list(name = "SURFACE6", c(4, 7, 6), c(7, 4, 5)),
+                            list(name = "SURFACE7", c(1, 3, 0), c(3, 1, 2)))
+    ), class = "r2ogs6_gml_prefab")
+}
+
+
+#'Constructor based on HydroMechanics/Verification/hm1_2Dsquare.gml
+new_r2ogs6_gml_hm1_2Dsquare <- function() {
+    structure(list(
+        gml_geometry_name = "geometry",
+        gml_points = tibble::tibble(x = c(0, 1, 1, 0, 0, 0, 1, 1),
+                                    y = c(0, 0, 1, 1, 1, 0, 0, 1),
+                                    z = c(0, 0, 0, 0, 0.10000000000000001, 0.10000000000000001,
+                                          0.10000000000000001, 0.10000000000000001),
+                                    name = c("POINT0", "POINT1", "POINT2", "POINT3",
+                                             "POINT6", "POINT7", "POINT17", "POINT18")),
+
+        gml_polylines = list(list(name = "POLYLINE1", c(0, 1, 2, 3, 0)),
+                             list(name = "POLYLINE2", c(0, 3, 4, 5, 0)),
+                             list(name = "POLYLINE3", c(0, 1, 6, 5, 0)),
+                             list(name = "POLYLINE4", c(0, 1, 2, 3, 0)),
+                             list(name = "POLYLINE5", c(5, 6, 7, 4, 5))),
+
+        gml_surfaces = list(list(name = "SURFACE1", c(0, 2, 1), c(2, 0, 3)),
+                            list(name = "SURFACE2", c(0, 4, 3), c(4, 0, 5)),
+                            list(name = "SURFACE3", c(0, 6, 1), c(6, 0, 5)),
+                            list(name = "SURFACE4", c(0, 2, 1), c(2, 0, 3)),
+                            list(name = "SURFACE5", c(5, 7, 6), c(7, 5, 4)))
+    ), class = "r2ogs6_gml_prefab")
+}
+
+
+#'Constructor based on HydroMechanics/Verification/hm1_3Dcube.gml
+new_r2ogs6_gml_hm1_3Dcube <- function() {
+    structure(list(
+        gml_geometry_name = "geometry",
+        gml_points = tibble::tibble(x = c(0, 1, 1, 0, 1, 0, 0, 1),
+                                    y = c(0, 0, 1, 1, 0, 0, 1, 1),
+                                    z = c(0, 0, 0, 0, 1, 1, 1, 1),
+                                    name = c("POINT0", "POINT1", "POINT2", "POINT3",
+                                             "POINT6", "POINT7", "POINT15", "POINT14")),
+
+        gml_polylines = list(list(name = "POLYLINE1", c(0, 1, 2, 3, 0)),
+                             list(name = "POLYLINE2", c(0, 1, 4, 5, 0)),
+                             list(name = "POLYLINE3", c(0, 3, 6, 5, 0)),
+                             list(name = "POLYLINE4", c(5, 4, 7, 6, 5))),
+
+        gml_surfaces = list(list(name = "SURFACE1", c(0, 2, 1), c(2, 0, 3)),
+                            list(name = "SURFACE2", c(0, 4, 1), c(4, 0, 5)),
+                            list(name = "SURFACE3", c(0, 6, 3), c(6, 0, 5)),
+                            list(name = "SURFACE4", c(5, 7, 4), c(7, 5, 6)))
+    ), class = "r2ogs6_gml_prefab")
+}
+
+#'Constructor based on HydroMechanics/Verification/hm1_3Dgravity.gml
+new_r2ogs6_gml_hm1_3Dgravity <- function() {
+    structure(list(
+        gml_geometry_name = "geometry",
+        gml_points = tibble::tibble(x = c(0, 60, 60, 0, 0, 60, 60, 0),
+                                    y = c(0, 0, 30, 30, 0, 0, 30, 30),
+                                    z = c(30, 30, 30, 30, 0, 0, 0, 0),
+                                    name = c("POINT0", "POINT1", "POINT2", "POINT3",
+                                             "POINT4", "POINT5", "POINT6", "POINT7")),
+
+        gml_polylines = list(list(name = "POLYLINE1", c(0, 1, 2, 3, 0)),
+                             list(name = "POLYLINE2", c(4, 5, 6, 7, 4)),
+                             list(name = "POLYLINE3", c(4, 5, 1, 0, 4)),
+                             list(name = "POLYLINE4", c(7, 6, 2, 3, 7)),
+                             list(name = "POLYLINE5", c(4, 7, 3, 0, 4)),
+                             list(name = "POLYLINE6", c(5, 6, 2, 1, 5))),
+
+        gml_surfaces = list(list(name = "SURFACE1", c(0, 2, 1), c(2, 0, 3)),
+                            list(name = "SURFACE2", c(4, 6, 5), c(6, 4, 7)),
+                            list(name = "SURFACE3", c(4, 1, 5), c(1, 4, 0)),
+                            list(name = "SURFACE4", c(7, 2, 6), c(2, 7, 3)),
+                            list(name = "SURFACE5", c(4, 3, 7), c(3, 4, 0)),
+                            list(name = "SURFACE6", c(5, 2, 6), c(2, 5, 1)))
+    ), class = "r2ogs6_gml_prefab")
+}
+
+
+#.../hm2_1d1bt.gml (add the hm2-gml prefabs)
+
+
+#'Constructor based on ogs6_benchmarks/Elliptic/circle_radius_1/circle_1_axi.gml
+#'
+#'
+new_r2ogs6_gml_circle_1_axi <- function() {
+    structure(list(
+        gml_geometry_name = "geometry",
+        gml_points = tibble::tibble(x = c(0, 1),
+                                    y = c(0, 0),
+                                    z = c(0, 0),
+                                    name = c("inner", "outer"))
+    ), class = "r2ogs6_gml_prefab")
+}
+
diff --git a/R/read_input_from_file.R b/R/read_input_from_file.R
new file mode 100644
index 0000000000000000000000000000000000000000..fec13607487636e618055d259affb7ae876bffd9
--- /dev/null
+++ b/R/read_input_from_file.R
@@ -0,0 +1,26 @@
+#This script contains various functions to read user input from a file (WIP)
+
+
+read_gml_input_from_file <- function(point_file_name, polyline_file_name, surface_file_name) {
+    read_gml_points(point_file_name)
+}
+
+#'Reads in a .tsv file and parses it as a tibble.
+#'A column specifying the name of the points is optional, but the number of columns must be consistent.
+#'If one row has an entry for the point name, all rows need an entry for the point name.
+#'If only some points need to be named, put "" (double quotes) for the points that don't have a name.
+#'@param file_name The name of the .tsv file containing the point tibble
+# @example
+#
+#
+read_gml_points <- function(file_name) {
+    point_tibble <- readr::read_tsv(file_name, col_names = FALSE, quoted_na = FALSE, comment = "#")
+    if(length(point_tibble) == 4) {
+        names(point_tibble) <- c("x", "y", "z", "name")
+    }else if(length(point_tibble) == 3) {
+        names(point_tibble) <- c("x", "y", "z")
+    }else{
+        stop(paste("Invalid number of columns detected while reading file containing .gml points"), call. = FALSE)
+    }
+    return(point_tibble)
+}
diff --git a/R/user_input.R b/R/user_input.R
index 86c868945205b835820b0bb51e1df695d4c34d00..2acfb732eca53f089217751a101e2e1a60472fc2 100644
--- a/R/user_input.R
+++ b/R/user_input.R
@@ -8,13 +8,13 @@
 geo_name <- "cube_1x1x1_geometry"
 
 #Some points (will become second XML child under root)
-my_points <- tibble::tibble(x = c(0, 0, 0, 0, 1, 1, 1, 1),
+my_gml_points <- tibble::tibble(x = c(0, 0, 0, 0, 1, 1, 1, 1),
                             y = c(0, 0, 1, 1, 0, 0, 1, 1),
                             z = c(0, 1, 1, 0, 0, 1, 1, 0),
-                            is_origin = c(TRUE, rep(FALSE,7)))
+                            name = c("origin", rep("", 7)))
 
 #Some polylines (will become third XML child under root)
-my_polylines <- list(list(name = "front_left", c(0, 1)),
+my_gml_polylines <- list(list(name = "front_left", c(0, 1)),
                      list(name = "front_right", c(4, 5)),
                      list(name = "front_bottom", c(0, 4)),
                      list(name = "front_top", c(1, 5)),
@@ -28,19 +28,38 @@ my_polylines <- list(list(name = "front_left", c(0, 1)),
                      list(name = "back_top", c(2, 6)))
 
 #Some surfaces (will become fourth XML child under root)
-my_surfaces <- list(list(name = "left", c(0, 1, 2), c(0, 3, 2)),
+my_gml_surfaces <- list(list(name = "left", c(0, 1, 2), c(0, 3, 2)),
                     list(name = "right", c(4, 6, 5), c(4, 6, 7)),
                     list(name = "top", c(1, 2, 5), c(5, 2, 6)),
                     list(name = "bottom", c(0, 3, 4), c(4, 3, 7)),
                     list(name = "front", c(0, 1, 4), c(4, 1, 5)),
                     list(name = "back", c(2, 3, 6), c(6, 3, 7)))
 
-#Turn those into XML structure
-my_gml <- data_to_xml(geo_name, my_points, my_polylines, my_surfaces)
+
+#============================== Data for VTU file ================================
+
+
+my_vtk_specification <- list(type="UnstructuredGrid",
+                             version="0.1",
+                             byte_order="LittleEndian",
+                             header_type="UInt32",
+                             compressor="vtkZLibDataCompressor")
+
+
+
+my_vtk_pieces <- list()
+
+
+
+
+#============================== Execution ================================
+
+#Turn GML data into XML structure
+my_gml_xml <- gml_data_to_xml(geo_name, my_gml_points, my_gml_polylines, my_gml_surfaces)
 
 
 #Let's export that as a new .gml XML file
-export_gml_to_file(my_gml, "my_experimental_cube.gml")
+export_xml_to_file(my_gml_xml, "my_experimental_cube.gml")
 
 #Now the .gml XML file should be in the r2ogs6 folder (not a perfect location but ok for a little test)
 
diff --git a/R/validate_gml.R b/R/validate_gml.R
new file mode 100644
index 0000000000000000000000000000000000000000..33111c96de89496e2c392d797630d05b15287137
--- /dev/null
+++ b/R/validate_gml.R
@@ -0,0 +1,66 @@
+#This script contains various functions to verify data for a .gml file (WIP)
+
+# @export
+validate_gml_data <- function(gml_geometry_name, gml_points = NA, gml_polylines = NA, gml_surfaces = NA) {
+
+    validate_points(gml_points)
+    validate_polylines(gml_polylines)
+    validate_surfaces(gml_surfaces)
+}
+
+#'Checks if the input is a tibble, if this tibble has the right number of elements,
+#' if those elements are named correctly, if the lists in the tibble are of the same length and
+#' if there is any overlapping points or duplicate point names
+#' @param gml_points A tibble with 3 vectors named 'x', 'y' and 'z' (and an optional 'name' vector)
+validate_points <- function(gml_points) {
+
+    if(!(inherits(gml_points, "tbl_df") && inherits(gml_points, "tbl") &&
+         inherits(gml_points, "data.frame"))){
+        stop(paste(gml_points, " is not of class 'tbl_df' "), call. = FALSE)
+    }
+
+    names <- names(gml_points)
+
+    if (!((length(gml_points) == 4 && names[[1]] == "x" && names[[2]] == "y" &&
+           names[[3]] == "z" && names[[4]] == "name") ||
+          (length(gml_points) == 3 && names[[1]] == "x" && names[[2]] == "y" &&
+           names[[3]] == "z"))){
+        stop(paste(gml_points, " column names do not fit to 'x, y, z, (name)' "), call. = FALSE)
+    }
+
+    has_names <- (length(gml_points) == 4)
+
+    #Find overlapping points and duplicate names
+    for(i in 1:(length(gml_points[[1]])-1)){
+        for(j in (i+1):length(gml_points[[1]])){
+            if(gml_points[[1]][[i]] == gml_points[[1]][[j]] &&
+               gml_points[[2]][[i]] == gml_points[[2]][[j]] &&
+               gml_points[[3]][[i]] == gml_points[[3]][[j]]){
+                stop("Overlapping .gml points (with the same coordinates) detected", call. = FALSE)
+            }
+
+            if(has_names){
+                if(gml_points[[4]][[i]] == gml_points[[4]][[j]] &&
+                   gml_points[[4]][[i]] != ""){
+                    warning("Duplicate .gml point names detected", call. = FALSE)
+                }
+            }
+        }
+    }
+
+    return(invisible(gml_points))
+}
+
+
+validate_polylines <- function(gml_polylines) {
+
+    if(!(inherits(gml_polylines, "list"))){
+        stop(paste(gml_polylines, " is not of class 'list' "), call. = FALSE)
+    }
+
+}
+
+
+validate_surfaces <- function(gml_surfaces) {
+
+}
\ No newline at end of file
diff --git a/inst/extdata/gml_points_with_names.tsv b/inst/extdata/gml_points_with_names.tsv
new file mode 100644
index 0000000000000000000000000000000000000000..b226882853265b510297708b8112ffa61768cd2b
--- /dev/null
+++ b/inst/extdata/gml_points_with_names.tsv
@@ -0,0 +1,8 @@
+0	0	0	origin
+0	0	1	""
+0	1	1	""
+0	1	0	""
+1	0	0	""
+1	0	1	""
+1	1	1	""
+1	1	0	""
\ No newline at end of file
diff --git a/inst/extdata/gml_points_without_names.tsv b/inst/extdata/gml_points_without_names.tsv
new file mode 100644
index 0000000000000000000000000000000000000000..53f58ddad60060f0ce584bfafdf311a539337866
--- /dev/null
+++ b/inst/extdata/gml_points_without_names.tsv
@@ -0,0 +1,8 @@
+0	0	0
+0	0	1
+0	1	1
+0	1	0
+1	0	0
+1	0	1
+1	1	1
+1	1	0
\ No newline at end of file
diff --git a/man/GMLPoints.Rd b/man/GMLPoints.Rd
index 38b26db0952884aa28fd0ae893c00b226e3b4153..2b84e1ac081a24e653d8ca5b81c7f40acc4300df 100644
--- a/man/GMLPoints.Rd
+++ b/man/GMLPoints.Rd
@@ -15,7 +15,6 @@ Class describing a set of points
 \item \href{#method-validate_ids}{\code{GMLPoints$validate_ids()}}
 \item \href{#method-validate_dim}{\code{GMLPoints$validate_dim()}}
 \item \href{#method-validate_coordinates}{\code{GMLPoints$validate_coordinates()}}
-\item \href{#method-validate_origin}{\code{GMLPoints$validate_origin()}}
 \item \href{#method-as_list}{\code{GMLPoints$as_list()}}
 \item \href{#method-as_node}{\code{GMLPoints$as_node()}}
 \item \href{#method-print}{\code{GMLPoints$print()}}
@@ -57,15 +56,6 @@ Class describing a set of points
 \if{html}{\out{<div class="r">}}\preformatted{GMLPoints$validate_coordinates(gml_points)}\if{html}{\out{</div>}}
 }
 
-}
-\if{html}{\out{<hr>}}
-\if{html}{\out{<a id="method-validate_origin"></a>}}
-\if{latex}{\out{\hypertarget{method-validate_origin}{}}}
-\subsection{Method \code{validate_origin()}}{
-\subsection{Usage}{
-\if{html}{\out{<div class="r">}}\preformatted{GMLPoints$validate_origin(gml_points)}\if{html}{\out{</div>}}
-}
-
 }
 \if{html}{\out{<hr>}}
 \if{html}{\out{<a id="method-as_list"></a>}}
diff --git a/man/export_gml_to_file.Rd b/man/export_xml_to_file.Rd
similarity index 59%
rename from man/export_gml_to_file.Rd
rename to man/export_xml_to_file.Rd
index 1170afc902968ad7bfa007452af3d07e214d5faa..8aa952faca91bfea7617cf221aa91226407bd1db 100644
--- a/man/export_gml_to_file.Rd
+++ b/man/export_xml_to_file.Rd
@@ -1,13 +1,13 @@
 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/create_gml.R
-\name{export_gml_to_file}
-\alias{export_gml_to_file}
+\name{export_xml_to_file}
+\alias{export_xml_to_file}
 \title{Export function}
 \usage{
-export_gml_to_file(gml_data, file_name)
+export_xml_to_file(xml_data, file_name)
 }
 \arguments{
-\item{gml_data}{The .gml data (already in XML friendly format)}
+\item{xml_data}{The .gml data (already in XML friendly format)}
 
 \item{file_name}{The name of the .gml file to be written}
 }
diff --git a/man/data_to_xml.Rd b/man/gml_data_to_xml.Rd
similarity index 82%
rename from man/data_to_xml.Rd
rename to man/gml_data_to_xml.Rd
index 2edaee904a70c109292cd3f4e2a61555b059db8c..919404a797f323c77c899b08c87ef6edfebbaf9c 100644
--- a/man/data_to_xml.Rd
+++ b/man/gml_data_to_xml.Rd
@@ -1,10 +1,10 @@
 % Generated by roxygen2: do not edit by hand
 % Please edit documentation in R/create_gml.R
-\name{data_to_xml}
-\alias{data_to_xml}
+\name{gml_data_to_xml}
+\alias{gml_data_to_xml}
 \title{Wrapper function to create a XML document based on the user input data}
 \usage{
-data_to_xml(geo_name, points_tibble, polylines_list, surfaces_list)
+gml_data_to_xml(geo_name, points_tibble, polylines_list, surfaces_list)
 }
 \arguments{
 \item{geo_name}{The name of the geometry specified by the user}
diff --git a/man/new_r2ogs6_gml_bar.Rd b/man/new_r2ogs6_gml_bar.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..73bf557e5dae263025597418371996ce5966cce9
--- /dev/null
+++ b/man/new_r2ogs6_gml_bar.Rd
@@ -0,0 +1,11 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gml_prefabs.R
+\name{new_r2ogs6_gml_bar}
+\alias{new_r2ogs6_gml_bar}
+\title{Constructor based on HydroMechanics/StaggeredScheme/InjectionProduction1D/bar.gml}
+\usage{
+new_r2ogs6_gml_bar()
+}
+\description{
+Constructor based on HydroMechanics/StaggeredScheme/InjectionProduction1D/bar.gml
+}
diff --git a/man/new_r2ogs6_gml_circle_1_axi.Rd b/man/new_r2ogs6_gml_circle_1_axi.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..14a8386869b320d366302a76be1942accec9c61c
--- /dev/null
+++ b/man/new_r2ogs6_gml_circle_1_axi.Rd
@@ -0,0 +1,11 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gml_prefabs.R
+\name{new_r2ogs6_gml_circle_1_axi}
+\alias{new_r2ogs6_gml_circle_1_axi}
+\title{Constructor based on ogs6_benchmarks/Elliptic/circle_radius_1/circle_1_axi.gml}
+\usage{
+new_r2ogs6_gml_circle_1_axi()
+}
+\description{
+Constructor based on ogs6_benchmarks/Elliptic/circle_radius_1/circle_1_axi.gml
+}
diff --git a/man/new_r2ogs6_gml_cube_1x1x1.Rd b/man/new_r2ogs6_gml_cube_1x1x1.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..903f5d5fa4f150570e12414fa2f5561b0daf17cb
--- /dev/null
+++ b/man/new_r2ogs6_gml_cube_1x1x1.Rd
@@ -0,0 +1,11 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gml_prefabs.R
+\name{new_r2ogs6_gml_cube_1x1x1}
+\alias{new_r2ogs6_gml_cube_1x1x1}
+\title{Constructor based on HydroMechanics/IdealGas/flow_free_expansion/cube_1x1x1.gml}
+\usage{
+new_r2ogs6_gml_cube_1x1x1()
+}
+\description{
+Constructor based on HydroMechanics/IdealGas/flow_free_expansion/cube_1x1x1.gml
+}
diff --git a/man/new_r2ogs6_gml_hm1_1Dbeam.Rd b/man/new_r2ogs6_gml_hm1_1Dbeam.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..0baceff5fbca4b972bf2af35cb777b5c048e85c5
--- /dev/null
+++ b/man/new_r2ogs6_gml_hm1_1Dbeam.Rd
@@ -0,0 +1,11 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gml_prefabs.R
+\name{new_r2ogs6_gml_hm1_1Dbeam}
+\alias{new_r2ogs6_gml_hm1_1Dbeam}
+\title{Constructor based on HydroMechanics/Verification/hm1_1Dbeam.gml}
+\usage{
+new_r2ogs6_gml_hm1_1Dbeam()
+}
+\description{
+Constructor based on HydroMechanics/Verification/hm1_1Dbeam.gml
+}
diff --git a/man/new_r2ogs6_gml_hm1_2Dsquare.Rd b/man/new_r2ogs6_gml_hm1_2Dsquare.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..c0dfbcb4af479f653822781696dabda6d838945a
--- /dev/null
+++ b/man/new_r2ogs6_gml_hm1_2Dsquare.Rd
@@ -0,0 +1,11 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gml_prefabs.R
+\name{new_r2ogs6_gml_hm1_2Dsquare}
+\alias{new_r2ogs6_gml_hm1_2Dsquare}
+\title{Constructor based on HydroMechanics/Verification/hm1_2Dsquare.gml}
+\usage{
+new_r2ogs6_gml_hm1_2Dsquare()
+}
+\description{
+Constructor based on HydroMechanics/Verification/hm1_2Dsquare.gml
+}
diff --git a/man/new_r2ogs6_gml_hm1_3Dcube.Rd b/man/new_r2ogs6_gml_hm1_3Dcube.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..0fad335c88910c6b7e79362c38433bee8d994003
--- /dev/null
+++ b/man/new_r2ogs6_gml_hm1_3Dcube.Rd
@@ -0,0 +1,11 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gml_prefabs.R
+\name{new_r2ogs6_gml_hm1_3Dcube}
+\alias{new_r2ogs6_gml_hm1_3Dcube}
+\title{Constructor based on HydroMechanics/Verification/hm1_3Dcube.gml}
+\usage{
+new_r2ogs6_gml_hm1_3Dcube()
+}
+\description{
+Constructor based on HydroMechanics/Verification/hm1_3Dcube.gml
+}
diff --git a/man/new_r2ogs6_gml_hm1_3Dgravity.Rd b/man/new_r2ogs6_gml_hm1_3Dgravity.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..b227f4d17ea56e5447e14d66bffdca805360b07f
--- /dev/null
+++ b/man/new_r2ogs6_gml_hm1_3Dgravity.Rd
@@ -0,0 +1,11 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gml_prefabs.R
+\name{new_r2ogs6_gml_hm1_3Dgravity}
+\alias{new_r2ogs6_gml_hm1_3Dgravity}
+\title{Constructor based on HydroMechanics/Verification/hm1_3Dgravity.gml}
+\usage{
+new_r2ogs6_gml_hm1_3Dgravity()
+}
+\description{
+Constructor based on HydroMechanics/Verification/hm1_3Dgravity.gml
+}
diff --git a/man/new_r2ogs6_gml_quad_1x10.Rd b/man/new_r2ogs6_gml_quad_1x10.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..8754e2a1801255f076dac1209d9ef02033ec7d42
--- /dev/null
+++ b/man/new_r2ogs6_gml_quad_1x10.Rd
@@ -0,0 +1,11 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gml_prefabs.R
+\name{new_r2ogs6_gml_quad_1x10}
+\alias{new_r2ogs6_gml_quad_1x10}
+\title{Constructor based on HydroMechanics/IdealGas/flow_pressure_boundary/quad_1x10.gml}
+\usage{
+new_r2ogs6_gml_quad_1x10()
+}
+\description{
+Constructor based on HydroMechanics/IdealGas/flow_pressure_boundary/quad_1x10.gml
+}
diff --git a/man/new_r2ogs6_gml_square_1x1.Rd b/man/new_r2ogs6_gml_square_1x1.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..9c5bdd2aeda19bf6a2fa94d7fbd6797f99295cb8
--- /dev/null
+++ b/man/new_r2ogs6_gml_square_1x1.Rd
@@ -0,0 +1,11 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gml_prefabs.R
+\name{new_r2ogs6_gml_square_1x1}
+\alias{new_r2ogs6_gml_square_1x1}
+\title{Constructor based on HydroMechanics/IdealGas/flow_no_strain/square_1x1.gml}
+\usage{
+new_r2ogs6_gml_square_1x1()
+}
+\description{
+Constructor based on HydroMechanics/IdealGas/flow_no_strain/square_1x1.gml
+}
diff --git a/man/points_to_xml.Rd b/man/points_to_xml.Rd
index 403b7ad54f749411ceba313ecc6b1fa4d8924e30..3028ad0f44121cb6b084a15b4e556fb415e91467 100644
--- a/man/points_to_xml.Rd
+++ b/man/points_to_xml.Rd
@@ -16,6 +16,6 @@ An XML node containing the points
 Turns a tibble of points into an XML node
 }
 \examples{
-my_tibble <- tibble::tibble(x = c(0, 0), y = c(0, 0), z = c(0, 1),is_origin = c(TRUE, FALSE))
+my_tibble <- tibble::tibble(x = c(0, 0), y = c(0, 0), z = c(0, 1), name = c("origin", ""))
 point_node <- points_to_xml(my_tibble)
 }
diff --git a/man/read_gml_points.Rd b/man/read_gml_points.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..56da0a8dcaf9dd334a7658a3b39a97ec9307567a
--- /dev/null
+++ b/man/read_gml_points.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/read_input_from_file.R
+\name{read_gml_points}
+\alias{read_gml_points}
+\title{Reads in a .tsv file and parses it as a tibble.
+A column specifying the name of the points is optional, but the number of columns must be consistent.
+If one row has an entry for the point name, all rows need an entry for the point name.
+If only some points need to be named, put "" (double quotes) for the points that don't have a name.}
+\usage{
+read_gml_points(file_name)
+}
+\arguments{
+\item{file_name}{The name of the .tsv file containing the point tibble}
+}
+\description{
+Reads in a .tsv file and parses it as a tibble.
+A column specifying the name of the points is optional, but the number of columns must be consistent.
+If one row has an entry for the point name, all rows need an entry for the point name.
+If only some points need to be named, put "" (double quotes) for the points that don't have a name.
+}
diff --git a/man/validate_points.Rd b/man/validate_points.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..f7819198cb1510d6ab0511840598843d86c6dc0c
--- /dev/null
+++ b/man/validate_points.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/validate_gml.R
+\name{validate_points}
+\alias{validate_points}
+\title{Checks if the input is a tibble, if this tibble has the right number of elements,
+if those elements are named correctly, if the lists in the tibble are of the same length and
+if there is any overlapping points or duplicate point names}
+\usage{
+validate_points(gml_points)
+}
+\arguments{
+\item{gml_points}{A tibble with 3 vectors named 'x', 'y' and 'z' (and an optional 'name' vector)}
+}
+\description{
+Checks if the input is a tibble, if this tibble has the right number of elements,
+if those elements are named correctly, if the lists in the tibble are of the same length and
+if there is any overlapping points or duplicate point names
+}
diff --git a/my_experimental_cube.gml b/my_experimental_cube.gml
new file mode 100644
index 0000000000000000000000000000000000000000..939bdf775322bb0eb3881d4facaeb7d046d2f7d2
--- /dev/null
+++ b/my_experimental_cube.gml
@@ -0,0 +1,90 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysGLI xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:ogs="http://www.opengeosys.org">
+  <name>cube_1x1x1_geometry</name>
+  <points>
+    <point id="0" x="0" y="0" z="0" name="origin"/>
+    <point id="1" x="0" y="0" z="1"/>
+    <point id="2" x="0" y="1" z="1"/>
+    <point id="3" x="0" y="1" z="0"/>
+    <point id="4" x="1" y="0" z="0"/>
+    <point id="5" x="1" y="0" z="1"/>
+    <point id="6" x="1" y="1" z="1"/>
+    <point id="7" x="1" y="1" z="0"/>
+  </points>
+  <polylines>
+    <polyline id="0" name="front_left">
+      <pnt>0</pnt>
+      <pnt>1</pnt>
+    </polyline>
+    <polyline id="1" name="front_right">
+      <pnt>4</pnt>
+      <pnt>5</pnt>
+    </polyline>
+    <polyline id="2" name="front_bottom">
+      <pnt>0</pnt>
+      <pnt>4</pnt>
+    </polyline>
+    <polyline id="3" name="front_top">
+      <pnt>1</pnt>
+      <pnt>5</pnt>
+    </polyline>
+    <polyline id="4" name="bottom_left">
+      <pnt>0</pnt>
+      <pnt>3</pnt>
+    </polyline>
+    <polyline id="5" name="bottom_right">
+      <pnt>4</pnt>
+      <pnt>7</pnt>
+    </polyline>
+    <polyline id="6" name="top_left">
+      <pnt>1</pnt>
+      <pnt>2</pnt>
+    </polyline>
+    <polyline id="7" name="top_right">
+      <pnt>5</pnt>
+      <pnt>6</pnt>
+    </polyline>
+    <polyline id="8" name="back_left">
+      <pnt>2</pnt>
+      <pnt>3</pnt>
+    </polyline>
+    <polyline id="9" name="back_right">
+      <pnt>6</pnt>
+      <pnt>7</pnt>
+    </polyline>
+    <polyline id="10" name="back_bottom">
+      <pnt>3</pnt>
+      <pnt>7</pnt>
+    </polyline>
+    <polyline id="11" name="back_top">
+      <pnt>2</pnt>
+      <pnt>6</pnt>
+    </polyline>
+  </polylines>
+  <surfaces>
+    <surface id="0" name="left">
+      <element p1="0" p2="1" p3="2"/>
+      <element p1="0" p2="3" p3="2"/>
+    </surface>
+    <surface id="1" name="right">
+      <element p1="4" p2="6" p3="5"/>
+      <element p1="4" p2="6" p3="7"/>
+    </surface>
+    <surface id="2" name="top">
+      <element p1="1" p2="2" p3="5"/>
+      <element p1="5" p2="2" p3="6"/>
+    </surface>
+    <surface id="3" name="bottom">
+      <element p1="0" p2="3" p3="4"/>
+      <element p1="4" p2="3" p3="7"/>
+    </surface>
+    <surface id="4" name="front">
+      <element p1="0" p2="1" p3="4"/>
+      <element p1="4" p2="1" p3="5"/>
+    </surface>
+    <surface id="5" name="back">
+      <element p1="2" p2="3" p3="6"/>
+      <element p1="6" p2="3" p3="7"/>
+    </surface>
+  </surfaces>
+</OpenGeoSysGLI>
diff --git a/tests/testthat/test-read_input_from_file.R b/tests/testthat/test-read_input_from_file.R
new file mode 100644
index 0000000000000000000000000000000000000000..2abeecada4d17c21aae58372c7fc2079ad27dee9
--- /dev/null
+++ b/tests/testthat/test-read_input_from_file.R
@@ -0,0 +1,12 @@
+
+test_that("read_gml_points correctly reads point tibble from file", {
+
+    point_tibble_1 <- read_gml_points(system.file("extdata", "gml_points_with_names.tsv", package = "r2ogs6"))
+    expect_equal(length(point_tibble_1), 4)
+    expect_setequal(names(point_tibble_1), c("x", "y", "z", "name"))
+
+    point_tibble_2 <- read_gml_points(system.file("extdata", "gml_points_without_names.tsv", package = "r2ogs6"))
+    expect_equal(length(point_tibble_2), 3)
+    expect_setequal(names(point_tibble_2), c("x", "y", "z"))
+})
+
diff --git a/tests/testthat/test-validate_gml.R b/tests/testthat/test-validate_gml.R
new file mode 100644
index 0000000000000000000000000000000000000000..2ad65b7565cab5dc0c4bb9450eb2bfb67ae82ec9
--- /dev/null
+++ b/tests/testthat/test-validate_gml.R
@@ -0,0 +1,24 @@
+test_that("validate_points function checks ...", {
+
+    point_list <- list(x = c(0, 0), y = c(1, 1), z = c(0, 1))
+
+    #Check class (should expect a tibble, not a list)
+    expect_error(validate_points(point_list))
+
+    point_tibble_inv_0 <- tibble::tibble(a = c(0, 0), b = c(1, 1), c = c(0, 1))
+    point_tibble_inv_1 <- tibble::tibble(x = c(0, 0), y = c(1, 1), z = c(0))
+    point_tibble_inv_2 <- tibble::tibble(x = c(0, 0), y = c(1, 1), z = c(0, 1), name = c("oh", "oh"))
+
+    #Check column names
+    expect_error(validate_points(point_tibble_inv_0))
+
+    #Check if duplicate points and point names are detected
+    expect_error(validate_points(point_tibble_inv_1))
+    expect_warning(validate_points(point_tibble_inv_2))
+
+    point_tibble_val <- tibble::tibble(x = c(0, 0), y = c(1, 1), z = c(0, 1), name = c("this", "works"))
+    point_tibble_val_2 <- tibble::tibble(x = c(0, 0), y = c(1, 1), z = c(0, 1))
+
+    expect_invisible(validate_points(point_tibble_val))
+    expect_invisible(validate_points(point_tibble_val_2))
+})