Adding save and load functions for dolfin functions.

fem_sol.py

 # Customized FEM solver for FRAP theory
 from dolfin import *
+import fem_utils
 import matplotlib.pyplot as plt
 import numpy as np
 import time
@@ -49,6 +50,7 @@ class frap_solver:
 
         self.cs = []  # Intensity profiles over time
 
+
     def solve_frap(self):
         '''
         Solver routine 
@@ -87,7 +89,8 @@ class frap_solver:
         
         ti = time.time()
         for i in range(self.T):
-            self.cs.append([Phi_u([x, self.cent_poin[1], self.cent_poin[2]])
+            fem_utils.save_time_point(Phi_u0,self.name+'t_p_'+str(i)+'.h5')
+            self.cs.append([Phi_u0([x, self.cent_poin[1], self.cent_poin[2]])
                             for x in np.linspace(4, 4.5, 100)])
             solve(Form == 0, Phi_u)
             assign(Phi_u0, Phi_u)

fem_utils.py

@@ -2,6 +2,7 @@
 Utilities for simplifying FEM modelling, e.g. Meshes.
 '''
 
+from dolfin import *
 import meshio
 import numpy as np
 
@@ -107,3 +108,19 @@ def convert_msh_to_xdmf(input_path, output_path):
                             cell_data={"name_to_read":[line_data]})
 
     meshio.write(output_path, tetra_mesh)
+
+
+def save_time_point(f, name):
+    ''' Save Fenics function to h5 file.'''
+    fFile = HDF5File(MPI.comm_world,name,"w")
+    fFile.write(f,"/f")
+    fFile.close()
+
+def load_time_point(name, mesh):
+    '''Load Fenics function from h5 file.'''
+    V = FunctionSpace(mesh, 'CG', 1)
+    f = Function(V)
+    fFile = HDF5File(MPI.comm_world, name, "r")
+    fFile.read(f, "/f")
+    fFile.close()
+    return f