diff --git a/src/bvpy/domains/geometry.py b/src/bvpy/domains/geometry.py
index a55865bbea884cb415d9fcb387c64d36b4ae602a..bb230f16be2e7f44fce645b2e7b3d3e8da535fd1 100644
--- a/src/bvpy/domains/geometry.py
+++ b/src/bvpy/domains/geometry.py
@@ -103,14 +103,14 @@ class Init_orientation_3D(fe.UserExpression):
 
 
 def normal(domain, scale=1, interior=None, degree=1, name='normal'):
-    """Compute a displacement field along the normals to a mesh.
+    """Computes the normal vector field to a given domain.
 
     Parameters
     ----------
     domain : subclass of :class:`Domain<bvpy.domains.AbstracDomain>`
         The meshed surface on which we want to compute the normal field.
     scale : float, optional
-        The displacement intensity to consider (the default is 1).
+        The norm of the vector field (the default is 1).
     interior : list(float), optional
         A 3D vector to set the mesh orientation (the default is None).
     degree : int, optional
@@ -128,7 +128,7 @@ def normal(domain, scale=1, interior=None, degree=1, name='normal'):
     if domain.mesh is not None:
         mesh = domain.mesh
     else:
-        print('WARNING - (geometry.py - l.84): no mesh on the domain.')
+        print('WARNING - (geometry.normal): no mesh on the domain.')
 
     np.seterr(divide='ignore', invalid='ignore')
     V = fe.VectorFunctionSpace(mesh, 'P', degree, dim=3)
@@ -224,7 +224,29 @@ def get_boundary_vertex(mesh):
     return bnd_vertex
 
 
-def boundary_normal(mesh, scale=1):
+def boundary_normal(domain, scale=1, name='boundary_normal'):
+    """Computes along the boundary the normals within the tangent plane.
+
+    Parameters
+    ----------
+    domain : subclass of :class:`Domain<bvpy.domains.AbstracDomain>`
+        The meshed surface on which we want to compute the normal field.
+    scale : float, optional
+        The norm of the vector field (the default is 1).
+    name : str, optional
+        the label we want to attach to the normal field.
+
+    Returns
+    -------
+    :class:`dolfin.cpp.function.Function<dolfin.cpp.function.Function>`
+        The computed normal field.
+
+    """
+    if domain.mesh is not None:
+        mesh = domain.mesh
+    else:
+        print('WARNING - (geometry.boundary_normal): no mesh on the domain.')
+
     mesh.init()
     bnd_facets = get_boundary_facets(mesh)
     facets_index = [f.index() for f in bnd_facets]
@@ -257,5 +279,7 @@ def boundary_normal(mesh, scale=1):
     d2v = fe.dof_to_vertex_map(V)
     u.vector().set_local(normals.flatten()[d2v])
 
+    u.rename(name, 'label')
+
     return u
 # geometry.py ends here