add ip factor for weak bnd

fluid/fluid_problem.c

@@ -210,7 +210,7 @@ static void f_boundary(WeakBoundary *wbnd, FluidProblem *problem,const double *n
     f00[(U+id)*n_fields+P] += c*(pid<0?0:-n[id]);
   }
   double c_du_o_c[D][D];
-  double sigma = (1+D)/(D*h)*mu*N_N;
+  double sigma = problem->ip_factor*(1+D)/(D*h)*mu*N_N;
   for (int iD = 0; iD < D; ++iD) {
     for (int jD = 0; jD < D; ++jD)
       c_du_o_c[iD][jD] = ds[(U+iD)*D+jD] -u[iD]/c*dc[jD];
@@ -1248,7 +1248,7 @@ static void compute_porosity(Mesh *mesh, double *node_volume, double *porosity,
   }
 }
 
-FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, int temporal, int advection)
+FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor, int temporal, int advection)
 {
   if (n_fluids != 1 && n_fluids != 2) {
     printf("error : only 1 or 2 fluids are supported\n");
@@ -1257,6 +1257,7 @@ FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho
   problem->n_fluids = n_fluids;
   problem->sigma = sigma;
   problem->drag_coeff_factor = drag_coeff_factor;
+  problem->ip_factor = ip_factor;
   problem->stab_param = 0;
   problem->drag_in_stab = drag_in_stab;
   problem->reduced_gravity = 0;

fluid/fluid_problem.h

@@ -59,6 +59,7 @@ struct FluidProblem {
   double *g;
   double coeffStab;
   double drag_coeff_factor;
+  double ip_factor;
   double sigma;
   double kinetic_energy;
   Mesh *mesh;
@@ -107,7 +108,7 @@ void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, dou
 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);
-FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor,int temporal, int advection);
+FluidProblem *fluid_problem_new(double *g, int n_fluids, double *mu, double *rho, double sigma, double coeffStab, double volume_drag, double quadratic_drag, int drag_in_stab, double drag_coeff_factor, double ip_factor,int temporal, int advection);
 void fluid_problem_load_msh(FluidProblem *problem, const char *mesh_file_name);
 void fluid_problem_adapt_mesh(FluidProblem *problem, double lcmax, double lcmin, double n_el, double cmax, double cmin, int old_n_particles, double *old_particle_position, double *old_particle_volume);
 int *fluid_problem_particle_element_id(FluidProblem *problem);

python/fluid.py

@@ -77,7 +77,7 @@ def _np2c(a,dtype=float) :
 class FluidProblem :
     """Create numerical representation of the fluid."""
 
-    def __init__(self, dim, g, mu, rho, sigma=0, coeff_stab=0.01, volume_drag=0., quadratic_drag=0., petsc_solver_type="", drag_in_stab=1, drag_coefficient_factor=1, temporal=True, advection=True):
+    def __init__(self, dim, g, mu, rho, sigma=0, coeff_stab=0.01, volume_drag=0., quadratic_drag=0., petsc_solver_type="", drag_in_stab=1, drag_coefficient_factor=1, ip_factor=1, temporal=True, advection=True):
         """Build the fluid structure.
 
         Keyword arguments:
@@ -112,7 +112,7 @@ class FluidProblem :
         n_fluids = np.require(rho,"float","C").reshape([-1]).shape[0]
         self._n_fluids = n_fluids
         self._dim = dim
-        self._fp = c_void_p(self._lib.fluid_problem_new(_np2c(g), n_fluids, _np2c(mu), _np2c(rho), c_double(sigma), c_double(coeff_stab), c_double(volume_drag), c_double(quadratic_drag), c_int(drag_in_stab),c_double(drag_coefficient_factor),c_int(temporal), c_int(advection)))
+        self._fp = c_void_p(self._lib.fluid_problem_new(_np2c(g), n_fluids, _np2c(mu), _np2c(rho), c_double(sigma), c_double(coeff_stab), c_double(volume_drag), c_double(quadratic_drag), c_int(drag_in_stab),c_double(drag_coefficient_factor),c_double(ip_factor),c_int(temporal), c_int(advection)))
         if self._fp == None :
             raise NameError("Cannot create fluid problem.")