diff --git a/README.md b/README.md
index 2932c49801e6c9e3e9051541d472d4bfc5244362..a91f7e9d8afad29f4ecb3caba9c14d060a95402a 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-![OGS workflow overview](https://ogstools.opengeosys.org/_static/ogstools.png "Supporting complex workflows - from preprocessing to simulation to postprocessing")
+![OGSTools](https://ogstools.opengeosys.org/_static/ogstools.png "Supporting complex workflows - from preprocessing to simulation to postprocessing")
 
 is a:
 
diff --git a/docs/development/index.md b/docs/development/index.md
index f4fccd2d1506f8efce73f148d4b7cbb8c456e3c5..c1bf2ce92781c36e58fa2b10aab559efc1bd8416 100644
--- a/docs/development/index.md
+++ b/docs/development/index.md
@@ -1,5 +1,20 @@
 # Development setup
 
+# Cloning the Source Repository
+
+You can clone the source repository from https://gitlab.opengeosys.org/ogs/tools/ogstools and install the latest version by running:
+
+`git clone git@gitlab.opengeosys.org:ogs/tools/ogstools.git` or
+`git clone https://gitlab.opengeosys.org/ogs/tools/ogstools.git`
+
+# Development environment
+
+Change into the directory with the cloned ogstools sources. (`.` assumes that your working directory contains the ogstools sources, the directory where you find the `pyproject.toml`).
+
+```bash
+cd ogstools
+```
+
 Create a virtual environment, activate it and install required packages:
 
 ```bash
diff --git a/docs/examples/howto_conversions/plot_C_feflowlib_2layers_model.py b/docs/examples/howto_conversions/plot_C_feflowlib_2layers_model.py
index 351cbc6d74996025881f26a407ca32648b83f0da..2c8ee75ce2f3abc8e4a8337e901a269eda4605e3 100644
--- a/docs/examples/howto_conversions/plot_C_feflowlib_2layers_model.py
+++ b/docs/examples/howto_conversions/plot_C_feflowlib_2layers_model.py
@@ -11,7 +11,7 @@ In this example we show how a simple FEFLOW model consisting of two layers can b
 # 1. Let us convert only the points and cells at first.
 import ifm_contrib as ifm
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools.examples import feflow_model_2layers
 from ogstools.feflowlib import (
     convert_geometry_mesh,
@@ -36,4 +36,4 @@ print(pv_mesh)
 pv_mesh.save("2layers_model.vtu")
 # %%
 # 4. Use the ogstools plotting functionalities.
-fig = ogs.plot.contourf(pv_mesh.slice("z"), "P_HEAD")
+fig = ot.plot.contourf(pv_mesh.slice("z"), "P_HEAD")
diff --git a/docs/examples/howto_conversions/plot_D_feflowlib_CT_simulation.py b/docs/examples/howto_conversions/plot_D_feflowlib_CT_simulation.py
index eea96efceb3f0a71f976e351eff2c796eb2da7e6..6831acd3e081d01fe9c60a6ad45efbffe0b2bb1d 100644
--- a/docs/examples/howto_conversions/plot_D_feflowlib_CT_simulation.py
+++ b/docs/examples/howto_conversions/plot_D_feflowlib_CT_simulation.py
@@ -17,7 +17,7 @@ import ifm_contrib as ifm
 import matplotlib.pyplot as plt
 import numpy as np
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools.examples import feflow_model_2D_CT_t_560
 from ogstools.feflowlib import (
     component_transport,
@@ -29,7 +29,7 @@ from ogstools.feflowlib import (
 )
 from ogstools.meshlib import Mesh
 
-ogs.plot.setup.show_element_edges = True
+ot.plot.setup.show_element_edges = True
 # %%
 # 1. Load a FEFLOW model (.fem) as a FEFLOW document, convert and save it. More details on
 # how the conversion function works can be found here: :py:mod:`ogstools.feflowlib.convert_properties_mesh`.
@@ -40,7 +40,7 @@ temp_dir = Path(tempfile.mkdtemp("feflow_test_simulation"))
 feflow_mesh_file = temp_dir / "2D_CT_model.vtu"
 feflow_pv_mesh.save(feflow_mesh_file)
 
-feflow_concentration = ogs.variables.Scalar(
+feflow_concentration = ot.variables.Scalar(
     data_name="single_species_P_CONC",
     output_name="concentration",
     data_unit="mg/l",
@@ -48,7 +48,7 @@ feflow_concentration = ogs.variables.Scalar(
 )
 # The original mesh is clipped to focus on the relevant part of it, where concentration is larger
 # than 1e-9 mg/l. The rest of the mesh has concentration values of 0.
-ogs.plot.contourf(
+ot.plot.contourf(
     feflow_pv_mesh.clip_scalar(
         scalars="single_species_P_CONC", invert=False, value=1.0e-9
     ),
@@ -61,7 +61,7 @@ write_point_boundary_conditions(temp_dir, feflow_pv_mesh)
 # %%
 # 3. Setup a prj-file (see: :py:mod:`ogstools.feflowlib.setup_prj_file`) to run a OGS-simulation.
 path_prjfile = feflow_mesh_file.with_suffix(".prj")
-prj = ogs.Project(output_file=path_prjfile)
+prj = ot.Project(output_file=path_prjfile)
 species = get_species(feflow_pv_mesh)
 CT_model = component_transport(
     saving_path=temp_dir / "sim_2D_CT_model",
@@ -95,7 +95,7 @@ ET.dump(ET.parse(path_prjfile))
 model.run_model(logfile=temp_dir / "out.log")
 # %%
 # 5. Read the results along a line on the upper edge of the mesh parallel to the x-axis and plot them.
-ms = ogs.MeshSeries(temp_dir / "sim_2D_CT_model.pvd")
+ms = ot.MeshSeries(temp_dir / "sim_2D_CT_model.pvd")
 # Read the last timestep:
 ogs_sim_res = ms.mesh(ms.timesteps[-1])
 """
@@ -108,7 +108,7 @@ profile = np.array([[0.038 + 1.0e-8, 0.005, 0], [0.045, 0.005, 0]])
 fig, ax = plt.subplots(1, 1, figsize=(7, 5))
 ogs_sim_res.plot_linesample(
     "dist",
-    ogs.variables.Scalar(
+    ot.variables.Scalar(
         data_name="single_species",
         output_name="concentration",
         data_unit="mg/l",
@@ -144,7 +144,7 @@ fig.tight_layout()
 ogs_sim_res["concentration_difference"] = (
     feflow_pv_mesh["single_species_P_CONC"] - ogs_sim_res["single_species"]
 )
-concentration_difference = ogs.variables.Scalar(
+concentration_difference = ot.variables.Scalar(
     data_name="concentration_difference",
     output_name="concentration",
     data_unit="mg/l",
@@ -153,7 +153,7 @@ concentration_difference = ogs.variables.Scalar(
 
 bounds = [0.038, 0.045, 0, 0.01, 0, 0]
 
-ogs.plot.contourf(
+ot.plot.contourf(
     ogs_sim_res.clip_box(bounds, invert=False),
     concentration_difference,
 )
diff --git a/docs/examples/howto_conversions/plot_E_feflowlib_H_simulation.py b/docs/examples/howto_conversions/plot_E_feflowlib_H_simulation.py
index 16fb39ae36b847abb62fb4de7d5dfd9648fad8b4..b474a301c5dcc5f489184072557e9130f04bedc2 100644
--- a/docs/examples/howto_conversions/plot_E_feflowlib_H_simulation.py
+++ b/docs/examples/howto_conversions/plot_E_feflowlib_H_simulation.py
@@ -18,7 +18,7 @@ import ifm_contrib as ifm
 import numpy as np
 import pyvista as pv
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools.examples import feflow_model_box_Neumann
 from ogstools.feflowlib import (
     convert_properties_mesh,
@@ -69,7 +69,7 @@ topsurface.save(path_topsurface)
 # %%
 # 3. Setup a prj-file (see: :py:mod:`ogstools.feflowlib.setup_prj_file`) to run a OGS-simulation.
 path_prjfile = feflow_mesh_file.with_suffix(".prj")
-prj = ogs.Project(output_file=path_prjfile)
+prj = ot.Project(output_file=path_prjfile)
 # Get the template prj-file configurations for a steady state diffusion process
 ssd_model = steady_state_diffusion(temp_dir / "sim_boxNeumann", prj)
 # Include the mesh specific configurations to the template.
@@ -90,7 +90,7 @@ ET.dump(model_prjfile)
 model.run_model(logfile=temp_dir / "out.log")
 # %%
 # 5. Read the results and plot them.
-ms = ogs.MeshSeries(temp_dir / "sim_boxNeumann.pvd")
+ms = ot.MeshSeries(temp_dir / "sim_boxNeumann.pvd")
 # Read the last timestep:
 ogs_sim_res = ms.mesh(ms.timesteps[-1])
 """
@@ -106,10 +106,8 @@ ogs_sim_res.plot(
 )
 # %%
 # 5.1 Plot the hydraulic head simulated in OGS with :py:mod:`ogstools.plot.contourf`.
-head = ogs.variables.Scalar(
-    data_name="HEAD_OGS", data_unit="m", output_unit="m"
-)
-fig = ogs.plot.contourf(ogs_sim_res.slice(normal="z", origin=[50, 50, 0]), head)
+head = ot.variables.Scalar(data_name="HEAD_OGS", data_unit="m", output_unit="m")
+fig = ot.plot.contourf(ogs_sim_res.slice(normal="z", origin=[50, 50, 0]), head)
 
 
 # %%
@@ -126,17 +124,17 @@ pyvista_mesh.plot(
 # %%
 # 6.1 Plot the differences in the hydraulic head with :py:mod:`ogstools.plot.contourf`.
 # Slices are taken along the z-axis.
-diff_head = ogs.variables.Scalar(
+diff_head = ot.variables.Scalar(
     data_name="diff_HEAD", data_unit="m", output_unit="m"
 )
 slices = np.reshape(list(pyvista_mesh.slice_along_axis(n=4, axis="z")), (2, 2))
-fig = ogs.plot.contourf(slices, diff_head)
+fig = ot.plot.contourf(slices, diff_head)
 for ax, slice in zip(fig.axes, np.ravel(slices), strict=False):
     ax.set_title(f"z = {slice.center[2]:.1f} {ms.spatial_output_unit}")
 
 # %%
 # Slices are taken along the y-axis.
 slices = np.reshape(list(pyvista_mesh.slice_along_axis(n=4, axis="y")), (2, 2))
-fig = ogs.plot.contourf(slices, diff_head)
+fig = ot.plot.contourf(slices, diff_head)
 for ax, slice in zip(fig.axes, np.ravel(slices), strict=False):
     ax.set_title(f"y = {slice.center[1]:.1f} {ms.spatial_output_unit}")
diff --git a/docs/examples/howto_conversions/plot_F_feflowlib_HT_simulation.py b/docs/examples/howto_conversions/plot_F_feflowlib_HT_simulation.py
index 56e35feead29a80607e90c9c32310e4935616085..672352366c454db05d1f453b3dc9ffa02bd14bf0 100644
--- a/docs/examples/howto_conversions/plot_F_feflowlib_HT_simulation.py
+++ b/docs/examples/howto_conversions/plot_F_feflowlib_HT_simulation.py
@@ -16,7 +16,7 @@ from pathlib import Path
 import ifm_contrib as ifm
 import pyvista as pv
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools.examples import feflow_model_2D_HT_model
 from ogstools.feflowlib import (
     convert_properties_mesh,
@@ -33,8 +33,8 @@ from ogstools.feflowlib.tools import (
 # how the conversion function works can be found here: :py:mod:`ogstools.feflowlib.convert_properties_mesh`.
 feflow_model = ifm.loadDocument(str(feflow_model_2D_HT_model))
 feflow_pv_mesh = convert_properties_mesh(feflow_model)
-feflow_temperature = ogs.variables.temperature.replace(data_name="P_TEMP")
-ogs.plot.contourf(feflow_pv_mesh, feflow_temperature)
+feflow_temperature = ot.variables.temperature.replace(data_name="P_TEMP")
+ot.plot.contourf(feflow_pv_mesh, feflow_temperature)
 
 temp_dir = Path(tempfile.mkdtemp("feflow_test_simulation"))
 feflow_mesh_file = temp_dir / "2D_HT_model.vtu"
@@ -54,7 +54,7 @@ plotter.show()
 # %%
 # 3. Setup a prj-file (see: :py:mod:`ogstools.feflowlib.setup_prj_file`) to run a OGS-simulation.
 path_prjfile = feflow_mesh_file.with_suffix(".prj")
-prj = ogs.Project(output_file=path_prjfile)
+prj = ot.Project(output_file=path_prjfile)
 # Get the template prj-file configurations for a hydro thermal process.
 HT_model = hydro_thermal(temp_dir / "sim_2D_HT_model", prj, dimension=2)
 # Include the mesh specific configurations to the template.
@@ -75,7 +75,7 @@ ET.dump(model_prjfile)
 model.run_model(logfile=temp_dir / "out.log")
 # %%
 # 5. Read the results and plot them.
-ms = ogs.MeshSeries(temp_dir / "sim_2D_HT_model.pvd")
+ms = ot.MeshSeries(temp_dir / "sim_2D_HT_model.pvd")
 # Read the last timestep:
 ogs_sim_res = ms.mesh(ms.timesteps[-1])
 """
@@ -85,28 +85,28 @@ ogs_sim_res = pv.read(
 )
 """
 # Plot the hydraulic head/height, which was simulated in OGS.
-hydraulic_head = ogs.variables.Scalar(
+hydraulic_head = ot.variables.Scalar(
     data_name="HEAD_OGS", data_unit="m", output_unit="m"
 )
-ogs.plot.contourf(ogs_sim_res, hydraulic_head)
+ot.plot.contourf(ogs_sim_res, hydraulic_head)
 # %%
 # Plot the temperature, which was simulated in OGS.
-ogs.plot.contourf(ogs_sim_res, ogs.variables.temperature)
+ot.plot.contourf(ogs_sim_res, ot.variables.temperature)
 
 # %%
 # 6. Plot the difference between the FEFLOW and OGS simulation.
 feflow_pv_mesh["HEAD"] = feflow_pv_mesh["P_HEAD"]
 ogs_sim_res["HEAD"] = ogs_sim_res["HEAD_OGS"]
 # Plot differences in hydraulic head/height.
-diff_mesh = ogs.meshlib.difference(feflow_pv_mesh, ogs_sim_res, "HEAD")
-hydraulic_head_diff = ogs.variables.Scalar(
+diff_mesh = ot.meshlib.difference(feflow_pv_mesh, ogs_sim_res, "HEAD")
+hydraulic_head_diff = ot.variables.Scalar(
     data_name="HEAD_difference", data_unit="m", output_unit="m"
 )
-ogs.plot.contourf(diff_mesh, hydraulic_head_diff)
+ot.plot.contourf(diff_mesh, hydraulic_head_diff)
 # %%
 feflow_pv_mesh["temperature"] = feflow_pv_mesh["P_TEMP"]
 # Plot differences in temperature.
-diff_mesh = ogs.meshlib.difference(
-    feflow_pv_mesh, ogs_sim_res, ogs.variables.temperature
+diff_mesh = ot.meshlib.difference(
+    feflow_pv_mesh, ogs_sim_res, ot.variables.temperature
 )
-ogs.plot.contourf(diff_mesh, ogs.variables.temperature.difference)
+ot.plot.contourf(diff_mesh, ot.variables.temperature.difference)
diff --git a/docs/examples/howto_plot/plot_animation.py b/docs/examples/howto_plot/plot_animation.py
index a38a3ef0f5d8be07ed1f9a4b051383c6e98d2e77..62c4e6d68be1185a384ec7fcaa7be917587f7fb2 100644
--- a/docs/examples/howto_plot/plot_animation.py
+++ b/docs/examples/howto_plot/plot_animation.py
@@ -13,7 +13,7 @@ example from the ogs benchmark gallery
 import matplotlib.pyplot as plt
 import numpy as np
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples
 
 mesh_series = examples.load_meshseries_CT_2D_XDMF()
@@ -46,7 +46,7 @@ timevalues = np.linspace(
 # method of MeshSeries.
 
 
-def mesh_func(mesh: ogs.Mesh) -> ogs.Mesh:
+def mesh_func(mesh: ot.Mesh) -> ot.Mesh:
     "Clip the left half of the mesh."
     return mesh.clip("-x", [0, 0, 0])
 
@@ -58,7 +58,7 @@ def plot_func(ax: plt.Axes, timevalue: float) -> None:
 
 # %%
 anim = mesh_series.transform(mesh_func).animate(
-    ogs.variables.saturation, timevalues, plot_func, vmin=0, vmax=100, dpi=50
+    ot.variables.saturation, timevalues, plot_func, vmin=0, vmax=100, dpi=50
 )
 
 # %% [markdown]
@@ -66,7 +66,7 @@ anim = mesh_series.transform(mesh_func).animate(
 #
 # ..  code-block:: python
 #
-#   ogs.plot.utils.save_animation(anim, "Saturation", fps=5)
+#   ot.plot.utils.save_animation(anim, "Saturation", fps=5)
 #
 
 # sphinx_gallery_start_ignore
diff --git a/docs/examples/howto_plot/plot_aspect_ratios.py b/docs/examples/howto_plot/plot_aspect_ratios.py
index 638d66feddd70f1168522969bef922cba693c413..319ed1cf507031435e37986dad51f9131776d92c 100644
--- a/docs/examples/howto_plot/plot_aspect_ratios.py
+++ b/docs/examples/howto_plot/plot_aspect_ratios.py
@@ -15,10 +15,10 @@ behaviour.
 import numpy as np
 import pyvista as pv
 
-import ogstools as ogs
+import ogstools as ot
 
-print(f"{ogs.plot.setup.min_ax_aspect=}")
-print(f"{ogs.plot.setup.max_ax_aspect=}")
+print(f"{ot.plot.setup.min_ax_aspect=}")
+print(f"{ot.plot.setup.max_ax_aspect=}")
 
 
 # sphinx_gallery_start_ignore
@@ -44,19 +44,19 @@ def custom_mesh(dx: float, dy: float):
 # proportions.
 
 # %%
-fig = ogs.plot.contourf(custom_mesh(np.pi * 2, np.pi), "example")
+fig = ot.plot.contourf(custom_mesh(np.pi * 2, np.pi), "example")
 # %% [markdown]
 # This one would be too wide and thus and gets compressed to fit the maximum
 # aspect ratio.
 
 # %%
-fig = ogs.plot.contourf(custom_mesh(np.pi * 4, np.pi), "example")
+fig = ot.plot.contourf(custom_mesh(np.pi * 4, np.pi), "example")
 # %% [markdown]
 # When plotting multiple meshes together, this applies to each subplot.
 # So here each subplot has true proportions again since each one fits the limits.
 
 # %%
-fig = ogs.plot.contourf(
+fig = ot.plot.contourf(
     [custom_mesh(np.pi * 2, np.pi), custom_mesh(np.pi * 2, np.pi)], "example"
 )
 # %% [markdown]
@@ -64,12 +64,12 @@ fig = ogs.plot.contourf(
 # ratio.
 
 # %%
-fig = ogs.plot.contourf(custom_mesh(np.pi, np.pi * 3), "example")
+fig = ot.plot.contourf(custom_mesh(np.pi, np.pi * 3), "example")
 # %% [markdown]
 # The same is true here:
 
 # %%
-fig = ogs.plot.contourf(
+fig = ot.plot.contourf(
     [custom_mesh(np.pi, np.pi * 3), custom_mesh(np.pi, np.pi * 3)], "example"
 )
 
@@ -79,12 +79,12 @@ fig = ogs.plot.contourf(
 # very wide figure.
 
 # %%
-ogs.plot.setup.min_ax_aspect = None
-ogs.plot.setup.max_ax_aspect = None
-fig = ogs.plot.contourf(custom_mesh(np.pi * 3, np.pi), "example")
+ot.plot.setup.min_ax_aspect = None
+ot.plot.setup.max_ax_aspect = None
+fig = ot.plot.contourf(custom_mesh(np.pi * 3, np.pi), "example")
 
 # %% [markdown]
 # And in this case we get a very tall figure.
 
 # %%
-fig = ogs.plot.contourf(custom_mesh(np.pi, np.pi * 3), "example")
+fig = ot.plot.contourf(custom_mesh(np.pi, np.pi * 3), "example")
diff --git a/docs/examples/howto_plot/plot_contourf_2d.py b/docs/examples/howto_plot/plot_contourf_2d.py
index 3549cccee467e274e8971c7f390aac15ea2f5ac2..a990e866f819e531b1dd4beb04cbe9220cee7415 100644
--- a/docs/examples/howto_plot/plot_contourf_2d.py
+++ b/docs/examples/howto_plot/plot_contourf_2d.py
@@ -13,10 +13,10 @@ available as keyword arguments in the function call. Please see
 """
 
 # %%
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples
 
-ogs.plot.setup.material_names = {i + 1: f"Layer {i+1}" for i in range(26)}
+ot.plot.setup.material_names = {i + 1: f"Layer {i+1}" for i in range(26)}
 mesh = examples.load_meshseries_THM_2D_PVD().mesh(1)
 
 # %% [markdown]
@@ -24,12 +24,12 @@ mesh = examples.load_meshseries_THM_2D_PVD().mesh(1)
 #
 # ..  code-block:: python
 #
-#   mesh_series = ogs.MeshSeries("filepath/filename_pvd_or_xdmf")
+#   mesh_series = ot.MeshSeries("filepath/filename_pvd_or_xdmf")
 #
 
 # %% First, let's plot the material ids (cell_data). Per default in
 # the setup, this will automatically show the element edges.
-fig = mesh.plot_contourf(ogs.variables.material_id)
+fig = mesh.plot_contourf(ot.variables.material_id)
 
 # %% [markdown]
 # Now, let's plot the temperature field (point_data) at the first timestep.
@@ -37,35 +37,35 @@ fig = mesh.plot_contourf(ogs.variables.material_id)
 # data as Kelvin and converts them to degrees Celsius.
 
 # %%
-fig = mesh.plot_contourf(ogs.variables.temperature, show_max=True)
+fig = mesh.plot_contourf(ot.variables.temperature, show_max=True)
 
 # %% [markdown]
 # We can also plot components of vector variables:
 
 # %%
 fig = mesh.plot_contourf(
-    ogs.variables.displacement[0], show_min=True, show_max=True
+    ot.variables.displacement[0], show_min=True, show_max=True
 )
 
 # %%
 fig = mesh.plot_contourf(
-    ogs.variables.displacement[1], show_max=True, show_edges=True
+    ot.variables.displacement[1], show_max=True, show_edges=True
 )
 
 # %% [markdown]
 # This example has hydraulically deactivated subdomains:
 
 # %%
-fig = mesh.plot_contourf(ogs.variables.pressure.get_mask(), fontsize=40)
+fig = mesh.plot_contourf(ot.variables.pressure.get_mask(), fontsize=40)
 
 # %% [markdown]
 # Let's plot the fluid velocity field.
 
 # %%
-fig = mesh.plot_contourf(ogs.variables.velocity, show_region_bounds=False)
+fig = mesh.plot_contourf(ot.variables.velocity, show_region_bounds=False)
 
 # %% [markdown]
 # Let's plot it again, this time log-scaled.
 
 # %%
-fig = mesh.plot_contourf(ogs.variables.velocity, log_scaled=True, vmin=-8)
+fig = mesh.plot_contourf(ot.variables.velocity, log_scaled=True, vmin=-8)
diff --git a/docs/examples/howto_plot/plot_contourf_3d.py b/docs/examples/howto_plot/plot_contourf_3d.py
index f21ea7e0f19e2c71749ae5d857efc659ce222e20..6450744fa8830c6ced9c13235775a0a659f4dfaf 100644
--- a/docs/examples/howto_plot/plot_contourf_3d.py
+++ b/docs/examples/howto_plot/plot_contourf_3d.py
@@ -16,10 +16,10 @@ data: "facies" with a native pyvista plot.
 import numpy as np
 from pyvista import examples
 
-import ogstools as ogs
+import ogstools as ot
 
 mesh = examples.load_channels()
-data = ogs.variables.Scalar("facies", categoric=True)
+data = ot.variables.Scalar("facies", categoric=True)
 mesh.plot(cmap="bwr")
 
 
@@ -27,7 +27,7 @@ mesh.plot(cmap="bwr")
 # Let's create multiple slices along the z axis and plot them in a 2 by 2 grid.
 
 slices = np.reshape(list(mesh.slice_along_axis(n=4, axis="z")), (2, 2))
-fig = ogs.plot.contourf(slices, data)
+fig = ot.plot.contourf(slices, data)
 for ax, slice in zip(fig.axes, np.ravel(slices), strict=False):
     ax.set_title(f"z = {slice.center[2]:.1f}")
 
@@ -35,7 +35,7 @@ for ax, slice in zip(fig.axes, np.ravel(slices), strict=False):
 # We can also slice along the y-axis and plot the meshes in one row.
 
 slices = np.reshape(mesh.slice_along_axis(n=3, axis="y"), (1, -1))
-fig = ogs.plot.contourf(slices, data)
+fig = ot.plot.contourf(slices, data)
 for ax, slice in zip(fig.axes, np.ravel(slices), strict=False):
     ax.set_title(f"y = {slice.center[1]:.1f}")
 
@@ -43,4 +43,4 @@ for ax, slice in zip(fig.axes, np.ravel(slices), strict=False):
 # Arbitrary oriented slices are also possible. They get projected to the
 # cardinal plane, from which they have the least rotational offset.
 
-fig = ogs.plot.contourf(mesh.slice([1, -2, 0]), data)
+fig = ot.plot.contourf(mesh.slice([1, -2, 0]), data)
diff --git a/docs/examples/howto_plot/plot_observation_points.py b/docs/examples/howto_plot/plot_observation_points.py
index 1d0663da4359c674d819f70095a5c8fdeb0ed587..80931ef3cd0b5fc47b6aa9486356d5e2e99409b0 100644
--- a/docs/examples/howto_plot/plot_observation_points.py
+++ b/docs/examples/howto_plot/plot_observation_points.py
@@ -24,11 +24,11 @@ Here we use a component transport example from the ogs benchmark gallery
 import matplotlib.pyplot as plt
 import numpy as np
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples
 
 mesh_series = examples.load_meshseries_CT_2D_XDMF()
-si = ogs.variables.saturation
+si = ot.variables.saturation
 
 # %% [markdown]
 # To read your own data as a mesh series you can do:
@@ -59,7 +59,7 @@ for i, point in enumerate(points):
 # And now probe the points and the values over time:
 
 # %%,
-labels = [f"{i}: {label}" for i, label in enumerate(ogs.plot.utils.justified_labels(points))]
+labels = [f"{i}: {label}" for i, label in enumerate(ot.plot.utils.justified_labels(points))]
 fig = mesh_series.plot_probe(
     points=points[:4], variable=si, time_unit="a", labels=labels[:4]
 )
diff --git a/docs/examples/howto_plot/plot_shared_axes.py b/docs/examples/howto_plot/plot_shared_axes.py
index 8789bfb8dea30b385603c085ef3c2a9c9b1f45fa..9832cf340a95f5f516a71ffb0b2748256e4430b4 100644
--- a/docs/examples/howto_plot/plot_shared_axes.py
+++ b/docs/examples/howto_plot/plot_shared_axes.py
@@ -12,20 +12,20 @@ subplots having shared axes.
 # Import packages, load example data set and define often used variables.
 import matplotlib.pyplot as plt
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples
 
 meshseries = examples.load_meshseries_THM_2D_PVD()
 mesh_0 = meshseries.mesh(0)
 mesh_1 = meshseries.mesh(1)
-variable = ogs.variables.temperature
+variable = ot.variables.temperature
 
 # %%
 # If you pass multiple meshes to :py:func:`ogstools.plot.contourf`
 # by default both x and y axes will shared. Thus, only the outer axes get
 # axes labels and tick label.
 
-fig = ogs.plot.contourf([mesh_0, mesh_1], variable)
+fig = ot.plot.contourf([mesh_0, mesh_1], variable)
 
 # %%
 # On user defined figure and axis the axis belonging to specific subplot has to
@@ -34,10 +34,10 @@ fig = ogs.plot.contourf([mesh_0, mesh_1], variable)
 fig, axs = plt.subplots(2, 2, figsize=(40, 17), sharex=True, sharey=True)
 diff_a = mesh_0.difference(mesh_1, variable)
 diff_b = mesh_1.difference(mesh_0, variable)
-ogs.plot.contourf(mesh_0, variable, fig=fig, ax=axs[0][0])
-ogs.plot.contourf(mesh_1, variable, fig=fig, ax=axs[1][0])
-ogs.plot.contourf(diff_a, variable, fig=fig, ax=axs[0][1])
-ogs.plot.contourf(diff_b, variable, fig=fig, ax=axs[1][1])
+ot.plot.contourf(mesh_0, variable, fig=fig, ax=axs[0][0])
+ot.plot.contourf(mesh_1, variable, fig=fig, ax=axs[1][0])
+ot.plot.contourf(diff_a, variable, fig=fig, ax=axs[0][1])
+ot.plot.contourf(diff_b, variable, fig=fig, ax=axs[1][1])
 fig.tight_layout()
 plt.show()
 
diff --git a/docs/examples/howto_plot/plot_timeslice.py b/docs/examples/howto_plot/plot_timeslice.py
index 72ea6c487b5e2a5c5f17a115774ac42a907862f6..33f3c211d4847e65389a106a7b4fea588ccab041 100644
--- a/docs/examples/howto_plot/plot_timeslice.py
+++ b/docs/examples/howto_plot/plot_timeslice.py
@@ -18,11 +18,11 @@ To see this benchmark results over all timesteps have a look at
 # vertical, horizontal and diagonal.
 import numpy as np
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples
 
 mesh_series = examples.load_meshseries_CT_2D_XDMF()
-si = ogs.variables.saturation
+si = ot.variables.saturation
 points_vert = np.linspace([25, 0, -75], [25, 0, 75], num=100)
 points_hori = np.linspace([0, 0, 60], [150, 0, 60], num=100)
 points_diag = np.linspace([25, 0, 75], [100, 0, 0], num=100)
diff --git a/docs/examples/howto_plot/plot_with_custom_fig_ax.py b/docs/examples/howto_plot/plot_with_custom_fig_ax.py
index 3230a0683a8b60546e8c368ca52147671524d8e3..15706060e1d78b908f6d4d43d509cbef96820b61 100644
--- a/docs/examples/howto_plot/plot_with_custom_fig_ax.py
+++ b/docs/examples/howto_plot/plot_with_custom_fig_ax.py
@@ -12,12 +12,12 @@ ogstools are to be used on different subplots within the same figure.
 # %%
 import matplotlib.pyplot as plt
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples
 
 meshseries = examples.load_meshseries_THM_2D_PVD()
 
-ogs.plot.setup.combined_colorbar = False
+ot.plot.setup.combined_colorbar = False
 
 # %% [markdown]
 # Compare different variables
@@ -28,8 +28,8 @@ ogs.plot.setup.combined_colorbar = False
 # same figure:
 
 fig, ax = plt.subplots(2, 1, figsize=(15, 15))
-meshseries.mesh(0).plot_contourf(ogs.variables.temperature, fig=fig, ax=ax[0])
-meshseries.mesh(1).plot_contourf(ogs.variables.displacement, fig=fig, ax=ax[1])
+meshseries.mesh(0).plot_contourf(ot.variables.temperature, fig=fig, ax=ax[0])
+meshseries.mesh(1).plot_contourf(ot.variables.displacement, fig=fig, ax=ax[1])
 fig.tight_layout()
 
 # %% [markdown]
@@ -43,14 +43,14 @@ fig.tight_layout()
 
 fig, ax = plt.subplots(3, 1, figsize=(20, 30))
 
-meshseries.mesh(0).plot_contourf(ogs.variables.temperature, fig=fig, ax=ax[0])
-meshseries.mesh(1).plot_contourf(ogs.variables.temperature, fig=fig, ax=ax[1])
+meshseries.mesh(0).plot_contourf(ot.variables.temperature, fig=fig, ax=ax[0])
+meshseries.mesh(1).plot_contourf(ot.variables.temperature, fig=fig, ax=ax[1])
 diff_mesh = meshseries.mesh(1).difference(
-    meshseries.mesh(0), ogs.variables.temperature
+    meshseries.mesh(0), ot.variables.temperature
 )
 
 
-diff_mesh.plot_contourf(ogs.variables.temperature, fig=fig, ax=ax[2])
+diff_mesh.plot_contourf(ot.variables.temperature, fig=fig, ax=ax[2])
 ax[0].set_title(r"$T(\mathrm{t}_{0})$")
 ax[1].set_title(r"$T(\mathrm{t}_{end})$")
 ax[2].set_title(r"$T(\mathrm{t}_{end})$-$T(\mathrm{t}_{0})$")
diff --git a/docs/examples/howto_postprocessing/plot_aggregate.py b/docs/examples/howto_postprocessing/plot_aggregate.py
index 393c4266aad32dbc8da60b82ed483494be631b47..1094f11e2be15e8a8c91e5a34ade80b63d91f3fb 100644
--- a/docs/examples/howto_postprocessing/plot_aggregate.py
+++ b/docs/examples/howto_postprocessing/plot_aggregate.py
@@ -14,11 +14,11 @@ To see this benchmark results over all timesteps have a look at
 """
 
 # %%
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples
 
 mesh_series = examples.load_meshseries_CT_2D_XDMF()
-saturation = ogs.variables.saturation
+saturation = ot.variables.saturation
 
 # %% [markdown]
 # To read your own data as a mesh series you can do:
@@ -52,7 +52,7 @@ fig = mesh.plot_contourf(saturation)
 
 # %%
 mesh = mesh_series.time_of_max(saturation)
-fig = mesh.plot_contourf(ogs.variables.Scalar("max_Saturation_time", "s", "a"))
+fig = mesh.plot_contourf(ot.variables.Scalar("max_Saturation_time", "s", "a"))
 
 # %% [markdown]
 # Likewise we can calculate and visualize the variance of the saturation:
diff --git a/docs/examples/howto_postprocessing/plot_calculate_diff.py b/docs/examples/howto_postprocessing/plot_calculate_diff.py
index 4f59c68d847b99704bd6f6391bde5d34f98ff1e7..56f2019db1cbeb0822f5a7740d0e4711142b7f32 100644
--- a/docs/examples/howto_postprocessing/plot_calculate_diff.py
+++ b/docs/examples/howto_postprocessing/plot_calculate_diff.py
@@ -17,7 +17,7 @@ This example shows how to calculate differences between meshes.
 
 # sphinx_gallery_end_ignore
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples
 
 # sphinx_gallery_start_ignore
@@ -27,35 +27,35 @@ from pathlib import Path
 from tempfile import mkdtemp
 
 
-ogs.plot.setup.dpi = 75
-ogs.plot.setup.show_element_edges = True
+ot.plot.setup.dpi = 75
+ot.plot.setup.show_element_edges = True
 
 tmp_dir = Path(mkdtemp())
 mesh_path = tmp_dir / "mesh.msh"
 vtu_path = tmp_dir / "mesh_domain.vtu"
 
 
-def custom_mesh(lengths: int, element_order: int, quads: bool) -> ogs.Mesh:
+def custom_mesh(lengths: int, element_order: int, quads: bool) -> ot.Mesh:
     "Creates a custom mesh and runs a Mechanics simulation on it."
-    ogs.meshlib.rect(
+    ot.meshlib.rect(
         lengths=lengths,
         n_edge_cells=21,
         structured_grid=quads,
         order=element_order,
         out_name=mesh_path,
     )
-    meshes = ogs.meshes_from_gmsh(mesh_path, log=False)
+    meshes = ot.meshes_from_gmsh(mesh_path, log=False)
     for name, mesh in meshes.items():
         pv.save_meshio(Path(tmp_dir, name + ".vtu"), mesh)
 
-    model = ogs.Project(
+    model = ot.Project(
         output_file=tmp_dir / "default.prj", input_file=examples.prj_mechanics
     )
     if element_order == 2:
         model.replace_text(4, ".//integration_order")
     model.write_input()
     model.run_model(write_logs=True, args=f"-m {tmp_dir} -o {tmp_dir}")
-    return ogs.MeshSeries(tmp_dir / "mesh.pvd").mesh(-1)
+    return ot.MeshSeries(tmp_dir / "mesh.pvd").mesh(-1)
 
 
 # sphinx_gallery_end_ignore
@@ -78,8 +78,8 @@ mesh2 = mesh_series.mesh(-1)
 # variable between the two provided meshes. Then we plot the difference:
 
 # %%
-mesh_diff = mesh1.difference(mesh2, ogs.variables.temperature)
-fig = mesh_diff.plot_contourf(ogs.variables.temperature)
+mesh_diff = mesh1.difference(mesh2, ot.variables.temperature)
+fig = mesh_diff.plot_contourf(ot.variables.temperature)
 
 # %% [markdown]
 # Difference between two meshes of different topology
@@ -96,7 +96,7 @@ tri_mesh = custom_mesh(lengths=1, element_order=2, quads=False)
 # size:
 
 # %%
-fig = ogs.plot.contourf([quad_mesh, tri_mesh], ogs.variables.stress["xx"])
+fig = ot.plot.contourf([quad_mesh, tri_mesh], ot.variables.stress["xx"])
 
 # %% [markdown]
 # To better quantify it we form the difference and plot the result. The base
@@ -105,8 +105,8 @@ fig = ogs.plot.contourf([quad_mesh, tri_mesh], ogs.variables.stress["xx"])
 # topologies is inherently affected by interpolation.
 
 # %%
-diff_mesh = quad_mesh.difference(tri_mesh, ogs.variables.stress)
-fig = diff_mesh.plot_contourf(ogs.variables.stress.difference["xx"])
+diff_mesh = quad_mesh.difference(tri_mesh, ot.variables.stress)
+fig = diff_mesh.plot_contourf(ot.variables.stress.difference["xx"])
 
 # %% [markdown]
 # Doing it the other way around shows the difference on the tri-mesh. Here, we
@@ -114,8 +114,8 @@ fig = diff_mesh.plot_contourf(ogs.variables.stress.difference["xx"])
 # which are outside of the domain of the second mesh are masked.
 
 # %%
-diff_mesh = tri_mesh.difference(quad_mesh, ogs.variables.stress)
-fig = diff_mesh.plot_contourf(ogs.variables.stress.difference["xx"])
+diff_mesh = tri_mesh.difference(quad_mesh, ot.variables.stress)
+fig = diff_mesh.plot_contourf(ot.variables.stress.difference["xx"])
 
 # %% [markdown]
 # Differences between multiple meshes
@@ -151,8 +151,8 @@ fig = diff_mesh.plot_contourf(ogs.variables.stress.difference["xx"])
 meshes_1 = [mesh1] * 3
 meshes_2 = [mesh2] * 3
 
-mesh_diff_pair_wise = ogs.meshlib.difference_pairwise(
-    meshes_1, meshes_2, ogs.variables.temperature
+mesh_diff_pair_wise = ot.meshlib.difference_pairwise(
+    meshes_1, meshes_2, ot.variables.temperature
 )
 print(f"Length of mesh_list1: {len(meshes_1)}")
 print(f"Length of mesh_list2: {len(meshes_2)}")
@@ -182,8 +182,8 @@ print(f"Shape of mesh_diff_pair_wise: {mesh_diff_pair_wise.shape}")
 # %%
 mesh_list = [mesh1, mesh2, mesh1, mesh2]
 
-mesh_diff_matrix = ogs.meshlib.difference_matrix(
-    mesh_list, variable=ogs.variables.temperature
+mesh_diff_matrix = ot.meshlib.difference_matrix(
+    mesh_list, variable=ot.variables.temperature
 )
 print(f"Length of mesh_list1: {len(mesh_list)}")
 print(f"Shape of mesh_list1: {mesh_diff_matrix.shape}")
@@ -218,8 +218,8 @@ print(f"Shape of mesh_list1: {mesh_diff_matrix.shape}")
 mesh_list_matrix_1 = [mesh1, mesh2, mesh1]
 mesh_list_matrix_2 = [mesh2, mesh1]
 
-mesh_diff_matrix = ogs.meshlib.difference_matrix(
-    mesh_list_matrix_1, mesh_list_matrix_2, ogs.variables.temperature
+mesh_diff_matrix = ot.meshlib.difference_matrix(
+    mesh_list_matrix_1, mesh_list_matrix_2, ot.variables.temperature
 )
 print(f"Length of mesh_list_matrix_1: {len(mesh_list_matrix_1)}")
 print(f"Length of mesh_list_matrix_2: {len(mesh_list_matrix_2)}")
diff --git a/docs/examples/howto_postprocessing/plot_convergence_study_nuclear_decay.py b/docs/examples/howto_postprocessing/plot_convergence_study_nuclear_decay.py
index edcad7ca5b9444c9f07fa0b0810622a792e9c1cc..c87a310dfed6ead5446ecbacb6642458682b9a0e 100644
--- a/docs/examples/howto_postprocessing/plot_convergence_study_nuclear_decay.py
+++ b/docs/examples/howto_postprocessing/plot_convergence_study_nuclear_decay.py
@@ -38,7 +38,7 @@ import pyvista as pv
 from IPython.display import HTML
 from scipy.constants import Julian_year as sec_per_yr
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples, physics, studies, workflow
 
 temp_dir = Path(mkdtemp(prefix="nuclear_decay"))
@@ -65,13 +65,11 @@ edge_cells = [5 * 2**i for i in range(n_refinements)]
 
 # %%
 for dt, n_cells in zip(time_step_sizes, edge_cells, strict=False):
-    ogs.meshlib.rect(
-        lengths=100.0, n_edge_cells=[n_cells, 1], out_name=msh_path
-    )
-    for name, mesh in ogs.meshes_from_gmsh(msh_path, log=False).items():
+    ot.meshlib.rect(lengths=100.0, n_edge_cells=[n_cells, 1], out_name=msh_path)
+    for name, mesh in ot.meshes_from_gmsh(msh_path, log=False).items():
         pv.save_meshio(Path(temp_dir, name + ".vtu"), mesh)
 
-    prj = ogs.Project(output_file=temp_dir / "default.prj", input_file=prj_path)
+    prj = ot.Project(output_file=temp_dir / "default.prj", input_file=prj_path)
     prj.replace_text(str(dt * sec_per_yr), ".//delta_t")
     prj.replace_text(prefix.format(dt), ".//prefix")
     prj.write_input()
@@ -90,12 +88,12 @@ fig, (ax1, ax2) = plt.subplots(figsize=(8, 8), nrows=2, sharex=True)
 ax2.plot(time, heat, lw=2, label="reference", color="k")
 
 for sim_result, dt in zip(sim_results, time_step_sizes, strict=False):
-    mesh_series = ogs.MeshSeries(sim_result)
+    mesh_series = ot.MeshSeries(sim_result)
     results = {"heat_flux": [], "temperature": []}
     for ts in mesh_series.timesteps:
         mesh = mesh_series.mesh(ts)
         results["temperature"] += [np.max(mesh.point_data["temperature"])]
-    max_T = ogs.variables.temperature.transform(results["temperature"])
+    max_T = ot.variables.temperature.transform(results["temperature"])
     # times 2 due to symmetry, area of repo, to kW
     results["heat_flux"] += [np.max(mesh.point_data["heat_flux"][:, 0])]
     tv = np.asarray(mesh_series.timevalues("a"))
@@ -162,9 +160,9 @@ HTML(workflow.jupyter_to_html(report_name, show_input=False))
 # model behavior.
 
 # %%
-mesh_series = [ogs.MeshSeries(sim_result) for sim_result in sim_results]
+mesh_series = [ot.MeshSeries(sim_result) for sim_result in sim_results]
 evolution_metrics = studies.convergence.convergence_metrics_evolution(
-    mesh_series, ogs.variables.temperature, units=["s", "yrs"]
+    mesh_series, ot.variables.temperature, units=["s", "yrs"]
 )
 
 # %% [markdown]
diff --git a/docs/examples/howto_postprocessing/plot_convergence_study_steady_state_diffusion.py b/docs/examples/howto_postprocessing/plot_convergence_study_steady_state_diffusion.py
index 5c9dd841122dc98192f5fe0c28b09b4b338250e0..1647011ce7c141f0e6f0f1d93aeb2d07a3a92d87 100644
--- a/docs/examples/howto_postprocessing/plot_convergence_study_steady_state_diffusion.py
+++ b/docs/examples/howto_postprocessing/plot_convergence_study_steady_state_diffusion.py
@@ -37,7 +37,7 @@ from tempfile import mkdtemp
 import pyvista as pv
 from IPython.display import HTML
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples, workflow
 from ogstools.studies import convergence
 
@@ -56,14 +56,14 @@ refinements = 6
 edge_cells = [2**i for i in range(refinements)]
 for n_edge_cells in edge_cells:
     msh_path = temp_dir / "square.msh"
-    ogs.meshlib.rect(
+    ot.meshlib.rect(
         n_edge_cells=n_edge_cells, structured_grid=True, out_name=msh_path
     )
-    meshes = ogs.meshes_from_gmsh(filename=msh_path, log=False)
+    meshes = ot.meshes_from_gmsh(filename=msh_path, log=False)
     for name, mesh in meshes.items():
         pv.save_meshio(Path(temp_dir, name + ".vtu"), mesh)
 
-    model = ogs.Project(
+    model = ot.Project(
         output_file=temp_dir / "default.prj",
         input_file=examples.prj_steady_state_diffusion,
     )
@@ -80,10 +80,10 @@ for n_edge_cells in edge_cells:
 # %%
 analytical_solution_path = temp_dir / "analytical_solution.vtu"
 solution = examples.analytical_diffusion(
-    ogs.MeshSeries(result_paths[-1]).mesh(0)
+    ot.MeshSeries(result_paths[-1]).mesh(0)
 )
-ogs.plot.setup.show_element_edges = True
-fig = ogs.plot.contourf(solution, ogs.variables.hydraulic_head)
+ot.plot.setup.show_element_edges = True
+fig = ot.plot.contourf(solution, ot.variables.hydraulic_head)
 solution.save(analytical_solution_path)
 
 # %% [markdown]
diff --git a/docs/examples/howto_postprocessing/plot_ipdata.py b/docs/examples/howto_postprocessing/plot_ipdata.py
index c8f929dee9b9879d4c7c70dd95e1b5d769f7640b..6e41d4438c6ca63e9c4b5036b5840a693f06c8b3 100644
--- a/docs/examples/howto_postprocessing/plot_ipdata.py
+++ b/docs/examples/howto_postprocessing/plot_ipdata.py
@@ -22,14 +22,14 @@ from tempfile import mkdtemp
 
 import pyvista as pv
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples
 from ogstools.meshlib.gmsh_meshing import rect
 
-ogs.plot.setup.dpi = 75
-ogs.plot.setup.show_element_edges = True
+ot.plot.setup.dpi = 75
+ot.plot.setup.show_element_edges = True
 
-sigma_ip = ogs.variables.stress.replace(
+sigma_ip = ot.variables.stress.replace(
     data_name="sigma_ip", output_name="IP_stress"
 )
 
@@ -45,22 +45,22 @@ def simulate_and_plot(elem_order: int, quads: bool, intpt_order: int):
         order=elem_order,
         out_name=msh_path,
     )
-    meshes = ogs.meshes_from_gmsh(msh_path, log=False)
+    meshes = ot.meshes_from_gmsh(msh_path, log=False)
     for name, mesh in meshes.items():
         pv.save_meshio(Path(tmp_dir, name + ".vtu"), mesh)
 
-    model = ogs.Project(
+    model = ot.Project(
         output_file=tmp_dir / "default.prj",
         input_file=examples.prj_mechanics,
     )
     model.replace_text(intpt_order, xpath=".//integration_order")
     model.write_input()
     model.run_model(write_logs=True, args=f"-m {tmp_dir} -o {tmp_dir}")
-    mesh = ogs.MeshSeries(tmp_dir / "mesh.pvd").mesh(-1)
+    mesh = ot.MeshSeries(tmp_dir / "mesh.pvd").mesh(-1)
     int_pts = mesh.to_ip_point_cloud()
     ip_mesh = mesh.to_ip_mesh()
 
-    fig = mesh.plot_contourf(ogs.variables.stress)
+    fig = mesh.plot_contourf(ot.variables.stress)
     fig.axes[0].scatter(
         int_pts.points[:, 0], int_pts.points[:, 1], color="k", s=10
     )
diff --git a/docs/examples/howto_postprocessing/plot_sample_mesh_line.py b/docs/examples/howto_postprocessing/plot_sample_mesh_line.py
index e0864b20ffdf5f5513c87eb74680b292f40de6aa..21a59f950f5da737b4bb017d53f99275dd22b9f3 100644
--- a/docs/examples/howto_postprocessing/plot_sample_mesh_line.py
+++ b/docs/examples/howto_postprocessing/plot_sample_mesh_line.py
@@ -9,7 +9,7 @@ Extract a 1D profile from 2D and plot it
 import matplotlib.pyplot as plt
 import numpy as np
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples
 
 # %% [markdown]
@@ -24,7 +24,7 @@ mesh = examples.load_meshseries_HT_2D_XDMF().mesh(-1)
 profile_HT = np.array([[4, 2, 0], [4, 18, 0]])
 
 # %%
-mesh_sp, mesh_kp = ogs.meshlib.sample_polyline(
+mesh_sp, mesh_kp = ot.meshlib.sample_polyline(
     mesh, ["pressure", "temperature"], profile_HT
 )
 
@@ -64,9 +64,9 @@ ax_twinx = mesh.plot_linesample(
     ax=ax_twinx,
     fontsize=15,
 )
-ogs.plot.utils.color_twin_axes(
+ot.plot.utils.color_twin_axes(
     [ax, ax_twinx],
-    [ogs.variables.pressure.color, ogs.variables.temperature.color],
+    [ot.variables.pressure.color, ot.variables.temperature.color],
 )
 fig.tight_layout()
 
@@ -76,7 +76,7 @@ fig.tight_layout()
 # We can see it in the following example using the Darcy velocity:
 
 # %%
-mesh_sp, mesh_kp = ogs.meshlib.sample_polyline(
+mesh_sp, mesh_kp = ot.meshlib.sample_polyline(
     mesh, "darcy_velocity", profile_HT
 )
 
@@ -106,8 +106,8 @@ profile_CT = np.array([[47.0, 1.17, 72.0], [-4.5, 1.17, -59.0]])
 mesh = examples.load_meshseries_CT_2D_XDMF().mesh(11)
 
 # %%
-mesh_sp, mesh_kp = ogs.meshlib.sample_polyline(
-    mesh, ogs.variables.saturation, profile_CT
+mesh_sp, mesh_kp = ot.meshlib.sample_polyline(
+    mesh, ot.variables.saturation, profile_CT
 )
 
 # %% [markdown]
@@ -124,7 +124,7 @@ mesh_sp.head(5)
 
 # %%
 fig, ax = mesh.plot_linesample_contourf(
-    ogs.variables.saturation, profile_CT, resolution=100
+    ot.variables.saturation, profile_CT, resolution=100
 )
 
 # %% [markdown]
@@ -159,9 +159,9 @@ profile_THM = np.array(
 mesh = examples.load_meshseries_THM_2D_PVD().mesh(-1)
 
 # %%
-ms_THM_sp, dist_at_knot = ogs.meshlib.sample_polyline(
+ms_THM_sp, dist_at_knot = ot.meshlib.sample_polyline(
     mesh,
-    [ogs.variables.pressure, ogs.variables.temperature],
+    [ot.variables.pressure, ot.variables.temperature],
     profile_THM,
     resolution=100,
 )
@@ -202,7 +202,7 @@ fig.tight_layout()
 # %%
 # plt.rcdefaults()
 fig, ax = mesh.plot_linesample_contourf(
-    [ogs.variables.pressure, ogs.variables.temperature],
+    [ot.variables.pressure, ot.variables.temperature],
     profile_THM,
     resolution=100,
 )
diff --git a/docs/examples/howto_postprocessing/plot_variables.py b/docs/examples/howto_postprocessing/plot_variables.py
index 27397b084f60d391ffc7f49cb5043182aadd1589..74d6b1ba827ece42202eccdde0042779f7dd8a26 100644
--- a/docs/examples/howto_postprocessing/plot_variables.py
+++ b/docs/examples/howto_postprocessing/plot_variables.py
@@ -10,10 +10,10 @@ mesh data. There are also several predefined variables.
 """
 
 # %%
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples
 
-ogs.variables.get_dataframe()
+ot.variables.get_dataframe()
 
 # %% [markdown]
 # Scalar, Vector and Matrix inherit from the class Variable with its
@@ -22,21 +22,21 @@ ogs.variables.get_dataframe()
 # applies a function if specified. In this case we convert from K to °C:
 
 # %%
-ogs.variables.temperature.transform(273.15, strip_unit=False)
+ot.variables.temperature.transform(273.15, strip_unit=False)
 
 # %% [markdown]
 # You can also create your own variables by creating a Scalar, Vector or Matrix
 # variable. The following doesn't do any unit conversion.
 
 # %%
-custom_temperature = ogs.variables.Scalar(
+custom_temperature = ot.variables.Scalar(
     data_name="temperature", data_unit="K", output_unit="K"
 )
 custom_temperature.transform(273.15, strip_unit=False)
 
 # %% [markdown]
 # Or use existing presets as a template and replace some parameters:
-custom_temperature = ogs.variables.temperature.replace(output_unit="°F")
+custom_temperature = ot.variables.temperature.replace(output_unit="°F")
 custom_temperature.transform(273.15, strip_unit=False)
 
 # %% [markdown]
@@ -47,10 +47,10 @@ custom_temperature.transform(273.15, strip_unit=False)
 # length 4 [xx, yy, zz, xy] or 6 [xx, yy, zz, xy, yz, xz].
 
 # %%
-ogs.variables.displacement[1].transform([0.01, 0.02, 0.03], strip_unit=False)
+ot.variables.displacement[1].transform([0.01, 0.02, 0.03], strip_unit=False)
 
 # %%
-ogs.variables.strain["xx"].transform(
+ot.variables.strain["xx"].transform(
     [0.01, 0.02, 0.03, 0.04, 0.05, 0.06], strip_unit=False
 )
 
@@ -58,7 +58,7 @@ ogs.variables.strain["xx"].transform(
 # Magnitude of a 2D displacement vector:
 
 # %%
-ogs.variables.displacement.magnitude.transform([0.03, 0.04], strip_unit=False)
+ot.variables.displacement.magnitude.transform([0.03, 0.04], strip_unit=False)
 
 # %% [markdown]
 # We suggest specifying the variables and their transformations once.
@@ -67,8 +67,8 @@ ogs.variables.displacement.magnitude.transform([0.03, 0.04], strip_unit=False)
 # task of processing the data (e.g. calculate the von Mises stress):
 
 # %%
-fig = ogs.plot.contourf(
-    examples.load_mesh_mechanics_2D(), ogs.variables.stress.von_Mises
+fig = ot.plot.contourf(
+    examples.load_mesh_mechanics_2D(), ot.variables.stress.von_Mises
 )
 
 # %% [markdown]
diff --git a/docs/examples/howto_preprocessing/plot_remeshing.py b/docs/examples/howto_preprocessing/plot_remeshing.py
index 6be519dc1a5adae3bdb4b80290f850ff6861892b..f86f2fcc97f340ae7463dccb333806535ce065e2 100644
--- a/docs/examples/howto_preprocessing/plot_remeshing.py
+++ b/docs/examples/howto_preprocessing/plot_remeshing.py
@@ -21,14 +21,14 @@ to adapt.
 from pathlib import Path
 from tempfile import mkdtemp
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples
 
 mesh = examples.load_meshseries_THM_2D_PVD().mesh(1)
 
 # %% This is our example mesh which we want to discretize with triangle
 # elements.
-fig = mesh.plot_contourf(ogs.variables.material_id)
+fig = mesh.plot_contourf(ot.variables.material_id)
 
 # %% [markdown]
 # Here, we do the remeshing and convert the resulting msh file to an
@@ -38,6 +38,6 @@ fig = mesh.plot_contourf(ogs.variables.material_id)
 mesh = examples.load_meshseries_THM_2D_PVD().mesh(1)
 temp_dir = Path(mkdtemp())
 msh_path = temp_dir / "tri_mesh.msh"
-ogs.meshlib.gmsh_meshing.remesh_with_triangles(mesh, msh_path)
-meshes = ogs.meshes_from_gmsh(msh_path, reindex=False, log=False)
-fig = meshes["domain"].plot_contourf(ogs.variables.material_id)
+ot.meshlib.gmsh_meshing.remesh_with_triangles(mesh, msh_path)
+meshes = ot.meshes_from_gmsh(msh_path, reindex=False, log=False)
+fig = meshes["domain"].plot_contourf(ot.variables.material_id)
diff --git a/docs/examples/howto_prjfile/plot_creation.py b/docs/examples/howto_prjfile/plot_creation.py
index ffd1f590c17ebc2c0eb0127eddbc33f938d71852..9ef53820120b9f79448d778f60e1e9b51e168610 100644
--- a/docs/examples/howto_prjfile/plot_creation.py
+++ b/docs/examples/howto_prjfile/plot_creation.py
@@ -15,11 +15,11 @@ The names of the method calls are based on the corresponding XML tags.
 from pathlib import Path
 from tempfile import mkdtemp
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools.definitions import EXAMPLES_DIR
 
 output_dir = Path(mkdtemp())
-prj = ogs.Project(output_file=output_dir / "mechanics_new.prj")
+prj = ot.Project(output_file=output_dir / "mechanics_new.prj")
 
 # %%
 # Define geometry and/or meshes:
diff --git a/docs/examples/howto_prjfile/plot_manipulation.py b/docs/examples/howto_prjfile/plot_manipulation.py
index f8916197ab0bbf2051b411e4e8b52f4cf079bf67..618ab9cfa6c3945be12f12bdf79494f4e5907bce 100644
--- a/docs/examples/howto_prjfile/plot_manipulation.py
+++ b/docs/examples/howto_prjfile/plot_manipulation.py
@@ -12,14 +12,14 @@ E.g., to iterate over three Young's moduli one can use the replace parameter met
 from pathlib import Path
 from tempfile import mkdtemp
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools.definitions import EXAMPLES_DIR
 
 model_dir = Path(mkdtemp())
 youngs_moduli = [1, 2, 3]
 filename = EXAMPLES_DIR / "prj/simple_mechanics.prj"
 for youngs_modulus in youngs_moduli:
-    prj = ogs.Project(input_file=filename, output_file=model_dir / filename)
+    prj = ot.Project(input_file=filename, output_file=model_dir / filename)
     prj.replace_parameter_value(name="E", value=youngs_modulus)
     prj.replace_text(
         f"out_E={youngs_modulus}", xpath="./time_loop/output/prefix"
diff --git a/docs/examples/howto_quickstart/plot_solid_mechanics.py b/docs/examples/howto_quickstart/plot_solid_mechanics.py
index 920dfa3157f9cbfa0a9d5e48da199a1901d8f04a..12d7bad15428dee7b24c3a4f0ea89bdde83d6b0e 100644
--- a/docs/examples/howto_quickstart/plot_solid_mechanics.py
+++ b/docs/examples/howto_quickstart/plot_solid_mechanics.py
@@ -21,11 +21,11 @@ stress analysis:
 # sphinx_gallery_end_ignore
 
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples
 
 mesh = examples.load_mesh_mechanics_2D()
-fig = mesh.plot_contourf(ogs.variables.displacement)
+fig = mesh.plot_contourf(ot.variables.displacement)
 
 # %% [markdown]
 # Tensor components
@@ -33,8 +33,8 @@ fig = mesh.plot_contourf(ogs.variables.displacement)
 # We can inspect the stress (or strain) tensor components by indexing.
 
 # %%
-fig = mesh.plot_contourf(ogs.variables.stress["xx"])
-fig = mesh.plot_contourf(ogs.variables.stress["xy"])
+fig = mesh.plot_contourf(ot.variables.stress["xx"])
+fig = mesh.plot_contourf(ot.variables.stress["xy"])
 
 # %% [markdown]
 # Principal stresses
@@ -45,16 +45,16 @@ fig = mesh.plot_contourf(ogs.variables.stress["xy"])
 # compressive principal stress.
 
 # %%
-eigvecs = ogs.variables.stress.eigenvectors
-fig = mesh.plot_contourf(variable=ogs.variables.stress.eigenvalues[0])
+eigvecs = ot.variables.stress.eigenvectors
+fig = mesh.plot_contourf(variable=ot.variables.stress.eigenvalues[0])
 mesh.plot_quiver(ax=fig.axes[0], variable=eigvecs[0], glyph_type="line")
 
 # %%
-fig = mesh.plot_contourf(variable=ogs.variables.stress.eigenvalues[1])
+fig = mesh.plot_contourf(variable=ot.variables.stress.eigenvalues[1])
 mesh.plot_quiver(ax=fig.axes[0], variable=eigvecs[1], glyph_type="line")
 
 # %%
-fig = mesh.plot_contourf(variable=ogs.variables.stress.eigenvalues[2])
+fig = mesh.plot_contourf(variable=ot.variables.stress.eigenvalues[2])
 mesh.plot_quiver(ax=fig.axes[0], variable=eigvecs[2], glyph_type="line")
 
 # %% [markdown]
@@ -63,7 +63,7 @@ mesh.plot_quiver(ax=fig.axes[0], variable=eigvecs[2], glyph_type="line")
 # see: :py:func:`ogstools.variables.tensor_math.mean`
 
 # %%
-fig = mesh.plot_contourf(ogs.variables.stress.mean)
+fig = mesh.plot_contourf(ot.variables.stress.mean)
 
 # %% [markdown]
 # Von Mises stress
@@ -71,7 +71,7 @@ fig = mesh.plot_contourf(ogs.variables.stress.mean)
 # see: :py:func:`ogstools.variables.tensor_math.von_mises`
 
 # %%
-fig = mesh.plot_contourf(ogs.variables.stress.von_Mises)
+fig = mesh.plot_contourf(ot.variables.stress.von_Mises)
 
 # %% [markdown]
 # octahedral shear stress
@@ -79,7 +79,7 @@ fig = mesh.plot_contourf(ogs.variables.stress.von_Mises)
 # see: :py:func:`ogstools.variables.tensor_math.octahedral_shear`
 
 # %%
-fig = mesh.plot_contourf(ogs.variables.stress.octahedral_shear)
+fig = mesh.plot_contourf(ot.variables.stress.octahedral_shear)
 
 # %% [markdown]
 # Integrity criteria
@@ -96,7 +96,7 @@ fig = mesh.plot_contourf(ogs.variables.stress.octahedral_shear)
 
 # %%
 mesh["pressure"] = mesh.p_fluid()
-fig = mesh.plot_contourf(ogs.variables.pressure)
+fig = mesh.plot_contourf(ot.variables.pressure)
 
 # %% [markdown]
 # But since this assumes that the top of the model is equal to the ground
@@ -107,7 +107,7 @@ fig = mesh.plot_contourf(ogs.variables.pressure)
 mesh["depth"] = mesh.depth(use_coords=True)
 fig = mesh.plot_contourf("depth")
 mesh["pressure"] = mesh.p_fluid()
-fig = mesh.plot_contourf(ogs.variables.pressure)
+fig = mesh.plot_contourf(ot.variables.pressure)
 
 # %% [markdown]
 # Dilantancy criterion
@@ -115,8 +115,8 @@ fig = mesh.plot_contourf(ogs.variables.pressure)
 # see: :py:func:`ogstools.variables.mesh_dependent.dilatancy_critescu`
 
 # %%
-fig = mesh.plot_contourf(ogs.variables.dilatancy_critescu_tot)
-fig = mesh.plot_contourf(ogs.variables.dilatancy_critescu_eff)
+fig = mesh.plot_contourf(ot.variables.dilatancy_critescu_tot)
+fig = mesh.plot_contourf(ot.variables.dilatancy_critescu_eff)
 
 # %% [markdown]
 # Fluid pressure criterion
@@ -124,4 +124,4 @@ fig = mesh.plot_contourf(ogs.variables.dilatancy_critescu_eff)
 # see: :py:func:`ogstools.variables.mesh_dependent.fluid_pressure_criterion`
 
 # %%
-fig = mesh.plot_contourf(ogs.variables.fluid_pressure_crit)
+fig = mesh.plot_contourf(ot.variables.fluid_pressure_crit)
diff --git a/docs/releases/ogstools-0.4.0.md b/docs/releases/ogstools-0.4.0.md
index 7c94743102db5661a57f021ad41afbf8a0abbfe2..c63495b9352fbeb990f94cbf06fc687110f2e017 100644
--- a/docs/releases/ogstools-0.4.0.md
+++ b/docs/releases/ogstools-0.4.0.md
@@ -17,7 +17,7 @@ The recommended strategy is:
 - Back up your environment with `pip freeze > yourfile.txt`
 - Upgrade ogs6py to 0.403 first and resolve issues (see breaking changes).
 - Uninstall ogs6py and install ogstools
-- Remove `import ogs6py`. Add `import ogstools as ogs`. The former `OGS` becomes `Project` and function parameters of `OGS.__init__()` are now with lower case names. See example in Features
+- Remove `import ogs6py`. Add `import ogstools as ot`. The former `OGS` becomes `Project` and function parameters of `OGS.__init__()` are now with lower case names. See example in Features
 
 ## API changes
 
@@ -51,9 +51,9 @@ The recommended strategy is:
 before:
 
 ```python
-import ogstools as ogs
+import ogstools as ot
 
-prj = ogs.Project(input_file="mechanics.prj", output_file="old_parameter_add.prj")
+prj = ot.Project(input_file="mechanics.prj", output_file="old_parameter_add.prj")
 prj.add_block(
     blocktag="parameter",
     parent_xpath="./parameters",
@@ -66,9 +66,9 @@ prj.write_input()
 now:
 
 ```python
-import ogstools as ogs
+import ogstools as ot
 
-prj = ogs.Project(input_file="mechanics.prj", output_file="new_parameter_add.prj")
+prj = ot.Project(input_file="mechanics.prj", output_file="new_parameter_add.prj")
 prj.parameters.add_parameter(name="density", type="Constant", value="1900")
 prj.write_input()
 ```
@@ -96,7 +96,7 @@ prj.write_input()
 - Mesh can be created from a shapefile
 
 ```
-mesh = ogs.Mesh.read(test_shapefile)
+mesh = ot.Mesh.read(test_shapefile)
 ```
 
 - Mesh can be remeshed with triangle
@@ -105,9 +105,9 @@ mesh = ogs.Mesh.read(test_shapefile)
 - MeshSeries allows multidimensional indexing on ndarrays <https://numpy.org/doc/stable/user/basics.indexing.html>
 
 ```python
-import ogstools as ogs
+import ogstools as ot
 
-ms = ogs.MeshSeries("filepath/filename_pvd_or_xdmf")
+ms = ot.MeshSeries("filepath/filename_pvd_or_xdmf")
 ms.data("darcy_velocity")[-2:, 1:4, :]  # shape is(2, 3, 2)
 result_mesh = ms[-1]
 for mesh in ms:
diff --git a/ogstools/examples/ogs6py/example_THM.py b/ogstools/examples/ogs6py/example_THM.py
index 0295419dea68e2865d6ec8fa61f2012643ef11f9..11624d9bab8f30435f36a6b704a57dd964a24ca3 100644
--- a/ogstools/examples/ogs6py/example_THM.py
+++ b/ogstools/examples/ogs6py/example_THM.py
@@ -1,8 +1,8 @@
-import ogstools as ogs
+import ogstools as ot
 
 # if MKL set vars vars script should be executed before OGS
-# prj = ogs.Project(output_file="thm_test.prj", MKL=True, OMP_NUM_THREADS=4)
-prj = ogs.Project(output_file="thm_test.prj", OMP_NUM_THREADS=4)
+# prj = ot.Project(output_file="thm_test.prj", MKL=True, OMP_NUM_THREADS=4)
+prj = ot.Project(output_file="thm_test.prj", OMP_NUM_THREADS=4)
 prj.geometry.add_geometry(filename="square_1x1_thm.gml")
 prj.mesh.add_mesh(filename="quarter_002_2nd.vtu", axially_symmetric="true")
 prj.processes.set_process(
diff --git a/tests/test_meshlib.py b/tests/test_meshlib.py
index 435282e0383269dec74210d860d708395f55fc0f..108b410a3ec4ba65b5a9a2c4975dd16a069fad28 100644
--- a/tests/test_meshlib.py
+++ b/tests/test_meshlib.py
@@ -8,7 +8,7 @@ import pkg_resources
 import pytest
 import pyvista as pv
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples
 
 
@@ -37,12 +37,12 @@ class TestUtils:
         tmp_dir = Path(mkdtemp())
         mesh_path = tmp_dir / "mesh.msh"
         for quads in [True, False]:
-            ogs.meshlib.rect(
+            ot.meshlib.rect(
                 1, 1, structured_grid=quads, order=2, out_name=mesh_path
             )
-            ogs.meshes_from_gmsh(mesh_path, tmp_dir, log=False, save=True)
+            ot.meshes_from_gmsh(mesh_path, tmp_dir, log=False, save=True)
 
-            model = ogs.Project(
+            model = ot.Project(
                 output_file=tmp_dir / "default.prj",
                 input_file=examples.prj_mechanics,
             )
@@ -50,7 +50,7 @@ class TestUtils:
             model.replace_text(4, xpath=".//integration_order")
             model.write_input()
             model.run_model(write_logs=False, args=f"-m {tmp_dir} -o {tmp_dir}")
-            ogs.MeshSeries(tmp_dir / "mesh_mesh_domain.xdmf").mesh(0)
+            ot.MeshSeries(tmp_dir / "mesh_mesh_domain.xdmf").mesh(0)
 
     @pytest.mark.parametrize(
         "ht",
@@ -150,7 +150,7 @@ class TestUtils:
         pvd = examples.load_meshseries_THM_2D_PVD()
         xdmf = examples.load_meshseries_HT_2D_XDMF()
         with pytest.raises(TypeError, match="Can only read"):
-            ogs.MeshSeries(__file__)
+            ot.MeshSeries(__file__)
 
         for ms in [pvd, xdmf]:
             try:
@@ -207,7 +207,7 @@ class TestUtils:
     def test_time_aggregate_mesh_dependent(self):
         "Test aggregation of mesh_dependent variable on meshseries."
         mesh_series = examples.load_meshseries_THM_2D_PVD()
-        prop = ogs.variables.dilatancy_alkan
+        prop = ot.variables.dilatancy_alkan
         agg_mesh = mesh_series.aggregate_over_time(prop, "max")
         assert not np.any(np.isnan(agg_mesh[prop.output_name + "_max"]))
         agg_mesh = mesh_series.time_of_max(prop)
@@ -248,22 +248,22 @@ class TestUtils:
         """Test creation of probe plots."""
         meshseries = examples.load_meshseries_THM_2D_PVD()
         points = meshseries.mesh(0).center
-        meshseries.plot_probe(points, ogs.variables.temperature)
+        meshseries.plot_probe(points, ot.variables.temperature)
         points = meshseries.mesh(0).points[[0, -1]]
-        meshseries.plot_probe(points, ogs.variables.temperature)
-        meshseries.plot_probe(points, ogs.variables.velocity)
-        meshseries.plot_probe(points, ogs.variables.stress)
-        meshseries.plot_probe(points, ogs.variables.stress.von_Mises)
+        meshseries.plot_probe(points, ot.variables.temperature)
+        meshseries.plot_probe(points, ot.variables.velocity)
+        meshseries.plot_probe(points, ot.variables.stress)
+        meshseries.plot_probe(points, ot.variables.stress.von_Mises)
         mesh_series = examples.load_meshseries_HT_2D_XDMF()
         points = mesh_series.mesh(0).center
-        meshseries.plot_probe(points, ogs.variables.temperature)
-        meshseries.plot_probe(points, ogs.variables.velocity)
+        meshseries.plot_probe(points, ot.variables.temperature)
+        meshseries.plot_probe(points, ot.variables.velocity)
 
     def test_diff_two_meshes(self):
         meshseries = examples.load_meshseries_THM_2D_PVD()
         mesh1 = meshseries.mesh(0)
         mesh2 = meshseries.mesh(-1)
-        mesh_diff = ogs.meshlib.difference(mesh1, mesh2, "temperature")
+        mesh_diff = ot.meshlib.difference(mesh1, mesh2, "temperature")
         # test, that no sampling occurs for equal topology
         np.testing.assert_array_equal(
             mesh_diff["temperature_difference"],
@@ -272,41 +272,41 @@ class TestUtils:
         # test same/different topology and scalar / vector variable
         for scaling in [1.0, 2.0]:
             for variable in ["temperature", "velocity"]:
-                mesh_diff = ogs.meshlib.difference(
+                mesh_diff = ot.meshlib.difference(
                     mesh1.scale(scaling), mesh2, variable
                 )
 
-        quad_tri_diff = ogs.meshlib.difference(
+        quad_tri_diff = ot.meshlib.difference(
             mesh1.triangulate(), mesh1, "temperature"
         )
-        quad_tri_diff_vals = ogs.variables.temperature.difference.transform(
+        quad_tri_diff_vals = ot.variables.temperature.difference.transform(
             quad_tri_diff
         )
         np.testing.assert_allclose(quad_tri_diff_vals, 0.0, atol=1e-12)
-        mesh_diff = ogs.meshlib.difference(
-            mesh1, mesh2, ogs.variables.temperature
+        mesh_diff = ot.meshlib.difference(
+            mesh1, mesh2, ot.variables.temperature
         )
         assert isinstance(mesh_diff, pv.UnstructuredGrid)
-        mesh_diff = ogs.meshlib.difference(mesh1, mesh2)
+        mesh_diff = ot.meshlib.difference(mesh1, mesh2)
 
     def test_diff_pairwise(self):
         n = 5
         meshseries = examples.load_meshseries_THM_2D_PVD()
         meshes1 = [meshseries.mesh(0)] * n
         meshes2 = [meshseries.mesh(-1)] * n
-        meshes_diff = ogs.meshlib.difference_pairwise(
-            meshes1, meshes2, ogs.variables.temperature
+        meshes_diff = ot.meshlib.difference_pairwise(
+            meshes1, meshes2, ot.variables.temperature
         )
         assert isinstance(meshes_diff, np.ndarray)
         assert len(meshes_diff) == n
 
-        meshes_diff = ogs.meshlib.difference_pairwise(meshes1, meshes2)
+        meshes_diff = ot.meshlib.difference_pairwise(meshes1, meshes2)
 
     def test_diff_matrix_single(self):
         meshseries = examples.load_meshseries_THM_2D_PVD()
         meshes1 = [meshseries.mesh(0), meshseries.mesh(-1)]
-        meshes_diff = ogs.meshlib.difference_matrix(
-            meshes1, variable=ogs.variables.temperature
+        meshes_diff = ot.meshlib.difference_matrix(
+            meshes1, variable=ot.variables.temperature
         )
         assert isinstance(meshes_diff, np.ndarray)
 
@@ -315,21 +315,21 @@ class TestUtils:
             len(meshes1),
         )
 
-        meshes_diff = ogs.meshlib.difference_matrix(meshes1)
+        meshes_diff = ot.meshlib.difference_matrix(meshes1)
 
     def test_diff_matrix_unequal(self):
         meshseries = examples.load_meshseries_THM_2D_PVD()
         meshes1 = [meshseries.mesh(0), meshseries.mesh(-1)]
         meshes2 = [meshseries.mesh(0), meshseries.mesh(-1), meshseries.mesh(-1)]
-        meshes_diff = ogs.meshlib.difference_matrix(
-            meshes1, meshes2, ogs.variables.temperature
+        meshes_diff = ot.meshlib.difference_matrix(
+            meshes1, meshes2, ot.variables.temperature
         )
         assert isinstance(meshes_diff, np.ndarray)
         assert meshes_diff.shape == (
             len(meshes1),
             len(meshes2),
         )
-        meshes_diff = ogs.meshlib.difference_matrix(meshes1, meshes2)
+        meshes_diff = ot.meshlib.difference_matrix(meshes1, meshes2)
 
     def test_depth_2D(self):
         mesh = examples.load_mesh_mechanics_2D()
@@ -340,7 +340,7 @@ class TestUtils:
         assert np.all(mesh["depth"] < -mesh.points[..., 1])
 
     def test_depth_3D(self):
-        mesh = ogs.Mesh(pv.SolidSphere(100, center=(0, 0, -101)))
+        mesh = ot.Mesh(pv.SolidSphere(100, center=(0, 0, -101)))
         mesh["depth"] = mesh.depth(use_coords=True)
         assert np.all(mesh["depth"] == -mesh.points[..., -1])
         mesh["depth"] = mesh.depth()
@@ -354,7 +354,7 @@ class TestUtils:
                 [100, -300, 6700],
             ]
         )
-        profile_points = ogs.meshlib.interp_points(profile, resolution=100)
+        profile_points = ot.meshlib.interp_points(profile, resolution=100)
         assert isinstance(profile_points, np.ndarray)
         # Check first point
         assert (profile_points[0, :] == profile[0, :]).all()
@@ -365,15 +365,15 @@ class TestUtils:
 
     def test_distance_in_segments(self):
         profile = np.array([[0, 0, 0], [1, 0, 0], [2, 0, 0]])
-        profile_points = ogs.meshlib.interp_points(profile, resolution=9)
-        dist_in_seg = ogs.meshlib.distance_in_segments(profile, profile_points)
+        profile_points = ot.meshlib.interp_points(profile, resolution=9)
+        dist_in_seg = ot.meshlib.distance_in_segments(profile, profile_points)
         assert len(np.where(dist_in_seg == 0)[0]) == 2
         assert len(np.where(dist_in_seg == 1)[0]) == 1
 
     def test_distance_in_profile(self):
         profile = np.array([[0, 0, 0], [1, 0, 0], [2, 0, 0]])
-        profile_points = ogs.meshlib.interp_points(profile, resolution=9)
-        dist_in_seg = ogs.meshlib.distance_in_profile(profile_points)
+        profile_points = ot.meshlib.interp_points(profile, resolution=9)
+        dist_in_seg = ot.meshlib.distance_in_profile(profile_points)
         # Check if distance is increasing
         assert np.all(np.diff(dist_in_seg) > 0)
         # Check if distances at the beginning and end of profile are correct
@@ -383,7 +383,7 @@ class TestUtils:
     def test_sample_over_polyline_single_segment(self):
         ms = examples.load_meshseries_HT_2D_XDMF()
         profile = np.array([[4, 2, 0], [4, 18, 0]])
-        ms_sp, _ = ogs.meshlib.sample_polyline(
+        ms_sp, _ = ot.meshlib.sample_polyline(
             ms.mesh(-1),
             ["pressure", "temperature"],
             profile,
@@ -407,13 +407,13 @@ class TestUtils:
                 [910, -590, 6700],
             ]
         )
-        ms_sp, _ = ogs.meshlib.sample_polyline(
+        ms_sp, _ = ot.meshlib.sample_polyline(
             ms.mesh(1),
-            ogs.variables.temperature,
+            ot.variables.temperature,
             profile,
             resolution=10,
         )
-        data = ms_sp[ogs.variables.temperature.data_name].to_numpy()
+        data = ms_sp[ot.variables.temperature.data_name].to_numpy()
         assert not np.any(np.isnan(data))
         assert (np.abs(data) > np.zeros_like(data)).all()
         # output should be in Celsius
@@ -423,7 +423,7 @@ class TestUtils:
     def test_sample_over_polyline_single_segment_vec_prop(self):
         ms = examples.load_meshseries_HT_2D_XDMF()
         profile = np.array([[4, 2, 0], [4, 18, 0]])
-        ms_sp, _ = ogs.meshlib.sample_polyline(
+        ms_sp, _ = ot.meshlib.sample_polyline(
             ms.mesh(-1),
             "darcy_velocity",
             profile,
@@ -458,9 +458,9 @@ class TestUtils:
 
         tmp_path = Path(mkdtemp())
         mesh_path = Path(tmp_path) / "mesh.msh"
-        sigma_ip = ogs.variables.stress.replace(data_name="sigma_ip")
+        sigma_ip = ot.variables.stress.replace(data_name="sigma_ip")
 
-        ogs.meshlib.rect(
+        ot.meshlib.rect(
             n_edge_cells=6,
             n_layers=2,
             structured_grid=quads,
@@ -469,17 +469,17 @@ class TestUtils:
             mixed_elements=mixed,
             jiggle=0.01,
         )
-        meshes = ogs.meshes_from_gmsh(mesh_path, log=False)
+        meshes = ot.meshes_from_gmsh(mesh_path, log=False)
         for name, mesh in meshes.items():
             pv.save_meshio(Path(tmp_path, name + ".vtu"), mesh)
-        model = ogs.Project(
+        model = ot.Project(
             output_file=tmp_path / "default.prj",
             input_file=examples.prj_mechanics,
         )
         model.replace_text(intpt_order, xpath=".//integration_order")
         model.write_input()
         model.run_model(write_logs=True, args=f"-m {tmp_path} -o {tmp_path}")
-        meshseries = ogs.MeshSeries(tmp_path / "mesh.pvd")
+        meshseries = ot.MeshSeries(tmp_path / "mesh.pvd")
         int_pts = meshseries.mesh(-1).to_ip_point_cloud()
         ip_ms = meshseries.ip_tesselated()
         ip_mesh = ip_ms.mesh(-1)
@@ -494,23 +494,23 @@ class TestUtils:
         )
 
     def test_reader(self):
-        assert isinstance(examples.load_meshseries_THM_2D_PVD(), ogs.MeshSeries)
-        assert isinstance(ogs.MeshSeries(examples.elder_xdmf), ogs.MeshSeries)
-        assert isinstance(ogs.Mesh.read(examples.mechanics_vtu), ogs.Mesh)
-        assert isinstance(ogs.Mesh.read(examples.test_shapefile), ogs.Mesh)
+        assert isinstance(examples.load_meshseries_THM_2D_PVD(), ot.MeshSeries)
+        assert isinstance(ot.MeshSeries(examples.elder_xdmf), ot.MeshSeries)
+        assert isinstance(ot.Mesh.read(examples.mechanics_vtu), ot.Mesh)
+        assert isinstance(ot.Mesh.read(examples.test_shapefile), ot.Mesh)
 
     def test_xdmf_quadratic(self):
         "Test reading of quadratic elements in xdmf."
 
         tmp_path = Path(mkdtemp())
         msh_path = Path(tmp_path) / "mesh.msh"
-        ogs.meshlib.rect(
+        ot.meshlib.rect(
             n_edge_cells=6, structured_grid=False, order=2, out_name=msh_path
         )
-        meshes = ogs.meshes_from_gmsh(msh_path, log=False)
+        meshes = ot.meshes_from_gmsh(msh_path, log=False)
         for name, mesh in meshes.items():
             pv.save_meshio(Path(tmp_path, name + ".vtu"), mesh)
-        model = ogs.Project(
+        model = ot.Project(
             input_file=examples.prj_mechanics,
             output_file=tmp_path / "default.prj",
         )
@@ -518,22 +518,22 @@ class TestUtils:
         model.replace_text("XDMF", xpath="./time_loop/output/type")
         model.write_input()
         model.run_model(write_logs=True, args=f"-m {tmp_path} -o {tmp_path}")
-        mesh = ogs.MeshSeries(tmp_path / "mesh_domain.xdmf").mesh(-1)
-        assert not np.any(np.isnan(ogs.variables.stress.transform(mesh)))
+        mesh = ot.MeshSeries(tmp_path / "mesh_domain.xdmf").mesh(-1)
+        assert not np.any(np.isnan(ot.variables.stress.transform(mesh)))
 
     def test_remesh_with_tri(self):
         mesh = examples.load_meshseries_THM_2D_PVD().mesh(1)
         temp_dir = Path(mkdtemp())
         msh_path = temp_dir / "tri_mesh.msh"
-        ogs.meshlib.gmsh_meshing.remesh_with_triangles(mesh, msh_path)
+        ot.meshlib.gmsh_meshing.remesh_with_triangles(mesh, msh_path)
         assert len(
-            ogs.meshes_from_gmsh(msh_path, reindex=False, log=False).items()
+            ot.meshes_from_gmsh(msh_path, reindex=False, log=False).items()
         ) == 1 + len(np.unique(mesh["MaterialIDs"]))
         # boundaries are not assigned a physical tag in remesh_with_trinagles
 
     def test_indexing(self):
         ms = examples.load_meshseries_HT_2D_XDMF()
-        assert isinstance(ms[1], ogs.Mesh)
+        assert isinstance(ms[1], ot.Mesh)
 
     def test_slice(self):
         ms = examples.load_meshseries_HT_2D_XDMF()
diff --git a/tests/test_plot.py b/tests/test_plot.py
index 4ca2e3fc4e7798c804caf0321590447a71821ee3..58a0529880cd6758a41167039a4f2839e0bc97e6 100644
--- a/tests/test_plot.py
+++ b/tests/test_plot.py
@@ -8,7 +8,7 @@ import numpy as np
 import pytest
 from pyvista import examples as pv_examples
 
-import ogstools as ogs
+import ogstools as ot
 from ogstools import examples
 from ogstools.plot import utils
 
@@ -33,30 +33,30 @@ class TestPlotting:
     def test_levels(self):
         """Test levels calculation."""
         assert_allclose(
-            ogs.plot.compute_levels(0.5, 10.1, 10), [0.5, *range(1, 11), 10.1]
+            ot.plot.compute_levels(0.5, 10.1, 10), [0.5, *range(1, 11), 10.1]
         )
         assert_allclose(
-            ogs.plot.compute_levels(293, 350, 10), [293, *range(295, 355, 5)]
+            ot.plot.compute_levels(293, 350, 10), [293, *range(295, 355, 5)]
         )
         assert_allclose(
-            ogs.plot.compute_levels(1e-3, 1.2, 5),
+            ot.plot.compute_levels(1e-3, 1.2, 5),
             [1e-3, *np.arange(0.2, 1.4, 0.2)],
         )
         assert_allclose(
-            ogs.plot.compute_levels(1e5, 9e6, 20),
+            ot.plot.compute_levels(1e5, 9e6, 20),
             [1e5, *np.arange(5e5, 9.5e6, 5e5)],
         )
         assert_allclose(
-            ogs.plot.compute_levels(1, 40, 20), [1, *range(2, 42, 2)]
+            ot.plot.compute_levels(1, 40, 20), [1, *range(2, 42, 2)]
         )
-        assert_allclose(ogs.plot.compute_levels(0.0, 0.0, 10), [0.0, 0.0])
-        assert_allclose(ogs.plot.compute_levels(1e9, 1e9, 10), [1e9, 1e9])
+        assert_allclose(ot.plot.compute_levels(0.0, 0.0, 10), [0.0, 0.0])
+        assert_allclose(ot.plot.compute_levels(1e9, 1e9, 10), [1e9, 1e9])
 
     def test_ticklabels(self):
         def compare(lower, upper, precision, ref_labels, ref_offset=None):
-            labels, offset = ogs.plot.contourplots.get_ticklabels(
+            labels, offset = ot.plot.contourplots.get_ticklabels(
                 np.asarray(
-                    ogs.plot.compute_levels(lower, upper, n_ticks=precision)
+                    ot.plot.compute_levels(lower, upper, n_ticks=precision)
                 )
             )
             assert np.all(labels == ref_labels)
@@ -108,24 +108,24 @@ class TestPlotting:
     def test_missing_data(self):
         """Test missing data in mesh."""
         mesh = pv_examples.load_uniform()
-        pytest.raises(KeyError, ogs.plot.contourf, mesh, "missing_data")
+        pytest.raises(KeyError, ot.plot.contourf, mesh, "missing_data")
 
     def test_plot_2_d(self):
         """Test creation of 2D plots."""
-        ogs.plot.setup.reset()
-        ogs.plot.setup.material_names = {
+        ot.plot.setup.reset()
+        ot.plot.setup.material_names = {
             i + 1: f"Layer {i+1}" for i in range(26)
         }
         meshseries = examples.load_meshseries_THM_2D_PVD()
         mesh = meshseries.mesh(1)
-        mesh.plot_contourf(ogs.variables.material_id)
-        mesh.plot_contourf(ogs.variables.temperature)
-        mesh.plot_contourf(ogs.variables.Scalar("pressure_active"))
-        ogs.plot.contourf(
-            mesh.threshold((1, 3), "MaterialIDs"), ogs.variables.velocity
+        mesh.plot_contourf(ot.variables.material_id)
+        mesh.plot_contourf(ot.variables.temperature)
+        mesh.plot_contourf(ot.variables.Scalar("pressure_active"))
+        ot.plot.contourf(
+            mesh.threshold((1, 3), "MaterialIDs"), ot.variables.velocity
         )
-        fig = mesh.plot_contourf(ogs.variables.displacement[0])
-        ogs.plot.shape_on_top(
+        fig = mesh.plot_contourf(ot.variables.displacement[0])
+        ot.plot.shape_on_top(
             fig.axes[0], mesh, lambda x: min(max(0, 0.1 * (x - 3)), 100)
         )
         plt.close()
@@ -139,10 +139,10 @@ class TestPlotting:
             "temperature_difference"
         )
         for prop in [
-            ogs.variables.temperature,
-            ogs.variables.displacement,
-            ogs.variables.stress,
-            ogs.variables.stress.von_Mises,
+            ot.variables.temperature,
+            ot.variables.displacement,
+            ot.variables.stress,
+            ot.variables.stress.von_Mises,
         ]:
             mesh1.difference(mesh0, prop).plot_contourf(prop)
         plt.close()
@@ -152,27 +152,27 @@ class TestPlotting:
         meshseries = examples.load_meshseries_THM_2D_PVD()
         fig, ax = plt.subplots(3, 1, figsize=(40, 30))
         meshseries.mesh(0).plot_contourf(
-            ogs.variables.temperature, fig=fig, ax=ax[0]
+            ot.variables.temperature, fig=fig, ax=ax[0]
         )
         meshseries.mesh(1).plot_contourf(
-            ogs.variables.temperature, fig=fig, ax=ax[1]
+            ot.variables.temperature, fig=fig, ax=ax[1]
         )
         diff_mesh = meshseries.mesh(0).difference(
-            meshseries.mesh(1), ogs.variables.temperature
+            meshseries.mesh(1), ot.variables.temperature
         )
-        diff_mesh.plot_contourf(ogs.variables.temperature, fig=fig, ax=ax[2])
+        diff_mesh.plot_contourf(ot.variables.temperature, fig=fig, ax=ax[2])
         plt.close()
 
     def test_user_defined_ax_two_variables(self):
         """Test creating plot with subfigures and user provided ax with different values plotted"""
         meshseries = examples.load_meshseries_THM_2D_PVD()
-        ogs.plot.setup.combined_colorbar = False
+        ot.plot.setup.combined_colorbar = False
         fig, ax = plt.subplots(2, 1, figsize=(40, 20))
         meshseries.mesh(0).plot_contourf(
-            ogs.variables.temperature, fig=fig, ax=ax[0]
+            ot.variables.temperature, fig=fig, ax=ax[0]
         )
         meshseries.mesh(1).plot_contourf(
-            ogs.variables.displacement, fig=fig, ax=ax[1]
+            ot.variables.displacement, fig=fig, ax=ax[1]
         )
         fig.suptitle("Test user defined ax")
         plt.close()
@@ -180,11 +180,11 @@ class TestPlotting:
     def test_user_defined_fig(self):
         """Test creating plot with subfigures and user provided fig"""
         meshseries = examples.load_meshseries_THM_2D_PVD()
-        ogs.plot.setup.combined_colorbar = False
+        ot.plot.setup.combined_colorbar = False
         fig, ax = plt.subplots(2, 1, figsize=(40, 20))
-        ogs.plot.contourf(
+        ot.plot.contourf(
             [meshseries.mesh(0), meshseries.mesh(1)],
-            ogs.variables.temperature,
+            ot.variables.temperature,
             fig=fig,
         )
         fig.suptitle("Test user defined fig")
@@ -203,7 +203,7 @@ class TestPlotting:
         meshseries = examples.load_meshseries_THM_2D_PVD()
         timevalues = np.linspace(0, meshseries.timevalues()[-1], num=3)
         anim = meshseries.animate(
-            ogs.variables.temperature,
+            ot.variables.temperature,
             timevalues,
             mesh_func=lambda mesh: mesh.clip("x"),
             plot_func=lambda ax, t: ax.set_title(str(t)),
@@ -215,7 +215,7 @@ class TestPlotting:
         """Test saving of an animation."""
         meshseries = examples.load_meshseries_THM_2D_PVD()
         timevalues = np.linspace(0, meshseries.timevalues()[-1], num=3)
-        anim = meshseries.animate(ogs.variables.temperature, timevalues)
+        anim = meshseries.animate(ot.variables.temperature, timevalues)
         if not utils.save_animation(anim, mkstemp()[1], 5):
             pytest.skip("Saving animation failed.")
         plt.close()
@@ -223,10 +223,10 @@ class TestPlotting:
     def test_plot_3_d(self):
         """Test creation of slice plots for 3D mesh."""
         mesh = pv_examples.load_uniform()
-        ogs.plot.contourf(mesh.slice((1, 1, 0)), "Spatial Point Data")
+        ot.plot.contourf(mesh.slice((1, 1, 0)), "Spatial Point Data")
         meshes = np.reshape(mesh.slice_along_axis(4, "x"), (2, 2))
-        ogs.plot.contourf(meshes, "Spatial Point Data")
-        ogs.plot.contourf(mesh.slice([1, -2, 0]), "Spatial Point Data")
+        ot.plot.contourf(meshes, "Spatial Point Data")
+        ot.plot.contourf(mesh.slice([1, -2, 0]), "Spatial Point Data")
         plt.close()
 
     def test_streamlines(self):
@@ -238,8 +238,8 @@ class TestPlotting:
         for axis in ["x", "y", "z", [1, 1, 0]]:
             ax: plt.Axes
             fig, ax = plt.subplots()
-            i_grid, j_grid, u, v, lw = ogs.plot.vectorplots._vectorfield(
-                mesh.slice(axis), ogs.variables.velocity
+            i_grid, j_grid, u, v, lw = ot.plot.vectorplots._vectorfield(
+                mesh.slice(axis), ot.variables.velocity
             )
             for vals in [i_grid, j_grid, u, v, lw]:
                 assert not np.all(np.isnan(vals))
@@ -251,13 +251,13 @@ class TestPlotting:
     def test_xdmf(self):
         """Test creation of 2D plots from xdmf data."""
         mesh = examples.load_meshseries_CT_2D_XDMF().mesh(0)
-        mesh.plot_contourf(ogs.variables.saturation)
+        mesh.plot_contourf(ot.variables.saturation)
         plt.close()
 
     def test_xdmf_with_slices(self):
         """Test creation of 2D plots from xdmf data."""
         mesh = examples.load_meshseries_HT_2D_XDMF().mesh(0)
-        mesh.plot_contourf(ogs.variables.pressure)
+        mesh.plot_contourf(ot.variables.pressure)
         plt.close()
 
     def test_lineplot(self):
@@ -287,7 +287,7 @@ class TestPlotting:
         ms_CT = examples.load_meshseries_CT_2D_XDMF()
         profile_CT = np.array([[47.0, 1.17, 72.0], [-4.5, 1.17, -59.0]])
         fig, ax = ms_CT.mesh(11).plot_linesample_contourf(
-            ogs.variables.saturation,
+            ot.variables.saturation,
             profile_CT,
             resolution=100,
             plot_nodal_pts=True,