diff --git a/.github/workflows/ci_docker.yml b/.github/workflows/ci_docker.yml
index 1506b0031..ce787ccfd 100644
--- a/.github/workflows/ci_docker.yml
+++ b/.github/workflows/ci_docker.yml
@@ -11,19 +11,28 @@ jobs:
         container_version: [v0.9.0, v0.10.0, nightly]
     container: dolfinx/dolfinx:${{ matrix.container_version }}
     steps:
-        - name: Checkout code
-          uses: actions/checkout@v2
-      
-        - name: Install local package and dependencies
-          run: |
-            pip install .[test]
+      - name: Checkout code
+        uses: actions/checkout@v6
 
-        - name: Run tests
-          run: |
-            python3 -m pytest test/ --cov festim --cov-report xml --cov-report term
+      - name: Install local package and dependencies
+        run: |
+          pip install .[test]
 
-        - name: Upload to codecov
-          uses: codecov/codecov-action@v4
-          with:
-            token: ${{ secrets.CODECOV_TOKEN }}
-            files: ./coverage.xml
\ No newline at end of file
+      - name: Overload adios4dolfinx
+        if: ${{ matrix.container_version == 'nightly' }}
+        run: python3 -m pip install git+https://github.com/jorgensd/adios4dolfinx.git
+
+      - name: Run tests
+        run: |
+          python3 -m pytest test/ --cov festim --cov-report xml --cov-report term
+
+      - name: Install curl
+        run: |
+          apt-get update
+          apt-get install -y curl
+
+      - name: Upload to codecov
+        uses: codecov/codecov-action@v5
+        with:
+          token: ${{ secrets.CODECOV_TOKEN }}
+          files: ./coverage.xml
diff --git a/src/festim/exports/surface_flux.py b/src/festim/exports/surface_flux.py
index fb3cb8f78..fe70e8f8c 100644
--- a/src/festim/exports/surface_flux.py
+++ b/src/festim/exports/surface_flux.py
@@ -47,12 +47,8 @@ def compute(
             entity_maps: entity maps relating parent mesh and submesh
         """
 
-        # obtain mesh normal from field
-        # if case multispecies, solution is an index, use sub_function_space
-        if isinstance(u, ufl.indexed.Indexed):
-            mesh = self.field.sub_function_space.mesh
-        else:
-            mesh = u.function_space.mesh
+        # obtain mesh normal from integration domain
+        mesh = ds.ufl_domain()
         n = ufl.FacetNormal(mesh)
 
         self.value = assemble_scalar(
diff --git a/src/festim/hydrogen_transport_problem.py b/src/festim/hydrogen_transport_problem.py
index 805511d4a..d4b7e47aa 100644
--- a/src/festim/hydrogen_transport_problem.py
+++ b/src/festim/hydrogen_transport_problem.py
@@ -13,12 +13,9 @@
 import tqdm.autonotebook
 import ufl
 from dolfinx import fem
-from dolfinx.nls.petsc import NewtonSolver
 from packaging.version import Version
 from scifem import BlockedNewtonSolver
 
-import festim.boundary_conditions
-import festim.problem
 from festim import (
     boundary_conditions,
     exports,
@@ -37,12 +34,12 @@
 )
 from festim.advection import AdvectionTerm
 from festim.helpers import (
+    KSPMonitor,
+    SnesMonitor,
     as_fenics_constant,
+    convergenceTest,
     get_interpolation_points,
     is_it_time_to_export,
-    SnesMonitor,
-    KSPMonitor,
-    convergenceTest,
     nmm_interpolate,
 )
 from festim.material import SolubilityLaw
@@ -433,7 +430,7 @@ def initialise_exports(self):
                             self.settings.stepsize.milestones.append(time)
                     self.settings.stepsize.milestones.sort()
 
-                if isinstance(export, festim.VTXTemperatureExport):
+                if isinstance(export, exports.VTXTemperatureExport):
                     self._temperature_as_function = (
                         self._get_temperature_field_as_function()
                     )
@@ -819,7 +816,7 @@ def create_formulation(self):
 
         for reaction in self.reactions:
             for reactant in reaction.reactant:
-                if isinstance(reactant, festim.species.Species):
+                if isinstance(reactant, _species.Species):
                     self.formulation += (
                         reaction.reaction_term(self.temperature_fenics)
                         * reactant.test_function
@@ -976,7 +973,7 @@ def post_processing(self):
                     else:
                         export.writer.write(float(self.t))
                 elif (
-                    isinstance(export, festim.VTXTemperatureExport)
+                    isinstance(export, exports.VTXTemperatureExport)
                     and self.temperature_time_dependent
                 ):
                     self._temperature_as_function.interpolate(
@@ -1422,14 +1419,14 @@ def create_subdomain_formulation(self, subdomain: _subdomain.VolumeSubdomain):
             if reaction.volume != subdomain:
                 continue
             for species in reaction.reactant + reaction.product:
-                if isinstance(species, festim.species.Species):
+                if isinstance(species, _species.Species):
                     # TODO remove
                     # temporarily overide the solution to the one of the subdomain
                     species.solution = species.subdomain_to_solution[subdomain]
 
             # reactant
             for reactant in reaction.reactant:
-                if isinstance(reactant, festim.species.Species):
+                if isinstance(reactant, _species.Species):
                     form += (
                         reaction.reaction_term(self.temperature_fenics)
                         * reactant.subdomain_to_test_function[subdomain]
@@ -2093,7 +2090,7 @@ def get_concentrations(species_list) -> list:
         # add reaction term to formulation
         # reactant
         for reactant in reaction.reactant:
-            if isinstance(reactant, festim.species.Species):
+            if isinstance(reactant, _species.Species):
                 self.formulation += (
                     reaction_term * reactant.test_function * self.dx(reaction.volume.id)
                 )
@@ -2123,7 +2120,7 @@ def override_post_processing_solution(self):
                 K_S0.x.array[entities] = subdomain.material.get_K_S_0(spe)
                 E_KS.x.array[entities] = subdomain.material.get_E_K_S(spe)
 
-            K_S = K_S0 * ufl.exp(-E_KS / (festim.k_B * self.temperature_fenics))
+            K_S = K_S0 * ufl.exp(-E_KS / (k_B * self.temperature_fenics))
 
             theta = spe.solution
 
@@ -2145,7 +2142,7 @@ def update_post_processing_solutions(self):
                 continue
             spe.post_processing_solution.interpolate(spe.dg_expr)
 
-    def create_dirichletbc_form(self, bc: festim.FixedConcentrationBC):
+    def create_dirichletbc_form(self, bc: boundary_conditions.FixedConcentrationBC):
         """Creates a dirichlet boundary condition form
 
         Args:
@@ -2164,7 +2161,7 @@ def create_dirichletbc_form(self, bc: festim.FixedConcentrationBC):
             K_S0.x.array[entities] = subdomain.material.get_K_S_0(bc.species)
             E_KS.x.array[entities] = subdomain.material.get_E_K_S(bc.species)
 
-        K_S = K_S0 * ufl.exp(-E_KS / (festim.k_B * self.temperature_fenics))
+        K_S = K_S0 * ufl.exp(-E_KS / (k_B * self.temperature_fenics))
 
         # create value_fenics
         bc.create_value(
diff --git a/src/festim/initial_condition.py b/src/festim/initial_condition.py
index 06f6a9269..7153dffc4 100644
--- a/src/festim/initial_condition.py
+++ b/src/festim/initial_condition.py
@@ -240,7 +240,12 @@ def create_expr_fenics(
 
 
 def read_function_from_file(
-    filename: str, name: str, timestamp: int | float, family="P", order: int = 1
+    filename: str,
+    name: str,
+    timestamp: int | float,
+    family="P",
+    order: int = 1,
+    mesh: dolfinx.mesh.Mesh | None = None,
 ) -> fem.Function:
     """
     Read a function from a file
@@ -255,6 +260,7 @@ def read_function_from_file(
         timestamp: the timestamp of the function
         family: the family of the function space
         order: the order of the function space
+        mesh: Mesh to create input space on.
 
     Returns:
         the function
@@ -268,4 +274,14 @@ def read_function_from_file(
         name=name,
         time=timestamp,
     )
-    return u_in
+    if mesh is None:
+        return u_in
+    else:
+        V = fem.functionspace(mesh, (family, order))
+        u = fem.Function(V)
+        num_cells = mesh.topology.index_map(mesh.topology.dim).size_local
+        cells = np.arange(num_cells, dtype=np.int32)
+        padding = 1e2 * np.finfo(mesh.geometry.x.dtype).eps
+        idata = fem.create_interpolation_data(V, V_in, cells=cells, padding=padding)
+        u.interpolate_nonmatching(u_in, cells, idata)
+        return u
diff --git a/test/system_tests/test_io.py b/test/system_tests/test_io.py
index 0d17d3755..ffae401a8 100644
--- a/test/system_tests/test_io.py
+++ b/test/system_tests/test_io.py
@@ -74,6 +74,7 @@ def test_writing_and_reading_of_species_function_using_checkpoints(tmpdir):
                 filename=tmpdir + "/out_checkpoint.bp",
                 name="H",
                 timestamp=10,
+                mesh=mesh,
             ),
             species=H,
             volume=vol,
@@ -83,6 +84,7 @@ def test_writing_and_reading_of_species_function_using_checkpoints(tmpdir):
                 filename=tmpdir + "/out_checkpoint.bp",
                 name="D",
                 timestamp=10,
+                mesh=mesh,
             ),
             species=D,
             volume=vol,
diff --git a/test/system_tests/test_misc.py b/test/system_tests/test_misc.py
index 3170c91f9..42f4c630c 100644
--- a/test/system_tests/test_misc.py
+++ b/test/system_tests/test_misc.py
@@ -186,7 +186,7 @@ def test_temp_dependent_bc_mixed_domain_temperature_as_function():
     """Test to catch bug 986"""
     mesh, mt, ct = generate_mesh(8)
 
-    V = dolfinx.fem.functionspace(mesh, ("CG", 1))
+    V = dolfinx.fem.functionspace(mesh, ("Lagrange", 1))
     T = dolfinx.fem.Function(V)
     T.interpolate(lambda x: 800.0 + 100 * x[0] + x[1])
 
diff --git a/test/test_advection.py b/test/test_advection.py
index 8b5cca117..88ed2f857 100644
--- a/test/test_advection.py
+++ b/test/test_advection.py
@@ -9,7 +9,7 @@
 import festim as F
 
 test_mesh_1d = F.Mesh1D(vertices=np.linspace(0, 1, 100))
-test_functionspace = fem.functionspace(test_mesh_1d.mesh, ("CG", 1))
+test_functionspace = fem.functionspace(test_mesh_1d.mesh, ("Lagrange", 1))
 
 
 @pytest.mark.parametrize(
diff --git a/test/test_h_transport_problem.py b/test/test_h_transport_problem.py
index 71dbcda84..68ce04b48 100644
--- a/test/test_h_transport_problem.py
+++ b/test/test_h_transport_problem.py
@@ -1366,7 +1366,7 @@ def test_traps_with_CG_elements():
         species=[trap],
     )
 
-    my_model._element_for_traps = "CG"
+    my_model._element_for_traps = "Lagrange"
 
     my_model.create_species_from_traps()
     my_model.define_function_spaces(element_degree=1)
@@ -1486,7 +1486,7 @@ def test_temperature_as_function_in_discontinuous():
     H = F.Species("H", subdomains=my_model.volume_subdomains)
     my_model.species = [H]
 
-    V = fem.functionspace(my_model.mesh.mesh, ("CG", 1))
+    V = fem.functionspace(my_model.mesh.mesh, ("Lagrange", 1))
     u = fem.Function(V)
     u.interpolate(lambda x: 100 + x[0])
     u.x.array[:] = 100
diff --git a/test/test_initial_condition.py b/test/test_initial_condition.py
index c8154f12a..6ced0b703 100644
--- a/test/test_initial_condition.py
+++ b/test/test_initial_condition.py
@@ -173,7 +173,7 @@ def f(x):
     my_problem.subdomains = [vol]
 
     function_initial_value = F.read_function_from_file(
-        filename=filename, name="my_function", timestamp=0.2
+        filename=filename, name="my_function", timestamp=0.2, mesh=mesh
     )
     my_problem.initial_conditions = [
         F.InitialConcentration(value=function_initial_value, species=H, volume=vol)
@@ -233,18 +233,14 @@ def f2(x):
     my_problem.initial_conditions = [
         F.InitialConcentration(
             value=F.read_function_from_file(
-                filename=filename,
-                name="my_function1",
-                timestamp=0.2,
+                filename=filename, name="my_function1", timestamp=0.2, mesh=mesh
             ),
             species=H,
             volume=vol,
         ),
         F.InitialConcentration(
             value=F.read_function_from_file(
-                filename=filename,
-                name="my_function2",
-                timestamp=0.3,
+                filename=filename, name="my_function2", timestamp=0.3, mesh=mesh
             ),
             species=D,
             volume=vol,
@@ -460,7 +456,7 @@ def test_initial_condition_mixed_domain():
     my_model.temperature = 300
 
     intial_cond_value = 100
-    V = fem.functionspace(my_model.mesh.mesh, ("CG", 1))
+    V = fem.functionspace(my_model.mesh.mesh, ("Lagrange", 1))
     u = fem.Function(V)
     u.x.array[:] = intial_cond_value
 
@@ -507,11 +503,11 @@ def test_initial_condition_mixed_domain_multispecies():
 
     intial_cond_value_1 = 100
     intial_cond_value_2 = 200
-    V = fem.functionspace(my_model.mesh.mesh, ("CG", 1))
+    V = fem.functionspace(my_model.mesh.mesh, ("Lagrange", 1))
     u = fem.Function(V)
     u.x.array[:] = intial_cond_value_1
 
-    V2 = fem.functionspace(my_model.mesh.mesh, ("CG", 1))
+    V2 = fem.functionspace(my_model.mesh.mesh, ("Lagrange", 1))
     u2 = fem.Function(V2)
     u2.x.array[:] = intial_cond_value_2
 
diff --git a/test/test_mesh.py b/test/test_mesh.py
index ee04351a4..f00872ef6 100644
--- a/test/test_mesh.py
+++ b/test/test_mesh.py
@@ -3,7 +3,10 @@
 
 from mpi4py import MPI
 
-import ipyparallel as ipp
+import pytest
+
+ipp = pytest.importorskip("ipyparallel")
+
 import numpy as np
 import pytest
 from dolfinx import mesh as fenics_mesh