Add class functions to simplify kGrid creation

examples/at_circular_piston_3D/at_circular_piston_3D.py

@@ -48,20 +48,9 @@
 # GRID
 # --------------------
 
-# calculate the grid spacing based on the PPW and F0
-dx: float = c0 / (ppw * source_f0)  # [m]
-
-# compute the size of the grid
-# is round_even needed?
-Nx: int = round_even(axial_size / dx)
-Ny: int = round_even(lateral_size / dx)
-Nz: int = Ny
-
-grid_size_points = Vector([Nx, Ny, Nz])
-grid_spacing_meters = Vector([dx, dx, dx])
-
-# create the k-space grid
-kgrid = kWaveGrid(grid_size_points, grid_spacing_meters)
+kgrid = kWaveGrid.from_domain(
+    dimensions=np.array([axial_size, lateral_size, lateral_size]), frequency=source_f0, sound_speed_min=c0, points_per_wavelength=ppw
+)
 
 # compute points per temporal period
 ppp: int = round(ppw / cfl)
@@ -113,8 +102,8 @@
 sensor = kSensor()
 
 # set sensor mask to record central plane, not including the source point
-sensor.mask = np.zeros((Nx, Ny, Nz), dtype=bool)
-sensor.mask[1:, :, Nz // 2] = True
+sensor.mask = np.zeros(kgrid.N, dtype=bool)
+sensor.mask[1:, :, kgrid.Nz // 2] = True
 
 # record the pressure
 sensor.record = ["p"]
@@ -143,10 +132,10 @@
 amp, _, _ = extract_amp_phase(sensor_data["p"].T, 1.0 / kgrid.dt, source_f0, dim=1, fft_padding=1, window="Rectangular")
 
 # reshape data
-amp = np.reshape(amp, (Nx - 1, Ny), order="F")
+amp = np.reshape(amp, (kgrid.Nx - 1, kgrid.Ny), order="F")
 
 # extract pressure on axis
-amp_on_axis = amp[:, Ny // 2]
+amp_on_axis = amp[:, kgrid.Ny // 2]
 
 # define axis vectors for plotting
 x_vec = kgrid.x_vec[1:, :] - kgrid.x_vec[0]
@@ -161,7 +150,7 @@
 
 # define radius and axis
 a: float = source_diam / 2.0
-x_max: float = (Nx - 1) * dx
+x_max: float = (kgrid.Nx - 1) * kgrid.dx
 delta_x: float = x_max / 10000.0
 x_ref: float = np.arange(0.0, x_max + delta_x, delta_x, dtype=float)
 
@@ -194,7 +183,7 @@
 # plot the source mask (pml is outside the grid in this example)
 fig2, ax2 = plt.subplots(1, 1)
 ax2.pcolormesh(
-    1e3 * np.squeeze(kgrid.y_vec), 1e3 * np.squeeze(kgrid.x_vec), np.flip(source.p_mask[:, :, Nz // 2], axis=0), shading="nearest"
+    1e3 * np.squeeze(kgrid.y_vec), 1e3 * np.squeeze(kgrid.x_vec), np.flip(source.p_mask[:, :, kgrid.Nz // 2], axis=0), shading="nearest"
 )
 ax2.set(xlabel="y [mm]", ylabel="x [mm]", title="Source Mask")
 

examples/at_focused_bowl_3D/at_focused_bowl_3D.py

@@ -57,18 +57,9 @@
 # --------------------
 
 # calculate the grid spacing based on the PPW and F0
-dx: float = c0 / (ppw * source_f0)  # [m]
-
-# compute the size of the grid
-Nx: int = round_even(axial_size / dx) + source_x_offset
-Ny: int = round_even(lateral_size / dx)
-Nz: int = Ny
-
-grid_size_points = Vector([Nx, Ny, Nz])
-grid_spacing_meters = Vector([dx, dx, dx])
-
-# create the k-space grid
-kgrid = kWaveGrid(grid_size_points, grid_spacing_meters)
+kgrid = kWaveGrid.from_domain(
+    dimensions=[axial_size, lateral_size, lateral_size], frequency=source_f0, sound_speed_min=c0, points_per_wavelength=ppw
+)
 
 # compute points per temporal period
 ppp: int = round(ppw / cfl)
@@ -124,8 +115,8 @@
 sensor = kSensor()
 
 # set sensor mask to record central plane, not including the source point
-sensor.mask = np.zeros((Nx, Ny, Nz), dtype=bool)
-sensor.mask[(source_x_offset + 1) : -1, :, Nz // 2] = True
+sensor.mask = np.zeros(kgrid.N, dtype=bool)
+sensor.mask[(source_x_offset + 1) : -1, :, kgrid.Nz // 2] = True
 
 # record the pressure
 sensor.record = ["p"]
@@ -155,10 +146,10 @@
 amp, _, _ = extract_amp_phase(sensor_data["p"].T, 1.0 / kgrid.dt, source_f0, dim=1, fft_padding=1, window="Rectangular")
 
 # reshape data
-amp = np.reshape(amp, (Nx - (source_x_offset + 2), Ny), order="F")
+amp = np.reshape(amp, (kgrid.Nx - (source_x_offset + 2), kgrid.Ny), order="F")
 
 # extract pressure on axis
-amp_on_axis = amp[:, Ny // 2]
+amp_on_axis = amp[:, kgrid.Ny // 2]
 
 # define axis vectors for plotting
 x_vec = kgrid.x_vec[(source_x_offset + 1) : -1, :] - kgrid.x_vec[source_x_offset]
@@ -172,7 +163,7 @@
 knumber = 2.0 * np.pi * source_f0 / c0
 
 # define axis
-x_max = Nx * dx
+x_max = kgrid.x_max
 delta_x = x_max / 10000.0
 x_ref = np.arange(0.0, x_max + delta_x, delta_x)
 
@@ -210,14 +201,14 @@
 ax2a.pcolormesh(
     1e3 * np.squeeze(kgrid.y_vec),
     1e3 * np.squeeze(kgrid.x_vec),
-    np.flip(source.p_mask[:, :, Nz // 2], axis=0),
+    np.flip(source.p_mask[:, :, kgrid.Nz // 2], axis=0),
     shading="nearest",
 )
 ax2a.set(xlabel="y [mm]", ylabel="x [mm]", title="Source Mask")
 ax2b.pcolormesh(
     1e3 * np.squeeze(kgrid.y_vec),
     1e3 * np.squeeze(kgrid.x_vec),
-    np.flip(grid_weights[:, :, Nz // 2], axis=0),
+    np.flip(grid_weights[:, :, kgrid.Nz // 2], axis=0),
     shading="nearest",
 )
 ax2b.set(xlabel="y [mm]", ylabel="x [mm]", title="Off-Grid Source Weights")

kwave/kgrid.py

@@ -37,6 +37,8 @@ def __init__(self, N, spacing):
         N, spacing = np.atleast_1d(N), np.atleast_1d(spacing)  # if inputs are lists
         assert N.ndim == 1 and spacing.ndim == 1  # ensure no multidimensional lists
         assert (1 <= N.size <= 3) and (1 <= spacing.size <= 3)  # ensure valid dimensionality
+        if spacing.size == 1:
+            spacing = spacing * np.ones(N.size)
         assert N.size == spacing.size, "Size list N and spacing list do not have the same size."
 
         self.N = N.astype(int)  #: grid size in each dimension [grid points]
@@ -699,3 +701,99 @@ def k_dtt(self, dtt_type):  # Not tested for correctness!
             # define product of implied period
             M = Mx * My * Mz
             return k, M
+
+    @classmethod
+    def from_geometry(cls, dimensions, min_element_width, points_per_wavelength=10):
+        """
+        Create a kWaveGrid based on domain dimensions and the smallest resolvable geometry element.
+
+        Args:
+            dimensions: List or array of physical domain sizes [m]
+            min_element_width: Width of the smallest resolvable geometry element [m]
+            points_per_wavelength: Number of points per wavelength (default=10)
+            cfl: CFL number (default=cls.CFL_DEFAULT)
+
+        Returns:
+            kWaveGrid instance with appropriate grid size and spacing
+        """
+        # Validate input parameters
+        dimensions = np.atleast_1d(dimensions)
+        if not np.all(dimensions > 0):
+            raise ValueError("Domain dimensions must be positive")
+        if not min_element_width > 0:
+            raise ValueError("Minimum element width must be positive")
+
+        # Calculate grid spacing based on minimum element width
+        # Ensure at least points_per_wavelength points across the smallest element
+        grid_spacing = min_element_width / points_per_wavelength
+
+        # Create a list of grid spacings with the same length as domain_size
+        grid_spacing_list = [grid_spacing] * dimensions.size
+
+        # Calculate grid size
+        N = np.ceil(dimensions / grid_spacing).astype(int)
+
+        # Create grid instance
+        grid = cls(N=N, spacing=grid_spacing_list)
+
+        # Note: Time parameters are left as "auto"
+        # The user can set them later using makeTime method
+
+        return grid
+
+    @classmethod
+    def from_domain(cls, dimensions, frequency, sound_speed_min, sound_speed_max=None, points_per_wavelength=10, cfl=None):
+        """
+        Create a kWaveGrid based on physical dimensions and acoustic properties.
+
+        Args:
+            dimensions: List or array of physical domain sizes [m]
+            frequency: Transmit frequency [Hz]
+            sound_speed_min: Minimum sound speed in the medium [m/s]
+            sound_speed_max: Maximum sound speed in the medium [m/s] (default=sound_speed_min)
+            points_per_wavelength: Number of points per wavelength (default=10)
+            cfl: CFL number (default=cls.CFL_DEFAULT)
+
+        Returns:
+            kWaveGrid instance with appropriate grid size, spacing, and time parameters
+        """
+        # Validate input parameters
+        dimensions = np.atleast_1d(dimensions)
+        if not np.all(dimensions > 0):
+            raise ValueError("Dimensions must be positive")
+        if not frequency > 0:
+            raise ValueError("Frequency must be positive")
+        if not sound_speed_min > 0:
+            raise ValueError("Sound speed must be positive")
+
+        # Use sound_speed_min for sound_speed_max if not provided
+        if sound_speed_max is None:
+            sound_speed_max = sound_speed_min
+        if not sound_speed_max > 0:
+            raise ValueError("Sound speed must be positive")
+
+        # Calculate wavelength
+        wavelength = sound_speed_min / frequency
+
+        # Calculate grid spacing
+        grid_spacing = wavelength / points_per_wavelength
+
+        # Calculate grid size
+        N = np.ceil(dimensions / grid_spacing).astype(int)
+
+        # Create grid instance
+        grid = cls(N=N, spacing=grid_spacing)
+
+        return grid
+
+    @property
+    def x_max(self):
+        return self.x_vec[-1] - self.x_vec[0]
+
+    @property
+    def y_max(self):
+        return self.y_vec[-1] - self.y_vec[0]
+
+    @property
+    def z_max(self):
+        return self.z_vec[-1] - self.z_vec[0]