Hi all,
I’d like to access time dependent functions from an XDFM File as in:
from fenics import *
mesh = UnitSquareMesh(20,20)
Q = FunctionSpace(mesh, "CG", 1)
F = Function(Q)
filename = "u.xdmf"
ufile = XDMFFile(filename)
F.rename("F", "label")
ufile.write(F, 0)
ufile.write(F, 1)
F0 = Function(Q)
# Something like :
# F0 = ufile.read("F", 0)
Would you know anything that could replace the last line above ?
Thank you all for your help!
Rem
You can use the ‘write_checkpoint’ and ‘read_checkpoint’ function, see: https://bitbucket.org/fenics-project/dolfin/raw/8ccef27841c14f75ba5f3e05c8661bd19f88098c/python/test/unit/io/test_XDMF.py
Hi, I deeply apologise for the very late answer I’ve been swamped for the last few months.
I tried to use the function read_checkpoint but for some reason it doesn’t work.
from fenics import *
mesh = UnitSquareMesh(2, 2)
V = FunctionSpace(mesh, 'P', 1)
ini_u = Expression("0", degree=2)
u = interpolate(ini_u, V)
#write
u.rename("u", "label")
ufile = XDMFFile("u.xdmf")
ufile.write_checkpoint(u, "u_out", 0)
#read
u1 = Function(V)
with XDMFFile(mesh.mpi_comm(), "u.xdmf") as file:
file.read_checkpoint(u1, "u", 0)
It raises this error:
XDMF file "u.xdmf" will be overwritten.
HDF file "u.h5" will be overwritten.
Traceback (most recent call last):
File "test.py", line 17, in <module>
file.read_checkpoint(u1, "u", 0)
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
*** fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to reading data from XDMF file.
*** Reason: Storage format "" is unknown.
*** Where: This error was encountered inside XDMFFile.cpp.
*** Process: 0
***
*** DOLFIN version: 2018.1.0
*** Git changeset: 948dc42cc4e06ed9227d0201ad50f94ac94cbf9f
*** -------------------------------------------------------------------------
Does it have something to do with the h5 format ?
Thanks for your help,
Rem
No, this is because you call
file.read_checkpoint(u1, "u", 0)
You have saved the data with the label u_out, so it should be
file.read_checkpoint(u1, "u_out", 0)
Sorry that was silly of me. Thanks !
It works now, but apparently I can’t write several times as the file is “still open”. Is there a particular argument that need to be passed to the function ?
Thanks anyway for your help !
from fenics import *
mesh = UnitSquareMesh(2, 2)
V = FunctionSpace(mesh, 'P', 1)
ini_u = Expression("0", degree=2)
u = interpolate(ini_u, V)
#write
with XDMFFile("u.xdmf") as file:
file.write_checkpoint(u, "u_out", 0)
#read
u1 = Function(V)
with XDMFFile("u.xdmf") as file:
file.read_checkpoint(u1, "u_out", 0)
with XDMFFile("u.xdmf") as file:
file.write_checkpoint(u1, "u_out", 1)
file.write_checkpoint(u1, "u_out", 2)
Using the other signature of write_checkpoint, with arguments encoding and append:
#write
with XDMFFile("u.xdmf") as file:
file.write_checkpoint(u, "u_out", 0, XDMFFile.Encoding.HDF5, False)
#read
u1 = Function(V)
with XDMFFile("u.xdmf") as file:
file.read_checkpoint(u1, "u_out", 0)
with XDMFFile("u.xdmf") as file:
file.write_checkpoint(u1, "u_out", 1, XDMFFile.Encoding.HDF5, True)
file.write_checkpoint(u1, "u_out", 2, XDMFFile.Encoding.HDF5, True)
I see thanks I struggle finding docs on these functions.
This does the trick cheers!
Any differences between the function write and write_checkpoint in the file itself ?
Hi how do we do this if u is a vector? For example:
solve(F == 0, u, bcs)
velocity, pressure = u.split()
with XDMFFile("Velocity_A.xdmf") as file:
file.write_checkpoint(velocity, "V_A", 0, XDMFFile.Encoding.HDF5, False)
vtkfile_u_a << velocity
Here velocity is vector and I’m using a mixed elements defined as:
P1 = VectorElement('P', mesh.ufl_cell(), 2) #Velocity is a vector
P2 = FiniteElement('P', mesh.ufl_cell(), 1) #Pressure p is a scalar
element = MixedElement([P1,P2])
V = FunctionSpace(mesh, element)
I then read the file as follows:
W = VectorFunctionSpace(mesh, 'P', 2) #VectorSpace for velocity fields
velocity_b = Function(W)
with XDMFFile("Velocity_A.xdmf") as file:
file.read_checkpoint(velocity_b, "V_A", 0)
vtkfile_w << velocity_b
When I compare the results between the vtkfile_w and vtkfile_u_a, the results are very different. I believe the error has to do with velocity being a vector? In that case how do I write and read the velocity solution? Thanks
VTK will only save the nodal values at the vertices (equivalent of interpolating the solution to a CG1 space) you can verify thisby looking at the point values in paraview . Xdmf with checkpointing actually saved the CG 2 function.
Sorry I don’t quite follow, should’t writing velocity as an xdmf file and reading the file to place it into another variable produce the same vtkfile? Essentially, velocity should equal velocity_b.
I intend to use the velocity field from the first solution as a variable in a subsequent problem to be solved.
Essentially what I need to do is:
if count == 0:
with XDMFFile(output_folder_name+output_vtk_folder+"Velocity_angioplasty.xdmf") as file:
file.write_checkpoint(velocity, "V_Angioplasty", 0, XDMFFile.Encoding.HDF5, False)
vtkfile_u_angio << velocity
vtkfile_p_angio << pressure
velocity_angio = velocity
elif count == 1:
with XDMFFile(output_folder_name+output_vtk_folder+"Velocity_post_angioplasty.xdmf") as file:
file.write_checkpoint(velocity, "V_Post_Angioplasty", 0, XDMFFile.Encoding.HDF5, False)
vtkfile_u_post_angio << velocity
vtkfile_p_post_angio << pressure
velocity_post_angio = velocity
Unfortunately, velocity_angio seems to get updated to the same value as velocity_post_angio when it goes through the for loop. I think the assignment is to a memory location of velocity. How do I go about doing this? Thanks
Could you post a minimal example (full code, that produces your issue)?
It does not throw any error messages. This is the main block:
for count in range (2):
#Define Boundary Conditions
p_out = 5000
if count == 0:
p_inn = 100
elif count == 1:
p_in = 1000
boundary_conditions = {
1 : {"Dirichlet" : p_in},
2 : {"Dirichlet" : p_out}
}
bcs = []
for i in boundary_conditions:
if "Dirichlet" in boundary_conditions[i]:
#Here V.sub(1) refers to pressure while 0 refers to velocity
bc = DirichletBC(V.sub(1), boundary_conditions[i]["Dirichlet"], boundaries, i) #arguments: Function Space, value, subdomain marker, sudomain id
bcs.append(bc)
#Define variational form
F = dot(u1,v1)*dx + dot(v2,div(u1))*dx
solve(F == 0, u, bcs)
if count == 0:
vtkfile_u_angio << velocity
vtkfile_p_angio << pressure
velocity_angio = velocity # How do I do this?
elif count == 1:
vtkfile_u_post_angio << velocity
vtkfile_p_post_angio << pressure
velocity_post_angio = velocity # How do I do this?
The problem is after the for loop velocity_angio = velocity_post_angio
which should not be the case
Since this is a mixed space, I would use FunctionAssigner to assign the the velocity to the variable, see:
https://bitbucket.org/fenics-project/dolfin/raw/ec57db53f2b13f1768214829e8b4f80fc32205cf/python/demo/undocumented/sub-function-assignment/demo_sub-function-assignment.py
Thanks! That worked perfectly 