I need to store a solution of expensive FEM calculation to use it in further analysis. Browsing through tutorials I have so far discovered, that I may store my results like this:
from fenics import *
mesh = Mesh('mesh/UnitSquare8x8.xml')
V = FunctionSpace(mesh, 'P', 1)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)
def boundary(x, on_boundary):
return on_boundary
bc = DirichletBC(V, u_D, boundary)
A = assemble(a)
b = assemble(L)
bc.apply(A, b)
u = Function(V)
solver = KrylovSolver("cg", "ilu")
solver.solve(A, u.vector(), b)
File('solution.xml') << u.vector()
and later load them like this:
from fenics import *
mesh = Mesh('mesh/UnitSquare8x8.xml')
V = FunctionSpace(mesh, 'P', 1)
u = Function(V)
File('solution.xml') >> u.vector()
Unfortunately I hardly know what exactly I am doing here. Is that a proper way of storing and loading calculated results? Is order of elements in u.vector() (for the same mesh file) fixed within/between different FEniCS versions, or it is just an implementation detail which may change any time? If it is unsafe, then what is the proper way of doing so?
I have found another (possibly even more dangerous) solution. I may use VALUES = u.vector().get_local() and u.vector().set_local(VALUES) methods, as VALUES is a numpy array which I may easily store and load.
