diff --git a/.github/workflows/run-benchmark-examples.yml b/.github/workflows/run-benchmark-examples.yml
new file mode 100644
index 0000000..9ec96c6
--- /dev/null
+++ b/.github/workflows/run-benchmark-examples.yml
@@ -0,0 +1,54 @@
+name: Examples-CI
+on:
+  push:
+
+  pull_request:
+    branches: [ main ]
+
+  # Allows you to run this workflow manually from the Actions tab
+  workflow_dispatch:
+
+  # Runs the workflow once per day at 3:15am
+  schedule:
+    - cron: '3 16 * * *'
+
+env:
+  CACHE_NUMBER: 1  # increase to reset cache manually
+  SNAKEMAKE_RESULT_FILE: metadata4ing_provenance
+  PROVENANACE_FILE_NAME: element_size_vs_max_mises_stress.pdf
+
+jobs:
+  run-examples:
+    runs-on: ubuntu-latest
+    steps:
+      - name: checkout repo content
+        uses: actions/checkout@v2
+
+      - name: Setup Mambaforge
+        uses: conda-incubator/setup-miniconda@v3
+        with:
+            miniforge-version: latest
+            activate-environment: model-validation 
+            use-mamba: true
+
+      - name: Set strict channel priority
+        run: conda config --set channel_priority strict
+        
+      - name: Update environment
+        run: mamba env update -n model-validation -f environment_benchmarks.yml
+
+      - name: Run linear-elastic-plate-with-hole using fenics
+        shell: bash -l {0}
+        run: |
+          cd $GITHUB_WORKSPACE/examples/linear-elastic-plate-with-hole/fenics/
+          python run_benchmark.py
+
+      - name: Archive results of the fenics run of linear-elastic-plate-with-hole
+        uses: actions/upload-artifact@v4
+        with:
+            name: snakemake_results_linear-elastic-plate-with-hole
+            path: |
+              examples/linear-elastic-plate-with-hole/fenics/results/
+
+
+
diff --git a/examples/linear-elastic-plate-with-hole/fenics/environment_simulation.yml b/examples/linear-elastic-plate-with-hole/fenics/environment_simulation.yml
new file mode 100644
index 0000000..6ed60ef
--- /dev/null
+++ b/examples/linear-elastic-plate-with-hole/fenics/environment_simulation.yml
@@ -0,0 +1,16 @@
+name: fenics_simulation
+# Environment file for fenics simulation scripts. Called by fenics tool workflow.
+
+channels:
+  - conda-forge
+
+channel_priority: strict
+
+dependencies:
+  - python=3.12
+  - fenics-dolfinx=0.9.*
+  - mpich
+  - petsc4py
+  - pint
+  - python-gmsh
+  - sympy
diff --git a/examples/linear-elastic-plate-with-hole/fenics/run_benchmark.py b/examples/linear-elastic-plate-with-hole/fenics/run_benchmark.py
new file mode 100644
index 0000000..929fd72
--- /dev/null
+++ b/examples/linear-elastic-plate-with-hole/fenics/run_benchmark.py
@@ -0,0 +1,75 @@
+from pathlib import Path
+import zipfile
+import json
+import shutil
+import subprocess
+
+"""
+The script performs the following steps:
+
+1. Extracts the benchmark files from a zip archive (currently assuming that it is an RO-Crate of the benchmark).
+2. Iterates through the parameter configuration files, checks the "element-size" value, and if it meets the specified condition (>= 0.025)
+, it executes the Snakemake workflow for that configuration.
+
+The results of each run (and the files used by it) are stored in the directory with the configuration name.
+"""
+
+####################################################################################################
+####################################################################################################
+# Benchmark Extraction
+####################################################################################################
+####################################################################################################
+
+root_zipped_benchmark_dir = Path(__file__).resolve().parent.parent
+root_unzipped_benchmark_dir = Path(__file__).resolve().parent
+
+with zipfile.ZipFile(root_zipped_benchmark_dir / "linear-elastic-plate-with-hole.zip", 'r') as zip_ref:
+    # Extract all files
+    zip_ref.extractall(root_unzipped_benchmark_dir)
+    
+    
+#Creates a directory to store the conda environments. The environments are shared across different parameter configurations.
+#To avoid redundant creation of environments, this path will be passed to all snakemake files during execution.
+        
+shared_env_dir = root_unzipped_benchmark_dir / "conda_envs"
+shared_env_dir.mkdir(parents=True, exist_ok=True)  
+
+####################################################################################################
+####################################################################################################
+# Conditional execution of parameter configurations 
+####################################################################################################
+####################################################################################################
+  
+for file in root_unzipped_benchmark_dir.glob("parameters_*.json"):
+    with open(file, "r") as f:
+        data = json.load(f)
+        if data.get("element-size").get("value") >= 0.025:
+            
+            # Create output directory for the configuration
+            output_dir = root_unzipped_benchmark_dir / "results" / data.get("configuration")
+            output_dir.mkdir(parents=True, exist_ok=True) 
+            
+            # Copy the selected parameter file to the output directory with a standardised name
+            with open(output_dir / "parameters.json", "w") as outfile:
+                json.dump(data, outfile, indent=2)
+
+            # Copy files from benchmark_dir to output_dir, excluding non-matching parameter files.
+            for item in root_unzipped_benchmark_dir.iterdir():
+                if item.is_file():
+                    if item.name.startswith("parameters_") and item.suffix == ".json":
+                        continue
+                    else:
+                        shutil.copy(item, output_dir / item.name)
+                            
+            # Run the Snakemake workflow for the configuration
+            subprocess.run(["snakemake", "--use-conda", "--force", "--cores", "all", "--conda-prefix", str(shared_env_dir)], check=True, cwd=output_dir)
+            print("Workflow executed successfully.")
+            
+            # For the scenario where the snakemake workflow doesn't exist, one can directly run the simulation script using the subprocess module, e.g.:
+            #subprocess.run(["python", "run_fenics_simulation.py" \
+                            #"--input_parameter_file" "parameters.json" \
+                            #"--input_mesh_file" "mesh.msh" \
+                            #"--output_solution_file_zip" "solution_field_data.zip" \
+                            #"--output_metrics_file" "solution_metrics.json"], check=True, cwd=output_dir)
+                            
+            #Assuming the mesh.msh and parameters.json files are present/copied to the output_dir.
\ No newline at end of file
diff --git a/examples/linear-elastic-plate-with-hole/fenics/run_fenics_simulation.py b/examples/linear-elastic-plate-with-hole/fenics/run_fenics_simulation.py
new file mode 100644
index 0000000..6b5a114
--- /dev/null
+++ b/examples/linear-elastic-plate-with-hole/fenics/run_fenics_simulation.py
@@ -0,0 +1,357 @@
+import json
+import sys
+from argparse import ArgumentParser
+
+from pathlib import Path
+import dolfinx as df
+import basix.ufl
+import numpy as np
+import ufl
+from dolfinx.fem.petsc import LinearProblem
+from petsc4py.PETSc import ScalarType
+from mpi4py import MPI
+from pint import UnitRegistry
+
+# Add parent directory to sys.path
+sys.path.insert(0, str(Path(__file__).resolve().parent.parent))
+from analytical_solution import AnalyticalSolution
+
+
+def run_fenics_simulation(
+    parameter_file: str, mesh_file: str, solution_file_zip: str, metrics_file: str
+) -> None:
+    ureg = UnitRegistry()
+    with open(parameter_file) as f:
+        parameters = json.load(f)
+
+    mesh, cell_tags, facet_tags = df.io.gmshio.read_from_msh(
+        mesh_file,
+        comm=MPI.COMM_WORLD,
+        gdim=2,
+    )
+
+    V = df.fem.functionspace(mesh, ("CG", parameters["element-degree"], (2,)))
+
+    tags_left = facet_tags.find(1)
+    tags_bottom = facet_tags.find(2)
+    tags_right = facet_tags.find(3)
+    tags_top = facet_tags.find(4)
+
+    # Boundary conditions
+    dofs_left = df.fem.locate_dofs_topological(V.sub(0), 1, tags_left)
+    dofs_bottom = df.fem.locate_dofs_topological(V.sub(1), 1, tags_bottom)
+    dofs_right = df.fem.locate_dofs_topological(V, 1, tags_right)
+    dofs_top = df.fem.locate_dofs_topological(V, 1, tags_top)
+
+    bc_left = df.fem.dirichletbc(0.0, dofs_left, V.sub(0))
+    bc_bottom = df.fem.dirichletbc(0.0, dofs_bottom, V.sub(1))
+
+    E = (
+        ureg.Quantity(
+            parameters["young-modulus"]["value"], parameters["young-modulus"]["unit"]
+        )
+        .to_base_units()
+        .magnitude
+    )
+    nu = (
+        ureg.Quantity(
+            parameters["poisson-ratio"]["value"], parameters["poisson-ratio"]["unit"]
+        )
+        .to_base_units()
+        .magnitude
+    )
+    radius = (
+        ureg.Quantity(parameters["radius"]["value"], parameters["radius"]["unit"])
+        .to_base_units()
+        .magnitude
+    )
+    L = (
+        ureg.Quantity(parameters["length"]["value"], parameters["length"]["unit"])
+        .to_base_units()
+        .magnitude
+    )
+    load = (
+        ureg.Quantity(parameters["load"]["value"], parameters["load"]["unit"])
+        .to_base_units()
+        .magnitude
+    )
+    displacement_evaluation_point = parameters["displacement-evaluation-point"]
+    displacement_evaluation_x = (
+        ureg.Quantity(
+            displacement_evaluation_point["x"]["value"],
+            displacement_evaluation_point["x"]["unit"],
+        )
+        .to_base_units()
+        .magnitude
+    )
+    displacement_evaluation_y = (
+        ureg.Quantity(
+            displacement_evaluation_point["y"]["value"],
+            displacement_evaluation_point["y"]["unit"],
+        )
+        .to_base_units()
+        .magnitude
+    )
+
+    analytical_solution = AnalyticalSolution(
+        E=E,
+        nu=nu,
+        radius=radius,
+        L=L,
+        load=load,
+    )
+
+    def eps(v):
+        return ufl.sym(ufl.grad(v))
+
+    def sigma(v):
+        # plane stress
+        epsilon = eps(v)
+        return (
+            E
+            / (1.0 - nu**2)
+            * ((1.0 - nu) * epsilon + nu * ufl.tr(epsilon) * ufl.Identity(2))
+        )
+
+    def as_tensor(v):
+        return ufl.as_matrix([[v[0], v[2]], [v[2], v[1]]])
+
+    dx = ufl.Measure(
+        "dx",
+        metadata={
+            "quadrature_degree": parameters["quadrature-degree"],
+            "quadrature_scheme": parameters["quadrature-rule"],
+        },
+    )
+    ds = ufl.Measure(
+        "ds",
+        domain=mesh,
+        subdomain_data=facet_tags,
+    )
+    stress_space = df.fem.functionspace(
+        mesh, ("CG", parameters["element-degree"], (2, 2))
+    )
+    stress_function = df.fem.Function(stress_space)
+
+    u = df.fem.Function(V, name="u")
+    u_prescribed = df.fem.Function(V, name="u_prescribed")
+    u_prescribed.interpolate(lambda x: analytical_solution.displacement(x))
+    u_prescribed.x.scatter_forward()
+
+    u_ = ufl.TestFunction(V)
+    v_ = ufl.TrialFunction(V)
+    a = df.fem.form(ufl.inner(sigma(u_), eps(v_)) * dx)
+
+    # set rhs to zero
+    f = df.fem.form(ufl.inner(df.fem.Constant(mesh, np.array([0.0, 0.0])), u_) * ufl.ds)
+
+    bc_right = df.fem.dirichletbc(u_prescribed, dofs_right)
+    bc_top = df.fem.dirichletbc(u_prescribed, dofs_top)
+    solver = LinearProblem(
+        a,
+        f,
+        bcs=[bc_left, bc_bottom, bc_right, bc_top],
+        u=u,
+        petsc_options={
+            "ksp_type": "gmres",
+            "ksp_rtol": 1e-14,
+            "ksp_atol": 1e-14,
+        },
+    )
+    solver.solve()
+
+    # Evaluate displacement at the specified evaluation point
+    displacement_eval_point = np.array(
+        [[displacement_evaluation_x, displacement_evaluation_y, 0.0]],
+        dtype=np.float64,
+    )
+    tree = df.geometry.bb_tree(mesh, mesh.topology.dim)
+    cell_candidates = df.geometry.compute_collisions_points(
+        tree, displacement_eval_point
+    )
+    colliding_cells = df.geometry.compute_colliding_cells(
+        mesh, cell_candidates, displacement_eval_point
+    )
+    local_displacement_x = None
+    if len(colliding_cells.links(0)) > 0:
+        cell = colliding_cells.links(0)[0]
+        local_displacement_x = u.eval(
+            displacement_eval_point, np.array([cell], dtype=np.int32)
+        )[0]
+
+
+    def project(
+        v: df.fem.Function | ufl.core.expr.Expr,
+        V: df.fem.FunctionSpace,
+        dx: ufl.Measure = ufl.dx,
+    ) -> df.fem.Function:
+        """
+        Calculates an approximation of `v` on the space `V`
+
+        Args:
+            v: The expression that we want to evaluate.
+            V: The function space on which we want to evaluate.
+            dx: The measure that is used for the integration. This is important, if
+            either `V` is a quadrature space or `v` is a ufl expression containing a quadrature space.
+
+        Returns:
+            A function if `u` is None, otherwise `None`.
+
+        """
+        dv = ufl.TrialFunction(V)
+        v_ = ufl.TestFunction(V)
+        a_proj = ufl.inner(dv, v_) * dx
+        b_proj = ufl.inner(v, v_) * dx
+
+        solver = LinearProblem(a_proj, b_proj)
+        uh = solver.solve()
+        return uh
+
+    plot_space_stress = df.fem.functionspace(
+        mesh, ("DG", parameters["element-degree"] - 1, (2, 2))
+    )
+    plot_space_mises = df.fem.functionspace(
+        mesh, ("DG", parameters["element-degree"] - 1, (1,))
+    )
+    stress_nodes_red = project(sigma(u), plot_space_stress, dx)
+    stress_nodes_red.name = "stress"
+
+    def mises_stress(u):
+        stress = sigma(u)
+        p = ufl.tr(stress) / 3.0
+        s = stress - p * ufl.Identity(2)
+        return ufl.as_vector([(3.0 / 2.0) ** 0.5 * (ufl.inner(s, s) + p * p) ** 0.5])
+
+    mises_stress_nodes = project(mises_stress(u), plot_space_mises, dx)
+    mises_stress_nodes.name = "von_mises_stress"
+
+    # Write each function to its own VTK file on all ranks
+    output_dir = Path(solution_file_zip).parent
+    with df.io.VTKFile(
+        MPI.COMM_WORLD,
+        str(
+            output_dir
+            / f"solution_field_data_displacements_{parameters['configuration']}.vtk"
+        ),
+        "w",
+    ) as vtk:
+        vtk.write_function(u, 0.0)
+    with df.io.VTKFile(
+        MPI.COMM_WORLD,
+        str(
+            output_dir / f"solution_field_data_stress_{parameters['configuration']}.vtk"
+        ),
+        "w",
+    ) as vtk:
+        vtk.write_function(stress_nodes_red, 0.0)
+    with df.io.VTKFile(
+        MPI.COMM_WORLD,
+        str(
+            output_dir
+            / f"solution_field_data_mises_stress_{parameters['configuration']}.vtk"
+        ),
+        "w",
+    ) as vtk:
+        vtk.write_function(mises_stress_nodes, 0.0)
+
+    # extract maximum von Mises stress
+    max_mises_stress_nodes = np.max(mises_stress_nodes.x.array)
+
+    # Compute von Mises stress at quadrature (Gauss) points and extract maximum (global across MPI)
+    quad_element = basix.ufl.quadrature_element(
+        mesh.topology.cell_name(),
+        value_shape=(1,),
+        degree=parameters["quadrature-degree"],
+    )
+
+    Q_mises = df.fem.functionspace(mesh, quad_element)
+    mises_qp = df.fem.Function(Q_mises, name="von_mises_stress_qp")
+    expr_qp = df.fem.Expression(mises_stress(u), Q_mises.element.interpolation_points())
+    mises_qp.interpolate(expr_qp)
+    max_mises_stress_gauss_points = MPI.COMM_WORLD.allreduce(
+        np.max(mises_qp.x.array), op=MPI.MAX
+    )
+    
+    displacement_x_at_evaluation_point = None
+    if MPI.COMM_WORLD.rank == 0:
+        displacement_x_candidates = (
+            MPI.COMM_WORLD.gather(local_displacement_x, root=0) or []
+        )
+        for value in displacement_x_candidates:
+            if value is not None:
+                displacement_x_at_evaluation_point = value
+                break
+
+        if displacement_x_at_evaluation_point is None:
+            raise ValueError(
+                "Could not evaluate displacement at the configured evaluation point."
+            )
+    else:
+        MPI.COMM_WORLD.gather(local_displacement_x, root=0)
+
+    # Save metrics
+    metrics = {
+        "max_von_mises_stress_nodes": max_mises_stress_nodes,
+        "max_von_mises_stress_gauss_points": max_mises_stress_gauss_points,
+        f"displacement_x_at_evaluation_point (x={displacement_evaluation_x}, y={displacement_evaluation_y})": displacement_x_at_evaluation_point,
+    }
+
+    if MPI.COMM_WORLD.rank == 0:
+        with open(metrics_file, "w") as f:
+            json.dump(metrics, f, indent=4)
+        # store all .vtu, .pvtu and .vtk files for this configuration in the zip file
+        import zipfile
+
+        config = parameters["configuration"]
+        file_patterns = [
+            str(output_dir / f"solution_field_data_displacements_{config}*"),
+            str(output_dir / f"solution_field_data_stress_{config}*"),
+            str(output_dir / f"solution_field_data_mises_stress_{config}*"),
+        ]
+
+        files_to_store = []
+        for pattern in file_patterns:
+            files_to_store.extend(
+                filter(
+                    # filter for all file endings because this is not possible with glob
+                    lambda path: path.suffix in [".vtk", ".vtu", ".pvtu"],
+                    Path().glob(pattern),
+                )
+            )
+            # files_to_store.extend(Path().glob(pattern))
+        with zipfile.ZipFile(solution_file_zip, "w") as zipf:
+            for filepath in files_to_store:
+                zipf.write(filepath, arcname=filepath.name)
+
+
+if __name__ == "__main__":
+    parser = ArgumentParser(
+        description="Run FEniCS simulation for a plate with a hole.\n"
+        "Inputs: --input_parameter_file, --input_mesh_file\n"
+        "Outputs: --output_solution_file_hdf5, --output_metrics_file"
+    )
+    parser.add_argument(
+        "--input_parameter_file",
+        required=True,
+        help="JSON file containing simulation parameters (input)",
+    )
+    parser.add_argument(
+        "--input_mesh_file", required=True, help="Path to the mesh file (input)"
+    )
+    parser.add_argument(
+        "--output_solution_file_zip",
+        required=True,
+        help="Path to the zipped solution files (output)",
+    )
+    parser.add_argument(
+        "--output_metrics_file",
+        required=True,
+        help="Path to the output metrics JSON file (output)",
+    )
+    args, _ = parser.parse_known_args()
+    run_fenics_simulation(
+        args.input_parameter_file,
+        args.input_mesh_file,
+        args.output_solution_file_zip,
+        args.output_metrics_file,
+    )
diff --git a/examples/linear-elastic-plate-with-hole/linear-elastic-plate-with-hole.zip b/examples/linear-elastic-plate-with-hole/linear-elastic-plate-with-hole.zip
new file mode 100644
index 0000000..9644a81
Binary files /dev/null and b/examples/linear-elastic-plate-with-hole/linear-elastic-plate-with-hole.zip differ