Commit b953997a authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

close all files in vtk.py and mesh.py

parent e3e92320
Pipeline #5435 passed with stage
in 53 seconds
......@@ -60,18 +60,17 @@ def _read_array(a) :
def _write_index(idxname,filename,i,t,reset_index) :
if reset_index is None : reset_index = (i==0 or not os.path.isfile(idxname))
if reset_index :
f = open(idxname,"wb")
f.write(b'<?xml version="1.0"?>\n')
f.write(b'<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">\n')
f.write(b'<Collection>\n')
else:
f = open(idxname,"rb+")
f.seek(-25,os.SEEK_END)
f.write(b'<DataSet timestep="%.16g" group="" part="0" file="%s"/>\n'%(t,filename.encode()))
f.write(b'</Collection>\n')
f.write(b'</VTKFile>\n')
f.close()
with open(idxname,"wb" if reset_index else "rb+") as f :
if reset_index :
f.write(b'<?xml version="1.0"?>\n')
f.write(b'<VTKFile type="Collection" version="0.1" byte_order="LittleEndian">\n')
f.write(b'<Collection>\n')
else:
f.seek(-25,os.SEEK_END)
f.write(b'<DataSet timestep="%.16g" group="" part="0" file="%s"/>\n'%(t,filename.encode()))
f.write(b'</Collection>\n')
f.write(b'</VTKFile>\n')
f.close()
def write(basename,i,t,elements,x,data=None,field_data=None,cell_data=None,reset_index=None,vtktype="vtu") :
"""Write VTK UnstructuredMesh (vtu) or PolyData (vtp) files.
......@@ -105,43 +104,43 @@ def write(basename,i,t,elements,x,data=None,field_data=None,cell_data=None,reset
filename = "%s_%05d.%s"%(basename,i,vtktype)
if not os.path.exists(os.path.dirname(filename)):
os.makedirs(os.path.dirname(filename))
f = open(filename,"wb")
f.write(b'<?xml version="1.0"?>\n')
vtkdescr = b"UnstructuredGrid" if vtktype=="vtu" else b"PolyData"
f.write(b'<VTKFile type="%s" format="0.1" compressor="vtkZLibDataCompressor">\n' % vtkdescr)
f.write(b'<%s>\n'%vtkdescr)
nel = elements[0].shape[0] if elements is not None else 0
npt = x.shape[0]
f.write(b'<Piece NumberOfPoints="%i" NumberOfCells="%i">\n'%(npt,nel))
f.write(b'<Points>\n')
_write_array(f,x,b'NumberOfComponents="3"')
f.write(b'</Points>\n')
if elements is not None:
(types,connectivity,offsets) = elements
f.write(b'<Cells>\n')
_write_array(f,connectivity,b'Name="connectivity"')
_write_array(f,offsets,b'Name="offsets"')
_write_array(f,types,b'Name="types"')
f.write(b'</Cells>\n')
if data is not None :
f.write(b'<PointData>\n')
for name,v in data :
_write_array(f,v,b'Name="%s" NumberOfComponents="%i"' % (name.encode("utf8"),v.shape[-1]))
f.write(b'</PointData>\n')
if cell_data is not None :
f.write(b'<CellData>\n')
for name,v in cell_data :
_write_array(f,v,b'Name="%s" NumberOfComponents="%i"' % (name.encode("utf8"),v.shape[-1]))
f.write(b'</CellData>\n')
f.write(b'</Piece>\n')
if field_data is not None :
f.write(b'<FieldData>\n');
for name,v in field_data :
_write_array(f,v,b'Name="%s" NumberOfTuples="%i" NumberOfComponents="%i"' % (name,v.shape[0],v.shape[1]))
f.write(b'</FieldData>\n');
f.write(b'</%s>\n' % vtkdescr)
f.write(b'</VTKFile>\n')
f.close()
with open(filename,"wb") as f :
f.write(b'<?xml version="1.0"?>\n')
vtkdescr = b"UnstructuredGrid" if vtktype=="vtu" else b"PolyData"
f.write(b'<VTKFile type="%s" format="0.1" compressor="vtkZLibDataCompressor">\n' % vtkdescr)
f.write(b'<%s>\n'%vtkdescr)
nel = elements[0].shape[0] if elements is not None else 0
npt = x.shape[0]
f.write(b'<Piece NumberOfPoints="%i" NumberOfCells="%i">\n'%(npt,nel))
f.write(b'<Points>\n')
_write_array(f,x,b'NumberOfComponents="3"')
f.write(b'</Points>\n')
if elements is not None:
(types,connectivity,offsets) = elements
f.write(b'<Cells>\n')
_write_array(f,connectivity,b'Name="connectivity"')
_write_array(f,offsets,b'Name="offsets"')
_write_array(f,types,b'Name="types"')
f.write(b'</Cells>\n')
if data is not None :
f.write(b'<PointData>\n')
for name,v in data :
_write_array(f,v,b'Name="%s" NumberOfComponents="%i"' % (name.encode("utf8"),v.shape[-1]))
f.write(b'</PointData>\n')
if cell_data is not None :
f.write(b'<CellData>\n')
for name,v in cell_data :
_write_array(f,v,b'Name="%s" NumberOfComponents="%i"' % (name.encode("utf8"),v.shape[-1]))
f.write(b'</CellData>\n')
f.write(b'</Piece>\n')
if field_data is not None :
f.write(b'<FieldData>\n');
for name,v in field_data :
_write_array(f,v,b'Name="%s" NumberOfTuples="%i" NumberOfComponents="%i"' % (name,v.shape[0],v.shape[1]))
f.write(b'</FieldData>\n');
f.write(b'</%s>\n' % vtkdescr)
f.write(b'</VTKFile>\n')
f.close()
_write_index(basename+".pvd",os.path.basename(filename),i,t,reset_index)
def read(fname) :
......@@ -158,19 +157,20 @@ def read(fname) :
node_data,cell_data,field_data -- dictionaries whose keys are the name of the node,cell and field data.
"""
root = ET.parse(open(fname,"rb")).getroot()
ftype = root.attrib["type"]
grid = root.find(ftype)
piece = grid.find("Piece")
x = _read_array(piece.find("Points/DataArray"))
cells = piece.find("Cells")
fd = grid.find("FieldData")
cd = piece.find("CellData")
pd = piece.find("PointData")
c = dict((f.attrib["Name"],_read_array(f)) for f in cells) if cells else None
fdata = dict((f.attrib["Name"],_read_array(f)) for f in fd) if fd else None
cdata = dict((f.attrib["Name"],_read_array(f)) for f in cd) if cd else None
pdata = dict((f.attrib["Name"],_read_array(f)) for f in pd) if pd else None
with open(fname,"rb") as f :
root = ET.parse(f).getroot()
ftype = root.attrib["type"]
grid = root.find(ftype)
piece = grid.find("Piece")
x = _read_array(piece.find("Points/DataArray"))
cells = piece.find("Cells")
fd = grid.find("FieldData")
cd = piece.find("CellData")
pd = piece.find("PointData")
c = dict((f.attrib["Name"],_read_array(f)) for f in cells) if cells else None
fdata = dict((f.attrib["Name"],_read_array(f)) for f in fd) if fd else None
cdata = dict((f.attrib["Name"],_read_array(f)) for f in cd) if cd else None
pdata = dict((f.attrib["Name"],_read_array(f)) for f in pd) if pd else None
return x,c,pdata,cdata,fdata
def write_multipart(basename,parts,i,t,reset_index=None):
......@@ -188,15 +188,15 @@ def write_multipart(basename,parts,i,t,reset_index=None):
filename = basename+"_%05d.vtm"%i
if not os.path.exists(os.path.dirname(filename)):
os.makedirs(os.path.dirname(filename))
f = open(filename,"wb")
f.write(b'<?xml version="1.0"?>\n')
f.write(b'<VTKFile type="vtkMultiBlockDataSet" version="1.0" byte_order="LittleEndian">\n');
f.write(b'<vtkMultiBlockDataSet>\n');
for ip,n in enumerate(parts):
f.write(b'<DataSet index="%i" file="%s"/>\n'%(ip,n.encode()));
f.write(b'</vtkMultiBlockDataSet>\n');
f.write(b'</VTKFile>\n');
_write_index(basename+".pvd",os.path.basename(filename),i,t,reset_index)
with open(filename,"wb") as f :
f.write(b'<?xml version="1.0"?>\n')
f.write(b'<VTKFile type="vtkMultiBlockDataSet" version="1.0" byte_order="LittleEndian">\n');
f.write(b'<vtkMultiBlockDataSet>\n');
for ip,n in enumerate(parts):
f.write(b'<DataSet index="%i" file="%s"/>\n'%(ip,n.encode()));
f.write(b'</vtkMultiBlockDataSet>\n');
f.write(b'</VTKFile>\n');
_write_index(basename+".pvd",os.path.basename(filename),i,t,reset_index)
def string_array_encode(a) :
......
......@@ -57,112 +57,112 @@ class Mesh() :
self.read(filename)
def write(self, filename, format3=None) :
output = open(filename, "w")
if format3 == None :
format3 = self.useFormat3
if format3 :
output.write("$MeshFormat\n3 0 8\n$EndMeshFormat\n")
output.write("$PhysicalNames\n%i\n" % (len(self.physicals[0]) + len(self.physicals[1]) + len(self.physicals[2]) + len(self.physicals[3])))
for dim, plist in enumerate(self.physicals) :
for name, tag in plist.items() :
output.write("%i %i \"%s\"\n" % (dim, tag, name))
output.write("$EndPhysicalNames\n")
output.write("$Entities\n%i\n" % (len(self.entities)))
for e in self.entities :
output.write("%i %i %i %s\n" % (e.tag, e.dimension, len(e.physicals), " ".join(repr(ip) for ip in e.physicals)))
output.write("$EndEntities\n")
output.write("$Nodes\n%d\n" % len(self.vertices))
for (x, y, z, j) in self.vertices :
output.write("%i %.16g %.16g %.16g 0\n" % (j, x, y, z))
output.write("$EndNodes\n")
output.write("$Elements\n%d\n" % sum ([len(e.elements) for e in self.entities]))
for entity in self.entities :
for e in entity.elements :
ntag = len(e.vertices) + ((1 + len(e.partition)) if e.partition else 0)
spart = " ".join(repr(i) for i in [len(e.partition)] + e.partition) if e.partition else ""
output.write("%i %i %i %i %s %s\n" % (e.tag, e.etype.tag, entity.tag, ntag, " ".join(repr(v[3]) for v in e.vertices), spart))
output.write("$EndElements\n")
else :
output.write("$MeshFormat\n2.2 0 8\n$EndMeshFormat\n")
output.write("$PhysicalNames\n%i\n" % (len(self.physicals[0]) + len(self.physicals[1]) + len(self.physicals[2]) + len(self.physicals[3])))
for dim, plist in enumerate(self.physicals) :
for name, tag in plist.items() :
output.write("%i %i \"%s\"\n" % (dim, tag, name))
output.write("$EndPhysicalNames\n")
output.write("$Nodes\n%d\n" % len(self.vertices))
for (x, y, z, j) in self.vertices :
output.write("%i %.16g %.16g %.16g\n" % (j, x, y, z))
output.write("$EndNodes\n")
output.write("$Elements\n%d\n" % sum ([len(e.elements) for e in self.entities]))
for entity in self.entities :
for e in entity.elements :
ntag = 2 + ((1 + len(e.partition)) if e.partition else 0)
spart = " ".join(repr(i) for i in [len(e.partition)] + e.partition) if e.partition else ""
ptag = entity.physicals[0] if entity.physicals else 0
output.write("%i %i %i %i %i %s %s\n" % (e.tag, e.etype.tag, ntag, ptag, entity.tag, spart, " ".join(repr(v[3]) for v in e.vertices)))
output.write("$EndElements\n")
with open(filename, "w") as output :
if format3 == None :
format3 = self.useFormat3
if format3 :
output.write("$MeshFormat\n3 0 8\n$EndMeshFormat\n")
output.write("$PhysicalNames\n%i\n" % (len(self.physicals[0]) + len(self.physicals[1]) + len(self.physicals[2]) + len(self.physicals[3])))
for dim, plist in enumerate(self.physicals) :
for name, tag in plist.items() :
output.write("%i %i \"%s\"\n" % (dim, tag, name))
output.write("$EndPhysicalNames\n")
output.write("$Entities\n%i\n" % (len(self.entities)))
for e in self.entities :
output.write("%i %i %i %s\n" % (e.tag, e.dimension, len(e.physicals), " ".join(repr(ip) for ip in e.physicals)))
output.write("$EndEntities\n")
output.write("$Nodes\n%d\n" % len(self.vertices))
for (x, y, z, j) in self.vertices :
output.write("%i %.16g %.16g %.16g 0\n" % (j, x, y, z))
output.write("$EndNodes\n")
output.write("$Elements\n%d\n" % sum ([len(e.elements) for e in self.entities]))
for entity in self.entities :
for e in entity.elements :
ntag = len(e.vertices) + ((1 + len(e.partition)) if e.partition else 0)
spart = " ".join(repr(i) for i in [len(e.partition)] + e.partition) if e.partition else ""
output.write("%i %i %i %i %s %s\n" % (e.tag, e.etype.tag, entity.tag, ntag, " ".join(repr(v[3]) for v in e.vertices), spart))
output.write("$EndElements\n")
else :
output.write("$MeshFormat\n2.2 0 8\n$EndMeshFormat\n")
output.write("$PhysicalNames\n%i\n" % (len(self.physicals[0]) + len(self.physicals[1]) + len(self.physicals[2]) + len(self.physicals[3])))
for dim, plist in enumerate(self.physicals) :
for name, tag in plist.items() :
output.write("%i %i \"%s\"\n" % (dim, tag, name))
output.write("$EndPhysicalNames\n")
output.write("$Nodes\n%d\n" % len(self.vertices))
for (x, y, z, j) in self.vertices :
output.write("%i %.16g %.16g %.16g\n" % (j, x, y, z))
output.write("$EndNodes\n")
output.write("$Elements\n%d\n" % sum ([len(e.elements) for e in self.entities]))
for entity in self.entities :
for e in entity.elements :
ntag = 2 + ((1 + len(e.partition)) if e.partition else 0)
spart = " ".join(repr(i) for i in [len(e.partition)] + e.partition) if e.partition else ""
ptag = entity.physicals[0] if entity.physicals else 0
output.write("%i %i %i %i %i %s %s\n" % (e.tag, e.etype.tag, ntag, ptag, entity.tag, spart, " ".join(repr(v[3]) for v in e.vertices)))
output.write("$EndElements\n")
def read(self, filename):
fin = open(filename, "r");
l = fin.readline()
vmap = {}
entitymap = {}
while l != "" :
w = l.split()
if w[0] == "$MeshFormat":
l = fin.readline().split()
if float(l[0]) == 3.:
self.useFormat3 = True
elif int(float(l[0])) == 2 :
self.useFormat3 = False
else :
print("error : cannot read mesh format " + l[0])
with open(filename, "r") as fin :
l = fin.readline()
elif w[0] == "$PhysicalNames" :
n = int(fin.readline())
for i in range(n) :
dim, tag, name = fin.readline().split()
self.physicals[int(dim)][name[1:-1]] = int(tag)
fin.readline()
elif w[0] == "$Entities" and self.useFormat3:
n = int(fin.readline())
for i in range(n) :
l = fin.readline().split()
j, dim, nphys = int(l[0]), int(l[1]), int(l[2])
self.entities.append(Mesh.entity(j, dim, [int(ip) for ip in l[3:3+nphys]]))
entitymap[(dim, j)] = self.entities[-1]
fin.readline()
elif w[0] == "$Nodes" :
n = int(fin.readline())
for i in range(n) :
if self.useFormat3 :
(j, x, y, z, t) = fin.readline().split()
else :
(j, x, y, z) = fin.readline().split()
self.vertices.append([float(x), float(y), float(z), int(j)])
vmap[int(j)] = self.vertices[-1]
elif w[0] == "$Elements" :
n = int(fin.readline())
for i in range(n) :
l = fin.readline().split()
if self.useFormat3 :
j, t, e, nf = int(l[0]), int(l[1]), int(l[2]), int(l[3])
nv = gmshType.Type[t].numVertices
vertices = [vmap[int(i)] for i in l[4:4+nv]]
partition = [int(i) for i in l[5 + nv : 5 + nv + int(l[4 + nv])]] if nf > nv else []
else :
j, t, nf, p, e = int(l[0]), int(l[1]), int(l[2]), int(l[3]), int(l[4])
vertices = [vmap[int(i)] for i in l[3 + nf:]]
partition = [int(i) for i in l[6 : 6 + int(l[5])]] if nf > 2 else []
edim = gmshType.Type[t].baseType.dimension
entity = entitymap.get((edim, e), None)
if not self.useFormat3 and not entity:
entity = Mesh.entity(e, edim, [p])
self.entities.append(entity)
entitymap[(edim, e)] = entity
entity.elements.append(Mesh.element(gmshType.Type[t], vertices, partition, j))
self.maxeid = max(self.maxeid, j)
l = fin.readline()
vmap = {}
entitymap = {}
while l != "" :
w = l.split()
if w[0] == "$MeshFormat":
l = fin.readline().split()
if float(l[0]) == 3.:
self.useFormat3 = True
elif int(float(l[0])) == 2 :
self.useFormat3 = False
else :
print("error : cannot read mesh format " + l[0])
l = fin.readline()
elif w[0] == "$PhysicalNames" :
n = int(fin.readline())
for i in range(n) :
dim, tag, name = fin.readline().split()
self.physicals[int(dim)][name[1:-1]] = int(tag)
fin.readline()
elif w[0] == "$Entities" and self.useFormat3:
n = int(fin.readline())
for i in range(n) :
l = fin.readline().split()
j, dim, nphys = int(l[0]), int(l[1]), int(l[2])
self.entities.append(Mesh.entity(j, dim, [int(ip) for ip in l[3:3+nphys]]))
entitymap[(dim, j)] = self.entities[-1]
fin.readline()
elif w[0] == "$Nodes" :
n = int(fin.readline())
for i in range(n) :
if self.useFormat3 :
(j, x, y, z, t) = fin.readline().split()
else :
(j, x, y, z) = fin.readline().split()
self.vertices.append([float(x), float(y), float(z), int(j)])
vmap[int(j)] = self.vertices[-1]
elif w[0] == "$Elements" :
n = int(fin.readline())
for i in range(n) :
l = fin.readline().split()
if self.useFormat3 :
j, t, e, nf = int(l[0]), int(l[1]), int(l[2]), int(l[3])
nv = gmshType.Type[t].numVertices
vertices = [vmap[int(i)] for i in l[4:4+nv]]
partition = [int(i) for i in l[5 + nv : 5 + nv + int(l[4 + nv])]] if nf > nv else []
else :
j, t, nf, p, e = int(l[0]), int(l[1]), int(l[2]), int(l[3]), int(l[4])
vertices = [vmap[int(i)] for i in l[3 + nf:]]
partition = [int(i) for i in l[6 : 6 + int(l[5])]] if nf > 2 else []
edim = gmshType.Type[t].baseType.dimension
entity = entitymap.get((edim, e), None)
if not self.useFormat3 and not entity:
entity = Mesh.entity(e, edim, [p])
self.entities.append(entity)
entitymap[(edim, e)] = entity
entity.elements.append(Mesh.element(gmshType.Type[t], vertices, partition, j))
self.maxeid = max(self.maxeid, j)
l = fin.readline()
def getPhysicalNumber(self, dim, name) :
t = self.physicals[dim].get(name, None)
......
......@@ -85,7 +85,7 @@ class Poiseuille(unittest.TestCase) :
t = 0
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
fluid.export_vtk(outputdir,0,0)
......@@ -93,7 +93,7 @@ class Poiseuille(unittest.TestCase) :
tic = time.time()
while ii < 100 :
#Fluid solver
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
fluid.implicit_euler(dt)
forces = fluid.compute_node_force(dt)
vn = p.velocity() + forces * dt / p.mass()
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment