diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py index 5ca40d6ea7bdf2da1641d1dfa5f83e784c669407..4f7edcaabf950fe3e72babec3e0c4bfc2278de7c 100644 --- a/boundaries/boundaryhandling.py +++ b/boundaries/boundaryhandling.py @@ -30,6 +30,9 @@ class BoundaryHandling(object): if name is None: name = boundaryFunction.__name__ + if name in self._nameToIndex: + return 2 ** self._nameToIndex[name] + newIdx = len(self._boundaryFunctions) self._nameToIndex[name] = newIdx self._boundaryFunctions.append(boundaryFunction) @@ -45,7 +48,7 @@ class BoundaryHandling(object): def getFlag(self, name): return 2 ** self._nameToIndex[name] - def setBoundary(self, function, indexExpr, clearOtherBoundaries=True): + def setBoundary(self, function, indexExpr, maskArr=None): if hasattr(function, '__name__'): name = function.__name__ elif hasattr(function, 'name'): @@ -57,13 +60,12 @@ class BoundaryHandling(object): self.addBoundary(function, name) flag = self.getFlag(name) - if clearOtherBoundaries: + + if maskArr is None: self.flagField[indexExpr] = flag else: - # clear fluid flag - np.bitwise_and(self.flagField[indexExpr], np.invert(self._fluidFlag), self.flagField[indexExpr]) - # add new boundary flag - np.bitwise_or(self.flagField[indexExpr], flag, self.flagField[indexExpr]) + flagFieldView = self.flagField[indexExpr] + flagFieldView[maskArr] = flag self.invalidateIndexCache() diff --git a/boundaries/geometry.py b/boundaries/geometry.py new file mode 100644 index 0000000000000000000000000000000000000000..c586dec36e4369a1d8660fa51679e2ff1a9facf0 --- /dev/null +++ b/boundaries/geometry.py @@ -0,0 +1,23 @@ +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)) + diff --git a/serialscenario.py b/serialscenario.py index c30d443629820e0099a0613217138bb69cd6e6e4..1655ba1e14cac2c8345655c6ec33993de9b70f73 100644 --- a/serialscenario.py +++ b/serialscenario.py @@ -69,6 +69,9 @@ def createScenario(domainSize, boundarySetupFunction, methodParameters, optimiza pdfArrays[0], pdfArrays[1] = pdfArrays[1], pdfArrays[0] getMacroscopic(pdfs=pdfArrays[0], density=densityArr[0], velocity=velocityArr[0]) + #for vComp in range(velocityArr[0].shape[-1]): + # v = velocityArr[0][..., vComp] + # v[boundaryHandling.flagField != boundaryHandling._fluidFlag] = 0 return pdfArrays[0], densityArr[0], velocityArr[0] def gpuTimeLoop(timeSteps): @@ -112,7 +115,8 @@ def createLidDrivenCavity(domainSize, lidVelocity=0.005, optimizationParams={}, def createPressureGradientDrivenChannel(dim, pressureDifference, domainSize=None, radius=None, length=None, - lbmKernel=None, optimizationParams={}, kernelParams={}, **kwargs): + lbmKernel=None, optimizationParams={}, boundarySetupFunctions=[], + kernelParams={}, **kwargs): assert dim in (2, 3) if radius is not None: @@ -150,6 +154,9 @@ def createPressureGradientDrivenChannel(dim, pressureDifference, domainSize=None for direction in ('N', 'S', 'T', 'B'): boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, method.dim)) + for userFunction in boundarySetupFunctions: + userFunction(boundaryHandling, method, domainSize) + assert domainSize is not None return createScenario(domainSize, boundarySetupFunction, kwargs, optimizationParams, lbmKernel=lbmKernel, kernelParams=kernelParams) @@ -188,7 +195,7 @@ def createForceDrivenChannel(dim, force, domainSize=None, radius=None, length=No for direction in ('N', 'S', 'T', 'B'): boundaryHandling.setBoundary(noSlip, sliceFromDirection(direction, method.dim)) for userFunction in boundarySetupFunctions: - userFunction(boundaryHandling, method) + userFunction(boundaryHandling, method, domainSize) def periodicity(pdfArr): pdfArr[0, :, :] = pdfArr[-2, :, :]