new variable for set_particles in c file

fluid/fluid_problem.c

@@ -1602,7 +1602,7 @@ void fluid_problem_after_import(FluidProblem *problem)
 }
 
 
-void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *contact, long int *elid, int reload)
+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)
 {
   
   struct timespec tic, toc;
@@ -1651,6 +1651,11 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
     }
   }
   fluid_problem_compute_porosity(problem);
+  if (init){
+    for (int i=0; i<problem->mesh->n_nodes; ++i){
+      problem->oldporosity[i] = problem->porosity[i];
+    }
+  }
 }
 
 int *fluid_problem_particle_element_id(FluidProblem *problem)

fluid/fluid_problem.h

@@ -98,7 +98,7 @@ struct FluidProblem {
 // complete force on the particle (including gravity)
 void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt, double *particle_force);
 int fluid_problem_implicit_euler(FluidProblem *problem, double dt, double newton_atol, double newton_rtol, int newton_max_it);
-void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, double *contact, long int *elid, int reload);
+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_concentration_cg(FluidProblem *problem, double *concentration);
 void fluid_problem_free(FluidProblem *problem);
 FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, const char *petsc_solver_type);

python/fluid.py

@@ -330,9 +330,7 @@ class FluidProblem :
             r.tmp = tmp
             return r
         self.n_particles = mass.shape[0]
-        self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),np2c(mass),np2c(volume),np2c(position),np2c(velocity),np2c(contact),None,c_int(reload))
-        if init:
-            self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),np2c(mass),np2c(volume),np2c(position),np2c(velocity),np2c(contact),None,c_int(reload))
+        self._lib.fluid_problem_set_particles(self._fp,c_int(mass.shape[0]),np2c(mass),np2c(volume),np2c(position),np2c(velocity),np2c(contact),None,c_int(init),c_int(reload))
 
 
     def _get_matrix(self, f_name, nrow, ncol,typec=c_double) :