Commit 7030a013 authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

Merge branch 'master' into margaux

parents d47983db a0bec0a5
......@@ -183,6 +183,7 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
double u[D], du[D][D], dp[D], dun[D], dunt[D];
double s_c = 0;
double unold = 0;
double unext = 0;
double dcn = 0;
for (int iD = 0; iD < D; ++iD) {
u[iD] = s[U+iD];
......@@ -190,6 +191,7 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
dun[iD] = 0;
dunt[iD] = 0;
unold += sold[U+iD]*n[iD];
unext += data[vid+iD]*n[iD];
dcn += dc[iD]*n[iD];
s_c += vid<0?0:(u[iD]-data[vid+iD])*n[iD];
for (int jD = 0; jD < D; ++jD) du[iD][jD] = ds[(U+iD)*D+jD];
......@@ -220,12 +222,16 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
f00[(U+id)*n_fields+U+jd] += (vid<0?0:n[id]*n[jd]*sigma);
}
}
if (wbnd->type == BND_OPEN && unold<0) {
for (int id = 0; id < D; ++id) {
f0[U+id] += unold*((vid<0 ? 0 : data[vid+id])-u[id])*rho/c;
f00[(U+id)*n_fields+U+id] -= unold*rho/c;
}
}
double unbnd = unold;
if (wbnd->type == BND_WALL) unbnd = unold/2;
else if(vid>0) unbnd = (unold+unext)/2;
if (unbnd<0) {
for (int id = 0; id < D; ++id) {
f0[U+id] += unbnd*((vid<0 ? 0 : data[vid+id])-u[id])*rho/c;
f00[(U+id)*n_fields+U+id] -= unbnd*rho/c;
}
}
}
static double pow2(double a) {return a*a;}
......
......@@ -167,7 +167,7 @@ class ParticleProblem :
else :
return np.array([],np.int32)
def previousContactForces(self):
def contact_forces(self):
if self._fc is None:
self._fc = np.zeros_like(self.velocity())
return(self._fc)
......@@ -245,7 +245,7 @@ class ParticleProblem :
if self._dim == 2 :
v = np.insert(v,self._dim,0,axis=1)
x = np.insert(x,self._dim,0,axis=1)
data = [("Mass",self.mass()), ("Radius",self.r()),("Velocity",v),("Omega",omega),("Material",material),("PreviousContactForces",self.previousContactForces())]
data = [("Mass",self.mass()), ("Radius",self.r()),("Velocity",v),("Omega",omega),("Material",material),("Contacts",self.contact_forces())]
nmat = self._lib.particleProblemNMaterial(self._p)
self._lib.particleProblemGetMaterialTagName.restype = c_char_p
tags = list([self._lib.particleProblemGetMaterialTagName(self._p,i) for i in range(nmat)])
......@@ -297,7 +297,7 @@ class ParticleProblem :
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))
self._fc = d["PreviousContactForces"][:,:self._dim]
self._fc = d["Contacts"][:,:self._dim]
x,cells,d,cdata,fdata = VTK.read(dirname+"/boundaries_%05d.vtu"%i)
tagnames = VTK.string_array_decode(fdata["TagNames"])
......
......@@ -124,7 +124,7 @@ fluid.set_wall_boundary("Top", pressure=0)
fluid.set_wall_boundary("Bottom")
fluid.set_wall_boundary("Lateral")
# 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.previousContactForces())
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.previousContactForces())
# Initial fluid velocities are computed based on the forces applied by the particles without moving the particles
# The initial mesh is obtained by adapting for times the mesh on the converged fields
......
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