How to save solutions in a loop in xml files and name the xml files using the loop index or read solutions from vtk or pvd files?

General
1

I am solving a hyperelastic problem. I have a domain on which I am applying pressure to minimize an energy functional. I increase the pressure using a loop. I want to save all solutions as separate xml files named with the loop index. For example, for loop i, the xml file will be solution_i.xml. I tried File('lumsoln%(d).xml'%(i))<<u2 in the below MWE, which returns the error:

File('lumsoln%(d).xml'%(i))<<u2
TypeError: format requires a mapping

What “mapping” are we talking about here and how to fix this?

Also, is there a way to read solutions from pvd or vtk files?

Let’s consider the following MWE:

msh = meshio.read("mesh.msh")

triangle_mesh = create_mesh(msh, "triangle", True)
line_mesh = create_mesh(msh, "line", True)
meshio.write("mesh.xdmf", triangle_mesh)
meshio.write("mf.xdmf", line_mesh) 
from dolfin import *
mesh = Mesh()
xdmf = XDMFFile(mesh.mpi_comm(),"mesh.xdmf")
xdmf.read(mesh)
mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim())
with XDMFFile("mesh.xdmf") as infile:
   infile.read(mvc, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc)
xdmf.close()

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile("mf.xdmf") as infile:
    infile.read(mvc, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc)

ds_custom = Measure("dS", domain=mesh, subdomain_data=mf)                                                               
dx_custom = Measure("dx", domain=mesh, subdomain_data=cf)

#************************************************************************************

combined_subdomains_2 = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
combined_subdomains_2.array()[cf.array()==5] = 1  # Assigning same number for all the submeshes for compiling them into one submesh
combined_subdomains_2.array()[cf.array()==6] = 1
combined_subdomains_2.array()[cf.array()==7] = 1
mesh_ima_ref = SubMesh(mesh, combined_subdomains_2, 1)

boundary_marker_ima_ref = MeshFunction("size_t", mesh_ima_ref, mesh_ima_ref.topology().dim()-1, 0)
ncells_ref = MeshFunction("size_t", mesh_ima_ref, mesh_ima_ref.topology().dim())
surface_marker_ima_ref = MeshFunction("size_t", mesh_ima_ref, mesh_ima_ref.topology().dim(), 0)
vmap_ref = mesh_ima_ref.data().array("parent_vertex_indices", 0)
cmap_ref = mesh_ima_ref.data().array("parent_cell_indices", mesh_ima_ref.topology().dim())

        
n = 0
for c in cells(mesh_ima_ref):
  parent_cell_ima_ref = Cell(mesh, cmap_ref[c.index()])
  surface_marker_ima_ref.array()[c.index()] = cf.array()[parent_cell_ima_ref.index()]
  for f in facets(parent_cell_ima_ref):
    for g in facets(c):
      g_vertices = vmap_ref[g.entities(0)]
      if set(f.entities(0)) == set(g_vertices):
        boundary_marker_ima_ref.array()[g.index()] = mf.array()[f.index()]
    n=n+1

ds_ima_ref = Measure("ds", domain=mesh_ima_ref, subdomain_data=boundary_marker_ima_ref) # Defining global boundary measure, for global surface measure, plug in "dS" instead of "ds",                                                                                              which is used for subdomain boundary measures
dx_ima_ref = Measure("dx", domain=mesh_ima_ref, subdomain_data=surface_marker_ima_ref)

V_ima = VectorFunctionSpace(mesh_ima_ref, 'P', 1)      #   Changed to second order
W_ima_ref = FunctionSpace(mesh_ima_ref, 'P', 1)
T_ima = TensorFunctionSpace(mesh_ima_ref, 'P', 1)  

initial_int_area = assemble(Constant(1.0)*dx_ima_ref(5))
initial_med_area = assemble(Constant(1.0)*dx_ima_ref(6))
initial_ad_area = assemble(Constant(1.0)*dx_ima_ref(7))

print('initial intima area =', assemble(Constant(1.0)*dx_ima_ref(5)))
print('initial media area =', assemble(Constant(1.0)*dx_ima_ref(6)))
print('initial adventitia area =', assemble(Constant(1.0)*dx_ima_ref(7)))


#***********************************************************************************************
ug1 = Function(V_ima)
vg1 = TestFunction(V_ima)

# Elasticity parameters 
nu = 0.4999   
mu1, mu2, mu3 = 27.9, 1.27, 7.56 


ds = 0.001

NumStepg = 51

for j in range(0,NumStepg):          
 
 tj = j*ds 
 
 dg = mesh_ima_ref.geometry().dim()
 Ig = Identity(dg)             # Identity tensor
 Fg = Ig + grad(ug1)             # Deformation gradient
                     
 g_1 = Expression((('1.0+p1*t1','0.0','0.0'),
                     ('0.0','1.0+p2*t1','0.0'),
                     ('0.0','0.0','1.0')),degree=1,p1=1.0,p2=1.0,t1=tj) 

 ginv = Expression((('1.0/(1.0+p1*t1)','0.0','0.0'),
                     ('0.0','1.0/(1.0+p2*t1)','0.0'),
                     ('0.0','0.0','1.0')),degree=1,p1=1.0,p2=1.0,t1=tj) 
                     
                     
 g_i = project(det(g_1),W_ima_ref)

 Ce = (Fg*inv(g_1)).T*(Fg*inv(g_1))                   # Right Cauchy-Green tensor for the Intima
 
 Fe = Fg*ginv
 
 Mg = cofac(Fg)
   
 # Invariants of deformation tensors
 I1 = tr(Ce)     # Ic for the Intima
 Jg_e = det(Fg)

 psig1 = (mu1/2)*(I1 - 3) + ((nu)*(mu1)/(1-2*nu))*(Jg_e-1)**2 - mu1*(ln(Jg_e)) 
 psig2 = (mu2/2)*(I1 - 3) + ((nu)*(mu2)/(1-2*nu))*(Jg_e-1)**2 - mu2*(ln(Jg_e))
 psig3 = (mu3/2)*(I1 - 3) + ((nu)*(mu3)/(1-2*nu))*(Jg_e-1)**2 - mu3*(ln(Jg_e)) 
   
 Pig = (det(g_1))*psig1*dx_ima_ref(5) + psig2*dx_ima_ref(6) + psig3*dx_ima_ref(7) 
 
 # Compute first variation of Pi (directional derivative about u in the direction of v)
 Fg1 = derivative(Pig, ug1, vg1)
 
 bcgs_ima = []
 
 ug1.set_allow_extrapolation(True)
 
 # Solving variational problem
 
 print(f"Solving for IMA displacement field at time T={tj}", flush=True)
 
 solve(Fg1 == 0, ug1, bcgs_ima, solver_parameters={"newton_solver": {"relative_tolerance": 1e-8,
      "absolute_tolerance": 1e-8,"maximum_iterations": 200}})
  
 ug1.set_allow_extrapolation(True)
##########################################################################################

 ug_ima = Function(V_ima)
 
 vg1 = TestFunction(V_ima)
 
 ug1.set_allow_extrapolation(True)
 ug_ima.assign(ug1)
 ug1 = interpolate(ug_ima, V_ima)
 ug1.set_allow_extrapolation(True)
 
 File('lumsoln%(d).xml'%(j))<<u2
 File('imasoln%(d).xml'%(j))<<u1

Any suggestions would be extremely helpful and would be greatly appreciated.

Thank you.

2

Dear Avishmj,

Please note that your code is far from a minimal example.
A minimal example would be something along the lines of:

from dolfin import *

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "Lagrange", 1)
u = Function(V)

for i in range(3):
    u.vector()[:] = i
    File(f"solution{i}.xml") << u

which resolves your problem.

Note that by creating small, standalone examples that highlight the functionality you are missing, rather than taking a full example, vastly increases the chances for either finding the bug in your own code, or anyone else being able to help you.

No, you would have to use HDF5File.write or XDMFFile.write_checkpoint for this.

3

Thanks, that works.

I have a follow-up question:

Suppose I save the solutions in different xml files. Is there a way to extract the solutions from the xml files and write them into one PVD file? Here is the algorithm of what I mean:

First, I save the solutions into different xml files.

from dolfin import *

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "Lagrange", 1)
u = Function(V)

for i in range(3):
    u.vector()[:] = i
    File(f"solution{i}.xml") << u

Now I want to extract the solutions and combine them into one PVD file like

from dolfin import *

mesh = UnitSquareMesh(10, 10)
V = FunctionSpace(mesh, "Lagrange", 1)
u = Function(V)
v = Function(V)

for i in range(3):
    u.vector()[:] = i
    File(f"solution{i}.xml") << u

vtk = File("allsolutions.pvd")

for i in range(3):
    v = Function(V, "solution{i}.xml")
    vtk<<(v,i)

Will the above code generate a pvd file with all the solutions u?