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

reload contacts

parent 63e9c1a8
Pipeline #9338 failed with stages
in 2 minutes and 21 seconds
......@@ -421,7 +421,7 @@ class ParticleProblem :
VTK.write(odir+"/contacts",i,t,None,x,field_data=fdata,vtktype="vtp")
VTK.write_multipart(odir+"/particle_problem",["particles_%05d.vtp"%i,"boundaries_%05d.vtu"%i,"periodicBoundaries_%05d.vtu"%i,"contacts_%05d.vtp"%i],i,t)
def read_vtk(self,dirname,i,contact_factor=1) :
def read_vtk(self,dirname,i) :
"""Reads an existing output file to set the particle data.
Keyword arguments:
dirname -- Path to the directory of the file to read
......@@ -500,14 +500,15 @@ class ParticleProblem :
x2 = periodicX[pidx+2,:self._dim]
self.add_boundary_periodic_triangle(tuple(x0), tuple(x1), tuple(x2), entity_id)
return
_,_,_,_,fdata = VTK.read(dirname+"/contacts_%05d.vtp"%i)
ks = ("particle_particle","particle_disk","particle_segment","particle_triangle")
v = np.vstack(list([fdata[k] for k in ks]))*contact_factor
v = np.vstack(list([fdata[k] for k in ks]))
o = np.vstack([fdata[k+"_idx"] for k in ks])
t = np.repeat([0,1,2,3],[fdata[k].shape[0] for k in ks])
self._lib.particle_problem_set_contacts(self._p,c_int(t.shape[0]),_np2c(o,np.uint64),_np2c(v),_np2c(t,np.int32))
self._lib.particle_problem_allocate_contacts(self._p, c_size_t(t.size))
self.contacts()["reaction"] = v
self.contacts()["o"] = o
self.contacts()["type"] = t
def add_boundary_disk(self, x0, r, tag, material="default") :
"""Adds a boundary disk.
......
......@@ -1351,15 +1351,7 @@ int particle_problem_iterate(ParticleProblem *p, double alert, double dt, double
}
for (size_t i = 0; i < vectorSize(p->bodies); ++i){
Body *body = &p->bodies[i];
for (int d = 0; d < DIMENSION; ++d) {
body->position[d] += body->velocity[d]*dt;
}
#if DIMENSION ==2
body->theta += body->omega*dt;
#else
//TODO
#endif
body_advance(&p->bodies[i], dt);
}
free(oldVelocity);
......
......@@ -43,15 +43,15 @@ dt = 0.0025
outf=15
t = 0
while t<tend :
print(i)
if i%100 == 0:
add_piece(p,n_inserted_pieces)
n_inserted_pieces += 1
t += dt
f = g*(np.pi*p.r()**2*rho)
p.iterate(dt,f,tol=1e-8,force_motion=1)
p.iterate(dt,f,tol=1e-7,force_motion=1)
if i %outf == 0 :
ioutput = int(i/outf) + 1
print(ioutput)
p.write_vtk(odir, ioutput, t)
#p.read_vtk(odir, ioutput, t)
p.read_vtk(odir, ioutput)
i += 1
Supports Markdown
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