lbmpy boundary geometry: image loading

boundaries/geometry.py

-import numpy as np
-import scipy.misc
-from pystencils.slicing import makeSlice, normalizeSlice, shiftSlice
-
-
-class BlackAndWhiteImageBoundary:
-
-    def __init__(self, imagePath, boundaryFunction, targetSlice=makeSlice[:, :, :]):
-        self.imgArr = scipy.misc.imread(imagePath, flatten=True).astype(int)
-        self.imgArr = np.rot90(self.imgArr, 3)
-
-        self._boundaryFunction = boundaryFunction
-        self._targetSlice = targetSlice
-
-    def __call__(self, boundaryHandling, method, domainSize, **kwargs):
-        normalizedSlice = normalizeSlice(self._targetSlice, domainSize)
-        normalizedSlice = shiftSlice(normalizedSlice, 1)
-        targetSize = [s.stop - s.start for s in normalizedSlice]
-        img = scipy.misc.imresize(self.imgArr, size=targetSize)
-        img[img <= 254] = 0
-        img[img > 254] = 1
-        boundaryHandling.setBoundary(self._boundaryFunction, normalizedSlice, maskArr=(img == 0))
-

boundaries/handlinginterface.py

 import numpy as np
 
 from pystencils.field import Field
-from pystencils.slicing import normalizeSlice
+from pystencils.slicing import normalizeSlice, shiftSlice
 from lbmpy.boundaries.boundary_kernel import generateIndexBoundaryKernel
 from lbmpy.boundaries.createindexlist import createBoundaryIndexList
 from pystencils.cache import memorycache
@@ -210,22 +210,53 @@ class GenericBoundaryHandling(object):
     def reserveFlag(self, boundaryObject):
         self._flagFieldInterface.getFlag(boundaryObject)
 
-    def setBoundary(self, boundaryObject, indexExpr=None, maskCallback=None):
+    def setBoundary(self, boundaryObject, indexExpr=None, maskCallback=None, includeGhostLayers=True):
+        """
+        Sets boundary using either a rectangular slice, a boolean mask or a combination of both
+        
+        :param boundaryObject: instance of a boundary object that should be set 
+        :param indexExpr: a slice object (can be created with makeSlice[]) that selects a part of the domain where
+                          the boundary should be set. If none, the complete domain is selected which makes only sense
+                          if a maskCallback is passed. The slice can have ':' placeholders, which are interpreted 
+                          depending on the 'includeGhostLayers' parameter i.e. if it is True, the slice extends
+                          into the ghost layers
+        :param maskCallback: callback function getting x,y (z) parameters of the cell midpoints and returning a 
+                             boolean mask with True entries where boundary cells should be set. 
+                             The x, y, z arrays have 2D/3D shape such that they can be used directly 
+                             to create the boolean return array. i.e return x < 10 sets boundaries in cells with
+                             midpoint x coordinate smaller than 10.
+        :param includeGhostLayers: if this parameter is False, boundaries can not be set in the ghost
+                                   layer, because index 0 is the first inner layer and -1 is interpreted in the Python
+                                   way as maximum. If this parameter is True, the lower ghost layers have index 0, and
+                                   placeholders ':' in index expressions extend into the ghost layers.
+        """
+        if indexExpr is None:
+            indexExpr = [slice(None, None, None)] * len(self.flagField.shape)
+        if not includeGhostLayers:
+            domainSize = [i - 2 * self.ghostLayers for i in self._flagFieldInterface.array.shape]
+            indexExpr = normalizeSlice(indexExpr, domainSize)
+            indexExpr = shiftSlice(indexExpr, self.ghostLayers)
+        else:
+            indexExpr = normalizeSlice(indexExpr, self._flagFieldInterface.array.shape)
+
         mask = None
-        flagField = self._flagFieldInterface.array
         if maskCallback is not None:
-            gridParams = [offset + np.arange(-self.ghostLayers, s - self.ghostLayers) + 0.5
-                          for s, offset in zip(flagField.shape, self.offset)]
+            gridParams = []
+            for s, offset in zip(indexExpr, self.offset):
+                if isinstance(s, slice):
+                    gridParams.append(np.arange(s.start, s.stop) + offset + 0.5 - self.ghostLayers)
+                else:
+                    gridParams.append(s + offset + 0.5 - self.ghostLayers)
             indices = np.meshgrid(*gridParams, indexing='ij')
             mask = maskCallback(*indices)
-        return self.setBoundaryWithMaskArray(boundaryObject, indexExpr, mask)
+        return self._setBoundaryWithMaskArray(boundaryObject, indexExpr, mask)
 
-    def setBoundaryWithMaskArray(self, boundaryObject, indexExpr=None, maskArr=None):
+    def _setBoundaryWithMaskArray(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 indexExpr: slice expression, where boundary should be set. ghost layers are expected to have coord=0
         :param maskArr: optional boolean (masked) array specifying where the boundary should be set
         """
         if indexExpr is None:

geometry.py

 import numpy as np
-from lbmpy.boundaries import NoSlip
+from lbmpy.boundaries import NoSlip, UBB
+from pystencils.slicing import normalizeSlice, shiftSlice, sliceIntersection
+
+
+def addParabolicVelocityInflow(boundaryHandling, u_max, indexExpr, velCoord=0, diameter=None):
+    def velocityInfoCallback(boundaryData):
+        for i, name in enumerate(['vel_0', 'vel_1', 'vel_2']):
+            if i != velCoord:
+                boundaryData[name] = 0
+        if diameter is None:
+            radius = boundaryHandling.domainShape[velCoord] // 2
+        else:
+            radius = diameter // 2
+
+        y, z = boundaryData.linkPositions(1), boundaryData.linkPositions(2)
+        centeredY = y - radius
+        centeredZ = z - radius
+        distToCenter = np.sqrt(centeredY ** 2 + centeredZ ** 2)
+        boundaryData['vel_%d' % (velCoord,)] = u_max * (1 - distToCenter / radius)
+
+    inflow = UBB(velocityInfoCallback, dim=boundaryHandling.dim)
+    boundaryHandling.setBoundary(inflow, indexExpr=indexExpr, includeGhostLayers=False)
 
 
 def addPipe(boundaryHandling, diameter, boundary=NoSlip()):
@@ -8,7 +29,7 @@ def addPipe(boundaryHandling, diameter, boundary=NoSlip()):
 
     :param boundaryHandling: boundary handling object, works for serial and parallel versions 
     :param diameter: pipe diameter, can be either a constant value or a callback function.
-                     the callback function has the signature (xCoordArray, domainShapeInCells) and has to return
+                     the callback function has the signature (xCoordArray, domainShapeInCells) andp has to return
                      a array of same shape as the received xCoordArray, with the diameter for each x position
     :param boundary: boundary object that is set at the wall, defaults to NoSlip (bounce back)
     """
@@ -33,3 +54,61 @@ def addPipe(boundaryHandling, diameter, boundary=NoSlip()):
         return distance > radiusSq
 
     boundaryHandling.setBoundary(boundary, maskCallback=callback)
+
+
+def addBlackAndWhiteImage(boundaryHandling, imageFile, targetSlice=None, plane=(0, 1), boundary=NoSlip(),
+                          keepAspectRatio=False):
+    try:
+        from scipy.misc import imread
+        from scipy.ndimage import zoom
+    except ImportError:
+        raise ImportError("scipy image read could not be imported! Install 'scipy' and 'pillow'")
+
+    domainSize = boundaryHandling.domainShape
+    if targetSlice is None:
+        targetSlice = [slice(None, None, None)] * len(domainSize)
+
+    dim = len(boundaryHandling.domainShape)
+
+    imageSlice = normalizeSlice(targetSlice, domainSize)
+    targetSize = [imageSlice[i].stop - imageSlice[i].start for i in plane]
+
+    imgArr = imread(imageFile, flatten=True).astype(int)
+    imgArr = np.rot90(imgArr, 3)
+
+    zoomFactor = [targetSize[i] / imgArr.shape[i] for i in range(2)]
+    if keepAspectRatio:
+        zoomFactor = min(zoomFactor)
+    zoomedImage = zoom(imgArr, zoomFactor, order=0)
+
+    # binarize
+    zoomedImage[zoomedImage <= 254] = 0
+    zoomedImage[zoomedImage > 254] = 1
+    zoomedImage = np.logical_not(zoomedImage.astype(np.bool))
+
+    # resize necessary if aspect ratio should be constant
+    if zoomedImage.shape != targetSize:
+        resizedImage = np.zeros(targetSize, dtype=np.bool)
+        mid = [(ts - s)//2 for ts, s in zip(targetSize, zoomedImage.shape)]
+        resizedImage[mid[0]:zoomedImage.shape[0]+mid[0], mid[1]:zoomedImage.shape[1]+mid[1]] = zoomedImage
+        zoomedImage = resizedImage
+
+    def callback(*coordinates):
+        result = np.zeros_like(coordinates[0], dtype=np.bool)
+        maskStart = [int(coordinates[i][(0,) * dim] - 0.5) for i in range(dim)]
+        maskEnd = [int(coordinates[i][(-1,) * dim] + 1 - 0.5) for i in range(dim)]
+
+        maskSlice = [slice(start, stop) for start, stop in zip(maskStart, maskEnd)]
+        intersectionSlice = sliceIntersection(maskSlice, imageSlice)
+        if intersectionSlice is None:
+            return result
+        else:
+            maskTargetSlice = shiftSlice(intersectionSlice, [-e for e in maskStart])
+            imageTargetSlice = shiftSlice(intersectionSlice, [-s.start for s in imageSlice])
+            maskTargetSlice = [maskTargetSlice[i] if i in plane else slice(None, None) for i in range(dim)]
+            imageTargetSlice = [imageTargetSlice[i] if i in plane else np.newaxis for i in range(dim)]
+            result[maskTargetSlice] = zoomedImage[imageTargetSlice]
+            return result
+
+    boundaryHandling.setBoundary(boundary, indexExpr=imageSlice, maskCallback=callback, includeGhostLayers=False)
+

parallel/boundaryhandling.py

@@ -34,7 +34,8 @@ class BoundaryHandling(object):
         for block in self._blocks:
             block[self._boundaryId].triggerReinitializationOfBoundaryData()
 
-    def setBoundary(self, boundaryObject, indexExpr=None, maskCallback=None, sliceNormalizationGhostLayers=1):
+    def setBoundary(self, boundaryObject, indexExpr=None, maskCallback=None, includeGhostLayers=True):
+
         if indexExpr is None:
             indexExpr = [slice(None, None, None)] * 3
 
@@ -42,6 +43,8 @@ class BoundaryHandling(object):
             return
 
         gl = self._blocks[0][self._boundaryId].ghostLayers
+        sliceNormalizationGhostLayers = gl if includeGhostLayers else 0
+
         for block in self._blocks:
             block[self._boundaryId].reserveFlag(boundaryObject)
 

plot2d.py

@@ -2,7 +2,7 @@ from pystencils.plot2d import *
 from pystencils.slicing import normalizeSlice
 
 
-def plotBoundaryHandling(boundaryHandling, slice=None, boundaryNameToColor=None, legend=True):
+def plotBoundaryHandling(boundaryHandling, indexExpr=None, boundaryNameToColor=None, legend=True):
     """
     Shows boundary cells
 
@@ -14,7 +14,7 @@ def plotBoundaryHandling(boundaryHandling, slice=None, boundaryNameToColor=None,
 
     boundaryHandling.prepare()
 
-    if len(boundaryHandling.flagField.shape) != 2 and slice is None:
+    if len(boundaryHandling.flagField.shape) != 2 and indexExpr is None:
         raise ValueError("To plot 3D boundary handlings a slice has to be passed")
 
     if boundaryNameToColor:
@@ -42,8 +42,8 @@ def plotBoundaryHandling(boundaryHandling, slice=None, boundaryNameToColor=None,
     norm = matplotlib.colors.BoundaryNorm(bounds, cmap.N)
 
     flagField = boundaryHandling.flagField
-    if slice:
-        flagField = flagField[normalizeSlice(slice, flagField.shape)]
+    if indexExpr:
+        flagField = flagField[normalizeSlice(indexExpr, flagField.shape)]
     flagField = flagField.swapaxes(0, 1)
     plt.imshow(flagField, interpolation='none', origin='lower',
                cmap=cmap, norm=norm)