From d3049429ce823e756d7abf1da897d1a1e5bf1023 Mon Sep 17 00:00:00 2001
From: Martin Bauer <martin.bauer@fau.de>
Date: Fri, 25 Aug 2017 14:37:39 +0200
Subject: [PATCH] Bugfixes in lbmpy boundary handling

---
 boundaries/boundaryhandling.py    |  5 +++--
 boundaries/periodicityhandling.py |  9 ++++++---
 scenarios.py                      | 12 ++++++------
 3 files changed, 15 insertions(+), 11 deletions(-)

diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py
index ce481c2d..bf74fe73 100644
--- a/boundaries/boundaryhandling.py
+++ b/boundaries/boundaryhandling.py
@@ -3,14 +3,14 @@ from lbmpy.boundaries.handlinginterface import GenericBoundaryHandling, FlagFiel
 from lbmpy.boundaries.periodicityhandling import PeriodicityHandling
 
 
-class BoundaryHandling(GenericBoundaryHandling, PeriodicityHandling):
+class BoundaryHandling(PeriodicityHandling, GenericBoundaryHandling):  # important: first periodicity, then boundary
 
     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)
 
         GenericBoundaryHandling.__init__(self, flagFieldInterface, pdfField, lbMethod, ghostLayers, target, openMP)
-        PeriodicityHandling.__init__(self, domainShape)
+        PeriodicityHandling.__init__(self, list(domainShape) + [len(lbMethod.stencil)])
 
     def setBoundary(self, boundaryObject, indexExpr=None, maskCallback=None, sliceNormalizationGhostLayers=1):
         mask = None
@@ -65,3 +65,4 @@ class NumpyFlagFieldInterface(FlagFieldInterface):
         self._boundaryObjectToName = {}
         self._boundaryObjectToFlag = {'fluid': np.dtype(self._array.dtype)(2**0)}
         self._nextFreeExponent = 1
+
diff --git a/boundaries/periodicityhandling.py b/boundaries/periodicityhandling.py
index 6c1fcb21..eb51d31e 100644
--- a/boundaries/periodicityhandling.py
+++ b/boundaries/periodicityhandling.py
@@ -6,7 +6,7 @@ class PeriodicityHandling(object):
         self._spatialShape = fieldShape[:-1]
         self._indexShape = fieldShape[-1]
         self._periodicity = list(periodicity)
-        self._dirty = False
+        self._periodicityDirty = False
         self._periodicityKernels = []
 
     @property
@@ -20,13 +20,16 @@ class PeriodicityHandling(object):
             assert isinstance(d, bool)
 
         self._periodicity = [x, y, z]
-        self._dirty = True
+        self._periodicityDirty = True
 
     def __call__(self, **kwargs):
         for k in self._periodicityKernels:
             k(**kwargs)
 
     def prepare(self):
+        if not self._periodicityDirty:
+            return
+
         self._periodicityKernels = []
         dim = len(self.flagField.shape)
         if dim == 2:
@@ -60,5 +63,5 @@ class PeriodicityHandling(object):
             else:
                 assert False
 
-        self._dirty = False
+        self._periodicityDirty = False
 
diff --git a/scenarios.py b/scenarios.py
index f225d6e4..1ef8e06e 100644
--- a/scenarios.py
+++ b/scenarios.py
@@ -141,12 +141,12 @@ def createForceDrivenChannel(force=1e-6, domainSize=None, dim=2, radius=None, le
             boundaryHandling.setBoundary(wallBoundary, sliceFromDirection(direction, dim))
     elif dim == 3:
         if roundChannel:
-            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] = wallIdx
+            def circleMaskCb(x, y, z):
+                yMid = np.max(y) // 2
+                zMid = np.max(z) // 2
+                return (y - yMid) ** 2 + (z - zMid) ** 2 > radius ** 2
+
+            boundaryHandling.setBoundary(wallBoundary, maskCallback=circleMaskCb)
         else:
             for direction in ('N', 'S', 'T', 'B'):
                 boundaryHandling.setBoundary(wallBoundary, sliceFromDirection(direction, dim))
-- 
GitLab