diff --git a/src/io4dolfinx/readers.py b/src/io4dolfinx/readers.py
index 38150a3..54450a6 100644
--- a/src/io4dolfinx/readers.py
+++ b/src/io4dolfinx/readers.py
@@ -382,6 +382,18 @@ def create_geometry_function_space(mesh: dolfinx.mesh.Mesh, N: int) -> dolfinx.f
     return dolfinx.fem.FunctionSpace(mesh, ufl_el, cpp_space)
 
 
+def _function_dtype(data_dtype: np.dtype, mesh: dolfinx.mesh.Mesh) -> np.dtype:
+    """The dtype of a function holding data read from file.
+
+    Real data takes the precision of the mesh (a function space only allows its mesh's real type,
+    e.g. float32 data on a float64 mesh), complex data the complex type of that precision.
+    """
+    real = np.dtype(mesh.geometry.x.dtype)
+    if np.issubdtype(data_dtype, np.complexfloating):
+        return np.result_type(real, np.complex64)
+    return real
+
+
 def read_point_data(
     filename: Path | str,
     name: str,
@@ -421,7 +433,7 @@ def read_point_data(
 
     # Create appropriate function space (based on coordinate map)
     V = create_geometry_function_space(mesh, num_components)
-    uh = dolfinx.fem.Function(V, name=name, dtype=dataset.dtype)
+    uh = dolfinx.fem.Function(V, name=name, dtype=_function_dtype(dataset.dtype, mesh))
     # Assume that mesh is first order for now
     x_dofmap = fenicsx_compat.dofmap(mesh)
     igi = np.array(mesh.geometry.input_global_indices, dtype=np.int64)
@@ -491,7 +503,7 @@ def read_cell_data(
     else:
         shape = (num_components,)
     V = dolfinx.fem.functionspace(mesh, ("DG", 0, shape))
-    u = dolfinx.fem.Function(V, dtype=dofs.dtype)
+    u = dolfinx.fem.Function(V, dtype=_function_dtype(dofs.dtype, mesh))
     data_array = u.x.array.reshape(-1, num_components)
     for i in range(num_components):
         local_entities, local_values = dolfinx.io.distribute_entity_data(
diff --git a/tests/test_pyvista.py b/tests/test_pyvista.py
index 246e2e6..531781a 100644
--- a/tests/test_pyvista.py
+++ b/tests/test_pyvista.py
@@ -48,3 +48,31 @@ def test_read_mesh_and_cell_data(tmp_path):
             np.testing.assert_allclose(pd.x.array[:], grid.point_data[name][igi])
         else:
             raise RuntimeError(f"Could not find {name} in grid")
+
+
+def test_single_precision_data_on_double_precision_mesh(tmp_path):
+    # A file with float32 points and data, read into a float64 mesh (backend_args dtype): the point
+    # and cell data must come back as functions on that mesh, in its precision
+    tmp_path = MPI.COMM_WORLD.bcast(tmp_path, root=0)
+    filename = tmp_path / "grid32.vtu"
+    grid = pyvista.examples.load_hexbeam().triangulate()
+    grid.points = grid.points.astype(np.float32)
+    grid.point_data["height"] = grid.points[:, 2].astype(np.float32)
+    grid.cell_data["label"] = np.arange(grid.n_cells, dtype=np.float32) % 3
+    if MPI.COMM_WORLD.rank == 0:
+        grid.save(filename)
+    MPI.COMM_WORLD.barrier()
+
+    mesh = io4dolfinx.read_mesh(
+        filename, MPI.COMM_WORLD, backend="pyvista", backend_args={"dtype": np.float64}
+    )
+    assert mesh.geometry.x.dtype == np.float64
+
+    height = io4dolfinx.read_point_data(filename, "height", mesh, backend="pyvista")
+    assert height.x.array.dtype == np.float64
+    x = height.function_space.tabulate_dof_coordinates()
+    assert np.allclose(height.x.array, x[:, 2], atol=1e-6)
+
+    label = io4dolfinx.read_cell_data(filename, "label", mesh, backend="pyvista")
+    assert label.x.array.dtype == np.float64
+    assert set(np.unique(label.x.array)) <= {0.0, 1.0, 2.0}