diff --git a/boundaries/boundaryhandling.py b/boundaries/boundaryhandling.py index c4b9cc684740ae782171f5bde61080d08241c70a..7473695dea50090fa8dc245561ef8c639609ff66 100644 --- a/boundaries/boundaryhandling.py +++ b/boundaries/boundaryhandling.py @@ -72,7 +72,6 @@ class BoundaryHandling(object): for boundaryIdx, boundaryFunc in enumerate(self._boundaryFunctions): idxField = createBoundaryIndexList(self.flagField, self._lbMethod.stencil, 2 ** boundaryIdx, self._fluidFlag, self._ghostLayers) - idxField = transformIndexListToStruct(idxField) ast = generateBoundaryHandling(self._symbolicPdfField, idxField, self._lbMethod, boundaryFunc) if self._target == 'cpu': @@ -134,23 +133,13 @@ class LbmMethodInfo(CustomCppCode): super(LbmMethodInfo, self).__init__(code, symbolsRead=set(), symbolsDefined=symbolsDefined) -def transformIndexListToStruct(arr): - #TODO create in correct form right away - dim = arr.shape[-1] -1 - coordinateNames = ['x', 'y', 'z'][:dim] - dataTypeInfo = [(name, np.int) for name in coordinateNames] + [('dir', np.int)] - indexArrStruct = np.empty((arr.shape[0]), dtype=np.dtype(dataTypeInfo)) - for idx, name in enumerate(coordinateNames): - indexArrStruct[name] = arr[:, idx] - indexArrStruct['dir'] = arr[:, -1] - return indexArrStruct - - def generateBoundaryHandling(pdfField, indexArr, lbMethod, boundaryFunctor): indexField = Field.createFromNumpyArray("indexField", indexArr) elements = [LbmMethodInfo(lbMethod)] - boundaryEqList = boundaryFunctor(pdfField, indexField[0]('dir'), 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 diff --git a/boundaries/createindexlist.py b/boundaries/createindexlist.py index 83ce85756f0edf6126f524197fbd0a4fb8b6029c..f81e43f27c8e2418c1f4968616b5a392fd78e3b8 100644 --- a/boundaries/createindexlist.py +++ b/boundaries/createindexlist.py @@ -13,6 +13,9 @@ except Exception: def _createBoundaryIndexListPython(flagFieldArr, nrOfGhostLayers, boundaryMask, fluidMask, stencil): + coordinateNames = ("x", "y", "z")[:len(flagFieldArr.shape)] + indexArrDtype = np.dtype([(name, np.int32) for name in coordinateNames] + [('dir', np.int32)]) + result = [] gl = nrOfGhostLayers for cell in itertools.product(*[range(gl, i-gl) for i in flagFieldArr.shape]): @@ -21,20 +24,24 @@ def _createBoundaryIndexListPython(flagFieldArr, nrOfGhostLayers, boundaryMask, for dirIdx, direction in enumerate(stencil): neighborCell = tuple([cell_i + dir_i for cell_i, dir_i in zip(cell, direction)]) if flagFieldArr[neighborCell] & boundaryMask: - result.append(list(cell) + [dirIdx]) + result.append(cell + (dirIdx,)) - return np.array(result, dtype=np.int32) + return np.array(result, dtype=indexArrDtype) def createBoundaryIndexList(flagField, stencil, boundaryMask, fluidMask, nrOfGhostLayers=1): + dim = len(flagField.shape) + coordinateNames = ("x", "y", "z")[:dim] + indexArrDtype = np.dtype([(name, np.int32) for name in coordinateNames] + [('dir', np.int32)]) + if cythonFuncsAvailable: stencil = np.array(stencil, dtype=np.int32) - if len(flagField.shape) == 2: + if dim == 2: idxList = createBoundaryIndexList2D(flagField, nrOfGhostLayers, boundaryMask, fluidMask, stencil) - elif len(flagField.shape) == 3: + elif dim == 3: idxList = createBoundaryIndexList3D(flagField, nrOfGhostLayers, boundaryMask, fluidMask, stencil) else: raise ValueError("Flag field has to be a 2 or 3 dimensional numpy array") - return np.array(idxList, dtype=np.int32) + return np.array(idxList, dtype=indexArrDtype) else: return _createBoundaryIndexListPython(flagField, nrOfGhostLayers, boundaryMask, fluidMask, stencil)