Saving a function space

Hi there,

is it possible in fenics to save/load a function space?

What do you imply when you say saving/loading a function space? What data structures would you expect to write/load from file?

1 Like

Hello everyone,

I am using the HDF5 file format to save all relevant data (functions/mesh/meshfunctions) + metadata in a self-contained manner. For example,

from fenics import *

mesh = UnitSquareMesh(4,4)

def WriteH5Sample():
    V = FunctionSpace(mesh,"CG",1)
    u = project(Expression(("x[0]*x[1]"), degree=1), V)

    with HDF5File(MPI.comm_world, 'test.h5', file_mode='w') as hdf5file:
        hdf5file.write(u, 'f')    

WriteH5Sample() # Write out a simple H5 file

I can then read back in the function using

def ReadH5Func(V):
    with HDF5File(MPI.comm_world, 'test.h5', file_mode='r') as hdf5file:
     #  V = ???
        u = Function(V)
        hdf5file.read(u, 'f')
    return u

# Define FunctionSpace. 
# Can this be automatically interpreted from the data?
newV = FunctionSpace(mesh, 'CG', 1)
u = ReadH5Func(newV)

The problem is that we must decide on a FunctionSpace newV prior to reading the data. ReadH5(V) takes V as an argument, so it is not really self-contained.

How can I go about “saving the function space” so that the appropriate V is generated automatically? Perhaps I can just save information from str(V) and use it to create newV?

Thanks for any advice!

This is my attempt at solving my own problem. It is admittedly clunky, but I wanted it to be robust enough to recreate function spaces (scalar/vector/tensor) as well as mixed function spaces from a small piece of metadata.

If anyone would like to revise my code, point out any flaws, or suggest any improvements that would be great.

from fenics import *

mesh = UnitSquareMesh(4,4)

Q = FunctionSpace(mesh, "DG", 2)
V = VectorFunctionSpace(mesh, "CG", 1)

# We will save info about W and use it later to recreate the equivalent functionspace _W.
W = MixedFunctionSpace(Q,V)
u = Function(W)

def FunctionSpaceInfo(V):
    E = V.ufl_element()
    info = {'family':E.family(),'cell':E.cell(),'degree':E.degree(),
                'num_sub_elements':E.num_sub_elements(),'symmetry':E.symmetry(),
                'value_shape':E.value_shape(),'value_size':E.value_size()}
    return info

def SaveFunctionSpace(V):
    saveinfo = []
    if V.num_sub_spaces()>0:
        for V_i in (V.sub_spaces()):
            saveinfo.append(FunctionSpaceInfo(V_i))
    else:
        saveinfo.append(FunctionSpaceInfo(V))      
    return saveinfo

def RecreateFunctionSpace(info):
    FS_list = []
    for FS in info:
        if FS['value_size']==1: # scalar
            FS_list.append(      FunctionSpace(mesh,FS['family'],FS['degree']))
        if FS['value_size']==2: # vector
            FS_list.append(VectorFunctionSpace(mesh,FS['family'],FS['degree']))
        if FS['value_size']==4: # tensor
            FS_list.append(TensorFunctionSpace(mesh,FS['family'],FS['degree']))           
    return FS_list
# Extract info from FunctionSpace(s)
saveinfo = SaveFunctionSpace(W)

# Recreate FunctionSpace
FS_list = RecreateFunctionSpace(saveinfo)
_W = FS_list[0] if len(saveinfo)==1 else MixedFunctionSpace(*FS_list)
_u = Function(_W)

# This is just a reminder to split a Function defined in a MixedFunctionSpace
if _u.num_sub_spaces()>1: print(_u.sub(0).name())
else: print(_u.name())