168
views
0
4 months ago by
I have an overloaded Expression class which I have written in order to apply an interpolation of an image for material properties. I am trying to run a sensitivity analysis on the model, and therefore need to make my function solves as time effective as possible. It takes ~10 minutes to interpolate the image onto the mesh, whereas the FEM itself takes ~2-3 minutes. I need the scipy package for interpolation, but on the server, scipy is not accessible. Therefore, I need to do the interpolation on my personal machine, save the results, and then open the Expression and FunctionSpace on the server, where it is ready launch multiple jobs.

Is there a way to save a FunctionSpace and Expression variable, and then reopen them on another machine? I know how to do it with the mesh, but I have not had luck with the FunctionSpace and Expression.

from scipy.interpolate import RegularGridInterpolator

class InterpolatedParameter(Expression):
def __init__(self,X,Y,Z,param_map,**kwagrs):
self.X=X   #1D Coordinate map for 3D image
self.Y=Y
self.Z=Z
self.param_map=param_map #3D image

def eval_cell(selv,values,x,cell):
interp_fcn = RegularGridInterpolator((self.X,self.Y,self.Z),self.param_map)
values[0] = interp_fcn(x)

V = FunctionSpace(mesh,'P',1)
u = TrialFunction(V)
v = TestFunction(V)

#Instantiate the interpolator

#Perform the interpolation onto the FunctionSpace
#This takes a relatively long time

###Save the Expression and Functionspace for another python script to open up
File('my_mesh.xml') << mesh
#File('my_FunctionSpace') << V
_______________________________________________________________________________

###Load the Expression and FunctionSpace created on a different machine
File('my_mesh.xml') >> mesh
#File('my_FunctionSpace') >> V
#File('my_Expression') >> ADC​

Community: FEniCS Project
Is interpolate your own function or is it fenics.interpolate? Are you sure that it returns a fenics.Expression? In my project, I use fenics.interpolate to make a fenics.Function.
written 4 months ago by Alexander G. Zimmerman
interpolate is the fenics function. The class I wrote inherits from Expression, so that's why I thought after interpolating it would still be an Expression. I was wrong about the type it returned. Since its a Function, Ruben's solution seems to work just fine.
written 4 months ago by Ryan W

4
4 months ago by
Hi Ryan,

I think that the only thing missing in your code is the extension, for your purpose save all in .xml just as the mesh, try it an let us know
File('my_Expression.xml') << V

#and later

File('my_Expression.xml') >> V​
1
Really, Downvoted for trying to help?
written 4 months ago by Ruben Gonzalez
Please, post only solution that you have tried out and it worked. If you post your (wrong) guesses you might seriously confuse other users.
written 4 months ago by Michal Habera
So, when I have another script open the mesh file, that mesh cannot be used. When I try and run V=FunctionSpace(mesh) in a fresh script, I get the following error:

Illegal argument, not a mesh: <module 'dolfin.mesh' from '/usr/lib/python2.7/dist-packages/dolfin/mesh/__init__.pyc'>.

When I run File(''my_Expression.xml'') >> ADC, I get the following error:

NameError: name 'ADC' is not defined

written 4 months ago by Ryan W
1
The way that I've tried this approach is

#to load the mesh
mesh = Mesh("my_mesh.xml")

#I do not think you need to save and load the FunctionSpace
V = FunctionSpace(mesh,'P',1)

#to load the saved parameters
File("my_expression.xml")>>ADC​
written 4 months ago by Ruben Gonzalez
I always kick myself when its one line I'm missing. Thanks for the help!
written 4 months ago by Ryan W
2
4 months ago by
Based on my understanding of your requirements, you should only need to save the mesh and a fenics.Function where you have interpolated the field (with fenics.interpolate, which returns a fenics.Function).

If so, then you could use the fenics.XDMFFile class with XDMFFile.write_checkpoint and XDMFFile.read_checkpoint. These checkpoint methods might be a new feature, so here's a discussion about what to do if these methods aren't available: https://www.allanswered.com/post/aqox/#pzew
I accepted the other answer because it addressed the main problem I was having. This will be helpful as, I have a 24 hr run time. Does this also work in mpirun/parallel? Or do I need to try out the HDF5 format?
written 4 months ago by Ryan W