diff --git a/boundaries/__init__.py b/boundaries/__init__.py
index 327e1bb2ce0cb7690f8e51b1a8ed811a1f3c3e39..925579cb7c5f3c38f37a8581f318f3c1e7226b99 100644
--- a/boundaries/__init__.py
+++ b/boundaries/__init__.py
@@ -1,2 +1,2 @@
 from lbmpy.boundaries.boundaryconditions import NoSlip, NoSlipFullWay, UBB, FixedDensity
-from lbmpy.boundaries.boundaryhandling import BoundaryHandling
+from lbmpy.boundaries.boundaryhandling import BoundaryHandling
\ No newline at end of file
diff --git a/boundaries/boundary_kernel.py b/boundaries/boundary_kernel.py
new file mode 100644
index 0000000000000000000000000000000000000000..2c17bf0fc33cd73a3824a5f61dc3e6b429e9df4b
--- /dev/null
+++ b/boundaries/boundary_kernel.py
@@ -0,0 +1,73 @@
+import sympy as sp
+from pystencils.backends.cbackend import CustomCppCode
+from pystencils import TypedSymbol, Field
+
+INV_DIR_SYMBOL = TypedSymbol("invDir", "int")
+WEIGHTS_SYMBOL = TypedSymbol("weights", "double")
+
+
+# -------------------------------------- Helper Functions --------------------------------------------------------------
+
+
+def offsetSymbols(dim):
+    return [TypedSymbol("c_%d" % (d,), "int") for d in range(dim)]
+
+
+def offsetFromDir(dirIdx, dim):
+    return tuple([sp.IndexedBase(symbol, shape=(1,))[dirIdx] for symbol in offsetSymbols(dim)])
+
+
+def invDir(dirIdx):
+    return sp.IndexedBase(INV_DIR_SYMBOL, shape=(1,))[dirIdx]
+
+
+def weightOfDirection(dirIdx):
+    return sp.IndexedBase(WEIGHTS_SYMBOL, shape=(1,))[dirIdx]
+
+
+# ------------------------------------- Kernel Generation --------------------------------------------------------------
+
+class LbmMethodInfo(CustomCppCode):
+    def __init__(self, lbMethod):
+        stencil = lbMethod.stencil
+        symbolsDefined = set(offsetSymbols(lbMethod.dim) + [INV_DIR_SYMBOL, WEIGHTS_SYMBOL])
+
+        offsetSym = offsetSymbols(lbMethod.dim)
+        code = "\n"
+        for i in range(lbMethod.dim):
+            offsetStr = ", ".join([str(d[i]) for d in stencil])
+            code += "const int %s [] = { %s };\n" % (offsetSym[i].name, offsetStr)
+
+        invDirs = []
+        for direction in stencil:
+            inverseDir = tuple([-i for i in direction])
+            invDirs.append(str(stencil.index(inverseDir)))
+
+        code += "const int %s [] = { %s };\n" % (INV_DIR_SYMBOL.name, ", ".join(invDirs))
+        weights = [str(w.evalf()) for w in lbMethod.weights]
+        code += "const double %s [] = { %s };\n" % (WEIGHTS_SYMBOL.name, ",".join(weights))
+        super(LbmMethodInfo, self).__init__(code, symbolsRead=set(), symbolsDefined=symbolsDefined)
+
+
+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'))]
+    if createInitializationKernel:
+        boundaryEqList += boundaryFunctor.additionalDataInit(pdfField=pdfField, directionSymbol=dirSymbol,
+                                                             lbMethod=lbMethod, indexField=indexField)
+    else:
+        boundaryEqList += boundaryFunctor(pdfField=pdfField, directionSymbol=dirSymbol, lbMethod=lbMethod,
+                                          indexField=indexField)
+    elements += boundaryEqList
+
+    if target == 'cpu':
+        from pystencils.cpu import createIndexedKernel
+        return createIndexedKernel(elements, [indexField])
+    elif target == 'gpu':
+        from pystencils.gpucuda import createdIndexedCUDAKernel
+        return createdIndexedCUDAKernel(elements, [indexField])
+
diff --git a/boundaries/boundaryconditions.py b/boundaries/boundaryconditions.py
index 1fd31309a40af2581153e71cc1f15023d89ccfbc..245f988c2c32c479e0ad4b40c61116e6f6987c82 100644
--- a/boundaries/boundaryconditions.py
+++ b/boundaries/boundaryconditions.py
@@ -3,7 +3,7 @@ import sympy as sp
 from lbmpy.simplificationfactory import createSimplificationStrategy
 from pystencils.sympyextensions import getSymmetricPart
 from pystencils import Field
-from lbmpy.boundaries.boundaryhandling import offsetFromDir, weightOfDirection, invDir
+from lbmpy.boundaries.boundary_kernel import offsetFromDir, weightOfDirection, invDir
 from pystencils.types import createTypeFromString
 
 
diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py
index 5863b4b39aa1451608c7bb08629dc1807d5812cd..92269a161390acd8c2bef9c694b95d1504fe5e1f 100644
--- a/boundaries/boundaryhandling.py
+++ b/boundaries/boundaryhandling.py
@@ -1,294 +1,67 @@
-import sympy as sp
 import numpy as np
+from lbmpy.boundaries.handlinginterface import GenericBoundaryHandling, FlagFieldInterface
+from lbmpy.boundaries.periodicityhandling import PeriodicityHandling
 
-from lbmpy.stencils import getStencil
-from pystencils import TypedSymbol, Field
-from pystencils.backends.cbackend import CustomCppCode
-from lbmpy.boundaries.createindexlist import createBoundaryIndexList
-from pystencils.slicing import normalizeSlice, makeSlice
 
-INV_DIR_SYMBOL = TypedSymbol("invDir", "int")
-WEIGHTS_SYMBOL = TypedSymbol("weights", "double")
+class BoundaryHandling(GenericBoundaryHandling, PeriodicityHandling):
 
+    def __init__(self, pdfField, domainShape, lbMethod, ghostLayers=1, target='cpu', openMP=True, flagDtype=np.uint32):
+        shapeWithGl = [a + 2 * ghostLayers for a in domainShape]
+        flagFieldInterface = NumpyFlagFieldInterface(shapeWithGl, flagDtype)
 
-class BoundaryHandling(object):
-    class BoundaryInfo(object):
-        def __init__(self, flag, object, kernel, ast):
-            self.flag = flag
-            self.object = object
-            self.kernel = kernel
-            self.ast = ast
+        GenericBoundaryHandling.__init__(self, flagFieldInterface, pdfField, lbMethod, ghostLayers, target, openMP)
+        PeriodicityHandling.__init__(self, domainShape)
 
-    def __init__(self, pdfField, domainShape, lbMethod, ghostLayers=1, target='cpu', openMP=True):
-        """
-        Class for managing boundary kernels
+    def setBoundary(self, boundaryObject, indexExpr=None, maskCallback=None, sliceNormalizationGhostLayers=1):
+        mask = None
+        if maskCallback is not None:
+            meshgridParams = [np.arange(start=-sliceNormalizationGhostLayers, stop=s-sliceNormalizationGhostLayers)
+                              for s in self.flagField.shape]
+            indices = np.meshgrid(*meshgridParams, indexing='ij')
+            mask = maskCallback(*indices)
+        return GenericBoundaryHandling.setBoundary(self, boundaryObject, indexExpr, mask)
 
-        :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'
-        """
-        symbolicPdfField = Field.createFromNumpyArray('pdfs', pdfField, indexDimensions=1)
-        assert pdfField.shape[:-1] == tuple(d + 2*ghostLayers for d in domainShape)
+    def __call__(self, *args, **kwargs):
+        for cls in BoundaryHandling.__bases__:
+            cls.__call__(self, *args, **kwargs)
 
-        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
+    def prepare(self):
+        for cls in BoundaryHandling.__bases__:
+            cls.prepare(self)
 
-        self._boundaryInfos = {}
-        self._periodicityKernels = []
+# ----------------------------------------------------------------------------------------------------------------------
 
-        self._dirty = False
-        self._periodicity = [False, False, False]
-        self._target = target
-        self.openMP = openMP
-        if target not in ('cpu', 'gpu'):
-            raise ValueError("Invalid target '%s' . Allowed values: 'cpu' or 'gpu'" % (target,))
 
-    @property
-    def periodicity(self):
-        """List that indicates for x,y (z) coordinate if domain is periodic in that direction"""
-        return self._periodicity
+class NumpyFlagFieldInterface(FlagFieldInterface):
+
+    def __init__(self, shape, dtype=np.uint32):
+        self._array = np.ones(shape, dtype=dtype)
+        self._boundaryObjectToName = {}
+        self._boundaryObjectToFlag = {'fluid': dtype(2**0)}
+        self._nextFreeExponent = 1
 
     @property
-    def fluidFlag(self):
-        """Flag that is set where the lattice Boltzmann update should happen"""
-        return self._fluidFlag
+    def array(self):
+        return self._array
 
     def getFlag(self, boundaryObject):
-        """Flag that represents the given boundary."""
-        return self._boundaryInfos[boundaryObject].flag
-
-    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"""
-        for d in (x, y, z):
-            assert isinstance(d, bool)
-
-        self._periodicity = [x, y, z]
-        self._compilePeriodicityKernels()
-
-    def hasBoundary(self, boundaryObject):
-        """Returns boolean indicating if a boundary with that name exists"""
-        return boundaryObject in self._boundaryInfos
-
-    def setBoundary(self, boundaryObject, indexExpr=None, maskArr=None):
-        """
-        Sets boundary in a rectangular region (slice)
-
-        :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
-        """
-        if indexExpr is None:
-            indexExpr = [slice(None, None, None)] * len(self.flagField.shape)
-        if boundaryObject == 'fluid':
-            flag = self._fluidFlag
-        else:
-            flag = self.addBoundary(boundaryObject)
-            assert flag != self._fluidFlag
-
-        indexExpr = normalizeSlice(indexExpr, self._shapeWithGhostLayers)
-        if maskArr is None:
-            self.flagField[indexExpr] = flag
-        else:
-            flagFieldView = self.flagField[indexExpr]
-            flagFieldView[maskArr] = flag
-        self._dirty = True
-
-    def addBoundary(self, boundaryObject):
-        """
-        Adds a boundary condition, i.e. reserves a flog in the flag field and returns that flag
-        If the boundary already exists, the existing flag is returned.
-        This flag can be logicalled or'ed to the boundaryHandling.flagField
+        if boundaryObject not in self._boundaryObjectToFlag:
+            self._boundaryObjectToFlag[boundaryObject] = 2 ** self._nextFreeExponent
+            self._nextFreeExponent += 1
+            name = self._makeBoundaryName(boundaryObject, self._boundaryObjectToName.values())
+            self._boundaryObjectToName[boundaryObject] = name
 
-        :param boundaryObject: boundary condition object, see :mod:`lbmpy.boundaries.boundaryconditions`
-        """
-        if boundaryObject in self._boundaryInfos:
-            return self._boundaryInfos[boundaryObject].flag
+        return self._boundaryObjectToFlag[boundaryObject]
 
-        newIdx = len(self._boundaryInfos) + 1  # +1 because 2**0 is reserved for fluid flag
-        boundaryInfo = self.BoundaryInfo(2 ** newIdx, boundaryObject, None, None)
-        self._boundaryInfos[boundaryObject] = boundaryInfo
-        self._dirty = True
-        return boundaryInfo.flag
+    def getName(self, boundaryObject):
+        return self._boundaryObjectToName[boundaryObject]
 
-    def indices(self, dx=1.0, includeGhostLayers=False):
-        if not includeGhostLayers:
-            params = [np.arange(start=-1, stop=s-1) * dx for s in self.flagField.shape]
-        else:
-            params = [np.arange(s) * dx for s in self.flagField.shape]
-        return np.meshgrid(*params, indexing='ij')
-
-    def __call__(self, **kwargs):
-        """Run the boundary handling, all keyword args are passed through to the boundary sweeps"""
-        if self._dirty:
-            self.prepare()
-        for boundary in self._boundaryInfos.values():
-            boundary.kernel(**kwargs)
-        for k in self._periodicityKernels:
-            k(**kwargs)
+    @property
+    def boundaryObjects(self):
+        return self._boundaryObjectToName.keys()
 
     def clear(self):
-        """Removes all boundaries and fills the domain with fluid"""
-        self.flagField.fill(self._fluidFlag)
-        self._dirty = False
-        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.values():
-            assert boundary.flag != self._fluidFlag
-            idxArray = createBoundaryIndexList(self.flagField, self._lbMethod.stencil,
-                                               boundary.flag, self._fluidFlag, self._ghostLayers)
-
-            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': idxArray})
-            elif self._target == 'gpu':
-                from pystencils.gpucuda import makePythonFunction as makePythonGpuFunction
-                import pycuda.gpuarray as gpuarray
-                idxGpuField = gpuarray.to_gpu(idxArray)
-                boundary.kernel = makePythonGpuFunction(ast, {'indexField': idxGpuField})
-            else:
-                assert False
-        self._dirty = False
-
-    def invalidateIndexCache(self):
-        """Invalidates the cache for optimization data structures. When setting boundaries the cache is automatically
-        invalidated, so there is no need to call this function manually, as long as the flag field is not manually
-        modified."""
-        self._dirty = True
-
-    def _compilePeriodicityKernels(self):
-        self._periodicityKernels = []
-        dim = len(self.flagField.shape)
-        if dim == 2:
-            stencil = getStencil("D2Q9")
-        elif dim == 3:
-            stencil = getStencil("D3Q27")
-        else:
-            assert False
-
-        filteredStencil = []
-        for direction in stencil:
-            useDirection = True
-            if direction == (0, 0) or direction == (0, 0, 0):
-                useDirection = False
-            for component, periodicity in zip(direction, self._periodicity):
-                if not periodicity and component != 0:
-                    useDirection = False
-            if useDirection:
-                filteredStencil.append(direction)
-
-        if len(filteredStencil) > 0:
-            if self._target == 'cpu':
-                from pystencils.slicing import getPeriodicBoundaryFunctor
-                self._periodicityKernels.append(getPeriodicBoundaryFunctor(filteredStencil, ghostLayers=1))
-            elif self._target == 'gpu':
-                from pystencils.gpucuda.periodicity import getPeriodicBoundaryFunctor
-                self._periodicityKernels.append(getPeriodicBoundaryFunctor(filteredStencil, self.flagField.shape,
-                                                                           indexDimensions=1,
-                                                                           indexDimShape=len(self._lbMethod.stencil),
-                                                                           ghostLayers=1))
-            else:
-                assert False
-
-
-# -------------------------------------- Helper Functions --------------------------------------------------------------
-
-
-def offsetSymbols(dim):
-    return [TypedSymbol("c_%d" % (d,), "int") for d in range(dim)]
-
-
-def offsetFromDir(dirIdx, dim):
-    return tuple([sp.IndexedBase(symbol, shape=(1,))[dirIdx] for symbol in offsetSymbols(dim)])
-
-
-def invDir(dirIdx):
-    return sp.IndexedBase(INV_DIR_SYMBOL, shape=(1,))[dirIdx]
-
-
-def weightOfDirection(dirIdx):
-    return sp.IndexedBase(WEIGHTS_SYMBOL, shape=(1,))[dirIdx]
-
-
-# ------------------------------------- Kernel Generation --------------------------------------------------------------
-
-class LbmMethodInfo(CustomCppCode):
-    def __init__(self, lbMethod):
-        stencil = lbMethod.stencil
-        symbolsDefined = set(offsetSymbols(lbMethod.dim) + [INV_DIR_SYMBOL, WEIGHTS_SYMBOL])
-
-        offsetSym = offsetSymbols(lbMethod.dim)
-        code = "\n"
-        for i in range(lbMethod.dim):
-            offsetStr = ", ".join([str(d[i]) for d in stencil])
-            code += "const int %s [] = { %s };\n" % (offsetSym[i].name, offsetStr)
-
-        invDirs = []
-        for direction in stencil:
-            inverseDir = tuple([-i for i in direction])
-            invDirs.append(str(stencil.index(inverseDir)))
-
-        code += "const int %s [] = { %s };\n" % (INV_DIR_SYMBOL.name, ", ".join(invDirs))
-        weights = [str(w.evalf()) for w in lbMethod.weights]
-        code += "const double %s [] = { %s };\n" % (WEIGHTS_SYMBOL.name, ",".join(weights))
-        super(LbmMethodInfo, self).__init__(code, symbolsRead=set(), symbolsDefined=symbolsDefined)
-
-
-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'))]
-    if createInitializationKernel:
-        boundaryEqList += boundaryFunctor.additionalDataInit(pdfField=pdfField, directionSymbol=dirSymbol,
-                                                             lbMethod=lbMethod, indexField=indexField)
-    else:
-        boundaryEqList += boundaryFunctor(pdfField=pdfField, directionSymbol=dirSymbol, lbMethod=lbMethod,
-                                          indexField=indexField)
-    elements += boundaryEqList
-
-    if target == 'cpu':
-        from pystencils.cpu import createIndexedKernel
-        return createIndexedKernel(elements, [indexField])
-    elif target == 'gpu':
-        from pystencils.gpucuda import createdIndexedCUDAKernel
-        return createdIndexedCUDAKernel(elements, [indexField])
-
+        self._array.fill(0)
+        self._boundaryObjectToName = {}
+        self._boundaryObjectToFlag = {'fluid': np.dtype(self._array.dtype)(2**0)}
+        self._nextFreeExponent = 1
diff --git a/boundaries/handlinginterface.py b/boundaries/handlinginterface.py
new file mode 100644
index 0000000000000000000000000000000000000000..7de44c7239a04eee6d57d28d71fa59466599700e
--- /dev/null
+++ b/boundaries/handlinginterface.py
@@ -0,0 +1,180 @@
+import numpy as np
+
+from pystencils.field import Field
+from pystencils.slicing import normalizeSlice
+from lbmpy.boundaries.boundary_kernel import generateIndexBoundaryKernel
+from lbmpy.boundaries.createindexlist import createBoundaryIndexList
+
+
+class FlagFieldInterface(object):
+
+    def getFlag(self, boundaryObject):
+        raise NotImplementedError()
+
+    def getName(self, boundaryObject):
+        raise NotImplementedError()
+
+    @property
+    def array(self):
+        raise NotImplementedError()
+
+    @property
+    def boundaryObjects(self):
+        raise NotImplementedError()
+
+    def clear(self):
+        raise NotImplementedError()
+
+    @staticmethod
+    def _makeBoundaryName(boundaryObject, existingNames):
+        baseName = boundaryObject.name
+        name = baseName
+        counter = 1
+        while name in existingNames:
+            name = "%s_%d" % (baseName, counter)
+        return name
+
+
+class GenericBoundaryHandling(object):
+
+    class BoundaryInfo(object):
+        def __init__(self, kernel, ast, indexArray):
+            self.kernel = kernel
+            self.ast = ast
+            self.indexArray = indexArray
+
+    def __init__(self, flagFieldInterface, pdfField, lbMethod, ghostLayers=1, target='cpu', openMP=True):
+        """
+        :param flagFieldInterface: implementation of FlagFieldInterface
+        :param pdfField: numpy array
+        :param lbMethod:
+        :param target: 'cpu' or 'gpu'
+        :param openMP:
+        """
+        self._flagFieldInterface = flagFieldInterface
+        self._pdfField = pdfField
+        self._lbMethod = lbMethod
+        self._target = target
+        self._openMP = openMP
+        self.ghostLayers = ghostLayers
+        self._dirty = False
+
+        self._boundaryInfos = {}  # mapping of boundary object to boundary info
+        self._symbolicPdfField = Field.createFromNumpyArray('pdfs', pdfField, indexDimensions=1)
+
+        if target not in ('cpu', 'gpu'):
+            raise ValueError("Invalid target '%s' . Allowed values: 'cpu' or 'gpu'" % (target,))
+
+    def getBoundaryNameToFlagDict(self):
+        result = {self._flagFieldInterface.getName(o): self._flagFieldInterface.getFlag(o) for o in self._boundaryInfos}
+        result['fluid'] = self._flagFieldInterface.getFlag('fluid')
+        return result
+
+    def hasBoundary(self, boundaryObject):
+        """Returns boolean indicating if this boundary exists in that handling"""
+        return boundaryObject in self._boundaryInfos
+
+    def __call__(self, **kwargs):
+        """
+        Runs boundary handling
+        :param kwargs: keyword arguments passed to boundary kernel
+        :return:
+        """
+        if self._dirty:
+            self.prepare()
+        for boundary in self._boundaryInfos.values():
+            boundary.kernel(**kwargs)
+
+    def getBoundaryIndexArray(self, boundaryObject):
+        return self._boundaryInfos[boundaryObject].indexArray
+
+    def clear(self):
+        """Removes all boundaries and fills the domain with fluid"""
+        self._flagFieldInterface.clear()
+        self._dirty = False
+        self._boundaryInfos = {}
+
+    def _invalidateCache(self):
+        self._dirty = True
+
+    def prepare(self):
+        """Compiles boundary kernels according to flag field. When setting boundaries the cache is automatically
+        invalidated, so there is no need to call this function manually, as long as the flag field is not manually
+        modified."""
+        if not self._dirty:
+            return
+
+        dim = self._lbMethod.dim
+        fluidFlag = self._flagFieldInterface.getFlag('fluid')
+
+        for boundaryObject in self._flagFieldInterface.boundaryObjects:
+            boundaryFlag = self._flagFieldInterface.getFlag(boundaryObject)
+            idxArray = createBoundaryIndexList(self.flagField, self._lbMethod.stencil,
+                                               boundaryFlag, fluidFlag, self.ghostLayers)
+
+            if boundaryObject.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 boundaryObject.additionalData])
+                extendedIdxField = np.empty(len(idxArray), dtype=indexArrDtype)
+                for prop in coordinateNames + ['dir']:
+                    extendedIdxField[prop] = idxArray[prop]
+
+                idxArray = extendedIdxField
+                if boundaryObject.additionalDataInit:
+                    initKernelAst = generateIndexBoundaryKernel(self._symbolicPdfField, idxArray, self._lbMethod,
+                                                                boundaryObject, 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, boundaryObject,
+                                              target=self._target)
+            if self._target == 'cpu':
+                from pystencils.cpu import makePythonFunction as makePythonCpuFunction, addOpenMP
+                addOpenMP(ast, numThreads=self._openMP)
+                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(idxArray)
+                kernel = makePythonGpuFunction(ast, {'indexField': idxGpuField})
+            else:
+                assert False
+
+            boundaryInfo = GenericBoundaryHandling.BoundaryInfo(kernel, ast, idxArray)
+            self._boundaryInfos[boundaryObject] = boundaryInfo
+
+        self._dirty = False
+
+    def setBoundary(self, boundaryObject, indexExpr=None, maskArr=None):
+        """
+        Sets boundary in a rectangular region (slice)
+
+        :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
+        """
+        if indexExpr is None:
+            indexExpr = [slice(None, None, None)] * len(self.flagField.shape)
+
+        flag = self._flagFieldInterface.getFlag(boundaryObject)
+        flagField = self._flagFieldInterface.array
+        indexExpr = normalizeSlice(indexExpr, flagField.shape)
+
+        if maskArr is None:
+            flagField[indexExpr] = flag
+        else:
+            flagFieldView = flagField[indexExpr].squeeze()
+            maskArr = maskArr.squeeze()
+            flagFieldView[maskArr] = flag
+        self._dirty = True
+
+    def getMask(self, boundaryObject):
+        return np.logical_and(self._flagFieldInterface.array, self._flagFieldInterface.getFlag(boundaryObject))
+
+    @property
+    def flagField(self):
+        return self._flagFieldInterface.array
diff --git a/boundaries/periodicityhandling.py b/boundaries/periodicityhandling.py
new file mode 100644
index 0000000000000000000000000000000000000000..6c1fcb21bb9ee80e86fefbb7b58f415e76cc9f0c
--- /dev/null
+++ b/boundaries/periodicityhandling.py
@@ -0,0 +1,64 @@
+from lbmpy.stencils import getStencil
+
+
+class PeriodicityHandling(object):
+    def __init__(self, fieldShape, periodicity=(False, False, False)):
+        self._spatialShape = fieldShape[:-1]
+        self._indexShape = fieldShape[-1]
+        self._periodicity = list(periodicity)
+        self._dirty = False
+        self._periodicityKernels = []
+
+    @property
+    def periodicity(self):
+        """List that indicates for x,y (z) coordinate if domain is periodic in that direction"""
+        return self._periodicity
+
+    def setPeriodicity(self, x=False, y=False, z=False):
+        """Enable periodic boundary conditions at the border of the domain"""
+        for d in (x, y, z):
+            assert isinstance(d, bool)
+
+        self._periodicity = [x, y, z]
+        self._dirty = True
+
+    def __call__(self, **kwargs):
+        for k in self._periodicityKernels:
+            k(**kwargs)
+
+    def prepare(self):
+        self._periodicityKernels = []
+        dim = len(self.flagField.shape)
+        if dim == 2:
+            stencil = getStencil("D2Q9")
+        elif dim == 3:
+            stencil = getStencil("D3Q27")
+        else:
+            assert False
+
+        filteredStencil = []
+        for direction in stencil:
+            useDirection = True
+            if direction == (0, 0) or direction == (0, 0, 0):
+                useDirection = False
+            for component, periodicity in zip(direction, self._periodicity):
+                if not periodicity and component != 0:
+                    useDirection = False
+            if useDirection:
+                filteredStencil.append(direction)
+
+        if len(filteredStencil) > 0:
+            if self._target == 'cpu':
+                from pystencils.slicing import getPeriodicBoundaryFunctor
+                self._periodicityKernels.append(getPeriodicBoundaryFunctor(filteredStencil, ghostLayers=1))
+            elif self._target == 'gpu':
+                from pystencils.gpucuda.periodicity import getPeriodicBoundaryFunctor
+                self._periodicityKernels.append(getPeriodicBoundaryFunctor(filteredStencil, self._spatialShape,
+                                                                           indexDimensions=1,
+                                                                           indexDimShape=self._indexShape,
+                                                                           ghostLayers=1))
+            else:
+                assert False
+
+        self._dirty = False
+
diff --git a/chapman_enskog/chapman_enskog.py b/chapman_enskog/chapman_enskog.py
index 60cc31e7e757bc0efe5a00ce5106de544e860ca5..c3e0ed05aebe160d08be578b8be484e9eaadff49 100644
--- a/chapman_enskog/chapman_enskog.py
+++ b/chapman_enskog/chapman_enskog.py
@@ -100,6 +100,9 @@ class CeMoment(sp.Symbol):
             result += "}"
         return result
 
+    def __repr__(self):
+        return "%s_(%d)_%s" % (self.name, self.ceIdx, self.momentTuple)
+
 
 class LbMethodEqMoments:
     def __init__(self, lbMethod):
@@ -113,12 +116,14 @@ class LbMethodEqMoments:
         return self._momentCache[moment]
 
 
-def insertMoments(eqn, lbMethodMoments, useSolvabilityConditions=True):
+def insertMoments(eqn, lbMethodMoments, momentName="\\Pi", useSolvabilityConditions=True):
     subsDict = {}
     if useSolvabilityConditions:
-        subsDict.update({m: 0 for m in eqn.atoms(CeMoment) if m.ceIdx > 0 and sum(m.momentTuple) <= 1})
+        condition = lambda m:  m.ceIdx > 0 and sum(m.momentTuple) <= 1 and m.name == momentName
+        subsDict.update({m: 0 for m in eqn.atoms(CeMoment) if condition(m)})
 
-    subsDict.update({m: lbMethodMoments(m.momentTuple) for m in eqn.atoms(CeMoment) if m.ceIdx == 0})
+    condition = lambda m:  m.ceIdx == 0 and m.name == momentName
+    subsDict.update({m: lbMethodMoments(m.momentTuple) for m in eqn.atoms(CeMoment) if condition(m)})
     return eqn.subs(subsDict)
 
 
@@ -151,7 +156,7 @@ def substituteCollisionOperatorMoments(expr, lbMethod, collisionOpMomentName='\\
     return expr.subs(subsDict)
 
 
-def takeMoments(eqn, pdfToMomentName=(('f', '\Pi'), ('Cf', '\\Upsilon')), velocityName='c', maxExpansion=5):
+def takeMoments(eqn, pdfToMomentName=(('f', '\Pi'), ('\Omega f', '\\Upsilon')), velocityName='c', maxExpansion=5):
 
     pdfSymbols = [tuple(expandedSymbol(name, superscript=i) for i in range(maxExpansion))
                   for name, _ in pdfToMomentName]
@@ -285,7 +290,7 @@ def removeErrorTerms(expr):
 # ----------------------------------------------------------------------------------------------------------------------
 
 
-def getTaylorExpandedLbEquation(pdfSymbolName="f", pdfsAfterCollisionOperator="Cf", velocityName="c",
+def getTaylorExpandedLbEquation(pdfSymbolName="f", pdfsAfterCollisionOperator="\Omega f", velocityName="c",
                                 dim=3, taylorOrder=2):
     dimLabels = [sp.Rational(i, 1) for i in range(dim)]
 
@@ -315,7 +320,7 @@ def getTaylorExpandedLbEquation(pdfSymbolName="f", pdfsAfterCollisionOperator="C
 
 
 def useChapmanEnskogAnsatz(equation, timeDerivativeOrders=(1, 3), spatialDerivativeOrders=(1, 2),
-                           pdfs=(['f', 0, 3], ['Cf', 1, 3])):
+                           pdfs=(['f', 0, 3], ['\Omega f', 1, 3])):
 
     t, eps = sp.symbols("t epsilon")
 
@@ -442,6 +447,12 @@ class ChapmanEnskogAnalysis(object):
         momentsUntilOrder1 = [1] + list(c)
         momentsOrder2 = [c_i * c_j for c_i, c_j in productSymmetric(c, c)]
 
+        symbolicRelaxationRates = [rr for rr in method.relaxationRates if isinstance(rr, sp.Symbol)]
+        if constants is None:
+            constants = set(symbolicRelaxationRates)
+        else:
+            constants.update(symbolicRelaxationRates)
+
         oEpsMoments1 = [expandUsingLinearity(self._takeAndInsertMoments(self.equationsGroupedByOrder[1] * moment),
                                              constants=constants)
                         for moment in momentsUntilOrder1]
@@ -570,3 +581,7 @@ class ChapmanEnskogAnalysis(object):
         solveRes = sp.solve(kinematicViscosity - nu, kinematicViscosity.atoms(sp.Symbol), dict=True)
         return solveRes[0]
 
+if __name__ == '__main__':
+    from lbmpy.creationfunctions import createLatticeBoltzmannMethod
+    m = createLatticeBoltzmannMethod(stencil='D2Q9')
+    ce = ChapmanEnskogAnalysis(m)
\ No newline at end of file
diff --git a/creationfunctions.py b/creationfunctions.py
index 1407199d4d2464c95e7ceb34fac63715fbe503c5..3f0397834d3871c1cc0ae89a1d635d678efd3dbc 100644
--- a/creationfunctions.py
+++ b/creationfunctions.py
@@ -215,6 +215,25 @@ def updateWithDefaultParameters(params, optParams, failOnUnknownParameter=True):
     return paramsResult, optParamsResult
 
 
+def switchToSymbolicRelaxationRatesForEntropicMethods(methodParameters, kernelParams):
+    """
+    For entropic kernels the relaxation rate has to be a variable. If a constant was passed a
+    new dummy variable is inserted and the value of this variable is later on passed to the kernel
+    """
+    if methodParameters['entropic']:
+        newRelaxationRates = []
+        for rr in methodParameters['relaxationRates']:
+            if not isinstance(rr, sp.Symbol):
+                dummyVar = sp.Dummy()
+                newRelaxationRates.append(dummyVar)
+                kernelParams[dummyVar.name] = rr
+            else:
+                newRelaxationRates.append(rr)
+        if len(newRelaxationRates) < 2:
+            newRelaxationRates.append(sp.Dummy())
+        methodParameters['relaxationRates'] = newRelaxationRates
+
+
 def createLatticeBoltzmannFunction(ast=None, optimizationParams={}, **kwargs):
     params, optParams = updateWithDefaultParameters(kwargs, optimizationParams)
 
diff --git a/parallel/__init__.py b/parallel/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/parallel/boundaryhandling.py b/parallel/boundaryhandling.py
new file mode 100644
index 0000000000000000000000000000000000000000..1ae128a21aaeb6d5e320b73e35f86e3adac6ac95
--- /dev/null
+++ b/parallel/boundaryhandling.py
@@ -0,0 +1,114 @@
+import numpy as np
+import waLBerla as wlb
+from lbmpy.boundaries.handlinginterface import GenericBoundaryHandling, FlagFieldInterface
+from pystencils.slicing import normalizeSlice
+
+
+class BoundaryHandling(object):
+
+    def __init__(self, blocks, lbMethod, pdfFieldId, flagFieldId, boundaryId='boundary', target='cpu', openMP=True):
+
+        def addBoundaryHandling(block, *args, **kwargs):
+            flagFieldInterface = WalberlaFlagFieldInterface(block, flagFieldId)
+            pdfField = wlb.field.toArray(block[pdfFieldId], withGhostLayers=True)
+            ghostLayers = block[pdfFieldId].nrOfGhostLayers
+            return GenericBoundaryHandling(flagFieldInterface, pdfField, lbMethod, ghostLayers=ghostLayers,
+                                           target=target, openMP=openMP)
+
+        self._boundaryId = boundaryId
+        self._blocks = blocks
+        self.dim = lbMethod.dim
+        blocks.addBlockData(boundaryId, addBoundaryHandling)
+
+    def __call__(self, **kwargs):
+        for block in self._blocks:
+            block[self._boundaryId](**kwargs)
+
+    def prepare(self):
+        for block in self._blocks:
+            block[self._boundaryId].prepare()
+
+    def setBoundary(self, boundaryObject, indexExpr=None, maskCallback=None, sliceNormalizationGhostLayers=1):
+        if indexExpr is None:
+            indexExpr = [slice(None, None, None)] * self.dim
+
+        domainCellBB = self._blocks.getDomainCellBB()
+        domainExtent = [s + 2 * sliceNormalizationGhostLayers for s in domainCellBB.size]
+        indexExpr = normalizeSlice(indexExpr, domainExtent)
+        targetCellBB = wlb.CellInterval.fromSlice(indexExpr)
+        targetCellBB.shift(*[a - sliceNormalizationGhostLayers for a in domainCellBB.min])
+
+        for block in self._blocks:
+            boundaryHandling = block[self._boundaryId]
+            ghostLayers = boundaryHandling.ghostLayers
+            intersection = self._blocks.getBlockCellBB(block).getExpanded(ghostLayers)
+            intersection.intersect(targetCellBB)
+            if not intersection.empty():
+                if maskCallback is not None:
+                    meshGridParams = [offset + np.arange(width)
+                                      for offset, width in zip(intersection.min, intersection.size)]
+                    indexArr = np.meshgrid(*meshGridParams, indexing='ij')
+                    mask = maskCallback(*indexArr)
+                else:
+                    mask = None
+                localTargetBB = self._blocks.transformGlobalToLocal(block, intersection)
+                localTargetBB.shift(ghostLayers, ghostLayers, ghostLayers)
+                block[self._boundaryId].setBoundary(boundaryObject, indexExpr=localTargetBB.toSlice(), maskArr=mask)
+
+
+# ----------------------------------------------------------------------------------------------------------------------
+
+
+class WalberlaFlagFieldInterface(FlagFieldInterface):
+
+    def __init__(self, block, flagFieldId):
+        self._block = block
+        self._flagFieldId = flagFieldId
+        self._flagArray = wlb.field.toArray(block[self._flagFieldId], withGhostLayers=True)
+        assert self._flagArray.shape[3] == 1
+        self._flagArray = self._flagArray[..., 0]
+
+        fluidFlag = self._block[self._flagFieldId].registerFlag('fluid')
+        self._boundaryObjectToName = {'fluid': fluidFlag}
+
+    def getFlag(self, boundaryObject):
+        if boundaryObject not in self._boundaryObjectToName:
+            name = self._makeBoundaryName(boundaryObject, self._boundaryObjectToName.values())
+            self._boundaryObjectToName[boundaryObject] = name
+            return self._block[self._flagFieldId].registerFlag(name)
+        else:
+            return self._block[self._flagFieldId].flag(self._boundaryObjectToName[boundaryObject])
+
+    def getName(self, boundaryObject):
+        return self._boundaryObjectToName[boundaryObject]
+
+    @property
+    def array(self):
+        return self._flagArray
+
+    @property
+    def boundaryObjects(self):
+        return self._boundaryObjectToName.keys()
+
+    def clear(self):
+        raise NotImplementedError()
+
+if __name__ == '__main__':
+    from lbmpy.creationfunctions import createLatticeBoltzmannMethod
+    from lbmpy.boundaries.boundaryconditions import NoSlip
+    from pystencils.slicing import makeSlice
+
+    blocks = wlb.createUniformBlockGrid(cellsPerBlock=(3, 3, 3), blocks=(2, 2, 2), oneBlockPerProcess=False)
+
+    lbMethod = createLatticeBoltzmannMethod(stencil='D3Q19', method='srt', relaxationRate=1.8)
+
+    wlb.field.addFlagFieldToStorage(blocks, 'flagField', nrOfBits=8, ghostLayers=1)
+    wlb.field.addToStorage(blocks, 'pdfField', float, fSize=len(lbMethod.stencil), ghostLayers=1)
+
+    bh = BoundaryHandling(blocks, lbMethod, 'pdfField', 'flagField')
+
+
+    def maskCallback(x, y, z):
+        return x > -100
+
+    bh.setBoundary(NoSlip(), makeSlice[0, :, :], maskCallback)
diff --git a/parallel/scenarios.py b/parallel/scenarios.py
new file mode 100644
index 0000000000000000000000000000000000000000..7e909e832e1fd016d7bb346a36a3d708c5225b74
--- /dev/null
+++ b/parallel/scenarios.py
@@ -0,0 +1,77 @@
+import sympy as sp
+import waLBerla as wlb
+from lbmpy.creationfunctions import createLatticeBoltzmannFunction, updateWithDefaultParameters, \
+    switchToSymbolicRelaxationRatesForEntropicMethods
+from lbmpy.parallel.boundaryhandling import BoundaryHandling
+
+
+class Scenario(object):
+
+    def __init__(self, blocks, methodParameters, optimizationParams, lbmKernel=None,
+                 initialVelocityCallback=None, preUpdateFunctions=[], kernelParams={},
+                 blockDataPrefix='', directCommunication=False):
+
+        methodParameters, optimizationParams = updateWithDefaultParameters(methodParameters, optimizationParams)
+        target = optimizationParams['target']
+
+        # ----- Create LBM kernel
+        if lbmKernel is None:
+            switchToSymbolicRelaxationRatesForEntropicMethods(methodParameters, kernelParams)
+            methodParameters['optimizationParams'] = optimizationParams
+            self._lbmKernel = createLatticeBoltzmannFunction(**methodParameters)
+        else:
+            self._lbmKernel = lbmKernel
+
+        self._blockDataPrefix = blockDataPrefix
+
+        # ----- Add fields
+        numPdfs = len(self._lbmKernel.method.stencils)
+        wlb.field.addToStorage(blocks, blockDataPrefix + 'pdfs', float, fSize=numPdfs, ghostLayers=1)
+        wlb.field.addFlagFieldToStorage(blocks, blockDataPrefix + "flags", nrOfBits=16, ghostLayers=1)
+
+        if target == 'gpu':
+            wlb.cuda.addGpuFieldToStorage(blocks, blockDataPrefix + "gpuPdfs", float, fSize=numPdfs,
+                                          usePitchedMem=False)
+
+        # ----- Create communication scheme
+        communicationStencil = methodParameters['stencil']
+
+        communicationFunctions = {
+            ('cpu', True): (wlb.createUniformDirectScheme, wlb.field.createMPIDatatypeInfo),
+            ('cpu', False): (wlb.createUniformBufferedScheme, wlb.field.createPackInfo),
+            ('gpu', True): (wlb.createUniformDirectScheme, wlb.cuda.createMPIDatatypeInfo),
+            ('gpu', False): (wlb.createUniformBufferedScheme, wlb.cuda.createPackInfo),
+        }
+        createScheme, createCommunicationData = communicationFunctions[(target, directCommunication)]
+        self._communicationScheme = createScheme(blocks, communicationStencil)
+        self._communicationScheme.addDataToCommunicate(createCommunicationData(blocks, blockDataPrefix + 'pdfs'))
+
+        # ----- Create boundary handling
+
+        self._boundaryHandling = BoundaryHandling(blocks, self._lbmKernel.method, blockDataPrefix + 'pdfs',
+                                                  blockDataPrefix + "flags", target=target)
+
+if __name__ == '__main__':
+    class A(object):
+        def __call__(self, call_from):
+            print("foo from A, call from %s" % call_from)
+
+
+    class B(object):
+        def __call__(self, call_from):
+            print( "foo from B, call from %s" % call_from)
+
+
+    class C(object):
+        def __call__(self, call_from):
+            print( "foo from C, call from %s" % call_from )
+
+
+    class D(A, B, C):
+        def foo(self):
+            for cls in D.__bases__:
+                cls.__call__(self, "D")
+
+
+    d = D()
+    d.foo()
\ No newline at end of file
diff --git a/plot2d.py b/plot2d.py
index 320c891343042ac21096f6db582b483de371fec0..38cc00a8a755cbffeb2edca638f25534b9ea5ea0 100644
--- a/plot2d.py
+++ b/plot2d.py
@@ -65,10 +65,11 @@ def plotBoundaryHandling(boundaryHandling, boundaryNameToColor=None):
             'noSlip': '#000000'
         }
 
-    boundaryNames = ['fluid'] + [b.name for b in boundaryHandling.getBoundaries()]
-    flagValues = [boundaryHandling.fluidFlag] + [boundaryHandling.getFlag(b)
-                                                 for b in boundaryHandling.getBoundaries()]
-
+    boundaryNames = []
+    flagValues = []
+    for name, flag in sorted(boundaryHandling.getBoundaryNameToFlagDict().items(), key=lambda l: l[1]):
+        boundaryNames.append(name)
+        flagValues.append(flag)
     defaultCycler = matplotlib.rcParams['axes.prop_cycle']
     colorValues = [fixedColors[name] if name in fixedColors else cycle['color']
                    for cycle, name in zip(defaultCycler, boundaryNames)]
diff --git a/scenarios.py b/scenarios.py
index feddd150763bc64d00d53965e7e550c79c9d32b2..f225d6e4b6132390b785670afa495529230b3e1b 100644
--- a/scenarios.py
+++ b/scenarios.py
@@ -24,10 +24,10 @@ at :mod:`lbmpy.creationfunctions`. The only mandatory keyword parameter is ``rel
 that defines the viscosity of the fluid (valid values being between 0 and 2).
 """
 import numpy as np
-import sympy as sp
 from pystencils.field import getLayoutOfArray, createNumpyArrayWithLayout
 from pystencils.slicing import sliceFromDirection, addGhostLayers, removeGhostLayers, normalizeSlice, makeSlice
-from lbmpy.creationfunctions import createLatticeBoltzmannFunction, updateWithDefaultParameters
+from lbmpy.creationfunctions import createLatticeBoltzmannFunction, updateWithDefaultParameters, \
+    switchToSymbolicRelaxationRatesForEntropicMethods
 from lbmpy.macroscopic_value_kernels import compileMacroscopicValuesGetter, compileMacroscopicValuesSetter
 from lbmpy.boundaries import BoundaryHandling, NoSlip, NoSlipFullWay, UBB, FixedDensity
 from lbmpy.stencils import getStencil
@@ -273,19 +273,7 @@ class Scenario(object):
 
         # Create kernel
         if lbmKernel is None:
-            if methodParameters['entropic']:
-                newRelaxationRates = []
-                for rr in methodParameters['relaxationRates']:
-                    if not isinstance(rr, sp.Symbol):
-                        dummyVar = sp.Dummy()
-                        newRelaxationRates.append(dummyVar)
-                        kernelParams[dummyVar.name] = rr
-                    else:
-                        newRelaxationRates.append(rr)
-                if len(newRelaxationRates) < 2:
-                    newRelaxationRates.append(sp.Dummy())
-                methodParameters['relaxationRates'] = newRelaxationRates
-
+            switchToSymbolicRelaxationRatesForEntropicMethods(methodParameters, kernelParams)
             optimizationParams['pdfArr'] = self._pdfArrays[0]
             methodParameters['optimizationParams'] = optimizationParams
             self._lbmKernel = createLatticeBoltzmannFunction(**methodParameters)
@@ -376,7 +364,7 @@ class Scenario(object):
     @property
     def velocity(self):
         """Velocity as numpy array"""
-        mask = np.logical_not(np.bitwise_and(self.boundaryHandling.flagField, self.boundaryHandling._fluidFlag))
+        mask = np.logical_not(self._boundaryHandling.getMask('fluid'))
         mask = np.repeat(mask[..., np.newaxis], self.dim, axis=2)
         return removeGhostLayers(np.ma.masked_array(self._velocity, mask), indexDimensions=1)
 
@@ -392,7 +380,7 @@ class Scenario(object):
     @property
     def density(self):
         """Density as numpy array"""
-        mask = np.logical_not(np.bitwise_and(self.boundaryHandling.flagField, self.boundaryHandling._fluidFlag))
+        mask = np.logical_not(self._boundaryHandling.getMask('fluid'))
         return removeGhostLayers(np.ma.masked_array(self._density, mask))
 
     @property
@@ -494,5 +482,3 @@ class Scenario(object):
         self._getMacroscopic(pdfs=self._pdfArrays[0], density=self._density, velocity=self._velocity,
                              **self.kernelParams)
 
-
-import pickle