Skip to content
Snippets Groups Projects

Use closest normal for boundary index list with single_link

Merged Markus Holzer requested to merge holzer/pystencils:CreateIndexList into master
All threads resolved!
Files
3
import itertools
import warnings
import warnings
import numpy as np
import numpy as np
@@ -35,45 +34,59 @@ def numpy_data_type_for_boundary_object(boundary_object, dim):
@@ -35,45 +34,59 @@ def numpy_data_type_for_boundary_object(boundary_object, dim):
+ [(i[0], i[1].numpy_dtype) for i in boundary_object.additional_data], align=True)
+ [(i[0], i[1].numpy_dtype) for i in boundary_object.additional_data], align=True)
def _create_boundary_neighbor_index_list_python(flag_field_arr, nr_of_ghost_layers, boundary_mask,
def _create_index_list_python(flag_field_arr, boundary_mask,
fluid_mask, stencil, single_link):
fluid_mask, stencil, single_link, inner_or_boundary=False, nr_of_ghost_layers=None):
 
 
if inner_or_boundary and nr_of_ghost_layers is None:
 
raise ValueError("If inner_or_boundary is set True the number of ghost layers "
 
"around the inner domain has to be specified")
 
 
if nr_of_ghost_layers is None:
 
nr_of_ghost_layers = 0
 
coordinate_names = boundary_index_array_coordinate_names[:len(flag_field_arr.shape)]
coordinate_names = boundary_index_array_coordinate_names[:len(flag_field_arr.shape)]
index_arr_dtype = np.dtype([(name, np.int32) for name in coordinate_names] + [(direction_member_name, np.int32)])
index_arr_dtype = np.dtype([(name, np.int32) for name in coordinate_names] + [(direction_member_name, np.int32)])
result = []
boundary_cells = np.transpose(np.where(flag_field_arr == boundary_mask))
gl = nr_of_ghost_layers
for i in range(len(flag_field_arr.shape)):
for cell in itertools.product(*reversed([range(gl, i - gl) for i in flag_field_arr.shape])):
boundary_cells = boundary_cells[boundary_cells[:, i].argsort(kind='mergesort')]
cell = cell[::-1]
if not flag_field_arr[cell] & fluid_mask:
fluid_cells = []
continue
for cell in boundary_cells:
 
cell = tuple(cell)
for dir_idx, direction in enumerate(stencil):
for dir_idx, direction in enumerate(stencil):
neighbor_cell = tuple([cell_i + dir_i for cell_i, dir_i in zip(cell, direction)])
neighbor_cell = tuple([cell_i + dir_i for cell_i, dir_i in zip(cell, direction)])
if flag_field_arr[neighbor_cell] & boundary_mask:
if any(not 0 + nr_of_ghost_layers <= e < upper - nr_of_ghost_layers
result.append(cell + (dir_idx,))
for e, upper in zip(neighbor_cell, flag_field_arr.shape)):
if single_link:
continue
break
if flag_field_arr[neighbor_cell] & fluid_mask and fluid_cells.count(neighbor_cell) == 0:
fluid_cells.append(neighbor_cell)
return np.array(result, dtype=index_arr_dtype)
 
if len(flag_field_arr.shape) == 3:
 
fluid_cells.sort(key=lambda tup: (tup[-1], tup[-2], tup[0]))
 
else:
 
fluid_cells.sort(key=lambda tup: (tup[-1], tup[0]))
def _create_boundary_cell_index_list_python(flag_field_arr, boundary_mask,
cells_to_iterate = fluid_cells if inner_or_boundary else boundary_cells
fluid_mask, stencil, single_link):
checkmask = boundary_mask if inner_or_boundary else fluid_mask
coordinate_names = boundary_index_array_coordinate_names[:len(flag_field_arr.shape)]
index_arr_dtype = np.dtype([(name, np.int32) for name in coordinate_names] + [(direction_member_name, np.int32)])
result = []
result = []
for cell in itertools.product(*reversed([range(0, i) for i in flag_field_arr.shape])):
for cell in cells_to_iterate:
cell = cell[::-1]
cell = tuple(cell)
if not flag_field_arr[cell] & boundary_mask:
sum_cells = np.zeros(len(cell))
continue
for dir_idx, direction in enumerate(stencil):
for dir_idx, direction in enumerate(stencil):
neighbor_cell = tuple([cell_i + dir_i for cell_i, dir_i in zip(cell, direction)])
neighbor_cell = tuple([cell_i + dir_i for cell_i, dir_i in zip(cell, direction)])
if any(not 0 <= e < upper for e, upper in zip(neighbor_cell, flag_field_arr.shape)):
if any(not 0 <= e < upper for e, upper in zip(neighbor_cell, flag_field_arr.shape)):
continue
continue
if flag_field_arr[neighbor_cell] & fluid_mask:
if flag_field_arr[neighbor_cell] & checkmask:
result.append(cell + (dir_idx,))
if single_link:
if single_link:
break
sum_cells += np.array(direction)
 
else:
 
result.append(tuple(cell) + (dir_idx,))
 
 
if single_link and any(sum_cells != 0):
 
idx = np.argmax(np.inner(sum_cells, stencil))
 
result.append(tuple(cell) + (idx,))
return np.array(result, dtype=index_arr_dtype)
return np.array(result, dtype=index_arr_dtype)
@@ -119,10 +132,8 @@ def create_boundary_index_list(flag_field, stencil, boundary_mask, fluid_mask,
@@ -119,10 +132,8 @@ def create_boundary_index_list(flag_field, stencil, boundary_mask, fluid_mask,
else:
else:
if flag_field.size > 1e6:
if flag_field.size > 1e6:
warnings.warn("Boundary setup may take very long! Consider installing cython to speed it up")
warnings.warn("Boundary setup may take very long! Consider installing cython to speed it up")
if inner_or_boundary:
return _create_index_list_python(*args_no_gl, inner_or_boundary=inner_or_boundary,
return _create_boundary_neighbor_index_list_python(*args)
nr_of_ghost_layers=nr_of_ghost_layers)
else:
return _create_boundary_cell_index_list_python(*args_no_gl)
def create_boundary_index_array(flag_field, stencil, boundary_mask, fluid_mask, boundary_object,
def create_boundary_index_array(flag_field, stencil, boundary_mask, fluid_mask, boundary_object,
Loading