Commit 90385209 authored by Michel Henry's avatar Michel Henry
Browse files

uniformisation lecture/ecriture VTK

parent dd6df073
Pipeline #8778 passed with stages
in 5 minutes and 20 seconds
......@@ -326,13 +326,12 @@ class FluidProblem :
bnds[bname.value] = nodes
return bnds
def export_vtk(self, output_dir, t, it,stab=False) :
def write_vtk(self, output_dir, it ,t, stab=False):
"""Writes output file for post-visualization.
Keyword arguments:
output_dir -- Output directory
t -- Computation time
it -- Computation iteration
t -- Computation time
stab -- If True exports the stabilisation parametres in the output files
"""
v = self.solution()[:,:self._dim]
......@@ -359,12 +358,17 @@ class FluidProblem :
offsets = np.cumsum(np.repeat([self._dim+1],connectivity.shape[0])).astype(np.int32)
VTK.write(output_dir+"/fluid",it,t,(types,connectivity,offsets),self.coordinates(),data,field_data,cell_data)
def import_vtk(self, filename) :
def export_vtk(self, output_dir, t, it, stab=False) :
print("Careful this function is deprecated... \n\tUse write_vtk instead !")
self.write_vtk(output_dir, it, t, stab)
def read_vtk(self, dirname, i):
"""Reads output file to reload computation.
Keyword arguments:
filename -- Name of the file to read
"""
filename = dirname + "/fluid_%05d.vtu"%i
x,cells,data,cdata,fdata = VTK.read(filename)
mesh_boundaries = {name[9:]:nodes for name,nodes in fdata.items() if name.startswith("Boundary ")}
cbnd_names = (c_char_p*len(mesh_boundaries))(*(i.encode() for i in mesh_boundaries.keys()))
......@@ -390,6 +394,34 @@ class FluidProblem :
self.old_porosity()[:] = data["old_porosity"]
self.sys = None
def import_vtk(self, filename) :
print("Careful this function is deprecated... \n\tUse read_vtk instead !")
x,cells,data,cdata,fdata = VTK.read(filename)
mesh_boundaries = {name[9:]:nodes for name,nodes in fdata.items() if name.startswith("Boundary ")}
cbnd_names = (c_char_p*len(mesh_boundaries))(*(i.encode() for i in mesh_boundaries.keys()))
el = cells["connectivity"].reshape([-1,self._dim+1])
nbnd = len(mesh_boundaries)
bnds = np.vstack(list(mesh_boundaries.values()))
bnd_tags = np.repeat(list(range(nbnd)),list([v.shape[0] for v in mesh_boundaries.values()]))
bnd_tags = np.require(bnd_tags,np.int32,"C")
self._lib.mesh_new_from_elements.restype = c_void_p
_mesh = c_void_p(self._lib.mesh_new_from_elements(
c_int(x.shape[0]),_np2c(x,np.float64),
c_int(el.shape[0]),_np2c(el,np.int32),
c_int(bnds.shape[0]),_np2c(bnds,np.int32),
_np2c(bnd_tags,np.int32),c_int(len(cbnd_names)),cbnd_names,
_np2c(data["parent_node_id"],np.int32) if "parent_node_id" in data else None))
self._lib.fluid_problem_set_mesh(self._fp, _mesh)
sol = self.solution()
sol[:,:self._dim] = data["velocity"][:,:self._dim]
sol[:,[self._dim]] = data["pressure"]
if self._n_fluids == 2 :
self.concentration_dg()[:] = cdata["concentration"]
self.porosity()[:] = data["porosity"]
self.old_porosity()[:] = data["old_porosity"]
self.sys = None
def compute_node_force(self, dt) :
"""Computes the forces to apply on each grain resulting from the fluid motion.
......
......@@ -454,10 +454,9 @@ class ParticleProblem :
periodicX, periodicCells,_, periodicCellsData, periodicFieldsData = VTK.read(dirname+"/periodicBoundaries_%05d.vtu"%i)
entityTransformation = periodicFieldsData["Entity_transformation"]
entity_ids = np.hstack(periodicCellsData["Entity_id"])
entity_ids = periodicCellsData["Entity_id"].ravel()
periodicOffsets = np.hstack([[0], periodicCells["offsets"]])
periodicEl = periodicCells["connectivity"]
# add entity
for entity_id in range(np.shape(entityTransformation)[0]):
trans = tuple(entityTransformation[entity_id])
......
......@@ -85,7 +85,7 @@ class Poiseuille(unittest.TestCase) :
#s[:,3] = -rho*g*(1-x[:,1])
fluid.export_vtk(outputdir,0,0)
fluid.write_vtk(outputdir,0,0)
ii = 0
tic = time.time()
......@@ -96,7 +96,7 @@ class Poiseuille(unittest.TestCase) :
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
......
......@@ -71,7 +71,7 @@ class Darcy(unittest.TestCase) :
x[i] = x[i]+x_r
if y[i]+y_r<L-r and y[i]+y_r>r-L:
y[i] = y[i]+y_r
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop);
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop)
p.write_vtk(filename,0,0)
n_test = 10
......@@ -115,7 +115,7 @@ class Darcy(unittest.TestCase) :
fluid.set_wall_boundary("Top",velocity=[0,0])
fp = np.zeros_like(fluid.porosity())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
fluid.export_vtk(outputdir,0,0)
fluid.write_vtk(outputdir,0,0)
t = 0
ii = 0
......@@ -129,7 +129,7 @@ class Darcy(unittest.TestCase) :
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
......
......@@ -105,7 +105,7 @@ class DieSwell(unittest.TestCase):
#
fs = time_integration.FreeSurface(fluid,"TopRight","BottomRight",isCentered=False)
fluid.export_vtk(outputdir,0,0)
fluid.write_vtk(outputdir,0,0)
tic = time.time()
#
# COMPUTATION LOOP
......@@ -120,7 +120,7 @@ class DieSwell(unittest.TestCase):
# Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
swelling_experiment = 1.02 # Value from "Goren, S. L., & Wronski, S. (1966).
......
......@@ -87,7 +87,7 @@ class Poiseuille(unittest.TestCase) :
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
fluid.export_vtk(outputdir,0,0)
fluid.write_vtk(outputdir,0,0)
G = p.mass()*g
ii = 0
......@@ -98,7 +98,7 @@ class Poiseuille(unittest.TestCase) :
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
p.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
......
......@@ -69,8 +69,8 @@ class Periodic2d(unittest.TestCase):
p.add_particle((r,0.5),r,r**2*np.pi*rhop)
p.add_particle((L - 2*r,0.5),r,r**2*np.pi*rhop)
p.velocity()[1,0] = 0.5
p.write_vtk(outputdir,0,0)
p.read_vtk(outputdir,0)
#Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
fluid = mbfluid.FluidProblem(2,g,mu,rho,petsc_solver_type="-pc_type lu")
......@@ -87,7 +87,7 @@ class Periodic2d(unittest.TestCase):
#set initial_condition
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces())
fluid.export_vtk(outputdir,0,0)
fluid.write_vtk(outputdir,0,0)
#Initial time and iteration
t = 0
......@@ -99,10 +99,13 @@ class Periodic2d(unittest.TestCase):
#
# COMPUTATION LOOP
#
haveCrossed = 0
while t < tEnd :
# Time integration
time_integration.iterate(fluid,p,dt,external_particles_forces=G)
if np.max(p.position()[:,0]) > L:
haveCrossed = 1
p.position()[:,0] = np.fmod(p.position()[:,0],L)
ind = np.where(p.position()[:,0] <= 0)
p.position()[ind,0] += L
......@@ -111,10 +114,11 @@ class Periodic2d(unittest.TestCase):
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
p.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
error = np.min(np.linalg.norm(p.position()[1,:] - p.position()[0,:]))
self.assertGreater(error, r, "error is too large in periodic 2d")
self.assertEqual(haveCrossed, 1, "No periodic boundary in periodic 2d")
......@@ -88,7 +88,7 @@ class Periodic2d(unittest.TestCase):
fluid.set_strong_boundary("X",3,0)
#Set location of the particles in the mesh and compute the porosity in each computation cell
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(), p.contact_forces())
fluid.export_vtk(outputdir,0,0)
fluid.write_vtk(outputdir,0,0)
# Initial time and iteration
......@@ -101,9 +101,12 @@ class Periodic2d(unittest.TestCase):
#
# COMPUTATION LOOP
#
haveCrossed = 0
while t < tEnd :
time_integration.iterate(None, p, dt,min_nsub=5, contact_tol = 1e-8,external_particles_forces=G)
time_integration.iterate(fluid, p, dt,min_nsub=5, contact_tol = 1e-8,external_particles_forces=G)
if np.max(p.position()[:,0]) > L:
haveCrossed = 1
p.position()[:,0] = np.fmod(p.position()[:,0],L)
ind = np.where(p.position()[:,0] <= 0)
p.position()[ind,0] += L
......@@ -113,9 +116,10 @@ class Periodic2d(unittest.TestCase):
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
p.write_vtk(outputdir, ioutput, t)
fluid.export_vtk(outputdir, t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
error = np.min(np.linalg.norm(p.position()[1,:] - p.position()[0,:]))
self.assertGreater(error, r, "error is too large in periodic 3d")
self.assertEqual(haveCrossed, 1, "No periodic boundary in periodic 3d")
\ No newline at end of file
......@@ -107,7 +107,7 @@ class Poiseuille(unittest.TestCase) :
mu = nu0*rho
vI = 1/(120*mu)
fluid.export_vtk(outputdir,0,0)
fluid.write_vtk(outputdir,0,0)
ii = 0
tic = time.time()
......@@ -118,7 +118,7 @@ class Poiseuille(unittest.TestCase) :
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
......
......@@ -79,7 +79,7 @@ class Poiseuille(unittest.TestCase) :
#set initial_condition
fluid.export_vtk(outputdir,0,0)
fluid.write_vtk(outputdir,0,0)
ii = 0
tic = time.time()
......@@ -90,7 +90,7 @@ class Poiseuille(unittest.TestCase) :
#Output files writting
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
fluid.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
......
......@@ -122,7 +122,7 @@ class Weight(unittest.TestCase) :
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
#p.write_vtk(outputdir, ioutput, t)
p.write_vtk(outputdir, ioutput, t)
print("%i : %.2g/%.2g" % (ii, t, tEnd))
error = np.sqrt(sum(sum(forces) - bnd_forces)**2)
error /= np.sqrt(sum(sum(forces)**2))
......
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