change syntax for strong boundaries

python/marblesbag/fluid.py

@@ -31,33 +31,19 @@ signal.signal(signal.SIGINT, signal.SIG_DFL)
 
 
 BNDCB = CFUNCTYPE(None,c_int,POINTER(c_double),POINTER(c_double))
-class StrongBoundary(Structure):
-    _fields_ = [("tag",c_char_p),("row",c_int),("field",c_int),("apply",BNDCB)]
 
 class fluid_problem :
 
-    def __init__(self, g, mu, rho, strong_boundaries):
-        nsb = len(strong_boundaries)
-        class Bnd :
-            def __init__(self, b) :
-                self._b = b
-            def apply(self, n, xp, vp) :
-                x = numpy.frombuffer(cast(xp, POINTER(int(n)*2*c_double)).contents).reshape([n,2])
-                v = numpy.frombuffer(cast(vp, POINTER(int(n)*c_double)).contents)
-                if callable(self._b) :
-                    v[:] = self._b(x)
-                else :
-                    v[:] = self._b
+    def __init__(self, g, mu, rho):
+        self.strong_cb_ref = []
         def np2c(a) :
             tmp = numpy.require(a,"float","C")
             r = c_void_p(tmp.ctypes.data)
             r.tmp = tmp
             return r
-        asb = (StrongBoundary*nsb)(*[StrongBoundary(i[0].encode(),i[1],i[2],BNDCB(Bnd(i[3]).apply)) for i in strong_boundaries])
-        self.asb = asb
         lib.fluid_problem_new.restype = c_void_p
         n_fluids = numpy.require(rho,"float","C").reshape([-1]).shape[0]
-        self._fp = c_void_p(lib.fluid_problem_new(c_double(g), n_fluids, np2c(mu), np2c(rho), nsb, asb))
+        self._fp = c_void_p(lib.fluid_problem_new(c_double(g), n_fluids, np2c(mu), np2c(rho)))
         if self._fp == None :
             raise NameError("cannot create fluid problem")
 
@@ -71,6 +57,23 @@ class fluid_problem :
     def set_weak_boundary(self, tag, bndtype) :
         lib.fluid_problem_set_weak_boundary(self._fp, tag.encode(),bndtype.encode())
 
+    def set_strong_boundary(self, tag, field,callback_or_value, row=None) :
+        if row is None :
+            row = field
+        class Bnd :
+            def __init__(self, b) :
+                self._b = b
+            def apply(self, n, xp, vp) :
+                x = numpy.frombuffer(cast(xp, POINTER(int(n)*2*c_double)).contents).reshape([n,2])
+                v = numpy.frombuffer(cast(vp, POINTER(int(n)*c_double)).contents)
+                if callable(self._b) :
+                    v[:] = self._b(x)
+                else :
+                    v[:] = self._b
+        bndcb = BNDCB(Bnd(callback_or_value).apply)
+        self.strong_cb_ref.append(bndcb)
+        lib.fluid_problem_set_strong_boundary(self._fp, tag.encode(), field, row, bndcb)
+
     def adapt_mesh(self, lcmax, lcmin, n_el) :
         lib.fluid_problem_adapt_mesh(self._fp, c_double(lcmax), c_double(lcmin), c_double(n_el))
 

src/fluid_problem.c

@@ -1003,7 +1003,7 @@ static void fluid_problem_compute_porosity(FluidProblem *problem)
   }
 }
 
-FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho, int n_strong_boundaries, StrongBoundary *strong_boundaries)
+FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho)
 {
   if (n_fluids != 1 && n_fluids != 2) {
     printf("error : only 1 or 2 fluids are supported\n");
@@ -1019,6 +1019,8 @@ FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho,
   FluidProblem *problem = (FluidProblem*)malloc(sizeof(FluidProblem));
   problem->n_fluids = n_fluids;
   int n_fields = fluid_problem_n_fields(problem);
+  problem->strong_boundaries = NULL;
+  problem->n_strong_boundaries = 0;
   problem->weak_boundaries = NULL;
   problem->n_weak_boundaries = 0;
   problem->mu = (double*)malloc(sizeof(double)*problem->n_fluids);
@@ -1028,14 +1030,6 @@ FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho,
     problem->mu[i] = mu[i];
   }
   problem->g = g;
-  problem->strong_boundaries = (StrongBoundary*)malloc(sizeof(StrongBoundary)*n_strong_boundaries);
-  problem->n_strong_boundaries = n_strong_boundaries;
-  for (int i = 0; i < n_strong_boundaries; ++i) {
-    problem->strong_boundaries[i].tag = strdup(strong_boundaries[i].tag);
-    problem->strong_boundaries[i].apply = strong_boundaries[i].apply;
-    problem->strong_boundaries[i].field = strong_boundaries[i].field;
-    problem->strong_boundaries[i].row = strong_boundaries[i].row;
-  }
   problem->node_volume = NULL;
   problem->mesh_tree = NULL;
   problem->mesh = NULL;
@@ -1464,6 +1458,17 @@ 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)
+{
+  problem->n_strong_boundaries++;
+  problem->strong_boundaries = realloc(problem->strong_boundaries,sizeof(StrongBoundary)*problem->n_strong_boundaries);
+  int i = problem->n_strong_boundaries-1;
+  problem->strong_boundaries[i].tag = strdup(tag);
+  problem->strong_boundaries[i].apply = apply;
+  problem->strong_boundaries[i].field = field;
+  problem->strong_boundaries[i].row = row;
+}
+
 void fluid_problem_set_weak_boundary(FluidProblem *p, const char *tag, const char *bndtype)
 {
   for (int i = 0; i < p->n_weak_boundaries; ++i) {

src/fluid_problem.h

@@ -34,11 +34,12 @@
 #define N_N 4
 #endif
 
+typedef void StrongBoundaryCallback(int n, const double *x, double *bnd);
 typedef struct {
   char *tag;
   int row;
   int field;
-  void (*apply)(int n, const double *x, double *bnd);
+  StrongBoundaryCallback *apply;
 } StrongBoundary;
 
 typedef struct FluidProblem FluidProblem;
@@ -89,7 +90,7 @@ void fluid_problem_compute_node_particle_force(FluidProblem *problem, double dt,
 int fluid_problem_implicit_euler(FluidProblem *problem, double dt);
 void fluid_problem_set_particles(FluidProblem *problem, int n, double *mass, double *volume, double *position, double *velocity, long int *elid, int reload);
 void fluid_problem_free(FluidProblem *problem);
-FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho, int n_strong_boundaries, StrongBoundary *strong_boundaries);
+FluidProblem *fluid_problem_new(double g, int n_fluids, double *mu, double *rho);
 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);
 int *fluid_problem_particle_element_id(FluidProblem *problem);
@@ -100,4 +101,5 @@ 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);
+void fluid_problem_set_strong_boundary(FluidProblem *p, const char *tag, int field, int row, StrongBoundaryCallback *apply);
 #endif

testcases/injectionWGrains/inject2fGrains.py

@@ -83,17 +83,19 @@ t = 0
 def outerBndV(x) :
     return max(1*t,0.1)
 
-strong_boundaries = [("Top",2,2,0.),("Injection",1,1,outerBndV),("Injection",3,3,1),("Injection",0,0,0)
-#        ("Top",3,3,1),("Lateral",3,3,1),("Bottom",3,3,1)
-        ]#,("Lateral",0,0,0),("Bottom",1,1,0)]
-fluid = fluid.fluid_problem(g,[nu0*rho0,nu1*rho1],[rho0,rho1],strong_boundaries)
+fluid = fluid.fluid_problem(g,[nu0*rho0,nu1*rho1],[rho0,rho1])
 fluid.load_msh("mesh.msh")
+
+fluid.set_strong_boundary("Top",2,2,0)
+fluid.set_strong_boundary("Injection",1,1,outerBndV)
+fluid.set_strong_boundary("Injection",3,3,1)
+fluid.set_strong_boundary("Injection",0,0,0)
+
 fluid.set_weak_boundary("Bottom","Wall")
 fluid.set_weak_boundary("Lateral","Wall")
 fluid.set_weak_boundary("Top","Outflow")
 fluid.set_weak_boundary("Injection","Inflow")
 
-
 fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity())
 
 #set initial_condition

validation/cavity/cavity.py

@@ -69,19 +69,17 @@ class Poiseuille(unittest.TestCase) :
 
         #Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
         #Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
-        strong_boundaries = [
-                ("Top",0,0,vUp),
-                ("Top",1,1,0),
-                ("Top",2,2,0),
-                ("Bottom",0,0,0),
-                ("Bottom",1,1,0),
-                ("Left",0,0,0),
-                ("Left",1,1,0),
-                ("Right",0,0,0),
-                ("Right",1,1,0),
-                ]
-        fluid = mbfluid.fluid_problem(g,nu*rho,rho,strong_boundaries)
+        fluid = mbfluid.fluid_problem(g,nu*rho,rho)
         fluid.load_msh("mesh.msh")
+        fluid.set_strong_boundary("Top",0,vUp)
+        fluid.set_strong_boundary("Top",1,0.)
+        fluid.set_strong_boundary("Top",2,0.)
+        fluid.set_strong_boundary("Bottom",0,0.)
+        fluid.set_strong_boundary("Bottom",1,0.)
+        fluid.set_strong_boundary("Left",0,0.)
+        fluid.set_strong_boundary("Left",1,0.)
+        fluid.set_strong_boundary("Right",0,0.)
+        fluid.set_strong_boundary("Right",1,0.)
         fluid.set_weak_boundary("Left","Wall")
         fluid.set_weak_boundary("Bottom","Wall")
         fluid.set_weak_boundary("Top","Wall")

validation/oneGrain/oneGrain.py

@@ -76,9 +76,9 @@ class Poiseuille(unittest.TestCase) :
 
             outf = 1                                       # number of iterations between output files
 
-            strong_boundaries = [("Top",2,2,0.)]#,("Lateral",0,0,0),("Bottom",1,1,0)]
-            fluid = mbfluid.fluid_problem(g,nu*rho,rho,strong_boundaries)
+            fluid = mbfluid.fluid_problem(g,nu*rho,rho)
             fluid.load_msh("mesh.msh")
+            fluid.set_strong_boundary("Top",2,0.)
             fluid.set_weak_boundary("Bottom","Wall")
             fluid.set_weak_boundary("Lateral","Wall")
             fluid.set_weak_boundary("Top","Wall")

validation/poiseuille/poiseuille.py

@@ -65,9 +65,17 @@ class Poiseuille(unittest.TestCase) :
         outf = 1                                       # number of iterations between output files
 
         #Object fluid creation + Boundary condition of the fluid (field 0 is horizontal velocity; field 1 is vertical velocity; field 2 is pressure)
-        #Format: strong_boundaries = [(Boundary tag, Fluid field, Value)
-        strong_boundaries = [("RightUp",2,2,0.),("RightDown",2,2,0.),("Top",0,0,0),("Bottom",0,0,0),("LeftUp",1,1,0),("RightUp",1,1,0),("LeftDown",1,1,0),("RightDown",1,1,0),("LeftUp",0,0,lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1])),("LeftDown",0,0,lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]))]#,("Lateral",0,0,0),("Bottom",1,1,0)]
-        fluid = mbfluid.fluid_problem(g,nu*rho,rho,strong_boundaries)
+        fluid = mbfluid.fluid_problem(g,nu*rho,rho)
+        fluid.set_strong_boundary("RightUp",2,0.)
+        fluid.set_strong_boundary("RightDown",2,0.)
+        fluid.set_strong_boundary("Top",0,0)
+        fluid.set_strong_boundary("Bottom",0,0)
+        fluid.set_strong_boundary("LeftUp",1,0)
+        fluid.set_strong_boundary("RightUp",1,0)
+        fluid.set_strong_boundary("LeftDown",1,0)
+        fluid.set_strong_boundary("RightDown",1,0)
+        fluid.set_strong_boundary("LeftUp",0,lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]))
+        fluid.set_strong_boundary("LeftDown",0,lambda x : 1/(20*mu)*x[:,1]*(1-x[:, 1]))
         fluid.load_msh("mesh.msh")
         fluid.set_weak_boundary("LeftUp","Inflow")
         fluid.set_weak_boundary("LeftDown","Inflow")