diff --git a/boundaries/__init__.py b/boundaries/__init__.py
index a7ab44b751d9fceabb2e4de7a8b8b47b0afcfcdc..0192bbc2d984d214a00da0d35d4fde50567a018a 100644
--- a/boundaries/__init__.py
+++ b/boundaries/__init__.py
@@ -1 +1,3 @@
 from pystencils.boundaries.boundaryhandling import BoundaryHandling
+from pystencils.boundaries.boundaryconditions import Neumann
+from pystencils.boundaries.inkernel import add_neumann_boundary
diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index 563940b30324271ce805d0656cee102ce1ee7a99..5032faff9aee6d7947d0395e8df8f52fff3f5885 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -29,7 +29,7 @@ class Boundary(object):
         """Return a list of (name, type) tuples for additional data items required in this boundary
         These data items can either be initialized in separate kernel see additional_data_kernel_init or by
         Python callbacks - see additional_data_callback """
-        return []
+        return ()
 
     @property
     def additional_data_init_callback(self):
diff --git a/boundaries/createindexlistcython.pyx b/boundaries/createindexlistcython.pyx
index fc3b2dc928437457d6efc22c4eff35cbb6eb1bb7..0ce36ae81a78e6788db3ff795f79777553904411 100644
--- a/boundaries/createindexlistcython.pyx
+++ b/boundaries/createindexlistcython.pyx
@@ -28,7 +28,7 @@ def create_boundary_index_list_2d(object[IntegerType, ndim=2] flag_field,
     for y in range(nr_of_ghost_layers, ys - nr_of_ghost_layers):
         for x in range(nr_of_ghost_layers, xs - nr_of_ghost_layers):
             if flag_field[x, y] & fluid_mask:
-                for dirIdx in range(1, num_directions):
+                for dirIdx in range(num_directions):
                     dx = stencil[dirIdx,0]
                     dy = stencil[dirIdx,1]
                     if flag_field[x + dx, y + dy] & boundary_mask:
@@ -38,25 +38,25 @@ def create_boundary_index_list_2d(object[IntegerType, ndim=2] flag_field,
 
 @cython.boundscheck(False) # turn off bounds-checking for entire function
 @cython.wraparound(False)  # turn off negative index wrapping for entire function
-def create_boundary_index_list_3d(object[IntegerType, ndim=3] flagField,
-                                  int nrOfGhostLayers, IntegerType boundaryMask, IntegerType fluidMask,
+def create_boundary_index_list_3d(object[IntegerType, ndim=3] flag_field,
+                                  int nr_of_ghost_layers, IntegerType boundary_mask, IntegerType fluid_mask,
                                   object[int, ndim=2] stencil):
     cdef int xs, ys, zs, x, y, z
     cdef int dirIdx, num_directions, dx, dy, dz
 
-    xs, ys, zs = flagField.shape
+    xs, ys, zs = flag_field.shape
     boundary_index_list = []
     num_directions = stencil.shape[0]
 
-    for z in range(nrOfGhostLayers, zs-nrOfGhostLayers):
-        for y in range(nrOfGhostLayers,ys-nrOfGhostLayers):
-            for x in range(nrOfGhostLayers,xs-nrOfGhostLayers):
-                if flagField[x, y, z] & fluidMask:
-                    for dirIdx in range(1, num_directions):
+    for z in range(nr_of_ghost_layers, zs - nr_of_ghost_layers):
+        for y in range(nr_of_ghost_layers, ys - nr_of_ghost_layers):
+            for x in range(nr_of_ghost_layers, xs - nr_of_ghost_layers):
+                if flag_field[x, y, z] & fluid_mask:
+                    for dirIdx in range(num_directions):
                         dx = stencil[dirIdx,0]
                         dy = stencil[dirIdx,1]
                         dz = stencil[dirIdx,2]
-                        if flagField[x + dx, y + dy, z + dz] & boundaryMask:
+                        if flag_field[x + dx, y + dy, z + dz] & boundary_mask:
                             boundary_index_list.append((x,y,z, dirIdx))
     return boundary_index_list