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

write/read VTK

parent 45ae898c
......@@ -66,7 +66,6 @@ class ParticleProblem :
if ptr is None :
return np.ndarray((0),dtype)
size = cast(ptr-(8 if is_64bits else 4), POINTER(c_size_t)).contents.value//dtype.itemsize
print(size)
buf = (size*np.ctypeslib.as_ctypes_type(dtype)).from_address(ptr)
return np.ctypeslib.as_array(buf)
......@@ -124,9 +123,9 @@ class ParticleProblem :
self._disktype = np.dtype(bndtype+[('x',np.float64,dim),('v',np.float64,dim),('r',np.float64)])
self._segmenttype = np.dtype(bndtype+[('p',np.float64,(2,dim)),('v',np.float64,(2,dim))])
self._triangletype = np.dtype(bndtype+[('p',np.float64,(3,dim)),('v',np.float64,(3,dim))])
self._periodicEntitytype = np.dtype([('etag', np.int32), ('edim', np.int32), ('periodic_transformation', np.float64, dim)])
self._periodicSegmenttype = np.dtype([('entity_id', np.int32),('p',np.float64,(2,dim))])
self._periodicTriangletype = np.dtype([('entity_id', np.int32),('p', np.float64,(3,dim))])
self._periodicEntitytype = np.dtype([('etag', np.int32),('edim', np.int32),('periodic_transformation', np.float64, dim)])
self._periodicSegmenttype = np.dtype([('entity_id', np.int64), ('p',np.float64,(2,dim))])
self._periodicTriangletype = np.dtype([('entity_id', np.int64),('p', np.float64,(3,dim))])
def __del__(self):
"""Deletes the particle solver structure."""
......@@ -221,7 +220,6 @@ class ParticleProblem :
def periodicSegments(self):
"""Returns the list of periodic segments."""
print("Periodic Segment")
return self._get_array("PeriodicSegment", self._periodicSegmenttype)
def periodicTriangles(self):
......@@ -365,6 +363,7 @@ class ParticleProblem :
self._lib.particleProblemGetMaterialTagName.restype = c_char_p
tags = list([self._lib.particleProblemGetMaterialTagName(self._p,i) for i in range(nmat)])
VTK.write(odir+"/particles",i,t,None,x,data,vtktype="vtp",field_data=[(b"MaterialNames",VTK.string_array_encode(tags))])
disks = self.disks()
segments = self.segments()
triangles = self.triangles()
......@@ -380,24 +379,21 @@ class ParticleProblem :
self._lib.particleProblemGetTagName.restype = c_char_p
tagnames = list([self._lib.particleProblemGetTagName(self._p,i) for i in range(ntagname)])
VTK.write(odir+"/boundaries",i,t,(types,connectivity,offsets),x,cell_data=[("Tag",tags.reshape([-1,1])),("Material",materials.reshape([-1,1]))],field_data=[(b"TagNames",VTK.string_array_encode(tagnames))])
periodicEntity = self.periodicEntity()
print("PEntity : \n", periodicEntity)
periodicSegments = self.periodicSegments()
print("periodicSegment itemsize :", self._periodicSegmenttype.itemsize)
print("PSegments : \n",periodicSegments)
periodicTriangles = self.periodicTriangles()
# print(periodicTriangles)
periodicX = np.vstack([periodicSegments["p"].reshape([-1,self._dim]), triangles["p"].reshape([-1,self._dim])])
periodicX = np.vstack([periodicSegments["p"].reshape([-1,self._dim]), periodicTriangles["p"].reshape([-1,self._dim])])
if self._dim == 2 :
periodicX = np.insert(x,self._dim, 0, axis=1)
periodicX = np.insert(periodicX,self._dim, 0, axis=1)
periodicConnectivity = np.arange(0,periodicX.shape[0], dtype = np.int32)
periodicTypes = np.repeat(np.array([3,5],dtype=np.int32),[periodicSegments.shape[0],periodicTriangles.shape[0]])
periodicOffsets = np.cumsum(np.repeat([2,3],[periodicSegments.shape[0],periodicTriangles.shape[0]]),dtype=np.int32)
periodicEntityIds = np.array(np.hstack([periodicSegments["entity_id"],periodicTriangles["entity_id"]]))
# entityTransformation = np.array(periodicEntity["periodic_transformation"])
periodicTypes = np.repeat(np.array([3,5],dtype=np.int32),[periodicSegments.shape[0], periodicTriangles.shape[0]])
periodicOffsets = np.cumsum(np.repeat([2,3],[periodicSegments.shape[0], periodicTriangles.shape[0]]),dtype=np.int32)
periodicEntityIds = np.array(np.hstack([periodicSegments["entity_id"], periodicTriangles["entity_id"]]))
entityTransformation = np.array(periodicEntity["periodic_transformation"])
VTK.write(odir+"/periodicBoundaries", i, t, (periodicTypes, periodicConnectivity, periodicOffsets),
periodicX, cell_data=[("Entity_id", periodicEntityIds.reshape([-1,1]))])
periodicX, cell_data=[("Entity_id", periodicEntityIds.reshape([-1,1]))], field_data=[(b"Entity_transformation", entityTransformation)])
#Contacts
self._lib.particleProblemContactN.restype = c_size_t
......@@ -433,12 +429,12 @@ class ParticleProblem :
self._lib.particleProblemClearParticles(self._p)
self._lib.particleProblemAddParticles(self._p,c_int(x.shape[0]),
_np2c(x[:,:self._dim]), _np2c(d["Mass"]), _np2c(d["Radius"]), _np2c(d["Velocity"][:,:self._dim]),_np2c(d["Omega"]),_np2c(partmat,dtype=np.int32),_np2c(forced,dtype=np.int32),_np2c(d["Contacts"][:,:self._dim]*contact_factor))
x,cells,d,cdata,fdata = VTK.read(dirname+"/boundaries_%05d.vtu"%i)
tagnames = VTK.string_array_decode(fdata["TagNames"])
tagmap = {}
for j,n in enumerate(tagnames) :
tagmap[j] = self._lib.particleProblemGetTagId(self._p,n)
offsets = np.hstack([[0],cells["offsets"]])
el = cells["connectivity"]
tags = np.vectorize(tagmap.get)(cdata["Tag"])
......@@ -455,6 +451,33 @@ class ParticleProblem :
x1 = x[el[offsets[idx]+1],:self._dim]
x2 = x[el[offsets[idx]+2],:self._dim]
self._lib.particleProblemAddBoundaryTriangleTagId(self._p, self._coord_type(*x0), self._coord_type(*x1), self._coord_type(*x2), c_int(tags[idx,0]), c_int(materials[idx,0]))
periodicX, periodicCells,_, periodicCellsData, periodicFieldsData = VTK.read(dirname+"/periodicBoundaries_%05d.vtu"%i)
entityTransformation = periodicFieldsData["Entity_transformation"]
entity_ids = np.hstack(periodicCellsData["Entity_id"])
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])
self.add_boundary_periodic_entity(entity_id, self._dim-1, trans)
# add periodic segment
for idx in np.where(periodicCells["types"] == 3)[0]:
pidx = periodicEl[periodicOffsets[idx]]
entity_id = entity_ids[pidx//2]
x0 = periodicX[pidx,:self._dim]
x1 = periodicX[pidx+1,:self._dim]
self.add_boundary_periodic_segment(tuple(x0), tuple(x1),entity_id)
# add periodic triangles
for idx in np.where(periodicCells["types"] == 5)[0]:
pidx = periodicEl[periodicOffsets[idx]]
entity_id = entity_ids[pidx//3]
x0 = periodicX[pidx,:self._dim]
x1 = periodicX[pidx+1,:self._dim]
x2 = periodicX[pidx+2,:self._dim]
self.add_boundary_periodic_triangle(tuple(x0), tuple(x1), tuple(x2), entity_id)
_,_,_,_,fdata = VTK.read(dirname+"/contacts_%05d.vtp"%i)
ks = ("particle_particle","particle_disk","particle_segment","particle_triangle")
v = []
......@@ -684,7 +707,6 @@ class ParticleProblem :
self.add_boundary_periodic_triangle(
x[pmap[l[0]]]+shift, x[pmap[l[1]]]+shift,
x[pmap[l[2]]]+shift, ptag)
# print(self.periodicSegments())
for dim, tag in gmsh.model.getEntities(0) :
stag = get_entity_name(dim,tag)
if (dim,tag) in periodic_entities : continue
......
......@@ -228,7 +228,7 @@ static int diskInitContact(size_t id, const Disk *d, size_t particle, const Part
return c->D < alert;
}
struct _PeriodicSegment{
int entity_id;
size_t entity_id;
double p[2][DIMENSION];
};
......@@ -330,7 +330,7 @@ struct _Triangle {
double v[3][DIMENSION];
};
struct _PeriodicTriangle {
int entity_id;
size_t entity_id;
double p[3][DIMENSION];
};
......@@ -357,7 +357,7 @@ size_t particleProblemAddBoundaryPeriodicTriangle(ParticleProblem *p, const dou
for(size_t ie = 0; ie < vectorSize(p->periodicEntities); ++ie){
PeriodicEntity *pe = &p->periodicEntities[ie];
if(etag == pe->etag){
t->entity_id = (int) ie;
t->entity_id = ie;
break;
}
}
......@@ -1387,7 +1387,7 @@ size_t particleProblemAddBoundaryPeriodicSegment(ParticleProblem *p, const doubl
for(size_t ie = 0; ie < vectorSize(p->periodicEntities); ++ie){
PeriodicEntity *pe = &p->periodicEntities[ie];
if (etag == pe->etag){
ps->entity_id = (int)ie;
ps->entity_id = ie;
break;
}
}
......@@ -1569,13 +1569,7 @@ double *particleProblemVelocity(ParticleProblem *p) {return &p->velocity[0];}
double *particleProblemDisk(ParticleProblem *p) { return (double*)&p->disks[0];}
double *particleProblemPeriodicEntity(ParticleProblem *p) {return (double*)&p->periodicEntities[0];}
double *particleProblemSegment(ParticleProblem *p) {return (double*)&p->segments[0];}
double *particleProblemPeriodicSegment(ParticleProblem *p){
for(size_t i = 0; i < vectorSize(p->periodicSegments); ++i){
const PeriodicSegment pe = p->periodicSegments[i];
printf("\ni : %lu\tetag : %d\tp1 = [%f %f],\tp2 = [%f %f]",i, pe.entity_id, pe.p[0][0], pe.p[0][1], pe.p[1][0], pe.p[1][1]);
}
printf("\n");
return (double*)&p->periodicSegments[0];}
double *particleProblemPeriodicSegment(ParticleProblem *p){return (double*)&p->periodicSegments[0];}
double *particleProblemParticle(ParticleProblem *p) {return (double*)&p->particles[0];}
double *particleProblemPosition(ParticleProblem *p){return &p->position[0];}
double *particleProblemContactForces(ParticleProblem *p){return &p->contactForces[0];}
......
......@@ -65,7 +65,7 @@ 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)
# p.read_vtk(outputdir,0)
#Initial time and iteration
t = 0
......@@ -76,7 +76,7 @@ tic = time.time()
#
# COMPUTATION LOOP
#
while t < 0 :
while t < tEnd :
# Time integration
time_integration.iterate(None,p,dt)
......
......@@ -22,6 +22,7 @@
#!/usr/bin/env python
from migflow import fluid as mbfluid
from migflow import time_integration
from migflow import scontact
import numpy as np
import os
import subprocess
......@@ -47,7 +48,7 @@ dir_path = os.path.dirname(os.path.realpath(__file__))
os.chdir(dir_path)
# Physical parameters for the drops are the ones presented by Metzger et al. (2007) "Falling clouds of particles in viscous fluids"
outputdir = "outputPoiseuilleVTK"
outputdir = "outputVTK"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
......@@ -62,43 +63,50 @@ rho = 1000 # fluid density
nu = 1e-3 # kinematic viscosity
mu = nu*rho # dynamic viscosity
tEnd = 1000 # final time
r = .25e-1 # particle radius
L = 1
#numerical parameters
dt = 10 # time step
outf = 10 # number of iterations between output files
tEnd = 2 # final time
dt = 0.02 # time step
outf = 5 # number of iterations between output files
# shutil.copy("mesh.msh", outputdir +"/mesh.msh")
#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,nu*rho,rho,petsc_solver_type="-pc_type lu")
# fluid.load_msh("mesh.msh") # a retirer
fluid.import_vtk("outputPoiseuille/fluid_%05d.vtu"%ii)
fluid.set_wall_boundary("Bottom",velocity=[0,0])
fluid.set_wall_boundary("Top",velocity=[0,0])
print("Import done !")
# fluid = mbfluid.FluidProblem(2,g,nu*rho,rho,petsc_solver_type="-pc_type lu")
# # fluid.load_msh("mesh.msh") # a retirer
# fluid.import_vtk("outputPoiseuille/fluid_%05d.vtu"%ii)
# fluid.set_wall_boundary("Bottom",velocity=[0,0])
# fluid.set_wall_boundary("Top",velocity=[0,0])
# print("Import done !")
p = scontact.ParticleProblem(2)
p.read_vtk("output",0)
p.write_vtk(outputdir,0,0)
ii = 0
t = ii*dt
#set initial_condition
fluid.export_vtk(outputdir,t,ii)
# fluid.export_vtk(outputdir,t,ii)
tic = time.time()
while t < tEnd:
#Fluid solver
time_integration.iterate(fluid,None,dt)
t += dt
while t < tEnd :
# Time integration
time_integration.iterate(None,p,dt)
p.position()[:,0] = np.fmod(p.position()[:,0],L)
ind = np.where(p.position()[:,0] <= 0)
p.position()[ind,0] += L
#Output files writting
t += dt
if ii %outf == 0 :
ioutput = int(ii/outf) + 1
fluid.export_vtk(outputdir, t, ioutput)
p.write_vtk(outputdir, ioutput, t)
ii += 1
print("%i : %.2g/%.2g (cpu %.6g)" % (ii, t, tEnd, time.time() - tic))
nodesRight = get_nodes(fluid, "Right")
nodesLeft = get_nodes(fluid, "Left")
s = fluid.solution()
error = np.abs(np.max(s[nodesRight,:],axis = 0) - np.max(s[nodesLeft,:], axis = 0))
eps = 1e-6
print(f"Periodic success ? {error} <? {eps} => {np.all(error < eps)}")
error = np.min(np.linalg.norm(p.position()[1,:] - p.position()[0,:]))
print(f"particle penetration ? {error < r}")
......@@ -36,7 +36,7 @@ os.chdir(dir_path)
outputdir = "output"
if not os.path.isdir(outputdir) :
os.makedirs(outputdir)
subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","1"])
subprocess.call(["gmsh", "-2", "mesh.geo","-clscale","3"])
# Physical parameters
rho = 1000 # fluid density
......@@ -58,11 +58,13 @@ outf = 5 # number of iterations between o
#
p = scontact.ParticleProblem(3)
# Initialise particles boundary
p.load_msh_boundaries("mesh.msh", ["Bottom", "Top", "X", "Z"])
p.add_particle((r,0.5,0.5), r, 4/3*r**3*np.pi*rhop)
p.add_particle((L-2*r,0.5,0.5), r, 4/3*r**3*np.pi*rhop)
p.velocity()[1,0] = 0.5
# p.load_msh_boundaries("mesh.msh", ["Bottom", "Top", "X", "Z"])
# p.add_particle((r,0.5,0.5), r, 4/3*r**3*np.pi*rhop)
# p.add_particle((L-2*r,0.5,0.5), r, 4/3*r**3*np.pi*rhop)
# p.velocity()[1,0] = 0.5
# p.write_vtk(outputdir,0,0)
p.read_vtk(outputdir,0)
p.write_vtk(outputdir,0,0)
#Initial time and iteration
......
......@@ -60,7 +60,7 @@ class Periodic2d(unittest.TestCase):
# Numerical parameters
outf = 1 # number of iterations between output files
dt = 0.1 # time step
tEnd = 2 # final time
tEnd = 2 # final time
#
# PARTICLE PROBLEM
......
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