diff --git a/ProcessLib/SmallDeformation/Tests.cmake b/ProcessLib/SmallDeformation/Tests.cmake
index e3bd7399531c78be6dfc027a1cea5937524ad70f..ef59fa740a51741f90557a3207348b62504baa64 100644
--- a/ProcessLib/SmallDeformation/Tests.cmake
+++ b/ProcessLib/SmallDeformation/Tests.cmake
@@ -352,7 +352,7 @@ AddTest(
 
 if(NOT OGS_USE_PETSC)
     NotebookTest(NOTEBOOKFILE Mechanics/CooksMembrane/CooksMembraneBbar.py RUNTIME 1)
-    NotebookTest(NOTEBOOKFILE Mechanics/Linear/SimpleMechanics.ipynb RUNTIME 5)
+    NotebookTest(NOTEBOOKFILE Mechanics/Linear/SimpleMechanics.py RUNTIME 5)
     NotebookTest(NOTEBOOKFILE Mechanics/EvaluatingBbarWithSimpleExamples/evaluating_bbbar_with_simple_examples.py RUNTIME 5)
     NotebookTest(NOTEBOOKFILE Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole.md RUNTIME 15)
     if(NOT WIN32)
diff --git a/Tests/Data/Mechanics/CooksMembrane/CooksMembraneBbar.py b/Tests/Data/Mechanics/CooksMembrane/CooksMembraneBbar.py
index 2f49571c4579f1caa5646afff5dfbe64bfde5ae6..b07f125089723441ae9b8c5579153a0028d2171b 100644
--- a/Tests/Data/Mechanics/CooksMembrane/CooksMembraneBbar.py
+++ b/Tests/Data/Mechanics/CooksMembrane/CooksMembraneBbar.py
@@ -66,9 +66,9 @@ from pathlib import Path
 import matplotlib.pyplot as plt
 import matplotlib.tri as tri
 import numpy as np
+import ogstools as ot
 import pyvista as pv
 import vtuIO
-from ogs6py.ogs import OGS
 
 out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
 if not out_dir.exists():
@@ -98,7 +98,9 @@ def get_top_uy(pvd_file_name):
 
 # %%
 def run_single_test(mesh_name, output_prefix, use_bbar="false"):
-    model = OGS(INPUT_FILE="CooksMembrane.prj", PROJECT_FILE=f"{out_dir}/modified.prj")
+    model = ot.Project(
+        input_file="CooksMembrane.prj", output_file=f"{out_dir}/modified.prj"
+    )
     model.replace_text(mesh_name, xpath="./mesh")
     model.replace_text(use_bbar, xpath="./processes/process/use_b_bar")
     model.replace_text(output_prefix, xpath="./time_loop/output/prefix")
@@ -273,258 +275,3 @@ for nedge, output_prefix in zip(nedges, output_prefices):
 
 # %% [markdown]
 # The contour plots show that even with the coarsest mesh, the B bar method still gives reasonable stress result.
-# ---
-# jupyter:
-#   jupytext:
-#     text_representation:
-#       extension: .py
-#       format_name: percent
-#       format_version: '1.3'
-#       jupytext_version: 1.14.5
-#   kernelspec:
-#     display_name: Python 3 (ipykernel)
-#     language: python
-#     name: python3
-# ---
-
-# %% [raw]
-# +++
-# title = "Cook's membrane example"
-# date = "2024-06-11"
-# author = "Wenqing Wang"
-# image = "figures/cooks_membrane.png"
-# web_subsection = "small-deformations"
-# weight = 3
-# +++
-
-# %% [markdown]
-# $$
-# \newcommand{\B}{\text{B}}
-# \newcommand{\F}{\text{F}}
-# \newcommand{\I}{\mathbf I}
-# \newcommand{\intD}[1]{\int_{\Omega_e}#1\mathrm{d}\Omega}
-# $$
-#
-# # Cook's membrane example for nearly icompressible solid
-#
-# ## B bar method
-#  Considering a strain decomposition: $\mathbf\epsilon = \underbrace{\mathbf\epsilon- \frac{1}{3}(\epsilon:\mathbf I)}_{\text{deviatoric}}\I +  \underbrace{\frac{1}{3}(\epsilon:\mathbf I)}_{\text{dilatational}} \I$.
-# 	   The idea of the B bar method is to use another quadrature rule to interpolate the dilatational part, which leads to a modified B matrix [1]:
-# $$
-# 	     \bar\B = \underbrace{\B - \B^{\text{dil}}}_{\text{original B elements}}+ \underbrace{{\bar\B}^{\text{dil}}}_{\text{by another quadrature rule} }
-# $$
-# There are several methods to form ${\bar\B}^{\text{dil}}$  such as selective integration, generalization of the mean-dilatation formulation. In the current OGS, we use the latter, which reads
-# $$
-# 		  {\bar\B}^{\text{dil}} = \frac{\intD{\B^{\text{dil}}(\xi)}}{\intD{}}
-# $$
-#
-# ## Example
-# To verify the implementation of the B bar method, the so called Cook's membrane is used as a benchmark.
-# Illustrated in the following figure, this example simulates a tapered
-# and swept panel of unit thickness. The left edge is clamped and the right edge is applied with  a distributed shearing load $F$ = 100 N/mm. The plane strain condition is considered. This numerical model is exactly the same as that is presented in the paper by T. Elguedj et al [1,2].
-#
-# <img src="figures/cooks_membrane.png" alt="Cook's membrane" width="320" height="320" />
-#
-# ## Reference
-#
-# [1] T.J.R. Hughes (1980). Generalization of selective integration procedures to anisotropic and nonlinear media. International Journal for Numerical Methods in Engineering, 15(9), 1413-1418.
-#
-# [2] T. Elguedj, Y. Bazilevs, V.M. Calo, T.J.R. Hughes (2008),
-#  $\bar\B$ and $\bar\F$ projection methods for nearly incompressible linear and non-linear elasticity and plasticity using higher-order NURBS elements, Computer Methods in Applied Mechanics and Engineering, 197(33--40), 2732-2762.
-#
-
-# %%
-import os
-import xml.etree.ElementTree as ET
-from pathlib import Path
-
-import matplotlib.pyplot as plt
-import matplotlib.tri as tri
-import numpy as np
-import ogstools as ot
-import pyvista as pv
-import vtuIO
-
-# %%
-out_dir = Path(os.environ.get("OT_TESTRUNNER_OUT_DIR", "_out"))
-if not out_dir.exists():
-    out_dir.mkdir(parents=True)
-
-
-# %%
-def get_last_vtu_file_name(pvd_file_name):
-    tree = ET.parse(Path(out_dir) / pvd_file_name)
-    root = tree.getroot()
-    # Get the last DataSet tag
-    last_dataset = root.findall(".//DataSet")[-1]
-
-    # Get the 'file' attribute of the last DataSet tag
-    file_attribute = last_dataset.attrib["file"]
-    return f"{out_dir}/" + file_attribute
-
-
-def get_top_uy(pvd_file_name):
-    top_point = (48.0e-3, 60.0e-3, 0)
-    file_name = get_last_vtu_file_name(pvd_file_name)
-    mesh = pv.read(file_name)
-    p_id = mesh.find_closest_point(top_point)
-    u = mesh.point_data["displacement"][p_id]
-    return u[1]
-
-
-# %%
-def run_single_test(mesh_name, output_prefix, use_bbar="false"):
-    model = ot.Project(
-        input_file="CooksMembrane.prj", output_file=f"{out_dir}/modified.prj"
-    )
-    model.replace_text(mesh_name, xpath="./mesh")
-    model.replace_text(use_bbar, xpath="./processes/process/use_b_bar")
-    model.replace_text(output_prefix, xpath="./time_loop/output/prefix")
-    model.replace_text(
-        "BiCGSTAB", xpath="./linear_solvers/linear_solver/eigen/solver_type"
-    )
-    model.replace_text("ILUT", xpath="./linear_solvers/linear_solver/eigen/precon_type")
-    vtu_file_name = output_prefix + "_ts_1_t_1.000000.vtu"
-    model.replace_text(vtu_file_name, xpath="./test_definition/vtkdiff[1]/file")
-    model.replace_text(vtu_file_name, xpath="./test_definition/vtkdiff[2]/file")
-    model.replace_text(vtu_file_name, xpath="./test_definition/vtkdiff[3]/file")
-
-    model.write_input()
-
-    # Run OGS
-    model.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir} -m .")
-
-    # Get uy at the top
-    return get_top_uy(output_prefix + ".pvd")
-
-
-# %%
-mesh_names = [
-    "mesh.vtu",
-    "mesh_n10.vtu",
-    "mesh_n15.vtu",
-    "mesh_n20.vtu",
-    "mesh_n25.vtu",
-    "mesh_n30.vtu",
-]
-output_prefices_non_bbar = [
-    "cooks_membrane_sd_edge_div_4_non_bbar",
-    "cooks_membrane_sd_refined_mesh_10_non_bbar",
-    "cooks_membrane_sd_refined_mesh_15_non_bbar",
-    "cooks_membrane_sd_refined_mesh_20_non_bbar",
-    "cooks_membrane_sd_refined_mesh_25_non_bbar",
-    "cooks_membrane_sd_refined_mesh_30_non_bbar",
-]
-
-uys_at_top_non_bbar = []
-for mesh_name, output_prefix in zip(mesh_names, output_prefices_non_bbar):
-    uy_at_top = run_single_test(mesh_name, output_prefix)
-    uys_at_top_non_bbar.append(uy_at_top)
-
-print(uys_at_top_non_bbar)
-
-# %%
-output_prefices = [
-    "cooks_membrane_sd_edge_div_4",
-    "cooks_membrane_sd_refined_mesh_10",
-    "cooks_membrane_sd_refined_mesh_15",
-    "cooks_membrane_sd_refined_mesh_20",
-    "cooks_membrane_sd_refined_mesh_25",
-    "cooks_membrane_sd_refined_mesh_30",
-]
-
-uys_at_top_bbar = []
-for mesh_name, output_prefix in zip(mesh_names, output_prefices):
-    uy_at_top = run_single_test(mesh_name, output_prefix, "true")
-    uys_at_top_bbar.append(uy_at_top)
-
-print(uys_at_top_bbar)
-
-# %%
-ne = [4, 10, 15, 20, 25, 30]
-
-
-def plot_data(ne, u_y_bbar, uy_non_bbar, file_name=""):
-    # Plotting
-    plt.rcParams["figure.figsize"] = [5, 5]
-
-    if len(u_y_bbar) != 0:
-        plt.plot(
-            ne, np.array(u_y_bbar) * 1e3, marker="o", linestyle="dashed", label="B bar"
-        )
-    if len(uy_non_bbar) != 0:
-        plt.plot(
-            ne,
-            np.array(uy_non_bbar) * 1e3,
-            marker="x",
-            linestyle="dashed",
-            label="non B bar",
-        )
-
-    plt.xlabel("Number of elements per side")
-    plt.ylabel("Top right corner displacement /mm")
-    plt.legend()
-
-    plt.tight_layout()
-    if file_name != "":
-        plt.savefig(file_name)
-    plt.show()
-
-
-# %% [markdown]
-# ## Result
-#
-# ### Vertical diplacement at the top point
-#
-# The following figure shows that the convergence of the solutions obtained by using the B bar method follows the one presented in the paper by T. Elguedj et al [1]. However, the results obtained without the B bar method are quit far from the converged solution with the finest mesh.
-
-# %%
-plot_data(ne, uys_at_top_bbar, uys_at_top_non_bbar, "b_bar_linear.png")
-
-# %% [markdown]
-# ### Contour plot
-
-# %%
-nedges = ["4", "10", "15", "20", "25", "30"]
-
-
-def contour_plot(pvd_file_name, title):
-    file_name = get_last_vtu_file_name(pvd_file_name)
-    m_plot = vtuIO.VTUIO(file_name, dim=2)
-    triang = tri.Triangulation(m_plot.points[:, 0], m_plot.points[:, 1])
-    triang = tri.Triangulation(m_plot.points[:, 0], m_plot.points[:, 1])
-    s_plot = m_plot.get_point_field("sigma")
-    s_trace = s_plot[:, 0] + s_plot[:, 1] + s_plot[:, 2]
-    u_plot = m_plot.get_point_field("displacement")
-
-    fig, ax = plt.subplots(ncols=2, figsize=(8, 3))
-    ax[0].set_title(title, loc="left", y=1.12)
-    plt.subplots_adjust(wspace=0.5)
-
-    contour_stress = ax[0].tricontourf(triang, s_trace, cmap="viridis")
-    contour_displacement = ax[1].tricontourf(triang, u_plot[:, 1], cmap="gist_rainbow")
-    fig.colorbar(contour_stress, ax=ax[0], label="Stress trace / MPa")
-    fig.colorbar(contour_displacement, ax=ax[1], label="Dispplacement / m")
-    fig.tight_layout()
-    plt.savefig(pvd_file_name + ".png")
-    plt.show()
-
-
-# %% [markdown]
-# #### Results obtained without the B bar method:
-
-# %%
-for nedge, output_prefix in zip(nedges, output_prefices_non_bbar):
-    contour_plot(output_prefix + ".pvd", "Number of elements per side: " + nedge)
-
-# %% [markdown]
-# #### Results obtained with the B bar method:
-
-# %%
-for nedge, output_prefix in zip(nedges, output_prefices):
-    contour_plot(output_prefix + ".pvd", "Number of elements per side: " + nedge)
-
-# %% [markdown]
-# The contour plots show that even with the coarsest mesh, the B bar method still gives reasonable stress result.
-
-# %%
diff --git a/Tests/Data/Mechanics/EvaluatingBbarWithSimpleExamples/evaluating_bbbar_with_simple_examples.py b/Tests/Data/Mechanics/EvaluatingBbarWithSimpleExamples/evaluating_bbbar_with_simple_examples.py
index d3b5fc8c171192ce933f87af6883e17cb19317e4..585deb3a6a13818d22a6e7eaad2d5f9f0f9d84f5 100644
--- a/Tests/Data/Mechanics/EvaluatingBbarWithSimpleExamples/evaluating_bbbar_with_simple_examples.py
+++ b/Tests/Data/Mechanics/EvaluatingBbarWithSimpleExamples/evaluating_bbbar_with_simple_examples.py
@@ -54,9 +54,9 @@ from pathlib import Path
 import matplotlib.pyplot as plt
 import matplotlib.tri as tri
 import numpy as np
+import ogstools as ot
 import pyvista as pv
 import vtuIO
-from ogs6py.ogs import OGS
 
 out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
 if not out_dir.exists():
@@ -138,8 +138,8 @@ class SingleSimulationModel:
         is_homogeneous_mode=False,
         poisson_ratio=0.2,
     ):
-        self.model = OGS(
-            INPUT_FILE="simple_b_bar_test.prj", PROJECT_FILE=f"{out_dir}/modified.prj"
+        self.model = ot.Project(
+            input_file="simple_b_bar_test.prj", output_file=f"{out_dir}/modified.prj"
         )
         self.model.replace_text(use_bbar, xpath="./processes/process/use_b_bar")
         self.model.replace_text(output_prefix, xpath="./time_loop/output/prefix")
diff --git a/Tests/Data/Notebooks/README.md b/Tests/Data/Notebooks/README.md
index e670d904838c1937c52574d9eef439b4453587a1..0614faf28b6cd74b763451477670ef37b63df746 100644
--- a/Tests/Data/Notebooks/README.md
+++ b/Tests/Data/Notebooks/README.md
@@ -21,7 +21,7 @@ nbdime config-git --enable --global
 ```bash
 # in virtualenv:
 pip install -r requirements.txt
-PATH=~/code/ogs/build/release/bin:$PATH python testrunner.py --out _out SimpleMechanics.ipynb
+PATH=~/code/ogs/build/release/bin:$PATH python testrunner.py --out _out SimpleMechanics.py
 ```
 
 In Notebook do checks with `assert False` or raise an exception, e.g. `raise SystemExit()`, on failure.
diff --git a/web/content/docs/devguide/documentation/jupyter-docs/index.md b/web/content/docs/devguide/documentation/jupyter-docs/index.md
index 7dd61455ff3ee0e1dc75f85da3d91266f29ecff7..c95021c5509e5c2fd3183d9150d619547b8fa63a 100644
--- a/web/content/docs/devguide/documentation/jupyter-docs/index.md
+++ b/web/content/docs/devguide/documentation/jupyter-docs/index.md
@@ -78,11 +78,11 @@ These notebooks are part of the regular CI testing. Please try to keep the noteb
 
 ### Create a new notebook
 
-Create a new notebook file in `Tests/Data` (if it should appear in the benchmark gallery) or in `web/content/docs` (e.g. for tutorials). Create it as a regular Markdown-file with Python code blocks. The notebook execution and conversion is done via [Jupytext](https://jupytext.readthedocs.io/en/latest). See examples:
+Create a new notebook file in `Tests/Data` (if it should appear in the benchmark gallery) or in `web/content/docs` (e.g. for tutorials). Create it as a regular Python-file with Python code blocks and Markdown blocks for text. The notebook execution and conversion is done via [Jupytext](https://jupytext.readthedocs.io/en/latest). See examples:
 
+- [SimpleMechanics.py](https://gitlab.opengeosys.org/ogs/ogs/-/blob/master/Tests/Data/Mechanics/Linear/SimpleMechanics.py) (notebooks written as regular `.py` files are preferred but Markdown-based or `.ipynb` files notebooks are also possible)
 - [Linear_Disc_with_hole.md](https://gitlab.opengeosys.org/ogs/ogs/-/blob/master/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole.md) (Jupytext-based benchmark notebook in Markdown)
 - [notebook-bhe_meshing.md](https://gitlab.opengeosys.org/ogs/ogs/-/blob/master/web/content/docs/tutorials/bhe_meshing/notebook-bhe_meshing.md) (Jupytext-based tutorial notebook in Markdown)
-- [SimpleMechanics.ipynb](https://gitlab.opengeosys.org/ogs/ogs/-/blob/master/Tests/Data/Mechanics/Linear/SimpleMechanics.ipynb) (regular `.ipynb`-notebook are also possible but Markdown-based notebooks are preferred)
 
 ### Add web meta information
 
@@ -103,7 +103,7 @@ web_subsection = "small-deformations" # required for notebooks in Tests/Data onl
 - Frontmatter needs to be in [TOML](https://toml.io)-format.
 - For notebooks describing benchmarks `web_subsection` needs to be set to a sub-folder in [web/content/docs/benchmarks](https://gitlab.opengeosys.org/ogs/ogs/-/tree/master/web/content/docs/benchmarks) (if not set the notebook page will not be linked from navigation bar / benchmark gallery on the web page).
 - If you edit a Markdown-based notebook with Jupyter and the Jupytext extension please don't add the two newlines but make sure that the frontmatter has its own cell (not mixed with markdown content).
-- For (deprecated) `.ipynb`-based notebooks the frontmatter needs to be given in the first cell. See existing notebooks (e.g. [SimpleMechanics.ipynb](https://gitlab.opengeosys.org/ogs/ogs/-/blob/master/Tests/Data/Mechanics/Linear/SimpleMechanics.ipynb)) for reference.
+- For (deprecated) `.ipynb`-based notebooks the frontmatter needs to be given in the first cell. See existing notebooks for reference.
 
 ### Notebook setup
 
@@ -171,7 +171,7 @@ Add the notebook to CTest ([example](https://gitlab.opengeosys.org/ogs/ogs/-/blo
 
 ```cmake
 if(NOT OGS_USE_PETSC)
-    NotebookTest(NOTEBOOKFILE Mechanics/Linear/SimpleMechanics.ipynb RUNTIME 10)
+    NotebookTest(NOTEBOOKFILE Mechanics/Linear/SimpleMechanics.py RUNTIME 10)
 
     # Notebooks in web/content need to be prefixed with 'notebook-'!
     NotebookTest(NOTEBOOKFILE ../../web/content/docs/tutorials/bhe_meshing/notebook-bhe_meshing.md