diff --git a/festim/hydrogen_transport_problem.py b/festim/hydrogen_transport_problem.py
index 3b9f2c549..97f547818 100644
--- a/festim/hydrogen_transport_problem.py
+++ b/festim/hydrogen_transport_problem.py
@@ -266,12 +266,10 @@ def run(self, final_time: float):
             mobile_xdmf.write_function(self.u, t)
 
             cm = self.species[0].solution
-            # TODO this should be a property of Mesh
-            n = FacetNormal(self.mesh.mesh)
             D = self.subdomains[0].material.get_diffusion_coefficient(
                 self.mesh.mesh, self.temperature
             )
-            surface_flux = form(D * dot(grad(cm), n) * self.ds(2))
+            surface_flux = form(D * dot(grad(cm), self.mesh.n) * self.ds(2))
             flux = assemble_scalar(surface_flux)
             flux_values.append(flux)
             times.append(t)
diff --git a/festim/mesh/mesh.py b/festim/mesh/mesh.py
index bb2aa0654..fc0980317 100644
--- a/festim/mesh/mesh.py
+++ b/festim/mesh/mesh.py
@@ -1,3 +1,6 @@
+import ufl
+
+
 class Mesh:
     """
     Mesh class
@@ -9,6 +12,7 @@ class Mesh:
         mesh (dolfinx.mesh.Mesh): the mesh
         vdim (int): the dimension of the mesh cells
         fdim (int): the dimension of the mesh facets
+        n (ufl.FacetNormal): the normal vector to the facets
     """
 
     def __init__(self, mesh=None):
@@ -32,3 +36,7 @@ def vdim(self):
     @property
     def fdim(self):
         return self.mesh.topology.dim - 1
+
+    @property
+    def n(self):
+        return ufl.FacetNormal(self.mesh)
diff --git a/test/test_permeation_problem.py b/test/test_permeation_problem.py
index 8ceddcda5..62b792b30 100644
--- a/test/test_permeation_problem.py
+++ b/test/test_permeation_problem.py
@@ -4,7 +4,7 @@
     dirichletbc,
     locate_dofs_topological,
 )
-from ufl import exp, FacetNormal
+from ufl import exp
 import numpy as np
 
 
@@ -44,11 +44,10 @@ def siverts_law(T, S_0, E_S, pressure):
         S = S_0 * exp(-E_S / F.k_B / T)
         return S * pressure**0.5
 
-    fdim = my_mesh.mesh.topology.dim - 1
     left_facets = my_model.facet_meshtags.find(1)
-    left_dofs = locate_dofs_topological(V, fdim, left_facets)
+    left_dofs = locate_dofs_topological(V, my_mesh.fdim, left_facets)
     right_facets = my_model.facet_meshtags.find(2)
-    right_dofs = locate_dofs_topological(V, fdim, right_facets)
+    right_dofs = locate_dofs_topological(V, my_mesh.fdim, right_facets)
 
     S_0 = 4.02e21
     E_S = 1.04
@@ -143,20 +142,15 @@ def test_permeation_problem_multi_volume():
     D = my_mat.get_diffusion_coefficient(my_mesh.mesh, temperature)
 
     V = my_model.function_space
-    u = mobile_H.solution
-
-    # TODO this should be a property of Mesh
-    n = FacetNormal(my_mesh.mesh)
 
     def siverts_law(T, S_0, E_S, pressure):
         S = S_0 * exp(-E_S / F.k_B / T)
         return S * pressure**0.5
 
-    fdim = my_mesh.mesh.topology.dim - 1
     left_facets = my_model.facet_meshtags.find(1)
-    left_dofs = locate_dofs_topological(V, fdim, left_facets)
+    left_dofs = locate_dofs_topological(V, my_mesh.fdim, left_facets)
     right_facets = my_model.facet_meshtags.find(2)
-    right_dofs = locate_dofs_topological(V, fdim, right_facets)
+    right_dofs = locate_dofs_topological(V, my_mesh.fdim, right_facets)
 
     S_0 = 4.02e21
     E_S = 1.04