diff --git a/Applications/Utils/ExtractBoundary.smk b/Applications/Utils/ExtractBoundary.smk
index 3809ec7a22442a450c649bb23e80b6c8142bd146..2d5837de3e3f2c31a937e22831292e18b98424f9 100644
--- a/Applications/Utils/ExtractBoundary.smk
+++ b/Applications/Utils/ExtractBoundary.smk
@@ -14,33 +14,26 @@ workdir: f"{config['Data_BINARY_DIR']}/{output_path}"
 elem_types = ['tri', 'quad']
 rule all:
     input:
-        expand("square_1x1_{type}_boundary_diff.out", type=elem_types)
+        expand("square_10_1x1_{type}_boundary_diff.out", type=elem_types)
 
-rule generate_meshes:
-    output:
-        "input_square_1x1_{type}.vtu"
-    shell:
-        """
-        generateStructuredMesh -e {wildcards.type} \
-            --lx 1 --ly 1 \
-            --nx 10 --ny 10 \
-            -o {output}
-        """
+include: f"{config['SOURCE_DIR']}/scripts/snakemake/modules/meshes.smk"
 
 rule extract_boundary:
     input:
-        "input_square_1x1_{type}.vtu"
+        rules.generate_square_mesh.output
     output:
-        "square_1x1_{type}_boundary.vtu"
+        "square_{size}_{lx}x{ly}_{type}_boundary.vtu"
     shell:
         "ExtractBoundary -i {input} -o {output}"
 
 rule vtkdiff:
     input:
-        "square_1x1_{type}_boundary.vtu"
+        a = rules.extract_boundary.output,
+        b = f"{config['Data_SOURCE_DIR']}/{output_path}/{rules.extract_boundary.output}"
     output:
-        "square_1x1_{type}_boundary_diff.out"
+        "square_{size}_{lx}x{ly}_{type}_boundary_diff.out"
     params:
+        check_mesh = True,
         fields = [
             # second field name can be omitted if identical
             ["bulk_node_ids", 0, 0],
diff --git a/Applications/Utils/VoxelGridFromLayers.smk b/Applications/Utils/VoxelGridFromLayers.smk
index 921c9e153094dd9de071961193e4c80525ed2c76..40e0b101222b804c62e829565419bf73b3c86bdd 100644
--- a/Applications/Utils/VoxelGridFromLayers.smk
+++ b/Applications/Utils/VoxelGridFromLayers.smk
@@ -11,8 +11,7 @@ out_dir = f"{config['Data_BINARY_DIR']}/MeshLib"
 rule all:
     input:
         f"{out_dir}/AREHS_test_diff_geometry.out",
-        f"{out_dir}/AREHS_test_fault_diff_material_ids.out",
-        f"{out_dir}/AREHS_test_fault_diff_geometry.out",
+        f"{out_dir}/AREHS_test_fault_diff.out",
         f"{out_dir}/AREHS_test_iso_diff_geometry.out"
 
 rule layers_to_grid:
@@ -41,23 +40,17 @@ rule add_fault_to_grid:
     shell:
         "AddFaultToVoxelGrid -i {input.grid} -f {input.fault} -o {output}"
 
-rule vtkdiff_material_ids:
+rule vtkdiff_fault:
     input:
-        out = rules.add_fault_to_grid.output,
-        ref = "AREHS_test_fault.vtu"
+        a = rules.add_fault_to_grid.output,
+        b = "AREHS_test_fault.vtu"
     output:
-        f"{out_dir}/AREHS_test_fault_diff_material_ids.out"
-    shell:
-        "vtkdiff -a MaterialIDs -b MaterialIDs {input.out} {input.ref} > {output}"
-
-rule vtkdiff_fault_geometry:
-    input:
-        out = rules.add_fault_to_grid.output,
-        ref = "AREHS_test_fault.vtu"
-    output:
-        f"{out_dir}/AREHS_test_fault_diff_geometry.out"
-    shell:
-        "vtkdiff -m {input.out} {input.ref} > {output}"
+        f"{out_dir}/AREHS_test_fault_diff.out"
+    params:
+        check_mesh = True,
+        fields = [["MaterialIDs", 0, 0]]
+    wrapper:
+        f"file://{config['SOURCE_DIR']}/scripts/snakemake/vtkdiff"
 
 rule layers_to_grid_iso:
     input:
diff --git a/Tests/Data/FileIO/square_1x1_quad_boundary.vtu b/Tests/Data/FileIO/square_10_1x1_quad_boundary.vtu
similarity index 100%
rename from Tests/Data/FileIO/square_1x1_quad_boundary.vtu
rename to Tests/Data/FileIO/square_10_1x1_quad_boundary.vtu
diff --git a/Tests/Data/FileIO/square_1x1_tri_boundary.vtu b/Tests/Data/FileIO/square_10_1x1_tri_boundary.vtu
similarity index 100%
rename from Tests/Data/FileIO/square_1x1_tri_boundary.vtu
rename to Tests/Data/FileIO/square_10_1x1_tri_boundary.vtu
diff --git a/scripts/snakemake/modules/meshes.smk b/scripts/snakemake/modules/meshes.smk
new file mode 100644
index 0000000000000000000000000000000000000000..08c0c6e247fbc575db39db305ce4261e16806b4e
--- /dev/null
+++ b/scripts/snakemake/modules/meshes.smk
@@ -0,0 +1,10 @@
+rule generate_square_mesh:
+    output:
+        "square_{size}_{lx}x{ly}_{type}.vtu"
+    shell:
+        """
+        generateStructuredMesh -e {wildcards.type} \
+            --lx {wildcards.lx} --ly {wildcards.ly} \
+            --nx {wildcards.size} --ny {wildcards.size} \
+            -o {output}
+        """
diff --git a/scripts/snakemake/vtkdiff/wrapper.py b/scripts/snakemake/vtkdiff/wrapper.py
index 564f9aeae03b8e3f363b81927f84f7bbec261fd8..74ed33ac140d8bae6f759e0173ad0cc63745101b 100644
--- a/scripts/snakemake/vtkdiff/wrapper.py
+++ b/scripts/snakemake/vtkdiff/wrapper.py
@@ -7,12 +7,12 @@ __license__ = "BSD"
 import os
 from snakemake.shell import shell
 
-output_path = (
-    os.getcwd().replace("\\", "/").replace(snakemake.config["Data_BINARY_DIR"], "")
-)
 if os.path.exists(snakemake.output[0]):
     os.remove(snakemake.output[0])
 
+if snakemake.params.check_mesh:
+    shell("vtkdiff {snakemake.input.a} {snakemake.input.b} -m > {snakemake.output[0]}")
+
 for field in snakemake.params.fields:
     field_a = field[0]
     offset = 0
@@ -24,8 +24,7 @@ for field in snakemake.params.fields:
 
     shell(
         """
-        vtkdiff {snakemake.input[0]} \
-          {snakemake.config[Data_SOURCE_DIR]}/{output_path}/{snakemake.input[0]} \
+        vtkdiff {snakemake.input.a} {snakemake.input.b} \
           -a {field_a} -b {field_b} \
           --abs {abs_tol} --rel {rel_tol} >> {snakemake.output[0]}
         """