From 551da07f3e8296fc3bbaa4c5138adb9a6eada77e Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Sun, 24 Oct 2021 18:59:44 -0700 Subject: [PATCH] Move `.grids` data into separate `cell_data` array. Also remove `Direction` enum --- gridlock/__init__.py | 1 - gridlock/draw.py | 70 +++++++++++++++++------------------- gridlock/grid.py | 86 ++++++++++++++++++++++---------------------- gridlock/read.py | 33 ++++++++--------- 4 files changed, 90 insertions(+), 100 deletions(-) diff --git a/gridlock/__init__.py b/gridlock/__init__.py index 171b4cc..d547794 100644 --- a/gridlock/__init__.py +++ b/gridlock/__init__.py @@ -16,7 +16,6 @@ Dependencies: - skimage [Grid.visualize_isosurface()] """ from .error import GridError -from .direction import Direction from .grid import Grid __author__ = 'Jan Petykiewicz' diff --git a/gridlock/draw.py b/gridlock/draw.py index 3a8ec8a..2dd18d7 100644 --- a/gridlock/draw.py +++ b/gridlock/draw.py @@ -6,15 +6,16 @@ from typing import List, Optional, Union, Sequence, Callable import numpy # type: ignore from float_raster import raster -from . import GridError, Direction from ._helpers import is_scalar +from . import GridError eps_callable_t = Callable[[numpy.ndarray, numpy.ndarray, numpy.ndarray], numpy.ndarray] def draw_polygons(self, - surface_normal: Union[Direction, int], + cell_data: numpy.ndarray, + surface_normal: int, center: numpy.ndarray, polygons: Sequence[numpy.ndarray], thickness: float, @@ -24,8 +25,8 @@ def draw_polygons(self, Draw polygons on an axis-aligned plane. Args: - surface_normal: Axis normal to the plane we're drawing on. Can be a `Direction` or - integer in `range(3)` + cell_data: Cell data to modify (e.g. created by `Grid.allocate()`) + surface_normal: Axis normal to the plane we're drawing on. Integer in `range(3)`. center: 3-element ndarray or list specifying an offset applied to all the polygons polygons: List of Nx2 or Nx3 ndarrays, each specifying the vertices of a polygon (non-closed, clockwise). If Nx3, the surface_normal coordinate is ignored. Each @@ -39,11 +40,6 @@ def draw_polygons(self, Raises: GridError """ - # Turn surface_normal into its integer representation - if isinstance(surface_normal, Direction): - surface_normal = surface_normal.value - assert(isinstance(surface_normal, int)) - if surface_normal not in range(3): raise GridError('Invalid surface_normal direction') @@ -66,8 +62,8 @@ def draw_polygons(self, % 'xyz'[surface_normal]) # Broadcast eps where necessary - if is_scalar(eps): - eps = [eps] * len(self.grids) + if numpy.size(eps) == 1: + eps = [eps] * len(cell_data) elif isinstance(eps, numpy.ndarray): raise GridError('ndarray not supported for eps') @@ -101,7 +97,7 @@ def draw_polygons(self, return numpy.insert(v_2d, surface_normal, (val,)) # iterate over grids - for i, grid in enumerate(self.grids): + for i, grid in enumerate(cell_data): # ## Evaluate or expand eps[i] if callable(eps[i]): # meshgrid over the (shifted) domain @@ -184,11 +180,12 @@ def draw_polygons(self, # ## Modify the grid g_slice = (i,) + tuple(numpy.s_[bdi_min[a]:bdi_max[a] + 1] for a in range(3)) - self.grids[g_slice] = (1 - w) * self.grids[g_slice] + w * eps_i + cell_data[g_slice] = (1 - w) * cell_data[g_slice] + w * eps_i def draw_polygon(self, - surface_normal: Union[Direction, int], + cell_data: numpy.ndarray, + surface_normal: int, center: numpy.ndarray, polygon: numpy.ndarray, thickness: float, @@ -198,8 +195,8 @@ def draw_polygon(self, Draw a polygon on an axis-aligned plane. Args: - surface_normal: Axis normal to the plane we're drawing on. Can be a Direction or - integer in range(3) + cell_data: Cell data to modify (e.g. created by `Grid.allocate()`) + surface_normal: Axis normal to the plane we're drawing on. Integer in `range(3)`. center: 3-element ndarray or list specifying an offset applied to the polygon polygon: Nx2 or Nx3 ndarray specifying the vertices of a polygon (non-closed, clockwise). If Nx3, the surface_normal coordinate is ignored. Must have at @@ -207,11 +204,12 @@ def draw_polygon(self, thickness: Thickness of the layer to draw eps: Value to draw with ('epsilon'). See `draw_polygons()` for details. """ - self.draw_polygons(surface_normal, center, [polygon], thickness, eps) + self.draw_polygons(cell_data, surface_normal, center, [polygon], thickness, eps) def draw_slab(self, - surface_normal: Union[Direction, int], + cell_data: numpy.ndarray, + surface_normal: int, center: numpy.ndarray, thickness: float, eps: Union[List[Union[float, eps_callable_t]], float, eps_callable_t], @@ -220,15 +218,13 @@ def draw_slab(self, Draw an axis-aligned infinite slab. Args: - surface_normal: Axis normal to the plane we're drawing on. Can be a `Direction` or - integer in `range(3)` + cell_data: Cell data to modify (e.g. created by `Grid.allocate()`) + surface_normal: Axis normal to the plane we're drawing on. Integer in `range(3)`. center: Surface_normal coordinate at the center of the slab thickness: Thickness of the layer to draw eps: Value to draw with ('epsilon'). See `draw_polygons()` for details. """ # Turn surface_normal into its integer representation - if isinstance(surface_normal, Direction): - surface_normal = surface_normal.value if surface_normal not in range(3): raise GridError('Invalid surface_normal direction') @@ -258,10 +254,11 @@ def draw_slab(self, [xyz_max[0], xyz_min[1]], [xyz_min[0], xyz_min[1]]], dtype=float) - self.draw_polygon(surface_normal, center_shift, p, thickness, eps) + self.draw_polygon(cell_data, surface_normal, center_shift, p, thickness, eps) def draw_cuboid(self, + cell_data: numpy.ndarray, center: numpy.ndarray, dimensions: numpy.ndarray, eps: Union[List[Union[float, eps_callable_t]], float, eps_callable_t], @@ -270,6 +267,7 @@ def draw_cuboid(self, Draw an axis-aligned cuboid Args: + cell_data: Cell data to modify (e.g. created by `Grid.allocate()`) center: 3-element ndarray or list specifying the cuboid's center dimensions: 3-element list or ndarray containing the x, y, and z edge-to-edge sizes of the cuboid @@ -280,11 +278,12 @@ def draw_cuboid(self, [+dimensions[0], -dimensions[1]], [-dimensions[0], -dimensions[1]]], dtype=float) / 2.0 thickness = dimensions[2] - self.draw_polygon(Direction.z, center, p, thickness, eps) + self.draw_polygon(cell_data, 2, center, p, thickness, eps) def draw_cylinder(self, - surface_normal: Union[Direction, int], + cell_data: numpy.ndarray, + surface_normal: int, center: numpy.ndarray, radius: float, thickness: float, @@ -295,8 +294,8 @@ def draw_cylinder(self, Draw an axis-aligned cylinder. Approximated by a num_points-gon Args: - surface_normal: Axis normal to the plane we're drawing on. Can be a `Direction` or - integer in `range(3)` + cell_data: Cell data to modify (e.g. created by `Grid.allocate()`) + surface_normal: Axis normal to the plane we're drawing on. Integer in `range(3)`. center: 3-element ndarray or list specifying the cylinder's center radius: cylinder radius thickness: Thickness of the layer to draw @@ -306,13 +305,14 @@ def draw_cylinder(self, theta = numpy.linspace(0, 2*numpy.pi, num_points, endpoint=False) x = radius * numpy.sin(theta) y = radius * numpy.cos(theta) - self.draw_polygon(surface_normal, center, polygon, thickness, eps) polygon = numpy.hstack((x[:, None], y[:, None])) + self.draw_polygon(cell_data, surface_normal, center, polygon, thickness, eps) def draw_extrude_rectangle(self, + cell_data: numpy.ndarray, rectangle: numpy.ndarray, - direction: Union[Direction, int], + direction: int, polarity: int, distance: float, ) -> None: @@ -320,16 +320,12 @@ def draw_extrude_rectangle(self, Extrude a rectangle of a previously-drawn structure along an axis. Args: + cell_data: Cell data to modify (e.g. created by `Grid.allocate()`) rectangle: 2x3 ndarray or list specifying the rectangle's corners - direction: Direction to extrude in. Direction enum or int in range(3) + direction: Direction to extrude in. Integer in `range(3)`. polarity: +1 or -1, direction along axis to extrude in distance: How far to extrude """ - # Turn extrude_direction into its integer representation - if isinstance(direction, Direction): - direction = direction.value - assert(isinstance(direction, int)) - s = numpy.sign(polarity) rectangle = numpy.array(rectangle, dtype=float) @@ -351,7 +347,7 @@ def draw_extrude_rectangle(self, thickness = distance eps_func = [] - for i, grid in enumerate(self.grids): + for i, grid in enumerate(cell_data): z = self.pos2ind(rectangle[0, :], i, round_ind=False, check_bounds=False)[direction] ind = [int(numpy.floor(z)) if i == direction else slice(None) for i in range(3)] @@ -373,5 +369,5 @@ def draw_extrude_rectangle(self, eps_func.append(f_eps) - self.draw_polygon(direction, center, p, thickness, eps_func) + self.draw_polygon(cell_data, direction, center, p, thickness, eps_func) diff --git a/gridlock/grid.py b/gridlock/grid.py index 3c1c4a8..3386823 100644 --- a/gridlock/grid.py +++ b/gridlock/grid.py @@ -7,8 +7,8 @@ import pickle import warnings import copy -from . import GridError, Direction from ._helpers import is_scalar +from . import GridError __author__ = 'Jan Petykiewicz' @@ -18,12 +18,18 @@ eps_callable_type = Callable[[numpy.ndarray, numpy.ndarray, numpy.ndarray], nump class Grid: """ - Simulation grid generator intended for electromagnetic simulations. - Can be used to generate non-uniform rectangular grids (the entire grid + Simulation grid metadata for finite-difference simulations. + + Can be used to generate non-uniform rectangular grids (the entire grid is generated based on the coordinates of the boundary points). Also does straightforward natural <-> grid unit conversion. - `self.grids[i][a,b,c]` contains the value of epsilon for the cell located around + This class handles data describing the grid, and should be paired with a + (separate) ndarray that contains the actual data in each cell. The `allocate()` + method can be used to create this ndarray. + + The resulting `cell_data[i, a, b, c]` should correspond to the value in the + `i`-th grid, in the cell centered around ``` (xyz[0][a] + dxyz[0][a] * shifts[i, 0], xyz[1][b] + dxyz[1][b] * shifts[i, 1], @@ -47,9 +53,6 @@ class Grid: exyz: List[numpy.ndarray] """Cell edges. Monotonically increasing without duplicates.""" - grids: numpy.ndarray - """epsilon (or mu, or whatever) grids. shape is (num_grids, X, Y, Z)""" - periodic: List[bool] """For each axis, determines how far the rightmost boundary gets shifted. """ @@ -103,6 +106,20 @@ class Grid: """ return numpy.array([coord.size - 1 for coord in self.exyz], dtype=int) + @property + def num_grids(self) -> int: + """ + The number of grids (number of shifts) + """ + return self.shifts.shape[0] + + @property + def cell_data_shape(self): + """ + The shape of the cell_data ndarray (num_grids, *self.shape). + """ + return numpy.hstack((self.num_grids, self.shape)) + @property def dxyz_with_ghost(self) -> List[numpy.ndarray]: """ @@ -218,16 +235,30 @@ class Grid: Returns: `[grid.shifted_dxyz(which_shifts=a)[a] for a in range(3)]` """ - if len(self.grids) != 3: - raise GridError('autoshifting requires exactly 3 grids') + if self.num_grids != 3: + raise GridError('Autoshifting requires exactly 3 grids') return [self.shifted_dxyz(which_shifts=a)[a] for a in range(3)] + def allocate(self, fill_value: Optional[float] = 1.0, dtype=numpy.float64) -> numpy.ndarray: + """ + Allocate an ndarray for storing grid data. + + Args: + fill_value: Value to initialize the grid to. If None, an + uninitialized array is returned. + dtype: Numpy dtype for the array. Default is `numpy.float64`. + + Returns: + The allocated array + """ + if fill_value is None: + return numpy.empty(self.cell_data_shape) + else: + return numpy.full(self.cell_data_shape, fill_value) def __init__(self, pixel_edge_coordinates: Sequence[numpy.ndarray], shifts: numpy.ndarray = Yee_Shifts_E, - initial: Union[float, numpy.ndarray] = 1.0, - num_grids: Optional[int] = None, periodic: Union[bool, Sequence[bool]] = False, ) -> None: """ @@ -238,12 +269,6 @@ class Grid: x=`x1`, the second has edges x=`x1` and x=`x2`, etc.) shifts: Nx3 array containing `[x, y, z]` offsets for each of N grids. E-field Yee shifts are used by default. - initial: Grids are initialized to this value. If scalar, all grids are initialized - with ndarrays full of the scalar. If a list of scalars, `grid[i]` is initialized to an - ndarray full of `initial[i]`. If a list of ndarrays of the same shape as the grids, `grid[i]` - is set to `initial[i]`. Default `1.0`. - num_grids: How many grids to create. Must be <= `shifts.shape[0]`. - Default is `shifts.shape[0]` periodic: Specifies how the sizes of edge cells are calculated; see main class documentation. List of 3 bool, or a single bool that gets broadcast. Default `False`. @@ -276,33 +301,6 @@ class Grid: # TODO: Test negative shifts warnings.warn('Negative shifts are still experimental and mostly untested, be careful!', stacklevel=2) - num_shifts = self.shifts.shape[0] - if num_grids is None: - num_grids = num_shifts - elif num_grids > num_shifts: - raise GridError('Number of grids exceeds number of shifts (%u)' % num_shifts) - - grids_shape = hstack((num_grids, self.shape)) - if isinstance(initial, (float, int)): - if isinstance(initial, int): - warnings.warn('Initial value is an int, grids will be integer-typed!', stacklevel=2) - self.grids = numpy.full(grids_shape, initial) - else: - if len(initial) < num_grids: - raise GridError('Too few initial grids specified!') - - self.grids = numpy.empty(grids_shape) - for i in range(num_grids): - if is_scalar(initial[i]): - if initial[i] is not None: - if isinstance(initial[i], int): - warnings.warn('Initial value is an int, grid {} will be integer-typed!'.format(i), stacklevel=2) - self.grids[i] = numpy.full(self.shape, initial[i]) - else: - if not numpy.array_equal(initial[i].shape, self.shape): - raise GridError('Initial grid sizes must match given coordinates') - self.grids[i] = initial[i] - @staticmethod def load(filename: str) -> 'Grid': """ diff --git a/gridlock/read.py b/gridlock/read.py index e049d77..472508e 100644 --- a/gridlock/read.py +++ b/gridlock/read.py @@ -5,8 +5,8 @@ from typing import Dict, Optional, Union, Any import numpy # type: ignore -from . import GridError, Direction from ._helpers import is_scalar +from . import GridError # .visualize_* uses matplotlib # .visualize_isosurface uses skimage @@ -14,7 +14,8 @@ from ._helpers import is_scalar def get_slice(self, - surface_normal: Union[Direction, int], + cell_data: numpy.ndarray, + surface_normal: int, center: float, which_shifts: int = 0, sample_period: int = 1 @@ -24,8 +25,8 @@ def get_slice(self, Interpolates if given a position between two planes. Args: - surface_normal: Axis normal to the plane we're displaying. Can be a `Direction` or - integer in `range(3)` + cell_data: Cell data to slice + surface_normal: Axis normal to the plane we're displaying. Integer in `range(3)`. center: Scalar specifying position along surface_normal axis. which_shifts: Which grid to display. Default is the first grid (0). sample_period: Period for down-sampling the image. Default 1 (disabled) @@ -43,9 +44,6 @@ def get_slice(self, if not is_scalar(which_shifts) or which_shifts < 0: raise GridError('Invalid which_shifts') - # Turn surface_normal into its integer representation - if isinstance(surface_normal, Direction): - surface_normal = surface_normal.value if surface_normal not in range(3): raise GridError('Invalid surface_normal direction') @@ -70,7 +68,7 @@ def get_slice(self, sliced_grid = numpy.zeros(self.shape[surface]) for ci, weight in zip(centers, w): s = tuple(ci if a == surface_normal else numpy.s_[::sp] for a in range(3)) - sliced_grid += weight * self.grids[which_shifts][tuple(s)] + sliced_grid += weight * cell_data[which_shifts][tuple(s)] # Remove extra dimensions sliced_grid = numpy.squeeze(sliced_grid) @@ -79,7 +77,8 @@ def get_slice(self, def visualize_slice(self, - surface_normal: Union[Direction, int], + cell_data: numpy.ndarray, + surface_normal: int, center: float, which_shifts: int = 0, sample_period: int = 1, @@ -91,8 +90,7 @@ def visualize_slice(self, Interpolates if given a position between two planes. Args: - surface_normal: Axis normal to the plane we're displaying. Can be a `Direction` or - integer in `range(3)` + surface_normal: Axis normal to the plane we're displaying. Integer in `range(3)`. center: Scalar specifying position along surface_normal axis. which_shifts: Which grid to display. Default is the first grid (0). sample_period: Period for down-sampling the image. Default 1 (disabled) @@ -100,14 +98,11 @@ def visualize_slice(self, """ from matplotlib import pyplot - # Set surface normal to its integer value - if isinstance(surface_normal, Direction): - surface_normal = surface_normal.value - if pcolormesh_args is None: pcolormesh_args = {} - grid_slice = self.get_slice(surface_normal=surface_normal, + grid_slice = self.get_slice(cell_data=cell_data, + surface_normal=surface_normal, center=center, which_shifts=which_shifts, sample_period=sample_period) @@ -129,6 +124,7 @@ def visualize_slice(self, def visualize_isosurface(self, + cell_data: numpy.ndarray, level: Optional[float] = None, which_shifts: int = 0, sample_period: int = 1, @@ -139,6 +135,7 @@ def visualize_isosurface(self, Draw an isosurface plot of the device. Args: + cell_data: Cell data to visualize level: Value at which to find isosurface. Default (None) uses mean value in grid. which_shifts: Which grid to display. Default is the first grid (0). sample_period: Period for down-sampling the image. Default 1 (disabled) @@ -150,8 +147,8 @@ def visualize_isosurface(self, # Claims to be unused, but needed for subplot(projection='3d') from mpl_toolkits.mplot3d import Axes3D - # Get data from self.grids - grid = self.grids[which_shifts][::sample_period, ::sample_period, ::sample_period] + # Get data from cell_data + grid = cell_data[which_shifts][::sample_period, ::sample_period, ::sample_period] if level is None: level = grid.mean()