Big refactoring of lbmpy boundary handling

- boundary conditions can now define their own state
- boundary classes instead of functions
- no boundary name any more
  (instead the identity is defined via instances)
- updated tutorials accordingly
- boundary handling now gets pdf numpy array instead of pystencils field

boundaries/__init__.py

-from lbmpy.boundaries.boundaryconditions import noSlip, ubb, fixedDensity
+from lbmpy.boundaries.boundaryconditions import NoSlip, NoSlipFullWay, UBB, FixedDensity
 from lbmpy.boundaries.boundaryhandling import BoundaryHandling

boundaries/boundaryconditions.py

@@ -4,81 +4,160 @@ from lbmpy.simplificationfactory import createSimplificationStrategy
 from pystencils.sympyextensions import getSymmetricPart
 from pystencils import Field
 from lbmpy.boundaries.boundaryhandling import offsetFromDir, weightOfDirection, invDir
+from pystencils.types import createTypeFromString
 
 
-def noSlip(pdfField, direction, lbMethod):
-    """No-Slip, simple bounce back boundary condition, enforcing zero velocity at obstacle"""
-    neighbor = offsetFromDir(direction, lbMethod.dim)
-    inverseDir = invDir(direction)
-    return [sp.Eq(pdfField[neighbor](inverseDir), pdfField(direction))]
+class Boundary(object):
+    """Base class all boundaries should derive from"""
 
+    def __call__(self, pdfField, directionSymbol, lbMethod, indexField):
+        """
+        This function defines the boundary behavior and must therefore be implemented by all boundaries.
+        Here the boundary is defined as a list of sympy equations, from which a boundary kernel is generated.
+        :param pdfField: pystencils field describing the pdf. The current cell is cell next to the boundary,
+                         which is influenced by the boundary cell i.e. has a link from the boundary cell to
+                         itself.
+        :param directionSymbol: a sympy symbol that can be used as index to the pdfField. It describes
+                                the direction pointing from the fluid to the boundary cell 
+        :param lbMethod: an instance of the LB method used. Use this to adapt the boundary to the method 
+                         (e.g. compressiblity)
+        :param indexField: the boundary index field that can be used to retrieve and update boundary data
+        :return: list of sympy equations
+        """
+        raise NotImplementedError("Boundary class has to overwrite __call__")
+
+    @property
+    def additionalData(self):
+        return []
+
+    @property
+    def name(self):
+        return type(self).__name__
+
+    def additionalDataInit(self, idxArray):
+        return []
+
+
+class NoSlip(Boundary):
+    """No-Slip, (half-way) simple bounce back boundary condition, enforcing zero velocity at obstacle"""
+    def __call__(self, pdfField, directionSymbol, lbMethod, **kwargs):
+        neighbor = offsetFromDir(directionSymbol, lbMethod.dim)
+        inverseDir = invDir(directionSymbol)
+        return [sp.Eq(pdfField[neighbor](inverseDir), pdfField(directionSymbol))]
+
+    def __hash__(self):
+        # All boundaries of these class behave equal -> should also be equal
+        return hash("NoSlip")
+
+    def __eq__(self, other):
+        return type(other) == NoSlip
+
+
+class NoSlipFullWay(Boundary):
+    """Full-way bounce back"""
+
+    @property
+    def additionalData(self):
+        return [('lastValue', createTypeFromString("double"))]
+
+    def additionalDataInit(self, pdfField, directionSymbol, indexField, **kwargs):
+        return [sp.Eq(indexField('lastValue'), pdfField(directionSymbol))]
+
+    def __call__(self, pdfField, directionSymbol, lbMethod, indexField, **kwargs):
+        neighbor = offsetFromDir(directionSymbol, lbMethod.dim)
+        inverseDir = invDir(directionSymbol)
+        return [sp.Eq(pdfField[neighbor](inverseDir), indexField('lastValue')),
+                sp.Eq(indexField('lastValue'), pdfField(directionSymbol))]
+
+    def __hash__(self):
+        # All boundaries of these class behave equal -> should also be equal
+        return hash("NoSlipFullWay")
+
+    def __eq__(self, other):
+        return type(other) == NoSlipFullWay
+
+
+class UBB(Boundary):
 
-def ubb(pdfField, direction, lbMethod, velocity, adaptVelocityToForce=False):
     """Velocity bounce back boundary condition, enforcing specified velocity at obstacle"""
 
-    assert len(velocity) == lbMethod.dim, \
-        "Dimension of velocity (%d) does not match dimension of LB method (%d)" % (len(velocity), lbMethod.dim)
-    neighbor = offsetFromDir(direction, lbMethod.dim)
-    inverseDir = invDir(direction)
+    def __init__(self, velocity, adaptVelocityToForce=False):
+        self._velocity = velocity
+        self._adaptVelocityToForce = adaptVelocityToForce
 
-    velocity = tuple(v_i.getShifted(*neighbor) if isinstance(v_i, Field.Access) else v_i for v_i in velocity)
+    def __call__(self, pdfField, directionSymbol, lbMethod, **kwargs):
+        vel = self._velocity
+        direction = directionSymbol
 
-    if adaptVelocityToForce:
-        cqc = lbMethod.conservedQuantityComputation
-        shiftedVelEqs = cqc.equilibriumInputEquationsFromInitValues(velocity=velocity)
-        velocity = [eq.rhs for eq in shiftedVelEqs.extract(cqc.firstOrderMomentSymbols).mainEquations]
+        assert len(vel) == lbMethod.dim, \
+            "Dimension of velocity (%d) does not match dimension of LB method (%d)" % (len(vel), lbMethod.dim)
+        neighbor = offsetFromDir(direction, lbMethod.dim)
+        inverseDir = invDir(direction)
+
+        velocity = tuple(v_i.getShifted(*neighbor) if isinstance(v_i, Field.Access) else v_i for v_i in vel)
+
+        if self._adaptVelocityToForce:
+            cqc = lbMethod.conservedQuantityComputation
+            shiftedVelEqs = cqc.equilibriumInputEquationsFromInitValues(velocity=velocity)
+            velocity = [eq.rhs for eq in shiftedVelEqs.extract(cqc.firstOrderMomentSymbols).mainEquations]
+
+        c_s_sq = sp.Rational(1, 3)
+        velTerm = 2 / c_s_sq * sum([d_i * v_i for d_i, v_i in zip(neighbor, velocity)]) * weightOfDirection(direction)
+
+        # Better alternative: in conserved value computation
+        # rename what is currently called density to "virtualDensity"
+        # provide a new quantity density, which is constant in case of incompressible models
+        if not lbMethod.conservedQuantityComputation.zeroCenteredPdfs:
+            cqc = lbMethod.conservedQuantityComputation
+            densitySymbol = sp.Symbol("rho")
+            pdfFieldAccesses = [pdfField(i) for i in range(len(lbMethod.stencil))]
+            densityEquations = cqc.outputEquationsFromPdfs(pdfFieldAccesses, {'density': densitySymbol})
+            densitySymbol = lbMethod.conservedQuantityComputation.definedSymbols()['density']
+            result = densityEquations.allEquations
+            result += [sp.Eq(pdfField[neighbor](inverseDir),
+                             pdfField(direction) - velTerm * densitySymbol)]
+            return result
+        else:
+            return [sp.Eq(pdfField[neighbor](inverseDir),
+                          pdfField(direction) - velTerm)]
+
+
+class FixedDensity(Boundary):
+
+    def __init__(self, density):
+        self._density = density
 
-    c_s_sq = sp.Rational(1, 3)
-    velTerm = 2 / c_s_sq * sum([d_i * v_i for d_i, v_i in zip(neighbor, velocity)]) * weightOfDirection(direction)
+    def __call__(self, pdfField, directionSymbol, lbMethod, **kwargs):
+        """Boundary condition that fixes the density/pressure at the obstacle"""
+
+        def removeAsymmetricPartOfMainEquations(eqColl, dofs):
+            newMainEquations = [sp.Eq(e.lhs, getSymmetricPart(e.rhs, dofs)) for e in eqColl.mainEquations]
+            return eqColl.copy(newMainEquations)
+
+        neighbor = offsetFromDir(directionSymbol, lbMethod.dim)
+        inverseDir = invDir(directionSymbol)
 
-    # Better alternative: in conserved value computation
-    # rename what is currently called density to "virtualDensity"
-    # provide a new quantity density, which is constant in case of incompressible models
-    if not lbMethod.conservedQuantityComputation.zeroCenteredPdfs:
         cqc = lbMethod.conservedQuantityComputation
-        densitySymbol = sp.Symbol("rho")
-        pdfFieldAccesses = [pdfField(i) for i in range(len(lbMethod.stencil))]
-        densityEquations = cqc.outputEquationsFromPdfs(pdfFieldAccesses, {'density': densitySymbol})
-        densitySymbol = lbMethod.conservedQuantityComputation.definedSymbols()['density']
-        result = densityEquations.allEquations
-        result += [sp.Eq(pdfField[neighbor](inverseDir),
-                         pdfField(direction) - velTerm * densitySymbol)]
-        return result
-    else:
-        return [sp.Eq(pdfField[neighbor](inverseDir),
-                      pdfField(direction) - velTerm)]
-
-
-def fixedDensity(pdfField, direction, lbMethod, density):
-    """Boundary condition that fixes the density/pressure at the obstacle"""
-
-    def removeAsymmetricPartOfMainEquations(eqColl, dofs):
-        newMainEquations = [sp.Eq(e.lhs, getSymmetricPart(e.rhs, dofs)) for e in eqColl.mainEquations]
-        return eqColl.copy(newMainEquations)
-
-    neighbor = offsetFromDir(direction, lbMethod.dim)
-    inverseDir = invDir(direction)
-
-    cqc = lbMethod.conservedQuantityComputation
-    velocity = cqc.definedSymbols()['velocity']
-    symmetricEq = removeAsymmetricPartOfMainEquations(lbMethod.getEquilibrium(), dofs=velocity)
-    substitutions = {sym: pdfField(i) for i, sym in enumerate(lbMethod.preCollisionPdfSymbols)}
-    symmetricEq = symmetricEq.copyWithSubstitutionsApplied(substitutions)
-
-    simplification = createSimplificationStrategy(lbMethod)
-    symmetricEq = simplification(symmetricEq)
-
-    densitySymbol = cqc.definedSymbols()['density']
-
-    densityEq = cqc.equilibriumInputEquationsFromInitValues(density=density).insertSubexpressions().mainEquations[0]
-    assert densityEq.lhs == densitySymbol
-    transformedDensity = densityEq.rhs
-
-    conditions = [(eq_i.rhs, sp.Equality(direction, i))
-                  for i, eq_i in enumerate(symmetricEq.mainEquations)] + [(0, True)]
-    eq_component = sp.Piecewise(*conditions)
-
-    subExprs = [sp.Eq(eq.lhs, transformedDensity if eq.lhs == densitySymbol else eq.rhs)
-                for eq in symmetricEq.subexpressions]
-    return subExprs + [sp.Eq(pdfField[neighbor](inverseDir), 2 * eq_component - pdfField(direction))]
+        velocity = cqc.definedSymbols()['velocity']
+        symmetricEq = removeAsymmetricPartOfMainEquations(lbMethod.getEquilibrium(), dofs=velocity)
+        substitutions = {sym: pdfField(i) for i, sym in enumerate(lbMethod.preCollisionPdfSymbols)}
+        symmetricEq = symmetricEq.copyWithSubstitutionsApplied(substitutions)
+
+        simplification = createSimplificationStrategy(lbMethod)
+        symmetricEq = simplification(symmetricEq)
+
+        densitySymbol = cqc.definedSymbols()['density']
+
+        density = self._density
+        densityEq = cqc.equilibriumInputEquationsFromInitValues(density=density).insertSubexpressions().mainEquations[0]
+        assert densityEq.lhs == densitySymbol
+        transformedDensity = densityEq.rhs
+
+        conditions = [(eq_i.rhs, sp.Equality(directionSymbol, i))
+                      for i, eq_i in enumerate(symmetricEq.mainEquations)] + [(0, True)]
+        eq_component = sp.Piecewise(*conditions)
+
+        subExprs = [sp.Eq(eq.lhs, transformedDensity if eq.lhs == densitySymbol else eq.rhs)
+                    for eq in symmetricEq.subexpressions]
+        return subExprs + [sp.Eq(pdfField[neighbor](inverseDir), 2 * eq_component - pdfField(directionSymbol))]
 

boundaries/boundaryhandling.py

@@ -13,10 +13,9 @@ WEIGHTS_SYMBOL = TypedSymbol("weights", "double")
 
 class BoundaryHandling(object):
     class BoundaryInfo(object):
-        def __init__(self, name, flag, function, kernel, ast):
-            self.name = name
+        def __init__(self, flag, object, kernel, ast):
             self.flag = flag
-            self.function = function
+            self.object = object
             self.kernel = kernel
             self.ast = ast
 
@@ -24,29 +23,24 @@ class BoundaryHandling(object):
         """
         Class for managing boundary kernels
 
-        :param pdfField: either pdf numpy array (including ghost layers), or pystencils.Field
+        :param pdfField: pdf numpy array including ghost layers
         :param domainShape: domain size without ghost layers
         :param lbMethod: lattice Boltzmann method
         :param ghostLayers: number of ghost layers
         :param target: either 'cpu' or 'gpu'
         """
-        if isinstance(pdfField, np.ndarray):
-            symbolicPdfField = Field.createFromNumpyArray('pdfs', pdfField, indexDimensions=1)
-            assert pdfField.shape[:-1] == tuple(d + 2*ghostLayers for d in domainShape)
-        elif isinstance(pdfField, Field):
-            symbolicPdfField = pdfField
-        else:
-            raise ValueError("pdfField has to be either a numpy array or a pystencils.Field")
+        symbolicPdfField = Field.createFromNumpyArray('pdfs', pdfField, indexDimensions=1)
+        assert pdfField.shape[:-1] == tuple(d + 2*ghostLayers for d in domainShape)
 
         self._symbolicPdfField = symbolicPdfField
+        self._pdfField = pdfField
         self._shapeWithGhostLayers = [d + 2 * ghostLayers for d in domainShape]
         self._fluidFlag = 2 ** 0
         self.flagField = np.full(self._shapeWithGhostLayers, self._fluidFlag, dtype=np.int32)
         self._ghostLayers = ghostLayers
         self._lbMethod = lbMethod
 
-        self._boundaryInfos = []
-        self._nameToBoundary = {}
+        self._boundaryInfos = {}
         self._periodicityKernels = []
 
         self._dirty = False
@@ -66,13 +60,12 @@ class BoundaryHandling(object):
         """Flag that is set where the lattice Boltzmann update should happen"""
         return self._fluidFlag
 
-    def getFlag(self, name):
-        """Flag that represents the boundary with given name. Raises KeyError if no such boundary exists."""
-        return self._nameToBoundary[name].flag
+    def getFlag(self, boundaryObject):
+        """Flag that represents the given boundary."""
+        return self._boundaryInfos[boundaryObject].flag
 
-    def getBoundaryNames(self):
-        """List of names of all registered boundary conditions"""
-        return [b.name for b in self._boundaryInfos]
+    def getBoundaries(self):
+        return [b.object for b in self._boundaryInfos.values()]
 
     def setPeriodicity(self, x=False, y=False, z=False):
         """Enable periodic boundary conditions at the border of the domain"""
@@ -82,37 +75,24 @@ class BoundaryHandling(object):
         self._periodicity = [x, y, z]
         self._compilePeriodicityKernels()
 
-    def hasBoundary(self, name):
+    def hasBoundary(self, boundaryObject):
         """Returns boolean indicating if a boundary with that name exists"""
-        return name in self._nameToBoundary
+        return boundaryObject in self._boundaryInfos
 
-    def setBoundary(self, function, indexExpr=None, maskArr=None, name=None):
+    def setBoundary(self, boundaryObject, indexExpr=None, maskArr=None):
         """
         Sets boundary in a rectangular region (slice)
 
-        :param function: boundary function or the string 'fluid' to remove boundary conditions
+        :param boundaryObject: boundary condition object or the string 'fluid' to remove boundary conditions
         :param indexExpr: slice expression, where boundary should be set, see :mod:`pystencils.slicing`
         :param maskArr: optional boolean (masked) array specifying where the boundary should be set
-        :param name: name of the boundary
         """
         if indexExpr is None:
             indexExpr = [slice(None, None, None)] * len(self.flagField.shape)
-        if function == 'fluid':
+        if boundaryObject == 'fluid':
             flag = self._fluidFlag
         else:
-            if name is None:
-                if hasattr(function, '__name__'):
-                    name = function.__name__
-                elif hasattr(function, 'name'):
-                    name = function.name
-                else:
-                    raise ValueError("Boundary function has to have a '__name__' or 'name' attribute "
-                                     "if name is not specified")
-
-            if not self.hasBoundary(name):
-                self.addBoundary(function, name)
-
-            flag = self.getFlag(name)
+            flag = self.addBoundary(boundaryObject)
             assert flag != self._fluidFlag
 
         indexExpr = normalizeSlice(indexExpr, self._shapeWithGhostLayers)
@@ -123,26 +103,20 @@ class BoundaryHandling(object):
             flagFieldView[maskArr] = flag
         self._dirty = True
 
-    def addBoundary(self, boundaryFunction, name=None):
+    def addBoundary(self, boundaryObject):
         """
         Adds a boundary condition, i.e. reserves a flog in the flag field and returns that flag
-        If a boundary with that name already exists, the existing flag is returned.
+        If the boundary already exists, the existing flag is returned.
         This flag can be logicalled or'ed to the boundaryHandling.flagField
 
-        :param boundaryFunction: boundary condition function, see :mod:`lbmpy.boundaries.boundaryconditions`
-        :param name: boundaries with different name are considered different. If not given
-                     ```boundaryFunction.__name__`` is used
+        :param boundaryObject: boundary condition object, see :mod:`lbmpy.boundaries.boundaryconditions`
         """
-        if name is None:
-            name = boundaryFunction.__name__
-
-        if self.hasBoundary(name):
-            return self._boundaryInfos[name].flag
+        if boundaryObject in self._boundaryInfos:
+            return self._boundaryInfos[boundaryObject].flag
 
         newIdx = len(self._boundaryInfos) + 1  # +1 because 2**0 is reserved for fluid flag
-        boundaryInfo = self.BoundaryInfo(name, 2 ** newIdx, boundaryFunction, None, None)
-        self._boundaryInfos.append(boundaryInfo)
-        self._nameToBoundary[name] = boundaryInfo
+        boundaryInfo = self.BoundaryInfo(2 ** newIdx, boundaryObject, None, None)
+        self._boundaryInfos[boundaryObject] = boundaryInfo
         self._dirty = True
         return boundaryInfo.flag
 
@@ -157,7 +131,7 @@ class BoundaryHandling(object):
         """Run the boundary handling, all keyword args are passed through to the boundary sweeps"""
         if self._dirty:
             self.prepare()
-        for boundary in self._boundaryInfos:
+        for boundary in self._boundaryInfos.values():
             boundary.kernel(**kwargs)
         for k in self._periodicityKernels:
             k(**kwargs)
@@ -166,27 +140,47 @@ class BoundaryHandling(object):
         """Removes all boundaries and fills the domain with fluid"""
         self.flagField.fill(self._fluidFlag)
         self._dirty = False
-        self._boundaryInfos = []
-        self._nameToBoundary = {}
+        self._boundaryInfos = {}
 
     def prepare(self):
         """Fills data structures to speed up the boundary handling and compiles all boundary kernels.
         This is automatically done when first called. With this function this can be triggered before running."""
-        for boundary in self._boundaryInfos:
+        for boundary in self._boundaryInfos.values():
             assert boundary.flag != self._fluidFlag
-            idxField = createBoundaryIndexList(self.flagField, self._lbMethod.stencil,
+            idxArray = createBoundaryIndexList(self.flagField, self._lbMethod.stencil,
                                                boundary.flag, self._fluidFlag, self._ghostLayers)
-            ast = generateBoundaryHandling(self._symbolicPdfField, idxField, self._lbMethod, boundary.function,
-                                           target=self._target)
+
+            dim = self._lbMethod.dim
+
+            if boundary.object.additionalData:
+                coordinateNames = ["x", "y", "z"][:dim]
+                indexArrDtype = np.dtype([(name, np.int32) for name in coordinateNames] +
+                                         [('dir', np.int32)] +
+                                         [(i[0], i[1].numpyDtype) for i in boundary.object.additionalData])
+                extendedIdxField = np.empty(len(idxArray), dtype=indexArrDtype)
+                for prop in coordinateNames + ['dir']:
+                    extendedIdxField[prop] = idxArray[prop]
+
+                idxArray = extendedIdxField
+                if boundary.object.additionalDataInit:
+                    initKernelAst = generateIndexBoundaryKernel(self._symbolicPdfField, idxArray, self._lbMethod,
+                                                                boundary.object, target='cpu',
+                                                                createInitializationKernel=True)
+                    from pystencils.cpu import makePythonFunction as makePythonCpuFunction
+                    initKernel = makePythonCpuFunction(initKernelAst, {'indexField': idxArray, 'pdfs': self._pdfField})
+                    initKernel()
+
+            ast = generateIndexBoundaryKernel(self._symbolicPdfField, idxArray, self._lbMethod, boundary.object,
+                                              target=self._target)
             boundary.ast = ast
             if self._target == 'cpu':
                 from pystencils.cpu import makePythonFunction as makePythonCpuFunction, addOpenMP
                 addOpenMP(ast, numThreads=self.openMP)
-                boundary.kernel = makePythonCpuFunction(ast, {'indexField': idxField})
+                boundary.kernel = makePythonCpuFunction(ast, {'indexField': idxArray})
             elif self._target == 'gpu':
                 from pystencils.gpucuda import makePythonFunction as makePythonGpuFunction
                 import pycuda.gpuarray as gpuarray
-                idxGpuField = gpuarray.to_gpu(idxField)
+                idxGpuField = gpuarray.to_gpu(idxArray)
                 boundary.kernel = makePythonGpuFunction(ast, {'indexField': idxGpuField})
             else:
                 assert False
@@ -276,19 +270,20 @@ class LbmMethodInfo(CustomCppCode):
         super(LbmMethodInfo, self).__init__(code, symbolsRead=set(), symbolsDefined=symbolsDefined)
 
 
-def generateBoundaryHandling(pdfField, indexArr, lbMethod, boundaryFunctor, target='cpu'):
+def generateIndexBoundaryKernel(pdfField, indexArr, lbMethod, boundaryFunctor, target='cpu',
+                                createInitializationKernel=False):
     indexField = Field.createFromNumpyArray("indexField", indexArr)
 
     elements = [LbmMethodInfo(lbMethod)]
     dirSymbol = TypedSymbol("dir", indexArr.dtype.fields['dir'][0])
     boundaryEqList = [sp.Eq(dirSymbol, indexField[0]('dir'))]
-    boundaryEqList += boundaryFunctor(pdfField, dirSymbol, lbMethod)
-    if type(boundaryEqList) is tuple:
-        boundaryEqList, additionalNodes = boundaryEqList
-        elements += boundaryEqList
-        elements += additionalNodes
+    if createInitializationKernel:
+        boundaryEqList += boundaryFunctor.additionalDataInit(pdfField=pdfField, directionSymbol=dirSymbol,
+                                                             lbMethod=lbMethod, indexField=indexField)
     else:
-        elements += boundaryEqList
+        boundaryEqList += boundaryFunctor(pdfField=pdfField, directionSymbol=dirSymbol, lbMethod=lbMethod,
+                                          indexField=indexField)
+    elements += boundaryEqList
 
     if target == 'cpu':
         from pystencils.cpu import createIndexedKernel

plot2d.py

@@ -65,9 +65,9 @@ def plotBoundaryHandling(boundaryHandling, boundaryNameToColor=None):
             'noSlip': '#000000'
         }
 
-    boundaryNames = ['fluid'] + boundaryHandling.getBoundaryNames()
-    flagValues = [boundaryHandling.fluidFlag] + [boundaryHandling.getFlag(n)
-                                                 for n in boundaryHandling.getBoundaryNames()]
+    boundaryNames = ['fluid'] + [b.name for b in boundaryHandling.getBoundaries()]
+    flagValues = [boundaryHandling.fluidFlag] + [boundaryHandling.getFlag(b)
+                                                 for b in boundaryHandling.getBoundaries()]
 
     defaultCycler = matplotlib.rcParams['axes.prop_cycle']
     colorValues = [fixedColors[name] if name in fixedColors else cycle['color']

scenarios.py

@@ -11,9 +11,9 @@ It is a good starting point if you are new to lbmpy.
 
 All scenarios can be modified, for example you can create a simple channel first, then place an object in it:
 
->>> from lbmpy.boundaries import noSlip
+>>> from lbmpy.boundaries import NoSlip
 >>> from pystencils.slicing import makeSlice
->>> scenario.boundaryHandling.setBoundary(noSlip, makeSlice[0.3:0.4, 0.0:0.3])
+>>> scenario.boundaryHandling.setBoundary(NoSlip(), makeSlice[0.3:0.4, 0.0:0.3])
 
 Functions for scenario setup:
 -----------------------------
@@ -25,12 +25,11 @@ that defines the viscosity of the fluid (valid values being between 0 and 2).
 """
 import numpy as np
 import sympy as sp
-from functools import partial
 from pystencils.field import getLayoutOfArray, createNumpyArrayWithLayout
 from pystencils.slicing import sliceFromDirection, addGhostLayers, removeGhostLayers, normalizeSlice, makeSlice
 from lbmpy.creationfunctions import createLatticeBoltzmannFunction, updateWithDefaultParameters
 from lbmpy.macroscopic_value_kernels import compileMacroscopicValuesGetter, compileMacroscopicValuesSetter
-from lbmpy.boundaries import BoundaryHandling, noSlip, ubb, fixedDensity
+from lbmpy.boundaries import BoundaryHandling, NoSlip, NoSlipFullWay, UBB, FixedDensity
 from lbmpy.stencils import getStencil
 from lbmpy.updatekernels import createPdfArray
 
@@ -78,18 +77,17 @@ def createLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParams={},
     """
     scenario = Scenario(domainSize, kwargs, optimizationParams, lbmKernel=lbmKernel, kernelParams=kernelParams)
 
-    myUbb = partial(ubb, velocity=[lidVelocity, 0, 0][:scenario.method.dim])
-    myUbb.name = 'ubb'
+    myUbb = UBB(velocity=[lidVelocity, 0, 0][:scenario.method.dim])
     dim = scenario.method.dim
     scenario.boundaryHandling.setBoundary(myUbb, sliceFromDirection('N', dim))
     for direction in ('W', 'E', 'S') if scenario.method.dim == 2 else ('W', 'E', 'S', 'T', 'B'):
-        scenario.boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, dim))
+        scenario.boundaryHandling.setBoundary(NoSlip(), sliceFromDirection(direction, dim))
 
     return scenario
 
 
 def createForceDrivenChannel(force=1e-6, domainSize=None, dim=2, radius=None, length=None, initialVelocity=None,
-                             optimizationParams={}, lbmKernel=None, kernelParams={}, **kwargs):
+                             optimizationParams={}, lbmKernel=None, kernelParams={}, halfWayBounceBack=True, **kwargs):
     """
     Creates a channel flow in x direction, which is driven by a constant force along the x axis
 
@@ -104,9 +102,12 @@ def createForceDrivenChannel(force=1e-6, domainSize=None, dim=2, radius=None, le
     :param optimizationParams: see :mod:`lbmpy.creationfunctions`
     :param lbmKernel: a LBM function, which would otherwise automatically created
     :param kernelParams: additional parameters passed to the sweep
+    :param halfWayBounceBack: determines wall boundary condition: if true half-way bounce back, if false full-way
     :param kwargs: other parameters are passed on to the method, see :mod:`lbmpy.creationfunctions`
     :return: instance of :class:`Scenario`
     """
+    wallBoundary = NoSlip() if halfWayBounceBack else NoSlipFullWay()
+
     if domainSize is not None:
         dim = len(domainSize)
     else:
@@ -137,18 +138,18 @@ def createForceDrivenChannel(force=1e-6, domainSize=None, dim=2, radius=None, le
     boundaryHandling.setPeriodicity(True, False, False)
     if dim == 2:
         for direction in ('N', 'S'):
-            boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, dim))
+            boundaryHandling.setBoundary(wallBoundary, sliceFromDirection(direction, dim))
     elif dim == 3:
         if roundChannel:
-            noSlipIdx = boundaryHandling.addBoundary(noSlip)
+            wallIdx = boundaryHandling.addBoundary(wallBoundary)
             ff = boundaryHandling.flagField
             yMid = ff.shape[1] // 2
             zMid = ff.shape[2] // 2
             y, z = np.meshgrid(range(ff.shape[1]), range(ff.shape[2]))
-            ff[(y - yMid) ** 2 + (z - zMid) ** 2 > radius ** 2] = noSlipIdx
+            ff[(y - yMid) ** 2 + (z - zMid) ** 2 > radius ** 2] = wallIdx
         else:
             for direction in ('N', 'S', 'T', 'B'):
-                boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, dim))
+                boundaryHandling.setBoundary(wallBoundary, sliceFromDirection(direction, dim))
 
     assert domainSize is not None
     if 'forceModel' not in kwargs:
@@ -159,7 +160,7 @@ def createForceDrivenChannel(force=1e-6, domainSize=None, dim=2, radius=None, le
 
 def createPressureGradientDrivenChannel(pressureDifference, domainSize=None, dim=2, radius=None, length=None,
                                         initialVelocity=None, optimizationParams={},
-                                        lbmKernel=None, kernelParams={}, **kwargs):
+                                        lbmKernel=None, kernelParams={}, halfWayBounceBack=True, **kwargs):
     """
     Creates a channel flow in x direction, which is driven by two pressure boundaries.
     Consider using :func:`createForceDrivenChannel` which does not have artifacts an inflow and outflow.
@@ -175,9 +176,12 @@ def createPressureGradientDrivenChannel(pressureDifference, domainSize=None, dim
     :param optimizationParams: see :mod:`lbmpy.creationfunctions`
     :param lbmKernel: a LBM function, which would otherwise automatically created
     :param kernelParams: additional parameters passed to the sweep
+    :param halfWayBounceBack: determines wall boundary condition: if true half-way bounce back, if false full-way
     :param kwargs: other parameters are passed on to the method, see :mod:`lbmpy.creationfunctions`
     :return: instance of :class:`Scenario`
     """
+    wallBoundary = NoSlip() if halfWayBounceBack else NoSlipFullWay()
+
     if domainSize is not None:
         dim = len(domainSize)
     else:
@@ -202,20 +206,18 @@ def createPressureGradientDrivenChannel(pressureDifference, domainSize=None, dim
     scenario = Scenario(domainSize, kwargs, optimizationParams, lbmKernel=lbmKernel,
                         initialVelocity=initialVelocity, kernelParams=kernelParams)
     boundaryHandling = scenario.boundaryHandling
-    pressureBoundaryInflow = partial(fixedDensity, density=1.0 + pressureDifference)
-    pressureBoundaryInflow.__name__ = "Inflow"
+    pressureBoundaryInflow = FixedDensity(density=1.0 + pressureDifference)
 
-    pressureBoundaryOutflow = partial(fixedDensity, density=1.0)
-    pressureBoundaryOutflow.__name__ = "Outflow"
+    pressureBoundaryOutflow = FixedDensity(density=1.0)
     boundaryHandling.setBoundary(pressureBoundaryInflow, sliceFromDirection('W', dim))
     boundaryHandling.setBoundary(pressureBoundaryOutflow, sliceFromDirection('E', dim))
 
     if dim == 2:
         for direction in ('N', 'S'):
-            boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, dim))
+            boundaryHandling.setBoundary(wallBoundary, sliceFromDirection(direction, dim))
     elif dim == 3:
         if roundChannel:
-            noSlipIdx = boundaryHandling.addBoundary(noSlip)
+            noSlipIdx = boundaryHandling.addBoundary(wallBoundary)
             ff = boundaryHandling.flagField
             yMid = ff.shape[1] // 2
             zMid = ff.shape[2] // 2
@@ -223,7 +225,7 @@ def createPressureGradientDrivenChannel(pressureDifference, domainSize=None, dim
             ff[(y - yMid) ** 2 + (z - zMid) ** 2 > radius ** 2] = noSlipIdx
         else:
             for direction in ('N', 'S', 'T', 'B'):
-                boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, dim))
+                boundaryHandling.setBoundary(wallBoundary, sliceFromDirection(direction, dim))
 
     return scenario