boundary callback

marblesbag/fluid.py

@@ -9,15 +9,24 @@ lib = CDLL(os.path.join(dir_path,"libmbfluid.so"))
 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),("field",c_int),("v",c_double)]
-
+    _fields_ = [("tag",c_char_p),("field",c_int),("apply",BNDCB)]
 
 class fluid_problem :
 
     def __init__(self, mesh_file_name, mu, rho, alpha, autoEpsilon, strong_boundaries):
         nsb = len(strong_boundaries)
-        asb = (StrongBoundary*nsb)(*[StrongBoundary(i[0].encode(),i[1],i[2]) for i in strong_boundaries])
+        def apply_bnd(n, xp, vp, b) :
+            x = numpy.frombuffer(cast(xp, POINTER(int(n)*2*c_double)).contents)
+            v = numpy.frombuffer(cast(vp, POINTER(int(n)*c_double)).contents)
+            if callable(b) :
+                v[:] = b(x)
+            else :
+                v[:] = b
+
+        asb = (StrongBoundary*nsb)(*[StrongBoundary(i[0].encode(),i[1],BNDCB(lambda n,x,v : apply_bnd(n,x,v,i[2]))) for i in strong_boundaries])
+        self.asb = asb
         lib.fluid_problem_new.restype = c_void_p
         self._fp = c_void_p(lib.fluid_problem_new(mesh_file_name.encode(), c_double(mu), c_double(rho), c_double(alpha), c_bool(autoEpsilon), nsb, asb))
         if self._fp == None :

marblesbag/src/fluid_problem.c

@@ -272,8 +272,18 @@ static void apply_boundary_conditions(const Mesh *mesh, double *solution, int nB
     if (iphys == mesh->n_phys){
       printf("Boundary tag \"%s\" not found.", bnd->tag);
     }
+    double *x = malloc(mesh->phys_n_nodes[iphys]*sizeof(double)*2);
+    double *v = malloc(mesh->phys_n_nodes[iphys]*sizeof(double));
+    for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i){
+      int j = mesh->phys_nodes[iphys][i];
+      x[i*2+0] = mesh->x[j*2+0];
+      x[i*2+1] = mesh->x[j*2+1];
+    }
+    bnd->apply(mesh->phys_n_nodes[iphys], x, v);
     for (int i = 0; i < mesh->phys_n_nodes[iphys]; ++i)
-      solution[mesh->phys_nodes[iphys][i]*n_fields+bnd->field] = bnd->v;
+      solution[mesh->phys_nodes[iphys][i]*n_fields+bnd->field] = v[i];
+    free(x);
+    free(v);
   }
 }
 
@@ -390,7 +400,7 @@ FluidProblem *fluid_problem_new(const char *mesh_file_name, double mu, double rh
   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].v = strong_boundaries[i].v;
+    problem->strong_boundaries[i].apply = strong_boundaries[i].apply;
     problem->strong_boundaries[i].field = strong_boundaries[i].field;
   }
   problem->porosity = malloc(mesh->n_nodes*sizeof(double));

marblesbag/src/fluid_problem.h

@@ -8,7 +8,7 @@
 typedef struct {
   char *tag;
   int field;
-  double v;
+  void (*apply)(int n, const double *x, double *bnd);
 }StrongBoundary;
 
 typedef struct {