diff --git a/docs/examples/howto_meshlib/plot_meshlib_vtu_input.py b/docs/examples/howto_meshlib/plot_meshlib_vtu_input.py
index e02669b19fb4e4f0773ff4981b30a40b6340747d..a2956693ca917b3d2cbeee87f221bcf561c9bb91 100644
--- a/docs/examples/howto_meshlib/plot_meshlib_vtu_input.py
+++ b/docs/examples/howto_meshlib/plot_meshlib_vtu_input.py
@@ -13,7 +13,7 @@ from pathlib import Path  # To get example vtu files
 import numpy as np  # For visualization only
 
 import ogstools.meshplotlib as mpl  # For visualization only
-from ogstools.definitions import ROOT_DIR  # To get example vtu files
+from ogstools.definitions import TESTS_DIR  # To get example vtu files
 from ogstools.meshlib.boundary import Layer
 from ogstools.meshlib.boundary_set import LayerSet
 from ogstools.meshlib.boundary_subset import Surface
@@ -32,7 +32,7 @@ mpl.setup.aspect_limits = [0.2, 5.0]
 # The loaded surfaces are defined within VTU files and adhere to properties such as non-intersecting boundaries with consistent x and y bounds. Alternatively, surfaces can also be created using PyVista with the same properties.
 
 
-surface_dir = ROOT_DIR / "meshlib/tests/data/mesh1/surface_data"
+surface_dir = TESTS_DIR / "data" / "meshlib" / "mesh1" / "surface_data"
 
 surface1 = Surface(Path(surface_dir / "00_KB.vtu"), material_id=0)
 surface2 = Surface(Path(surface_dir / "01_q.vtu"), material_id=5)
diff --git a/ogstools/definitions.py b/ogstools/definitions.py
index 6b2605693db1ba9ba92daafe849f61d7f3207c9d..d79a75c893573c6a70f6c1ac77f4f17bbac3b770 100644
--- a/ogstools/definitions.py
+++ b/ogstools/definitions.py
@@ -1,3 +1,4 @@
 import pathlib
 
 ROOT_DIR = pathlib.Path(__file__).parent.resolve()
+TESTS_DIR = pathlib.Path(__file__).parent.parent.resolve() / "tests"
diff --git a/ogstools/meshlib/boundary_subset.py b/ogstools/meshlib/boundary_subset.py
index e8c5c01cc874c985a0aa6c9842126eb8bf1b198d..6a7f7b6a2c5d2bdd9569c227d61bdd7bcd596826 100644
--- a/ogstools/meshlib/boundary_subset.py
+++ b/ogstools/meshlib/boundary_subset.py
@@ -30,15 +30,18 @@ class Surface:
         if isinstance(input, Path):
             self.filename = input
             if self.filename.exists() is False:
-                print(self.filename, "does not exist.")
-                raise ValueError
+                msg = f"{self.filename} does not exist."
+                raise ValueError(msg)
             self.mesh = pv.get_reader(self.filename).read()
         elif isinstance(input, pv.DataSet):
             self.mesh = input
             self.filename = Path(tempfile.mkstemp(".vtu", "surface")[1])
             pv.save_meshio(self.filename, self.mesh, file_format="vtu")
         else:
-            raise ValueError
+            msg = "{} given, must be either Path or pyvista Dataset.".format(
+                self.filename
+            )
+            raise ValueError(msg)
 
         self.mesh.cell_data["MaterialIDs"] = (
             np.ones(self.mesh.n_cells) * self.material_id
diff --git a/ogstools/meshlib/tests/__init__.py b/ogstools/meshlib/tests/__init__.py
deleted file mode 100644
index fdae684c33504f57a79b5467824b6948d4361243..0000000000000000000000000000000000000000
--- a/ogstools/meshlib/tests/__init__.py
+++ /dev/null
@@ -1,3 +0,0 @@
-from .test_utils import MeshPath, dataframe_from_csv
-
-__all__ = ["MeshPath", "dataframe_from_csv"]
diff --git a/ogstools/meshlib/tests/data/compose_geomodel/layersets.csv b/tests/data/meshlib/compose_geomodel/layersets.csv
similarity index 100%
rename from ogstools/meshlib/tests/data/compose_geomodel/layersets.csv
rename to tests/data/meshlib/compose_geomodel/layersets.csv
diff --git a/ogstools/meshlib/tests/data/compose_geomodel/materialset.csv b/tests/data/meshlib/compose_geomodel/materialset.csv
similarity index 100%
rename from ogstools/meshlib/tests/data/compose_geomodel/materialset.csv
rename to tests/data/meshlib/compose_geomodel/materialset.csv
diff --git a/ogstools/meshlib/tests/data/mesh1/set1.csv b/tests/data/meshlib/mesh1/set1.csv
similarity index 100%
rename from ogstools/meshlib/tests/data/mesh1/set1.csv
rename to tests/data/meshlib/mesh1/set1.csv
diff --git a/ogstools/meshlib/tests/data/mesh1/surface_data/00_KB.vtu b/tests/data/meshlib/mesh1/surface_data/00_KB.vtu
similarity index 100%
rename from ogstools/meshlib/tests/data/mesh1/surface_data/00_KB.vtu
rename to tests/data/meshlib/mesh1/surface_data/00_KB.vtu
diff --git a/ogstools/meshlib/tests/data/mesh1/surface_data/01_q.vtu b/tests/data/meshlib/mesh1/surface_data/01_q.vtu
similarity index 100%
rename from ogstools/meshlib/tests/data/mesh1/surface_data/01_q.vtu
rename to tests/data/meshlib/mesh1/surface_data/01_q.vtu
diff --git a/ogstools/meshlib/tests/data/mesh1/surface_data/02_krl.vtu b/tests/data/meshlib/mesh1/surface_data/02_krl.vtu
similarity index 100%
rename from ogstools/meshlib/tests/data/mesh1/surface_data/02_krl.vtu
rename to tests/data/meshlib/mesh1/surface_data/02_krl.vtu
diff --git a/ogstools/meshlib/tests/data/mesh1/surface_data/03_S3.vtu b/tests/data/meshlib/mesh1/surface_data/03_S3.vtu
similarity index 100%
rename from ogstools/meshlib/tests/data/mesh1/surface_data/03_S3.vtu
rename to tests/data/meshlib/mesh1/surface_data/03_S3.vtu
diff --git a/ogstools/meshlib/tests/data/mesh1/surface_data/04_krp.vtu b/tests/data/meshlib/mesh1/surface_data/04_krp.vtu
similarity index 100%
rename from ogstools/meshlib/tests/data/mesh1/surface_data/04_krp.vtu
rename to tests/data/meshlib/mesh1/surface_data/04_krp.vtu
diff --git a/tests/meshlib/__init__.py b/tests/meshlib/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..09d11fc4e887cb4359ceb6fb21d8b11f1b3fd620
--- /dev/null
+++ b/tests/meshlib/__init__.py
@@ -0,0 +1 @@
+__all__ = ["dataframe_from_csv"]
diff --git a/ogstools/meshlib/tests/test_demo.py b/tests/meshlib/test_demo.py
similarity index 88%
rename from ogstools/meshlib/tests/test_demo.py
rename to tests/meshlib/test_demo.py
index c77ec558e5a6491f71c1baef7af4ebc9f7f85b11..6de2d9cb5b298766a5658337fc0bf963f13b20c8 100644
--- a/ogstools/meshlib/tests/test_demo.py
+++ b/tests/meshlib/test_demo.py
@@ -1,6 +1,7 @@
 import unittest
 from collections import Counter
 
+from ogstools.definitions import TESTS_DIR
 from ogstools.meshlib.boundary import Layer
 from ogstools.meshlib.boundary_set import LayerSet
 from ogstools.meshlib.boundary_subset import Gaussian2D, Surface
@@ -10,26 +11,27 @@ from ogstools.meshlib.region import (
     to_region_tetraeder,
     to_region_voxel,
 )
-from ogstools.meshlib.tests import MeshPath
+
+meshpath = TESTS_DIR / "data" / "meshlib"
 
 
 class DemoTest(unittest.TestCase):
     def test_allcompare(self):
         # To define a mesh with 3 layers from example input, create 4 surfaces (3 bottom surface + 1 top surface)
         surface1 = Surface(
-            MeshPath("data/mesh1/surface_data/00_KB.vtu"),
+            meshpath / "mesh1/surface_data/00_KB.vtu",
             material_id=0,
         )
         surface2 = Surface(
-            MeshPath("data/mesh1/surface_data/01_q.vtu"),
+            meshpath / "mesh1/surface_data/01_q.vtu",
             material_id=5,
         )
         surface3 = Surface(
-            MeshPath("data/mesh1/surface_data/02_krl.vtu"),
+            meshpath / "mesh1/surface_data/02_krl.vtu",
             material_id=2,
         )
         surface4 = Surface(
-            MeshPath("data/mesh1/surface_data/03_S3.vtu"),
+            meshpath / "mesh1/surface_data/03_S3.vtu",
             material_id=3,
         )
         layer1 = Layer(top=surface1, bottom=surface2, num_subdivisions=2)
diff --git a/ogstools/meshlib/tests/test_layer.py b/tests/meshlib/test_layer.py
similarity index 94%
rename from ogstools/meshlib/tests/test_layer.py
rename to tests/meshlib/test_layer.py
index a551b8eb5663c201469a452964977dcdd72ad8ad..4fd1d364f0f66eece5d86d0b957c9730ea679802 100644
--- a/ogstools/meshlib/tests/test_layer.py
+++ b/tests/meshlib/test_layer.py
@@ -6,6 +6,7 @@ from pathlib import Path
 import numpy as np
 import pyvista as pv
 
+from ogstools.definitions import TESTS_DIR
 from ogstools.meshlib.boundary import Layer, LocationFrame, Raster
 from ogstools.meshlib.boundary_set import LayerSet
 from ogstools.meshlib.boundary_subset import Gaussian2D, Surface
@@ -15,7 +16,8 @@ from ogstools.meshlib.region import (
     to_region_tetraeder,
     to_region_voxel,
 )
-from ogstools.meshlib.tests import MeshPath
+
+meshpath = TESTS_DIR / "data" / "meshlib"
 
 
 class LayerTest(unittest.TestCase):
@@ -25,11 +27,11 @@ class LayerTest(unittest.TestCase):
         """
         layer1 = Layer(
             top=Surface(
-                MeshPath("data/mesh1/surface_data/00_KB.vtu"),
+                meshpath / "mesh1/surface_data/00_KB.vtu",
                 0,
             ),
             bottom=Surface(
-                MeshPath("data/mesh1/surface_data/01_q.vtu"),
+                meshpath / "mesh1/surface_data/01_q.vtu",
                 1,
             ),
             material_id=0,
diff --git a/ogstools/meshlib/tests/test_layer_set.py b/tests/meshlib/test_layer_set.py
similarity index 84%
rename from ogstools/meshlib/tests/test_layer_set.py
rename to tests/meshlib/test_layer_set.py
index bffed486c80dfccdee86933888db40937d71d2ee..1b06a40eba37ad0c90a1ef855c4171adc7e6295a 100644
--- a/ogstools/meshlib/tests/test_layer_set.py
+++ b/tests/meshlib/test_layer_set.py
@@ -1,16 +1,18 @@
 import unittest
 
+from ogstools.definitions import TESTS_DIR
 from ogstools.meshlib._utils import dataframe_from_csv
 from ogstools.meshlib.boundary import Layer
 from ogstools.meshlib.boundary_set import LayerSet
 from ogstools.meshlib.boundary_subset import Surface
-from ogstools.meshlib.tests import MeshPath
+
+meshpath = TESTS_DIR / "data" / "meshlib"
 
 
 class LayerSetTest(unittest.TestCase):
-    layerset = MeshPath("data/compose_geomodel/layersets.csv")
-    materialset = MeshPath("data/compose_geomodel/materialset.csv")
-    surfacedata = MeshPath("data/mesh1/surface_data/")
+    layerset = meshpath / "compose_geomodel/layersets.csv"
+    materialset = meshpath / "compose_geomodel/materialset.csv"
+    surfacedata = meshpath / "mesh1/surface_data/"
 
     def test_create_with_3_intermediate(self):
         mesh3_df = dataframe_from_csv(
@@ -36,15 +38,15 @@ class LayerSetTest(unittest.TestCase):
     def test_create_with_no_intermediate(self):
         resolution = 300
         surface1 = Surface(
-            MeshPath("data/mesh1/surface_data/00_KB.vtu"),
+            meshpath / "mesh1/surface_data/00_KB.vtu",
             0,
         )
         surface2 = Surface(
-            MeshPath("data/mesh1/surface_data/01_q.vtu"),
+            meshpath / "mesh1/surface_data/01_q.vtu",
             1,
         )
         surface3 = Surface(
-            MeshPath("data/mesh1/surface_data/02_krl.vtu"),
+            meshpath / "mesh1/surface_data/02_krl.vtu",
             1,
         )
 
diff --git a/ogstools/meshlib/tests/test_prism_mesh.py b/tests/meshlib/test_prism_mesh.py
similarity index 83%
rename from ogstools/meshlib/tests/test_prism_mesh.py
rename to tests/meshlib/test_prism_mesh.py
index 04c74b5063f7f4af4e9eee5d877519cf277da0dc..e6fab2b39f14a8e64ae89bd876b7d617b8f811e2 100644
--- a/ogstools/meshlib/tests/test_prism_mesh.py
+++ b/tests/meshlib/test_prism_mesh.py
@@ -1,15 +1,17 @@
 import unittest
 
+from ogstools.definitions import TESTS_DIR
 from ogstools.meshlib._utils import dataframe_from_csv
 from ogstools.meshlib.boundary_set import LayerSet
 from ogstools.meshlib.region import to_region_prism
-from ogstools.meshlib.tests import MeshPath
+
+meshpath = TESTS_DIR / "data" / "meshlib"
 
 
 class PrismMeshTest(unittest.TestCase):
-    layerset = MeshPath("data/compose_geomodel/layersets.csv")
-    materialset = MeshPath("data/compose_geomodel/materialset.csv")
-    surfacedata = MeshPath("data/mesh1/surface_data/")
+    layerset = meshpath / "compose_geomodel/layersets.csv"
+    materialset = meshpath / "compose_geomodel/materialset.csv"
+    surfacedata = meshpath / "mesh1/surface_data/"
 
     def test_Mesh_fineXY_coarseZ(self):
         mesh1_df = dataframe_from_csv(
diff --git a/ogstools/meshlib/tests/test_read_xdmf.py b/tests/meshlib/test_read_xdmf.py
similarity index 100%
rename from ogstools/meshlib/tests/test_read_xdmf.py
rename to tests/meshlib/test_read_xdmf.py
diff --git a/ogstools/meshlib/tests/test_region.py b/tests/meshlib/test_region.py
similarity index 65%
rename from ogstools/meshlib/tests/test_region.py
rename to tests/meshlib/test_region.py
index 560396cc8175e05f2a7ab2192dc11d9573255db9..3814c0e91b73fe7af6ab7693f21125541d960d6c 100644
--- a/ogstools/meshlib/tests/test_region.py
+++ b/tests/meshlib/test_region.py
@@ -1,16 +1,19 @@
 import unittest
 
+from ogstools.definitions import TESTS_DIR
+from ogstools.meshlib._utils import dataframe_from_csv
 from ogstools.meshlib.boundary_set import LayerSet
 from ogstools.meshlib.region import (
     to_region_simplified,
 )
-from ogstools.meshlib.tests import MeshPath, dataframe_from_csv
+
+meshpath = TESTS_DIR / "data" / "meshlib"
 
 
 class RegionTest(unittest.TestCase):
-    layerset = MeshPath("data/compose_geomodel/layersets.csv")
-    materialset = MeshPath("data/compose_geomodel/materialset.csv")
-    surfacedata = MeshPath("data/mesh1/surface_data/")
+    layerset = meshpath / "compose_geomodel/layersets.csv"
+    materialset = meshpath / "compose_geomodel/materialset.csv"
+    surfacedata = meshpath / "mesh1/surface_data/"
 
     def test_top_boundary(self):
         mesh_df_coarseZ = dataframe_from_csv(
diff --git a/ogstools/meshlib/tests/test_simplified_mesh.py b/tests/meshlib/test_simplified_mesh.py
similarity index 92%
rename from ogstools/meshlib/tests/test_simplified_mesh.py
rename to tests/meshlib/test_simplified_mesh.py
index a5d0476a75f0dbbff3cb67be51a2dd5cb91fb287..c2cbc5c0c2cba86dad9be60598450de2d1bfd7d4 100644
--- a/ogstools/meshlib/tests/test_simplified_mesh.py
+++ b/tests/meshlib/test_simplified_mesh.py
@@ -1,15 +1,17 @@
 import unittest
 
+from ogstools.definitions import TESTS_DIR
 from ogstools.meshlib._utils import dataframe_from_csv
 from ogstools.meshlib.boundary_set import LayerSet
 from ogstools.meshlib.region import to_region_simplified
-from ogstools.meshlib.tests import MeshPath
+
+meshpath = TESTS_DIR / "data" / "meshlib"
 
 
 class SimplifiedMeshTest(unittest.TestCase):
-    layerset = MeshPath("data/compose_geomodel/layersets.csv")
-    materialset = MeshPath("data/compose_geomodel/materialset.csv")
-    surfacedata = MeshPath("data/mesh1/surface_data/")
+    layerset = meshpath / "compose_geomodel/layersets.csv"
+    materialset = meshpath / "compose_geomodel/materialset.csv"
+    surfacedata = meshpath / "mesh1/surface_data/"
 
     def test_3D_points_and_cells_of_coarse_to_fine_meshes(self):
         mesh_df_coarseZ = dataframe_from_csv(
diff --git a/ogstools/meshlib/tests/test_surface.py b/tests/meshlib/test_surface.py
similarity index 85%
rename from ogstools/meshlib/tests/test_surface.py
rename to tests/meshlib/test_surface.py
index b39dcddb593b7db9a8d3383d25811dbeb964dc05..10273aa8f4ec89b7d7b322f2ffcb7b6b2bfbc842 100644
--- a/ogstools/meshlib/tests/test_surface.py
+++ b/tests/meshlib/test_surface.py
@@ -3,8 +3,10 @@ import unittest
 import numpy as np
 import pyvista as pv
 
+from ogstools.definitions import TESTS_DIR
 from ogstools.meshlib.boundary_subset import Surface
-from ogstools.meshlib.tests import MeshPath
+
+meshpath = TESTS_DIR / "data" / "meshlib"
 
 
 class SurfaceTest(unittest.TestCase):
@@ -16,7 +18,7 @@ class SurfaceTest(unittest.TestCase):
         """
         OK if file can be loaded - if not it raises an exception
         """
-        filename = MeshPath("data/mesh1/surface_data/00_KB.vtu")
+        filename = meshpath / "mesh1/surface_data/00_KB.vtu"
         s = Surface(filename, 0)  # checks
         self.assertEqual(s.filename, filename)
 
@@ -44,14 +46,14 @@ class SurfaceTest(unittest.TestCase):
         self.assertRaises(
             Exception,
             Surface.__init__,
-            MeshPath("data/mesh1/surface_data/notexisting.vtu"),
+            meshpath / "mesh1/surface_data/notexisting.vtu",
             0,
             0,
         )
 
     def testsurfaceToRaster(self):
         s1 = Surface(
-            MeshPath("data/mesh1/surface_data/00_KB.vtu"),
+            meshpath / "mesh1/surface_data/00_KB.vtu",
             0,
         )
         outfile = s1.create_raster_file(10)
diff --git a/ogstools/meshlib/tests/test_tetraeder_mesh.py b/tests/meshlib/test_tetraeder_mesh.py
similarity index 68%
rename from ogstools/meshlib/tests/test_tetraeder_mesh.py
rename to tests/meshlib/test_tetraeder_mesh.py
index 26ef05a82ff3e9a29ce2006a3424029af0ed3932..82b76e230b7180cb8dac4929b45975d984374b73 100644
--- a/ogstools/meshlib/tests/test_tetraeder_mesh.py
+++ b/tests/meshlib/test_tetraeder_mesh.py
@@ -1,14 +1,17 @@
 import unittest
 
+from ogstools.definitions import TESTS_DIR
+from ogstools.meshlib._utils import dataframe_from_csv
 from ogstools.meshlib.boundary_set import LayerSet
 from ogstools.meshlib.region import to_region_tetraeder
-from ogstools.meshlib.tests import MeshPath, dataframe_from_csv
+
+meshpath = TESTS_DIR / "data" / "meshlib"
 
 
 class TetraederTest(unittest.TestCase):
-    layerset = MeshPath("data/compose_geomodel/layersets.csv")
-    materialset = MeshPath("data/compose_geomodel/materialset.csv")
-    surfacedata = MeshPath("data/mesh1/surface_data/")
+    layerset = meshpath / "compose_geomodel/layersets.csv"
+    materialset = meshpath / "compose_geomodel/materialset.csv"
+    surfacedata = meshpath / "mesh1/surface_data/"
 
     def test_Mesh_coarseXYZ(self):
         mesh1_df = dataframe_from_csv(
diff --git a/ogstools/meshlib/tests/test_utils.py b/tests/meshlib/test_utils.py
similarity index 50%
rename from ogstools/meshlib/tests/test_utils.py
rename to tests/meshlib/test_utils.py
index 493a842e6a162ad7917d83eff2993c994e559914..2cfc7e055cf4027c6e38c6d124b16688e5c5c9e4 100644
--- a/ogstools/meshlib/tests/test_utils.py
+++ b/tests/meshlib/test_utils.py
@@ -1,16 +1,9 @@
 import unittest
-from pathlib import Path
 
+from ogstools.definitions import TESTS_DIR
 from ogstools.meshlib._utils import dataframe_from_csv
 
-
-def MeshPath(filenamepath: str):
-    """
-    Never use MeshPath in your projects.
-    It is supposed to be used in tests, only.
-    """
-    current_dir = Path(__file__).parent
-    return current_dir / Path(filenamepath)
+meshpath = TESTS_DIR / "data" / "meshlib"
 
 
 class GeoModelComposeExampleTest(unittest.TestCase):
@@ -22,9 +15,9 @@ class GeoModelComposeExampleTest(unittest.TestCase):
     def test_compose1(self):
         layerset3_df = dataframe_from_csv(
             3,
-            MeshPath("data/compose_geomodel/layersets.csv"),
-            MeshPath("data/compose_geomodel/materialset.csv"),
-            MeshPath("data/mesh1/surface_data/"),
+            meshpath / "compose_geomodel/layersets.csv",
+            meshpath / "compose_geomodel/materialset.csv",
+            meshpath / "mesh1/surface_data/",
         )
         self.assertEqual(len(layerset3_df), 5)
         # self.assertEqual(df["layer_id"][3], 3)
@@ -34,7 +27,7 @@ class GeoModelComposeExampleTest(unittest.TestCase):
             Exception,
             dataframe_from_csv,
             20,
-            MeshPath("data/compose_geomodel/layersets.csv"),
-            MeshPath("data/compose_geomodel/materialset.csv"),
-            MeshPath("data/mesh1/surface_data/"),
+            meshpath / "compose_geomodel/layersets.csv",
+            meshpath / "compose_geomodel/materialset.csv",
+            meshpath / "mesh1/surface_data/",
         )
diff --git a/ogstools/meshlib/tests/test_voxel_mesh.py b/tests/meshlib/test_voxel_mesh.py
similarity index 77%
rename from ogstools/meshlib/tests/test_voxel_mesh.py
rename to tests/meshlib/test_voxel_mesh.py
index 9471aade94387e216d4e01e5a45241bc090557ae..f78f9750d60e6de40f728be53ef3fe782647ca34 100644
--- a/ogstools/meshlib/tests/test_voxel_mesh.py
+++ b/tests/meshlib/test_voxel_mesh.py
@@ -1,15 +1,17 @@
 import unittest
 
+from ogstools.definitions import TESTS_DIR
 from ogstools.meshlib._utils import dataframe_from_csv
 from ogstools.meshlib.boundary_set import LayerSet
 from ogstools.meshlib.region import to_region_voxel
-from ogstools.meshlib.tests import MeshPath
+
+meshpath = TESTS_DIR / "data" / "meshlib"
 
 
 class VoxelMeshTest(unittest.TestCase):
-    layerset = MeshPath("data/compose_geomodel/layersets.csv")
-    materialset = MeshPath("data/compose_geomodel/materialset.csv")
-    surfacedata = MeshPath("data/mesh1/surface_data/")
+    layerset = meshpath / "compose_geomodel/layersets.csv"
+    materialset = meshpath / "compose_geomodel/materialset.csv"
+    surfacedata = meshpath / "mesh1/surface_data/"
 
     def test_create(self):
         mesh1_df = dataframe_from_csv(