weak bnd : change python interface

fluid/fluid_problem.c

@@ -1554,7 +1554,7 @@ int fluid_problem_n_nodes(const FluidProblem *p)
   return p->mesh->n_nodes;
 }
 
-void fluid_problem_set_strong_boundary(FluidProblem *problem, const char *tag, int field, int row, StrongBoundaryCallback *apply)
+void fluid_problem_set_strong_boundary(FluidProblem *problem, const char *tag, int field, int row, BoundaryCallback *apply)
 {
   problem->n_strong_boundaries++;
   problem->strong_boundaries = realloc(problem->strong_boundaries,sizeof(StrongBoundary)*problem->n_strong_boundaries);
@@ -1565,7 +1565,7 @@ void fluid_problem_set_strong_boundary(FluidProblem *problem, const char *tag, i
   problem->strong_boundaries[i].row = row;
 }
 
-void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const char *bndtype, StrongBoundaryCallback *apply)
+void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const char *bndtype, BoundaryCallback *apply)
 {
   for (int i = 0; i < p->n_weak_boundaries; ++i) {
     if (strcmp(p->weak_boundaries[i].tag, tag) == 0){

fluid/fluid_problem.h

@@ -34,12 +34,12 @@
 #define N_N 4
 #endif
 
-typedef void StrongBoundaryCallback(int n, const double *x, double *bnd);
+typedef void BoundaryCallback(int n, const double *x, double *bnd);
 typedef struct {
   char *tag;
   int row;
   int field;
-  StrongBoundaryCallback *apply;
+  BoundaryCallback *apply;
 } StrongBoundary;
 
 typedef struct FluidProblem FluidProblem;
@@ -48,7 +48,7 @@ typedef struct {
   char *tag;
   f_cb *cb;
   int n_value;
-  StrongBoundaryCallback *apply;
+  BoundaryCallback *apply;
 } WeakBoundary;
 
 
@@ -100,8 +100,8 @@ int fluid_problem_n_fields(const FluidProblem *p);
 double *fluid_problem_solution(const FluidProblem *p);
 double *fluid_problem_coordinates(const FluidProblem *p);
 int fluid_problem_n_nodes(const FluidProblem *p);
-void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const char *bndtype, StrongBoundaryCallback *apply);
-void fluid_problem_set_strong_boundary(FluidProblem *p, const char *tag, int field, int row, StrongBoundaryCallback *apply);
+void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const char *bndtype, BoundaryCallback *apply);
+void fluid_problem_set_strong_boundary(FluidProblem *p, const char *tag, int field, int row, BoundaryCallback *apply);
 int fluid_problem_n_elements(const FluidProblem *p);
 int *fluid_problem_elements(const FluidProblem *p);
 double *fluid_problem_old_porosity(const FluidProblem *p);

python/fluid.py

@@ -94,30 +94,30 @@ class FluidProblem :
         """
         self._lib.fluid_problem_load_msh(self._fp, mesh_file_name.encode())
 
-    def set_weak_boundary(self, tag, bnd_type, callback_or_value=[], wall_pressure = False) :
+    def set_wall_boundary(self, tag, pressure=None) :
         """Set the weak boundaries (=normal fluxes) for the fluid problem.
 
         Keyword arguments:
         tag -- physical tag (set in the mesh.geo file) of the boundary on which the flux type is added
         bnd_type -- type of the boundary flux: wall, pressure, velocity
-        callback_or_value -- list of values (functions) to specify the weak boundary values
 
         Raises:
         ValueError -- if the number of values set in the list callback_or_value do not match what is expected for the boundary type : wall do not take any values; pressure takes a value for the pressure and concentration if any; velocity takes values for each velocity dimension and the concentration if any.
         """
-        if bnd_type.lower() == 'wall' :
-            if wall_pressure :
-                self._n_value = 1
-            else :
-                self._n_value = 0
-        elif bnd_type.lower() == 'velocity' :
-            self._n_value = self._dim + self._n_fluids-1                      # u, v, w, p, a
-        elif bnd_type.lower() == 'pressure' :
-            self._n_value = self._n_fluids                                    # p, a
-        
-        if len(callback_or_value) != self._n_value :
-            raise ValueError("Number of weak boundary values don't match the weak boundary type.")
+        if pressure :
+            raise NameError("Wall pressure not yet implemented")
+        self._lib.fluid_problem_set_weak_boundary(self._fp,tag.encode(),b"wall",None)
 
+    def set_open_boundary(self, tag, velocity=None, pressure=None, porosity=1, concentration=1):
+        """Set the weak boundaries (=normal fluxes) for the fluid problem.
+
+        Keyword arguments:
+        tag -- physical tag (set in the mesh.geo file) of the boundary on which the flux type is added
+        bnd_type -- type of the boundary flux: wall, pressure, velocity
+
+        Raises:
+        ValueError -- if the number of values set in the list callback_or_value do not match what is expected for the boundary type : wall do not take any values; pressure takes a value for the pressure and concentration if any; velocity takes values for each velocity dimension and the concentration if any.
+        """
         class Bnd :
             def __init__(self, b, dim, n_value) :
                 self._b = b
@@ -131,7 +131,19 @@ class FluidProblem :
                         v[:,i] = self._b[i](x)
                     else :
                         v[:,i] = self._b[i]
-        bndcb = BNDCB(Bnd(callback_or_value, self._dim, self._n_value).apply)
+        if porosity != 1 :
+            raise NameError("Only porosity=1 is implemented")
+        if (velocity is None and pressure is None) or (velocity is not None and pressure is not None)  :
+            raise NameError("Pressure or Velocity (but not both) should be specified at open boundaries")
+        if velocity is not None :
+            n_value = self._dim+self._n_fluids-1    # u, v, w, a
+            cb_or_value = [*velocity,concentration] if self._n_fluids == 2 else velocity
+            bnd_type = "velocity"
+        else :
+            n_value = self._n_fluids                # p, a
+            cb_or_value = [pressure,concentration] if self._n_fluids == 2 else [pressure]
+            bnd_type = "pressure"
+        bndcb = BNDCB(Bnd(cb_or_value, self._dim, n_value).apply)
         self.weak_cb_ref.append(bndcb)
         self._lib.fluid_problem_set_weak_boundary(self._fp, tag.encode(),bnd_type.encode(),bndcb)
 

testcases/injectionWGrains/Q0.05Lyon/inject2fGrains.py

@@ -30,9 +30,6 @@ import time
 import shutil
 import random
 
-
-
-
 outputdir = "outputVid1"
 outputdir1 = "outputVid"
 if not os.path.isdir(outputdir) :
@@ -73,27 +70,28 @@ rhop = 1059
 dt = .001                                  # time step
 
 shutil.copy("mesh.msh", outputdir +"/mesh.msh")
-genInitialPosition(outputdir, N, r, lx, ly, rhop)
+#genInitialPosition(outputdir, N, r, lx, ly, rhop)
 friction=0.3
-lmgc90Interface.scontactTolmgc90(outputdir, 2, 0, friction)
-p = lmgc90Interface.ParticleProblem(2)
-p.write_vtk(outputdir,0,0)
-#p = scontact.ParticleProblem(2)
-#p.read_vtk(outputdir,0)
+#lmgc90Interface.scontactTolmgc90(outputdir, 2, 0, friction)
+#p = lmgc90Interface.ParticleProblem(2)
+#p.write_vtk(outputdir,0,0)
+p = scontact.ParticleProblem(2)
+p.read_vtk(outputdir,0)
 # number of iterations between output files
 outf = 25
 ii = 0
 t = 0
 def outerBndV(x) :
+    return 0.3
     #print(0.265258*min((6*t**5-15*t**4+10*t**3),1))
     return 0.265258*min((6*t**5-15*t**4+10*t**3),1)
 
 fluid = fluid.FluidProblem(2,g,[nu0*rho0,nu1*rho1],[rho0,rho1],coeff_stab=0.001)
 fluid.load_msh("mesh.msh")
-fluid.set_weak_boundary("Bottom","Wall")
-fluid.set_weak_boundary("Lateral","Wall")
-fluid.set_weak_boundary("Top","pressure",[0,1])
-fluid.set_weak_boundary("Injection","velocity", [0,outerBndV,1])
+fluid.set_wall_boundary("Bottom",pressure=None)
+fluid.set_wall_boundary("Lateral",pressure=None)
+fluid.set_open_boundary("Top",pressure=0,porosity=1,concentration=1)
+fluid.set_open_boundary("Injection",velocity=[0,outerBndV],porosity=1,concentration=1)
 
 fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())