particle state

fluid/fluid_problem.c

@@ -1728,7 +1728,7 @@ void fluid_problem_after_import(FluidProblem *problem)
   compute_porosity(problem->mesh, problem->node_volume, problem->porosity, problem->n_particles, problem->particle_position, problem->particle_volume, problem->particle_element_id, problem->particle_uvw);
 }
 
-void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *contact, long int *elid, int init, int reload)
+void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *state, double *contact, long int *elid, int init, int reload)
 {
   
   struct timespec tic, toc;
@@ -1761,16 +1761,17 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
       problem->particle_element_id[i]=elid[i];
     }
   }
-  mesh_tree_particle_to_mesh(problem->mesh_tree, n, position, problem->particle_element_id, problem->particle_uvw);
+  int rowl = (D==2 ? 5 : 9);
   for (int i = 0; i < n; ++i) {
     problem->particle_mass[i] = mass[i];
     problem->particle_volume[i] = volume[i];
     for (int k = 0; k < D; ++k) {
-      problem->particle_position[i*D+k] = position[i*D+k];
-      problem->particle_velocity[i*D+k] = velocity[i*D+k];
+      problem->particle_position[i*D+k] = state[i*rowl+k];
+      problem->particle_velocity[i*D+k] = state[i*rowl+k];
       problem->particle_contact[i*D+k] = contact[i*D+k];
     }
   }
+  mesh_tree_particle_to_mesh(problem->mesh_tree, n, problem->particle_position, problem->particle_element_id, problem->particle_uvw);
   if (!reload){
     for (int i=0; i<problem->mesh->n_nodes; ++i){
       problem->oldporosity[i] = problem->porosity[i];

fluid/fluid_problem.h

@@ -104,7 +104,7 @@ struct FluidProblem {
 
 // complete force on the particle (including gravity)
 void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force);
-void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *contact, long int *elid, int init, int reload);
+void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *state, double *contact, long int *elid, int init, int reload);
 void fluid_problem_set_bulk_force(FluidProblem *problem, double *force);
 void fluid_problem_set_concentration_cg(FluidProblem *problem, double *concentration);
 void fluid_problem_free(FluidProblem *problem);

python/scontact.py

@@ -58,6 +58,28 @@ def _np2c(a,dtype=np.float64) :
     r.tmp = tmp
     return r
 
+class ParticleState(np.ndarray) :
+
+    @property
+    def x(self) :
+        dim = 2 if self.shape[1] == 5 else 3
+        return self[:,:dim]
+
+    @property
+    def v(self) :
+        dim = 2 if self.shape[1] == 5 else 3
+        return self[:,dim:2*dim]
+
+    @property
+    def omega(self) :
+        dim = 2 if self.shape[1] == 5 else 3
+        return self[:,dim:2*dim]
+
+    @classmethod
+    def new(cls, nparticles,dim) :
+        state = cls((nparticles, 5 if dim == 2 else 9))
+        return state
+
 
 class ParticleProblem :
     """Create the numerical structure containing all the physical particles that appear in the problem"""
@@ -152,18 +174,21 @@ class ParticleProblem :
         """Return the list of particle masses."""
         return self._get_matrix("Particle", 2)[:, 1][:,None]
 
-    def position(self):
-        """Return the list of particle centre position vectors."""
-        return self._get_matrix("Position", self._dim)
-
-    def velocity(self):
-        """Return the list of particle velocity vectors."""
-        return self._get_matrix("Velocity", self._dim)
+    def state(self) :
+        state = ParticleState.new(self.n_particles,self._dim)
+        state.x[:] = self._get_matrix("Position", self._dim)
+        state.vel[:] = self._get_matrix("Velocity", self._dim)
+        if self._friction_enabled and self._rotation_enabled :
+            state.omega[:] = self._get_matrix("Omega", 1 if self._dim==2 else 3))
+        else :
+            omega = np.zeros(vel.sahpe[0],(1 if self._dim==2 else 3))
+        return state
 
-    def omega(self):
-        """Return the list of particle angular velocity."""
-        assert self._friction_enabled and self._rotation_enabled
-        return (self._get_matrix("Omega", 1)) if self._dim==2 else self._get_matrix("Omega", 3)
+    def set_state(self, state) :
+        self._get_matrix("Position", self._dim)[:] = state.x
+        self._get_matrix("Velocity", self._dim)[:] = state.v
+        if self._friction_enabled and self._rotation_enabled :
+            elf._get_matrix("Omega", 1 if self._dim==2 else 3))[:] = state.omega
 
     def particle_material(self):
         return self._get_idx("ParticleMaterial")
@@ -249,7 +274,6 @@ class ParticleProblem :
         disks["vt"][disks["tag"]==tag] = vt
         segments["vt"][segments["tag"]==tag] = vt
 
-    
     def set_friction_coefficient(self, mu, mat0="default",mat1="default") :
         """ Sets the friction coefficient between two materials.
 

python/time_integration.py

@@ -34,9 +34,10 @@ def _advance_particles(particles, f, dt, min_nsub,contact_tol,iteration=0,after_
     max_split -- maximum number of times the time step can be further split if conergence is not reached
     """
     # Compute free velocities of the grains
+    state = particles.state()
     if f is None:
-        f = np.zeros_like(particles.velocity())
-    vn = particles.velocity() + f*dt / particles.mass()
+        f = np.zeros_like(state.v)
+    vn = state.v + f*dt / particles.mass()
     # Estimation of the solid sub-time steps
     if particles.r() is not None:
       vmax = np.max(np.linalg.norm(vn,axis=1))
@@ -76,13 +77,12 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
     after_sub_iter -- callback to execute once a sub iteration has been made
     max_split -- maximum number of times the time step can be further split if conergence is not reached
     """
+    state0 = particles.state()
     if external_particles_forces is None:
-        external_particles_forces = np.zeros_like(particles.velocity())
+        external_particles_forces = np.zeros_like(state0.v)
     # Copy the fluid solution of the previous time step before computing the prediction 
     sf0 = np.copy(fluid.solution())
     # Copy the solid solution of the previous time step before computing the prediction
-    vp0 = np.copy(particles.velocity())
-    pp0 = np.copy(particles.position())
     cp0 = np.copy(particles.contact_forces())
 
     # Predictor
@@ -92,10 +92,9 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
     sf1 = np.copy(fluid.solution())
     f0 = np.copy(fluid.compute_node_force(dt))+external_particles_forces
     _advance_particles(particles,f0,dt,min_nsub,contact_tol,max_split=max_split)
+    state1 = particles.state()
     # Keep same contact forces and positions while giving to the fluid the predicted solid velocities to compute the correction
-    particles.position()[:,:] = pp0
-    particles.contact_forces()[:,:] = cp0
-    fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(), particles.contact_forces(), reload=True)
+    fluid.set_particles(particles.mass(), particles.volume(), state0.x, state1.v, cp0, reload=True)
 
     # Corrector
     #
@@ -111,10 +110,11 @@ def predictor_corrector_iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e
     if fluid.n_fluids() == 2 :
         fluid.advance_concentration(dt)
     # Reset solid velocities
-    particles.velocity()[:,:] = vp0 
+    particles.set_state(state0)
     _advance_particles(particles,(alpha*f0+(1-alpha)*f1),dt,min_nsub,contact_tol,after_sub_iter=after_sub_iter,max_split=max_split)
     # Give to the fluid the solid information
-    fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(), particles.contact_forces()) 
+    state2 = particles.state()
+    fluid.set_particles(particles.mass(), particles.volume(), state2.x, state2.v, particles.contact_forces()) 
 
 def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, fixed_grains=False, after_sub_iter=None,max_split=1) :  
     """Migflow solver: the solver type depends on the fluid and particles given as arguments.
@@ -137,7 +137,7 @@ def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particl
     if particles is None and fluid is None:
         raise ValueError("fluid and particles are both None: not possible to iterate!")
     if particles is not None and external_particles_forces is not None:
-        if external_particles_forces.shape!=particles.velocity().shape:
+        if external_particles_forces.shape!=(particles.n_particles,particles.dim())
             raise ValueError("external_particles_forces must have shape (number of particles,dimension!")
     if (particles is None or particles.r() is None) and fluid is not None:
         fluid.implicit_euler(dt)
@@ -150,14 +150,16 @@ def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particl
     if fluid is not None and particles is not None and particles.r() is not None:
         if fixed_grains:
             fluid.implicit_euler(dt)
-            fluid.set_particles(particles.mass(), particles.volume(), particles.position(), particles.velocity(),-fluid.compute_node_force(dt))
+            state = particles.state()
+            fluid.set_particles(particles.mass(), particles.volume(), state.x, state.v,-fluid.compute_node_force(dt))
         else:
             predictor_corrector_iterate(fluid,particles,dt,min_nsub,contact_tol,external_particles_forces,after_sub_iter=after_sub_iter,max_split=max_split)
 
 def particle_changed(fluid,p) :
     co = fluid.old_porosity().copy()
     c = fluid.porosity().copy()
-    fluid.set_particles(p.mass(),p.volume(), p.position(), p.velocity(), p.contact_forces())
+    state = p.state()
+    fluid.set_particles(p.mass(),p.volume(), p.x, p.v, p.contact_forces())
     c2 = fluid.porosity()
     fluid.old_porosity()[:] = co+c2-c
     fluid.velocity()[:] *= c2/c