diff --git a/src/foam2dolfinx/open_foam_reader.py b/src/foam2dolfinx/open_foam_reader.py
index a17b938..66e9cd9 100644
--- a/src/foam2dolfinx/open_foam_reader.py
+++ b/src/foam2dolfinx/open_foam_reader.py
@@ -176,7 +176,9 @@ def _create_dolfinx_mesh(self, subdomain: str | None = "default"):
         )
 
         # Create dolfinx Mesh
-        mesh_ufl = ufl.Mesh(self.mesh_vector_element)
+        mesh_ufl = ufl.Mesh(
+            basix.ufl.element("Lagrange", cell.cellname(), degree, shape=(3,))
+        )
         self.dolfinx_meshes_dict[subdomain] = create_mesh(
             comm=MPI.COMM_WORLD,
             cells=self.connectivities_dict[subdomain],
@@ -184,7 +186,49 @@ def _create_dolfinx_mesh(self, subdomain: str | None = "default"):
             e=mesh_ufl,
         )
 
-    def create_dolfinx_function(
+    def create_dolfinx_function_with_cell_data(
+        self, t: float, name: str = "U", subdomain: str | None = "default"
+    ) -> dolfinx.fem.Function:
+        """Creates a dolfinx.fem.Function from the OpenFOAM file.
+
+        Args:
+            t: timestamp of the data to read
+            name: Name of the field in the OpenFOAM file, defaults to "U" for velocity
+            subdomain: Name of the subdmain in the OpenFOAM file, from which a field is
+                extracted
+
+        Returns:
+            the dolfinx function
+        """
+
+        # read the OpenFOAM data in the filename provided
+        self._read_with_pyvista(t=t, subdomain=subdomain)
+
+        # create the dolfinx mesh
+        if subdomain not in self.dolfinx_meshes_dict:
+            self._create_dolfinx_mesh(subdomain=subdomain)
+
+        mesh = self.dolfinx_meshes_dict[subdomain]
+
+        if name == "U":
+            element = basix.ufl.element("DG", mesh.topology.cell_name(), 0, shape=(3,))
+        else:
+            element = basix.ufl.element("DG", mesh.topology.cell_name(), 0, shape=())
+
+        function_space = dolfinx.fem.functionspace(mesh, element)
+        u = dolfinx.fem.Function(function_space)
+
+        # Assign values in OF_mesh to dolfinx_mesh
+        assert hasattr(self.OF_meshes_dict[subdomain], "cell_data")
+        u.x.array[:] = (
+            self.OF_meshes_dict[subdomain]
+            .cell_data[name][mesh.topology.original_cell_index]
+            .flatten()
+        )
+
+        return u
+
+    def create_dolfinx_function_with_point_data(
         self, t: float, name: str = "U", subdomain: str | None = "default"
     ) -> dolfinx.fem.Function:
         """Creates a dolfinx.fem.Function from the OpenFOAM file.
diff --git a/test/test_example.py b/test/test_example.py
index cd4f24d..1231ab5 100644
--- a/test/test_example.py
+++ b/test/test_example.py
@@ -1,13 +1,15 @@
-from foam2dolfinx import OpenFOAMReader
-import dolfinx
-from pyvista import examples
 import zipfile
 from pathlib import Path
 
+import dolfinx
+from pyvista import examples
+
+from foam2dolfinx import OpenFOAMReader
+
 
 def test_reading_and_writing_cavity_example():
     my_of_reader = OpenFOAMReader(filename=examples.download_cavity(load=False))
-    vel = my_of_reader.create_dolfinx_function(t=2.5)
+    vel = my_of_reader.create_dolfinx_function_with_point_data(t=2.5)
 
     assert isinstance(vel, dolfinx.fem.Function)
 
@@ -28,8 +30,8 @@ def test_baby_example(tmpdir):
     # read the .foam file
     my_of_reader = OpenFOAMReader(filename=str(foam_file), cell_type=10)
 
-    vel = my_of_reader.create_dolfinx_function(t=time, name="U")
-    T = my_of_reader.create_dolfinx_function(t=time, name="T")
+    vel = my_of_reader.create_dolfinx_function_with_point_data(t=time, name="U")
+    T = my_of_reader.create_dolfinx_function_with_point_data(t=time, name="T")
 
     assert isinstance(vel, dolfinx.fem.Function)
     assert isinstance(T, dolfinx.fem.Function)
@@ -51,9 +53,9 @@ def test_hot_room(tmpdir):
     # read the .foam file
     my_of_reader = OpenFOAMReader(filename=str(foam_file), cell_type=12)
 
-    vel = my_of_reader.create_dolfinx_function(t=time, name="U")
-    T = my_of_reader.create_dolfinx_function(t=time, name="T")
-    nut = my_of_reader.create_dolfinx_function(t=time, name="nut")
+    vel = my_of_reader.create_dolfinx_function_with_point_data(t=time, name="U")
+    T = my_of_reader.create_dolfinx_function_with_point_data(t=time, name="T")
+    nut = my_of_reader.create_dolfinx_function_with_point_data(t=time, name="nut")
 
     assert isinstance(vel, dolfinx.fem.Function)
     assert isinstance(T, dolfinx.fem.Function)