Various fixes for proper reading on input data

Author: jorgensd

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(

src/festim/initial_condition.py

@@ -240,7 +240,7 @@ 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 +255,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 +269,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

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,

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)
@@ -236,6 +236,7 @@ def f2(x):
                 filename=filename,
                 name="my_function1",
                 timestamp=0.2,
+                mesh=mesh
             ),
             species=H,
             volume=vol,
@@ -245,6 +246,7 @@ def f2(x):
                 filename=filename,
                 name="my_function2",
                 timestamp=0.3,
+                mesh=mesh
             ),
             species=D,
             volume=vol,

test/test_mesh.py

@@ -1,9 +1,11 @@
 import logging
 import os
-
+import pytest
 from mpi4py import MPI
 
-import ipyparallel as ipp
+
+ipp = pytest.importorskip("ipyparallel")
+
 import numpy as np
 import pytest
 from dolfinx import mesh as fenics_mesh