diff --git a/Applications/Utils/GMSH2OGS_ExtractBoundary.smk b/Applications/Utils/GMSH2OGS_ExtractBoundary.smk
new file mode 100644
index 0000000000000000000000000000000000000000..08c084083dd2b84755035e631f8f135e45d3aa79
--- /dev/null
+++ b/Applications/Utils/GMSH2OGS_ExtractBoundary.smk
@@ -0,0 +1,84 @@
+# Usage, e.g.:
+#   snakemake -s GMSH2OGS_ExtractBoundary.smk -j 1 --configfile $HOME/code/ogs6/build/buildinfo.yaml
+#
+# buildinfo.yaml contains variables such as Data_BINARY_DIR
+
+import os
+os.environ["PATH"] += os.pathsep + os.pathsep.join([config['BIN_DIR']])
+workdir: f"{config['Data_BINARY_DIR']}/MeshLib"
+
+#import ipdb
+#ipdb.set_trace()
+
+VTKDIFF = f"file://{config['SOURCE_DIR']}/scripts/snakemake/vtkdiff"
+SOURCE_DIR = f"{config['Data_SOURCE_DIR']}/MeshLib"
+
+MESH_NAME = "A2"
+BULK_MESH = f"{MESH_NAME}.vtu"
+GMSH_MESH = f"{MESH_NAME}-gmsh.msh"
+INNER_BOUNDARIES = "[12]"
+OUTER_BOUNDARIES = "[034567]"
+
+ALL_REFERENCE_FILES = glob_wildcards(f"{SOURCE_DIR}/{{file, {MESH_NAME}.*}}.vtu").file
+
+rule all:
+    input:
+        expand("{file}.vtkdiff.out", file=ALL_REFERENCE_FILES)
+
+rule gmsh2ogs:
+    input:
+        f"{SOURCE_DIR}/{GMSH_MESH}"
+    output:
+        expand("{file}.vtu", file=ALL_REFERENCE_FILES)
+    shell:
+        "GMSH2OGS -i {input} -o {BULK_MESH} -e -b -v"
+
+rule vtkdiff_geometry_bulk:
+    input:
+        a = f"{SOURCE_DIR}/{BULK_MESH}",
+        b = BULK_MESH
+    output:
+        f"{MESH_NAME}.vtkdiff.out"
+    params:
+        check_mesh = True,
+        fields = [["MaterialIDs", 0, 0]]
+    wrapper:
+        VTKDIFF
+
+rule vtkdiff_geometry_boundary:
+    input:
+        a = f"{SOURCE_DIR}/{MESH_NAME}_{{x}}.vtu",
+        b = f"{MESH_NAME}_{{x}}.vtu"
+    output:
+        f"{MESH_NAME}_{{x}}.vtkdiff.out-bulk_node_ids"
+    params:
+        check_mesh = True,
+        fields = [["bulk_node_ids", 0, 0]]
+    wrapper:
+        VTKDIFF
+
+rule vtkdiff_geometry_inner:
+    input:
+        f"{MESH_NAME}_{{x}}.vtkdiff.out-bulk_node_ids",
+        a = f"{SOURCE_DIR}/{MESH_NAME}_{{x}}.vtu",
+        b = f"{MESH_NAME}_{{x}}.vtu"
+    output:
+        f"{MESH_NAME}_{{x, {INNER_BOUNDARIES}}}.vtkdiff.out"
+    params:
+        check_mesh = True,
+        fields = [["number_bulk_elements", 0, 0]]
+    wrapper:
+        VTKDIFF
+
+rule vtkdiff_geometry_outer:
+    input:
+        f"{MESH_NAME}_{{x}}.vtkdiff.out-bulk_node_ids",
+        a = f"{SOURCE_DIR}/{MESH_NAME}_{{x}}.vtu",
+        b = f"{MESH_NAME}_{{x}}.vtu"
+    output:
+        f"{MESH_NAME}_{{x, {OUTER_BOUNDARIES}}}.vtkdiff.out"
+    params:
+        check_mesh = True,
+        fields = [["bulk_element_ids", 0, 0]]
+    wrapper:
+        VTKDIFF
diff --git a/Applications/Utils/Tests.cmake b/Applications/Utils/Tests.cmake
index a1ac37c19da430f3a2922a6909c03def3687d3b3..cb08f5a929c00381be6b689935adb1548cf9100a 100644
--- a/Applications/Utils/Tests.cmake
+++ b/Applications/Utils/Tests.cmake
@@ -481,6 +481,17 @@ if(TARGET VerticalSliceFromLayers AND GMSH)
         REQUIREMENTS NOT OGS_USE_MPI
         DIFF_DATA AmmerSlice.vtu AmmerSlice.vtu 1e-16
     )
+
+endif()
+
+
+if(TARGET GMSH2OGS AND SNAKEMAKE AND NOT OGS_USE_MPI)
+    add_test(NAME snakemake_GMSH2OGS_ExtractBoundary
+        COMMAND ${SNAKEMAKE} --cores all
+        --configfile ${PROJECT_BINARY_DIR}/buildinfo.yaml
+        -s ${CMAKE_CURRENT_SOURCE_DIR}/GMSH2OGS_ExtractBoundary.smk
+        )
+    add_dependencies(ctest GMSH2OGS)
 endif()
 
 foreach(criterion ElementSize EdgeRatio EquiAngleSkew RadiusEdgeRatio SizeDifference)