diff --git a/src/foam2dolfinx/open_foam_reader.py b/src/foam2dolfinx/open_foam_reader.py
index a86d50f..cea4dd5 100644
--- a/src/foam2dolfinx/open_foam_reader.py
+++ b/src/foam2dolfinx/open_foam_reader.py
@@ -30,6 +30,10 @@ class OpenFOAMReader:
             connectivities_dict: dictionary of the OpenFOAM mesh cell connectivity with
                 vertices reordered in a sorted order for mapping with the dolfinx mesh.
             dolfinx_meshes_dict: dictionary of dolfinx meshes
+            vertex_maps_dict: dictionary of vertex maps per subdomain, mapping dolfinx
+                vertex indices to the original OpenFOAM point indices. Built once per
+                subdomain on the first call to create_dolfinx_function_with_point_data
+                and cached for subsequent calls.
 
         Notes:
             The cell type refers to the VTK cell type, a full list of cells and their
@@ -49,6 +53,7 @@ class OpenFOAMReader:
     OF_cells_dict: dict[str, np.ndarray]
     connectivities_dict: dict[str, np.ndarray]
     dolfinx_meshes_dict: dict[str, dolfinx.mesh.Mesh]
+    vertex_maps_dict: dict[str, np.ndarray]
 
     def __init__(self, filename, cell_type: int = 12):
         self.filename = filename
@@ -61,6 +66,7 @@ def __init__(self, filename, cell_type: int = 12):
         self.OF_cells_dict = {}
         self.connectivities_dict = {}
         self.dolfinx_meshes_dict = {}
+        self.vertex_maps_dict = {}
 
     @property
     def cell_type(self):
@@ -104,9 +110,6 @@ def _read_with_pyvista(self, t: float, subdomain: str | None = "default"):
         else:
             self.OF_meshes_dict[subdomain] = OF_multiblock["internalMesh"]
 
-        # Ensure the mesh has cell data
-        assert hasattr(self.OF_meshes_dict[subdomain], "cells_dict")
-
         # obtain dictionary of cell types in OF_mesh
         OF_cell_type_dict = self.OF_meshes_dict[subdomain].cells_dict
 
@@ -156,16 +159,6 @@ def _create_dolfinx_mesh(self, subdomain: str | None = "default"):
             rows, args_conn
         ]
 
-        # Define mesh element
-        if self.cell_type == 12:
-            shape = "hexahedron"
-        elif self.cell_type == 10:
-            shape = "tetrahedron"
-        else:
-            raise ValueError(
-                f"Cell type: {self.cell_type}, not supported, please use"
-                " either 12 (hexahedron) or 10 (tetrahedron) cells in OF mesh"
-            )
         degree = 1  # Set polynomial degree
         cell = ufl.Cell(shape)
         # ufl.Cell.cellname became a property after dolfinx v0.10
@@ -188,25 +181,22 @@ def _create_dolfinx_mesh(self, subdomain: str | None = "default"):
             e=mesh_ufl,
         )
 
-    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 using cell data.
+    def _get_mesh(
+        self, t: float, name: str, subdomain: str | None
+    ) -> dolfinx.mesh.Mesh:
+        """Reads pyvista data for the given time and subdomain, validates that the
+        requested field exists, builds the dolfinx mesh on first call, and returns it.
 
         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
+            name: name of the field to validate against the available cell data keys
+            subdomain: subdomain key in the OpenFOAM multiblock dataset
 
         Returns:
-            the dolfinx function
+            the dolfinx mesh for this subdomain
         """
-
-        # read the OpenFOAM data in the filename provided
         self._read_with_pyvista(t=t, subdomain=subdomain)
 
-        # test if name is in the cell data of the OF mesh
         if name not in self.OF_meshes_dict[subdomain].cell_data.keys():
             raise ValueError(
                 f"Function name: {name} not found in the subdomain: {subdomain}, "
@@ -215,28 +205,38 @@ def create_dolfinx_function_with_cell_data(
                 f"{self.OF_meshes_dict[subdomain].cell_data.keys()}"
             )
 
-        # create the dolfinx mesh
         if subdomain not in self.dolfinx_meshes_dict:
             self._create_dolfinx_mesh(subdomain=subdomain)
 
-        mesh = self.dolfinx_meshes_dict[subdomain]
+        return 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=())
+    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 using cell data.
 
-        function_space = dolfinx.fem.functionspace(mesh, element)
-        u = dolfinx.fem.Function(function_space)
+        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
 
-        # 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()
+        Returns:
+            the dolfinx function
+        """
+        mesh = self._get_mesh(t, name, subdomain)
+
+        data = self.OF_meshes_dict[subdomain].cell_data[name]
+        # () for scalars, (3,) for vectors
+        data_shape = data.shape[1:] if data.ndim > 1 else ()
+        element = basix.ufl.element(
+            "DG", mesh.topology.cell_name(), 0, shape=data_shape
         )
 
+        function_space = dolfinx.fem.functionspace(mesh, element)
+        u = dolfinx.fem.Function(function_space)
+        u.x.array[:] = data[mesh.topology.original_cell_index].flatten()
+
         return u
 
     def create_dolfinx_function_with_point_data(
@@ -253,59 +253,40 @@ def create_dolfinx_function_with_point_data(
         Returns:
             the dolfinx function
         """
+        mesh = self._get_mesh(t, name, subdomain)
 
-        # read the OpenFOAM data in the filename provided
-        self._read_with_pyvista(t=t, subdomain=subdomain)
+        # select scalar or vector element based on data dimensionality
+        data = self.OF_meshes_dict[subdomain].point_data[name]
+        element = (
+            self.mesh_vector_element if data.ndim > 1 else self.mesh_scalar_element
+        )
 
-        # test if name is in the cell data of the OF mesh
-        if name not in self.OF_meshes_dict[subdomain].cell_data.keys():
-            raise ValueError(
-                f"Function name: {name} not found in the subdomain: {subdomain}, "
-                "in the OpenFOAM file. "
-                f"Available functions in subdomain: {subdomain} : "
-                f"{self.OF_meshes_dict[subdomain].cell_data.keys()}"
+        if subdomain not in self.vertex_maps_dict:
+            num_vertices = (
+                mesh.topology.index_map(0).size_local
+                + mesh.topology.index_map(0).num_ghosts
             )
-
-        # 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 = self.mesh_vector_element
-        else:
-            element = self.mesh_scalar_element
+            vertex_map = np.empty(num_vertices, dtype=np.int32)
+            c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
+            num_cells = (
+                mesh.topology.index_map(mesh.topology.dim).size_local
+                + mesh.topology.index_map(mesh.topology.dim).num_ghosts
+            )
+            vertices = np.array([c_to_v.links(cell) for cell in range(num_cells)])
+            flat_vertices = np.concatenate(vertices)
+            cell_indices = np.repeat(np.arange(num_cells), [len(v) for v in vertices])
+            vertex_positions = np.concatenate([np.arange(len(v)) for v in vertices])
+            vertex_map[flat_vertices] = self.connectivities_dict[subdomain][
+                mesh.topology.original_cell_index
+            ][cell_indices, vertex_positions]
+            self.vertex_maps_dict[subdomain] = vertex_map
 
         function_space = dolfinx.fem.functionspace(mesh, element)
         u = dolfinx.fem.Function(function_space)
-
-        num_vertices = (
-            mesh.topology.index_map(0).size_local
-            + mesh.topology.index_map(0).num_ghosts
-        )
-        vertex_map = np.empty(num_vertices, dtype=np.int32)
-
-        # Get cell-to-vertex connectivity
-        c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
-        # Map the OF_mesh vertices to dolfinx_mesh vertices
-        num_cells = (
-            mesh.topology.index_map(mesh.topology.dim).size_local
-            + mesh.topology.index_map(mesh.topology.dim).num_ghosts
-        )
-        vertices = np.array([c_to_v.links(cell) for cell in range(num_cells)])
-        flat_vertices = np.concatenate(vertices)
-        cell_indices = np.repeat(np.arange(num_cells), [len(v) for v in vertices])
-        vertex_positions = np.concatenate([np.arange(len(v)) for v in vertices])
-
-        vertex_map[flat_vertices] = self.connectivities_dict[subdomain][
-            mesh.topology.original_cell_index
-        ][cell_indices, vertex_positions]
-
-        # Assign values in OF_mesh to dolfinx_mesh
-        assert hasattr(self.OF_meshes_dict[subdomain], "point_data")
         u.x.array[:] = (
-            self.OF_meshes_dict[subdomain].point_data[name][vertex_map].flatten()
+            self.OF_meshes_dict[subdomain]
+            .point_data[name][self.vertex_maps_dict[subdomain]]
+            .flatten()
         )
 
         return u