Merge branch '26-mpi-geometry-objects' into mpi

gprMax/output_controllers/grid_view.py

@@ -98,9 +98,12 @@ class GridView(Generic[GridType]):
     def nz(self) -> int:
         return self.size[2]
 
-    def get_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice:
+    def getter_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice:
         """Create a slice object for the specified dimension.
 
+        This is used to slice and get a view of arrays owned by the
+        grid.
+
         Args:
             dimension: Dimension to create the slice object for. Values
                 0, 1, and 2 map to the x, y, and z dimensions
@@ -119,7 +122,150 @@ class GridView(Generic[GridType]):
 
         return slice(self.start[dimension], stop, self.step[dimension])
 
-    def slice_array(self, array: npt.NDArray, upper_bound_exclusive: bool = True) -> npt.NDArray:
+    def setter_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice:
+        """Create a slice object for the specified dimension.
+
+        This is used to slice arrays owned by the grid in order to set
+        their value.
+
+        Args:
+            dimension: Dimension to create the slice object for. Values
+                0, 1, and 2 map to the x, y, and z dimensions
+                respectively.
+            upper_bound_exclusive: Optionally specify if the upper bound
+                of the slice should be exclusive or inclusive. Defaults
+                to True.
+
+        Returns:
+            slice: Slice object
+        """
+        return self.getter_slice(dimension, upper_bound_exclusive)
+
+    def get_output_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice:
+        """Create an output slice object for the specified dimension.
+
+        This provides a slice of the grid view for the section of the
+        grid view managed by this process. This can be used when writing
+        out arrays provided by the grid view as part of a collective
+        operation.
+
+        For example:
+        ```
+        dset_slice = (
+                grid_view.get_output_slice(0),
+                grid_view.get_output_slice(1),
+                grid_view.get_output_slice(2),
+            )
+
+        dset[dset_slice] = grid_view.get_solid()
+        ```
+
+        Args:
+            dimension: Dimension to create the slice object for. Values
+                0, 1, and 2 map to the x, y, and z dimensions
+                respectively.
+            upper_bound_exclusive: Optionally specify if the upper bound
+                of the slice should be exclusive or inclusive. Defaults
+                to True.
+
+        Returns:
+            slice: Slice object
+        """
+        if upper_bound_exclusive:
+            size = self.size[dimension]
+        else:
+            size = self.size[dimension] + 1
+
+        return slice(0, size)
+
+    def get_3d_output_slice(self, upper_bound_exclusive: bool = True) -> Tuple[slice, slice, slice]:
+        """Create a 3D output slice object.
+
+        This provides a slice of the grid view for the section of the
+        grid view managed by this process. This can be used when writing
+        out arrays provided by the grid view as part of a collective
+        operation.
+
+        For example:
+        `dset[grid_view.get_3d_output_slice()] = grid_view.get_solid()`
+
+        Args:
+            upper_bound_exclusive: Optionally specify if the upper bound
+                of the slice should be exclusive or inclusive. Defaults
+                to True.
+
+        Returns:
+            slice: 3D Slice object
+        """
+        return (
+            self.get_output_slice(0, upper_bound_exclusive),
+            self.get_output_slice(1, upper_bound_exclusive),
+            self.get_output_slice(2, upper_bound_exclusive),
+        )
+
+    def get_read_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice:
+        """Create a read slice object for the specified dimension.
+
+        This provides a slice of the grid view for the section of the
+        grid view managed by this rank. This can be used when reading
+        arrays provided by the grid view as part of a collective
+        operation.
+
+        For example:
+        ```
+        dset_slice = (
+                grid_view.get_read_slice(0),
+                grid_view.get_read_slice(1),
+                grid_view.get_read_slice(2),
+            )
+
+        grid_view.set_solid(dset[dset_slice])
+        ```
+
+        Args:
+            dimension: Dimension to create the slice object for. Values
+                0, 1, and 2 map to the x, y, and z dimensions
+                respectively.
+            upper_bound_exclusive: Optionally specify if the upper bound
+                of the slice should be exclusive or inclusive. Defaults
+                to True.
+
+        Returns:
+            slice: Slice object
+        """
+        return self.get_output_slice(dimension, upper_bound_exclusive)
+
+    def get_3d_read_slice(self, upper_bound_exclusive: bool = True) -> Tuple[slice, slice, slice]:
+        """Create a 3D read slice object.
+
+        This provides a slice of the grid view for the section of the
+        grid view managed by this rank. This can be used when reading
+        arrays provided by the grid view as part of a collective
+        operation.
+
+        For example:
+        ```
+        solid = dset[grid_view.get_3d_read_slice()]
+        grid_view.set_solid(solid)
+        ```
+
+        Args:
+            upper_bound_exclusive: Optionally specify if the upper bound
+                of the slice should be exclusive or inclusive. Defaults
+                to True.
+
+        Returns:
+            slice: 3D Slice object
+        """
+        return (
+            self.get_read_slice(0, upper_bound_exclusive),
+            self.get_read_slice(1, upper_bound_exclusive),
+            self.get_read_slice(2, upper_bound_exclusive),
+        )
+
+    def get_array_slice(
+        self, array: npt.NDArray, upper_bound_exclusive: bool = True
+    ) -> npt.NDArray:
         """Slice an array according to the dimensions of the grid view.
 
         It is assumed the last 3 dimensions of the provided array
@@ -142,12 +288,41 @@ class GridView(Generic[GridType]):
         return np.ascontiguousarray(
             array[
                 ...,
-                self.get_slice(0, upper_bound_exclusive),
-                self.get_slice(1, upper_bound_exclusive),
-                self.get_slice(2, upper_bound_exclusive),
+                self.getter_slice(0, upper_bound_exclusive),
+                self.getter_slice(1, upper_bound_exclusive),
+                self.getter_slice(2, upper_bound_exclusive),
             ]
         )
 
+    def set_array_slice(
+        self, array: npt.NDArray, value: npt.NDArray, upper_bound_exclusive: bool = True
+    ):
+        """Set value of an array according to the dimensions of the grid view.
+
+        It is assumed the last 3 dimensions of the array represent the
+        x, y, z spacial information. Other dimensions will not be
+        sliced.
+
+        E.g. If setting the value of an array of shape (10, 100, 50, 50)
+        the new values should have shape (10, x, y, z) where x, y, and z
+        are specified by the size/shape of the grid view.
+
+        Args:
+            array: Array to set the values of. Must have at least 3
+                dimensions.
+            value: New values. Its shape must match 'array' after
+                'array' has been sliced.
+            upper_bound_exclusive: Optionally specify if the upper bound
+                of the slice should be exclusive or inclusive. Defaults
+                to True.
+        """
+        array[
+            ...,
+            self.setter_slice(0, upper_bound_exclusive),
+            self.setter_slice(1, upper_bound_exclusive),
+            self.setter_slice(2, upper_bound_exclusive),
+        ] = value
+
     def initialise_materials(self, filter_materials: bool = True):
         """Create a new ID map for materials in the grid view.
 
@@ -209,32 +384,64 @@ class GridView(Generic[GridType]):
             ID: View of the ID array.
         """
         if self._ID is None or force_refresh:
-            self._ID = self.slice_array(self.grid.ID, upper_bound_exclusive=False)
+            self._ID = self.get_array_slice(self.grid.ID, upper_bound_exclusive=False)
         return self._ID
 
+    def set_ID(self, value: npt.NDArray[np.uint32]):
+        """Set the value of the ID array.
+
+        Args:
+            value: Array of new values.
+        """
+        self.set_array_slice(self.grid.ID, value, upper_bound_exclusive=False)
+
     def get_solid(self) -> npt.NDArray[np.uint32]:
         """Get a view of the solid array.
 
         Returns:
-            solid: View of the solid array
+            solid: View of the solid array.
         """
-        return self.slice_array(self.grid.solid)
+        return self.get_array_slice(self.grid.solid)
+
+    def set_solid(self, value: npt.NDArray[np.uint32]):
+        """Set the value of the solid array.
+
+        Args:
+            value: Array of new values.
+        """
+        self.set_array_slice(self.grid.solid, value)
 
     def get_rigidE(self) -> npt.NDArray[np.int8]:
         """Get a view of the rigidE array.
 
         Returns:
-            rigidE: View of the rigidE array
+            rigidE: View of the rigidE array.
         """
-        return self.slice_array(self.grid.rigidE)
+        return self.get_array_slice(self.grid.rigidE)
+
+    def set_rigidE(self, value: npt.NDArray[np.uint32]):
+        """Set the value of the rigidE array.
+
+        Args:
+            value: Array of new values.
+        """
+        self.set_array_slice(self.grid.rigidE, value)
 
     def get_rigidH(self) -> npt.NDArray[np.int8]:
         """Get a view of the rigidH array.
 
         Returns:
-            rigidH: View of the rigidH array
+            rigidH: View of the rigidH array.
         """
-        return self.slice_array(self.grid.rigidH)
+        return self.get_array_slice(self.grid.rigidH)
+
+    def set_rigidH(self, value: npt.NDArray[np.uint32]):
+        """Set the value of the rigidH array.
+
+        Args:
+            value: Array of new values.
+        """
+        self.set_array_slice(self.grid.rigidH, value)
 
     def get_Ex(self) -> npt.NDArray[np.float32]:
         """Get a view of the Ex array.
@@ -242,7 +449,7 @@ class GridView(Generic[GridType]):
         Returns:
             Ex: View of the Ex array
         """
-        return self.slice_array(self.grid.Ex, upper_bound_exclusive=False)
+        return self.get_array_slice(self.grid.Ex, upper_bound_exclusive=False)
 
     def get_Ey(self) -> npt.NDArray[np.float32]:
         """Get a view of the Ey array.
@@ -250,7 +457,7 @@ class GridView(Generic[GridType]):
         Returns:
             Ey: View of the Ey array
         """
-        return self.slice_array(self.grid.Ey, upper_bound_exclusive=False)
+        return self.get_array_slice(self.grid.Ey, upper_bound_exclusive=False)
 
     def get_Ez(self) -> npt.NDArray[np.float32]:
         """Get a view of the Ez array.
@@ -258,7 +465,7 @@ class GridView(Generic[GridType]):
         Returns:
             Ez: View of the Ez array
         """
-        return self.slice_array(self.grid.Ez, upper_bound_exclusive=False)
+        return self.get_array_slice(self.grid.Ez, upper_bound_exclusive=False)
 
     def get_Hx(self) -> npt.NDArray[np.float32]:
         """Get a view of the Hx array.
@@ -266,7 +473,7 @@ class GridView(Generic[GridType]):
         Returns:
             Hx: View of the Hx array
         """
-        return self.slice_array(self.grid.Hx, upper_bound_exclusive=False)
+        return self.get_array_slice(self.grid.Hx, upper_bound_exclusive=False)
 
     def get_Hy(self) -> npt.NDArray[np.float32]:
         """Get a view of the Hy array.
@@ -274,7 +481,7 @@ class GridView(Generic[GridType]):
         Returns:
             Hy: View of the Hy array
         """
-        return self.slice_array(self.grid.Hy, upper_bound_exclusive=False)
+        return self.get_array_slice(self.grid.Hy, upper_bound_exclusive=False)
 
     def get_Hz(self) -> npt.NDArray[np.float32]:
         """Get a view of the Hz array.
@@ -282,7 +489,7 @@ class GridView(Generic[GridType]):
         Returns:
             Hz: View of the Hz array
         """
-        return self.slice_array(self.grid.Hz, upper_bound_exclusive=False)
+        return self.get_array_slice(self.grid.Hz, upper_bound_exclusive=False)
 
 
 class MPIGridView(GridView[MPIGrid]):
@@ -386,9 +593,12 @@ class MPIGridView(GridView[MPIGrid]):
     def gz(self) -> int:
         return self.global_size[2]
 
-    def get_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice:
+    def getter_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice:
         """Create a slice object for the specified dimension.
 
+        This is used to slice and get a view of arrays owned by the
+        grid.
+
         Args:
             dimension: Dimension to create the slice object for. Values
                 0, 1, and 2 map to the x, y, and z dimensions
@@ -409,6 +619,35 @@ class MPIGridView(GridView[MPIGrid]):
 
         return slice(self.start[dimension], stop, self.step[dimension])
 
+    def setter_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice:
+        """Create a slice object for the specified dimension.
+
+        This is used to slice arrays owned by the grid in order to set
+        their value.
+
+        Args:
+            dimension: Dimension to create the slice object for. Values
+                0, 1, and 2 map to the x, y, and z dimensions
+                respectively.
+            upper_bound_exclusive: Optionally specify if the upper bound
+                of the slice should be exclusive or inclusive. Defaults
+                to True.
+
+        Returns:
+            slice: Slice object
+        """
+        if upper_bound_exclusive:
+            stop = self.stop[dimension]
+        else:
+            stop = self.stop[dimension] + self.step[dimension]
+
+        if self.has_negative_neighbour[dimension]:
+            start = self.start[dimension] - self.step[dimension]
+        else:
+            start = self.start[dimension]
+
+        return slice(start, stop, self.step[dimension])
+
     def get_output_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice:
         """Create an output slice object for the specified dimension.
 
@@ -450,30 +689,48 @@ class MPIGridView(GridView[MPIGrid]):
 
         return slice(offset, offset + size)
 
-    def get_3d_output_slice(self, upper_bound_exclusive: bool = True) -> Tuple[slice, slice, slice]:
-        """Create a 3D output slice object.
+    def get_read_slice(self, dimension: int, upper_bound_exclusive: bool = True) -> slice:
+        """Create a read slice object for the specified dimension.
 
         This provides a slice of the grid view for the section of the
-        grid view managed by this rank. This can be used when writing
-        out arrays provided by the grid view as part of a collective
+        grid view managed by this rank. This can be used when reading
+        arrays provided by the grid view as part of a collective
         operation.
 
         For example:
-        `dset[grid_view.get_3d_output_slice()] = grid_view.get_solid()`
+        ```
+        dset_slice = (
+                grid_view.get_read_slice(0),
+                grid_view.get_read_slice(1),
+                grid_view.get_read_slice(2),
+            )
+
+        grid_view.get_solid()[:] = dset[dset_slice]
+        ```
 
         Args:
+            dimension: Dimension to create the slice object for. Values
+                0, 1, and 2 map to the x, y, and z dimensions
+                respectively.
             upper_bound_exclusive: Optionally specify if the upper bound
                 of the slice should be exclusive or inclusive. Defaults
                 to True.
 
         Returns:
-            slice: 3D Slice object
+            slice: Slice object
         """
-        return (
-            self.get_output_slice(0, upper_bound_exclusive),
-            self.get_output_slice(1, upper_bound_exclusive),
-            self.get_output_slice(2, upper_bound_exclusive),
-        )
+        if upper_bound_exclusive:
+            size = self.size[dimension]
+        else:
+            size = self.size[dimension] + 1
+
+        offset = self.offset[dimension] // self.step[dimension]
+
+        if self.has_negative_neighbour[dimension]:
+            offset -= 1
+            size += 1
+
+        return slice(offset, offset + size)
 
     def initialise_materials(self, filter_materials: bool = True):
         """Create a new ID map for materials in the grid view.

gprMax/output_controllers/read_geometry_object.py

@@ -0,0 +1,178 @@
+from contextlib import AbstractContextManager
+from os import PathLike
+from types import TracebackType
+from typing import Optional
+
+import h5py
+import numpy as np
+import numpy.typing as npt
+from mpi4py import MPI
+
+from gprMax.grid.fdtd_grid import FDTDGrid
+from gprMax.grid.mpi_grid import MPIGrid
+from gprMax.output_controllers.grid_view import GridView, MPIGridView
+
+
+def create_read_geometry_object(
+    filename: PathLike,
+    grid: FDTDGrid,
+    start: npt.NDArray[np.int32],
+    num_existing_materials: int,
+):
+    if isinstance(grid, MPIGrid) and not grid.local_bounds_overlap_grid(start, stop):
+        # The MPIGridView created by the ReadGeometryObject will
+        # create a new communicator using MPI_Split. Calling this
+        # here prevents deadlock if not all ranks create the new
+        # ReadGeometryObject.
+        grid.comm.Split(MPI.UNDEFINED)
+        return None
+    else:
+        return ReadGeometryObject(
+            filename,
+            grid,
+            start[0],
+            start[1],
+            start[2],
+            stop[0],
+            stop[1],
+            stop[2],
+            num_existing_materials,
+        )
+
+
+class ReadGeometryObject(AbstractContextManager):
+    def __init__(
+        self,
+        filename: PathLike,
+        grid: FDTDGrid,
+        start: npt.NDArray[np.int32],
+        num_existing_materials: int,
+    ) -> None:
+        self.file_handler = h5py.File(filename)
+
+        data = self.file_handler["/data"]
+        assert isinstance(data, h5py.Dataset)
+        stop = start + data.shape
+
+        if isinstance(grid, MPIGrid):
+            if grid.local_bounds_overlap_grid(start, stop):
+                self.grid_view = MPIGridView(
+                    grid, start[0], start[1], start[2], stop[0], stop[1], stop[2]
+                )
+            else:
+                # The MPIGridView will create a new communicator using
+                # MPI_Split. Calling this here prevents deadlock if not
+                # all ranks need to read the geometry object.
+                grid.comm.Split(MPI.UNDEFINED)
+                self.grid_view = None
+
+        else:
+            self.grid_view = GridView(grid, start[0], start[1], start[2], stop[0], stop[1], stop[2])
+
+        self.num_existing_materials = num_existing_materials
+
+    def __enter__(self):
+        return self
+
+    def __exit__(
+        self,
+        exc_type: Optional[type[BaseException]],
+        exc_value: Optional[BaseException],
+        traceback: Optional[TracebackType],
+    ) -> Optional[bool]:
+        """Close the file when the context is exited.
+
+        The parameters describe the exception that caused the context to
+        be exited. If the context was exited without an exception, all
+        three arguments will be None. Any exception will be
+        processed normally upon exit from this method.
+
+        Returns:
+            suppress_exception (optional): Returns True if the exception
+                should be suppressed (i.e. not propagated). Otherwise,
+                the exception will be processed normally upon exit from
+                this method.
+        """
+        self.close()
+
+    def close(self) -> None:
+        """Close the file handler"""
+        self.file_handler.close()
+
+    def has_valid_discritisation(self) -> bool:
+        if self.grid_view is None:
+            return True
+
+        dx_dy_dz = self.file_handler.attrs["dx_dy_dz"]
+        return not isinstance(dx_dy_dz, h5py.Empty) and all(dx_dy_dz == self.grid_view.grid.dl)
+
+    def has_ID_array(self) -> bool:
+        ID_class = self.file_handler.get("ID", getclass=True)
+        return ID_class == h5py.Dataset
+
+    def has_rigid_arrays(self) -> bool:
+        rigidE_class = self.file_handler.get("rigidE", getclass=True)
+        rigidH_class = self.file_handler.get("rigidH", getclass=True)
+        return rigidE_class == h5py.Dataset and rigidH_class == h5py.Dataset
+
+    def read_data(self):
+        if self.grid_view is None:
+            return
+
+        data = self.file_handler["/data"]
+        assert isinstance(data, h5py.Dataset)
+        data = data[self.grid_view.get_3d_read_slice()]
+
+        # Should be int16 to allow for -1 which indicates background, i.e.
+        # don't build anything, but AustinMan/Woman maybe uint16
+        if data.dtype != "int16":
+            data = data.astype("int16")
+
+        self.grid_view.set_solid(data + self.num_existing_materials)
+
+    def get_data(self) -> Optional[npt.NDArray[np.int16]]:
+        if self.grid_view is None:
+            return None
+
+        data = self.file_handler["/data"]
+        assert isinstance(data, h5py.Dataset)
+        data = data[self.grid_view.get_3d_read_slice()]
+
+        # Should be int16 to allow for -1 which indicates background, i.e.
+        # don't build anything, but AustinMan/Woman maybe uint16
+        if data.dtype != "int16":
+            data = data.astype("int16")
+
+        return data + self.num_existing_materials
+
+    def read_rigidE(self):
+        if self.grid_view is None:
+            return
+
+        rigidE = self.file_handler["/rigidE"]
+        assert isinstance(rigidE, h5py.Dataset)
+
+        dset_slice = self.grid_view.get_3d_read_slice()
+        self.grid_view.set_rigidE(rigidE[:, dset_slice[0], dset_slice[1], dset_slice[2]])
+
+    def read_rigidH(self):
+        if self.grid_view is None:
+            return
+
+        rigidH = self.file_handler["/rigidH"]
+        assert isinstance(rigidH, h5py.Dataset)
+
+        dset_slice = self.grid_view.get_3d_read_slice()
+        self.grid_view.set_rigidH(rigidH[:, dset_slice[0], dset_slice[1], dset_slice[2]])
+
+    def read_ID(self):
+        if self.grid_view is None:
+            return
+
+        ID = self.file_handler["/ID"]
+        assert isinstance(ID, h5py.Dataset)
+
+        dset_slice = self.grid_view.get_3d_read_slice(upper_bound_exclusive=False)
+        self.grid_view.set_ID(
+            ID[:, dset_slice[0], dset_slice[1], dset_slice[2]] + self.num_existing_materials
+        )

gprMax/user_objects/cmds_geometry/geometry_objects_read.py

@@ -25,8 +25,8 @@ import gprMax.config as config
 from gprMax.cython.geometry_primitives import build_voxels_from_array
 from gprMax.grid.fdtd_grid import FDTDGrid
 from gprMax.hash_cmds_file import get_user_objects
+from gprMax.output_controllers.read_geometry_object import ReadGeometryObject
 from gprMax.user_objects.user_objects import GeometryUserObject
-from gprMax.utilities.utilities import round_value
 
 logger = logging.getLogger(__name__)
 
@@ -49,12 +49,6 @@ class GeometryObjectsRead(GeometryUserObject):
             logger.exception(f"{self.__str__()} requires exactly five parameters")
             raise
 
-        # Discretise the point using uip object. This has different behaviour
-        # depending on the type of uip object. So we can use it for
-        # the main grid or the subgrid.
-        uip = self._create_uip(grid)
-        xs, ys, zs = uip.discretise_point(p1)
-
         # See if material file exists at specified path and if not try input
         # file directory
         matfile = Path(matfile)
@@ -75,9 +69,21 @@ class GeometryObjectsRead(GeometryUserObject):
                 if (line.startswith("#") and not line.startswith("##") and line.rstrip("\n"))
             ]
 
+        # Avoid redefining default builtin materials
+        pec = f"#material: 1 inf 1 0 pec{{{matstr}}}\n"
+        free_space = f"#material: 1 0 1 0 free_space{{{matstr}}}\n"
+        if materials[0] == pec and materials[1] == free_space:
+            materials.pop(0)
+            materials.pop(1)
+            numexistmaterials -= 2
+        elif materials[0] == pec or materials[0] == free_space:
+            materials.pop(0)
+            numexistmaterials -= 1
+
         # Build scene
         # API for multiple scenes / model runs
         scene = config.get_model_config().get_scene()
+        assert scene is not None
         material_objs = get_user_objects(materials, checkessential=False)
         for material_obj in material_objs:
             scene.add(material_obj)
@@ -99,69 +105,53 @@ class GeometryObjectsRead(GeometryUserObject):
         if not geofile.exists():
             geofile = Path(config.sim_config.input_file_path.parent, geofile)
 
-        # Open geometry object file and read/check spatial resolution attribute
-        f = h5py.File(geofile, "r")
-        dx_dy_dz = f.attrs["dx_dy_dz"]
-        if round_value(
-            (dx_dy_dz[0] / grid.dx) != 1
-            or round_value(dx_dy_dz[1] / grid.dy) != 1
-            or round_value(dx_dy_dz[2] / grid.dz) != 1
-        ):
-            logger.exception(
-                f"{self.__str__()} requires the spatial resolution "
-                "of the geometry objects file to match the spatial "
-                "resolution of the model"
-            )
-            raise ValueError
+        # Discretise the point using uip object. This has different behaviour
+        # depending on the type of uip object. So we can use it for
+        # the main grid, MPI grids or the subgrid.
+        uip = self._create_uip(grid)
+        discretised_p1 = uip.discretise_point(p1)
+        p2 = uip.round_to_grid_static_point(p1)
 
-        data = f["/data"][:]
+        with ReadGeometryObject(geofile, grid, discretised_p1, numexistmaterials) as f:
+            # Check spatial resolution attribute
+            if not f.has_valid_discritisation():
+                raise ValueError(
+                    f"{self.__str__()} requires the spatial resolution "
+                    "of the geometry objects file to match the spatial "
+                    "resolution of the model"
+                )
 
-        # Should be int16 to allow for -1 which indicates background, i.e.
-        # don't build anything, but AustinMan/Woman maybe uint16
-        if data.dtype != "int16":
-            data = data.astype("int16")
+            if f.has_rigid_arrays() and f.has_ID_array():
+                f.read_data()
+                f.read_ID()
+                f.read_rigidE()
+                f.read_rigidH()
 
-        # Look to see if rigid and ID arrays are present (these should be
-        # present if the original geometry objects were written from gprMax)
-        try:
-            rigidE = f["/rigidE"][:]
-            rigidH = f["/rigidH"][:]
-            ID = f["/ID"][:]
-            grid.solid[
-                xs : xs + data.shape[0], ys : ys + data.shape[1], zs : zs + data.shape[2]
-            ] = (data + numexistmaterials)
-            grid.rigidE[
-                :, xs : xs + rigidE.shape[1], ys : ys + rigidE.shape[2], zs : zs + rigidE.shape[3]
-            ] = rigidE
-            grid.rigidH[
-                :, xs : xs + rigidH.shape[1], ys : ys + rigidH.shape[2], zs : zs + rigidH.shape[3]
-            ] = rigidH
-            grid.ID[:, xs : xs + ID.shape[1], ys : ys + ID.shape[2], zs : zs + ID.shape[3]] = (
-                ID + numexistmaterials
-            )
-            logger.info(
-                f"{self.grid_name(grid)}Geometry objects from file {geofile} "
-                f"inserted at {xs * grid.dx:g}m, {ys * grid.dy:g}m, "
-                f"{zs * grid.dz:g}m, with corresponding materials file "
-                f"{matfile}."
-            )
-        except KeyError:
-            averaging = False
-            build_voxels_from_array(
-                xs,
-                ys,
-                zs,
-                numexistmaterials,
-                averaging,
-                data,
-                grid.solid,
-                grid.rigidE,
-                grid.rigidH,
-                grid.ID,
-            )
-            logger.info(
-                f"{self.grid_name(grid)}Geometry objects from file "
-                f"(voxels only){geofile} inserted at {xs * grid.dx:g}m, "
-                f"{ys * grid.dy:g}m, {zs * grid.dz:g}m, with corresponding "
-                f"materials file {matfile}."
-            )
+                logger.info(
+                    f"{self.grid_name(grid)}Geometry objects from file {geofile}"
+                    f" inserted at {p2[0]:g}m, {p2[1]:g}m, {p2[2]:g}m,"
+                    f" with corresponding materials file"
+                    f" {matfile}."
+                )
+            else:
+                data = f.get_data()
+                if data is not None:
+                    averaging = False
+                    build_voxels_from_array(
+                        discretised_p1[0],
+                        discretised_p1[1],
+                        discretised_p1[2],
+                        numexistmaterials,
+                        averaging,
+                        data,
+                        grid.solid,
+                        grid.rigidE,
+                        grid.rigidH,
+                        grid.ID,
+                    )
+                logger.info(
+                    f"{self.grid_name(grid)}Geometry objects from file "
+                    f"(voxels only){geofile} inserted at {p2[0]:g}m, "
+                    f"{p2[1]:g}m, {p2[2]:g}m, with corresponding "
+                    f"materials file {matfile}."
+                )

reframe_tests/tests/base_tests.py

@@ -161,9 +161,11 @@ class GprMaxBaseTest(RunOnlyRegressionTest):
         if self.test_dependency is None:
             return None
 
-        # Always filter by the model parameter, but allow child classes
+        # Always filter by the model parameter (unless the test
+        # dependency only runs a single model), but allow child classes
         # (or mixins) to override how models are filtered.
-        kwargs.setdefault("model", self.model)
+        if len(self.test_dependency.model.values) > 1:
+            kwargs.setdefault("model", self.model)
 
         variant_nums = self.test_dependency.get_variant_nums(**kwargs)
 

reframe_tests/tests/mixins.py

@@ -1,10 +1,12 @@
 from pathlib import Path
+from shutil import copyfile
 from typing import Optional
 
+import reframe.utility.sanity as sn
 import reframe.utility.typecheck as typ
 from numpy import prod
 from reframe import RegressionMixin
-from reframe.core.builtins import parameter, required, run_after, variable
+from reframe.core.builtins import parameter, run_after, variable
 from typing_extensions import TYPE_CHECKING
 
 from reframe_tests.tests.base_tests import GprMaxBaseTest
@@ -31,7 +33,11 @@ class ReceiverMixin(GprMaxMixin):
 
     @run_after("setup", always_last=True)
     def add_receiver_regression_checks(self):
-        reference_file = self.build_reference_filepath(self.output_file)
+        test_dependency = self.get_test_dependency()
+        if test_dependency is not None:
+            reference_file = self.build_reference_filepath(test_dependency.output_file)
+        else:
+            reference_file = self.build_reference_filepath(self.output_file)
 
         if self.number_of_receivers > 0:
             for i in range(self.number_of_receivers):
@@ -88,16 +94,7 @@ class GeometryOnlyMixin(GprMaxMixin):
         self.executable_opts += ["--geometry-only"]
 
 
-class GeometryObjectMixin(GprMaxMixin):
-    """Add regression tests for geometry objects.
-
-    Attributes:
-        geometry_objects (list[str]): List of geometry objects to run
-            regression checks on.
-    """
-
-    geometry_objects = variable(typ.List[str], value=[])
-
+class GeometryObjectMixinBase(GprMaxMixin):
     def build_geometry_object_filepath(self, geometry_object: str) -> Path:
         """Build filepath to the specified geometry object.
 
@@ -114,6 +111,63 @@ class GeometryObjectMixin(GprMaxMixin):
         """
         return Path(f"{geometry_object}_materials").with_suffix(".txt")
 
+
+class GeometryObjectsReadMixin(GeometryObjectMixinBase):
+    geometry_objects_read = variable(typ.Dict[str, str], value={})
+
+    @run_after("setup")
+    def copy_geometry_objects_from_test_dependency(self):
+        self.skip_if(
+            len(self.geometry_objects_read) < 0,
+            f"Must provide a list of geometry objects being read by the test.",
+        )
+
+        test_dependency = self.get_test_dependency()
+
+        self.skip_if(
+            test_dependency is None,
+            f"GeometryObjectsReadMixin must be used with a test dependency.",
+        )
+
+        for geometry_object_input, geometry_object_output in self.geometry_objects_read.items():
+            geometry_object_input_file = self.build_geometry_object_filepath(geometry_object_input)
+            geometry_object_input_file = Path(test_dependency.stagedir, geometry_object_input_file)
+
+            materials_input_file = self.build_materials_filepath(geometry_object_input)
+            materials_input_file = Path(test_dependency.stagedir, materials_input_file)
+
+            self.skip_if(
+                not sn.path_isfile(geometry_object_input_file),
+                f"Test dependency did not create geometry object file.",
+            )
+
+            self.skip_if(
+                not sn.path_isfile(materials_input_file),
+                f"Test dependency did not create geometry object materials file.",
+            )
+
+            geometry_object_output_file = self.build_geometry_object_filepath(
+                geometry_object_output
+            )
+            geometry_object_output_file = Path(self.stagedir, geometry_object_output_file)
+
+            materials_output_file = self.build_materials_filepath(geometry_object_output)
+            materials_output_file = Path(self.stagedir, materials_output_file)
+
+            copyfile(geometry_object_input_file, geometry_object_output_file)
+            copyfile(materials_input_file, materials_output_file)
+
+
+class GeometryObjectsWriteMixin(GeometryObjectMixinBase):
+    """Add regression tests for geometry objects.
+
+    Attributes:
+        geometry_objects_write (list[str]): List of geometry objects to
+            run regression checks on.
+    """
+
+    geometry_objects_write = variable(typ.List[str], value=[])
+
     @run_after("setup")
     def add_geometry_object_regression_checks(self):
         """Add a regression check for each geometry object.
@@ -121,11 +175,11 @@ class GeometryObjectMixin(GprMaxMixin):
         The test will be skipped if no geometry objects have been specified.
         """
         self.skip_if(
-            len(self.geometry_objects) < 0,
-            f"Must provide a list of geometry objects.",
+            len(self.geometry_objects_write) < 0,
+            f"Must provide a list of geometry objects being written by the test.",
         )
 
-        for geometry_object in self.geometry_objects:
+        for geometry_object in self.geometry_objects_write:
             # Add materials regression check first as if this fails,
             # checking the .h5 file will almost definitely fail.
             materials_file = self.build_materials_filepath(geometry_object)
@@ -138,7 +192,9 @@ class GeometryObjectMixin(GprMaxMixin):
             self.regression_checks.append(materials_regression_check)
 
             geometry_object_file = self.build_geometry_object_filepath(geometry_object)
-            reference_file = self.build_reference_filepath(geometry_object)
+            reference_file = self.build_reference_filepath(
+                geometry_object, suffix=geometry_object_file.suffix
+            )
             regression_check = GeometryObjectRegressionCheck(geometry_object_file, reference_file)
             self.regression_checks.append(regression_check)
 

reframe_tests/tests/src/geometry_object_tests/geometry_object_move.in

@@ -0,0 +1,12 @@
+#title: Hertzian dipole in free-space
+#domain: 0.100 0.100 0.100
+#dx_dy_dz: 0.001 0.001 0.001
+#time_window: 3e-9
+
+#waveform: gaussiandot 1 1e9 myWave
+#hertzian_dipole: z 0.050 0.050 0.050 myWave
+
+#rx: 0.08 0.08 0.08
+
+#geometry_objects_read: 0.03 0.03 0.03 full_volume_read.h5 full_volume_read_materials.txt
+#geometry_objects_write: 0.05 0.05 0.05 0.09 0.09 0.09 partial_volume

reframe_tests/tests/src/geometry_object_tests/geometry_object_read_full_volume.in

@@ -0,0 +1,11 @@
+#title: Hertzian dipole in free-space
+#domain: 0.100 0.100 0.100
+#dx_dy_dz: 0.001 0.001 0.001
+#time_window: 3e-9
+
+#waveform: gaussiandot 1 1e9 myWave
+#hertzian_dipole: z 0.050 0.050 0.050 myWave
+
+#rx: 0.08 0.08 0.08
+
+#geometry_objects_read: 0 0 0 full_volume_read.h5 full_volume_read_materials.txt

reframe_tests/tests/src/geometry_object_tests/geometry_object_read_write.in

@@ -0,0 +1,13 @@
+#title: Hertzian dipole in free-space
+#domain: 0.100 0.100 0.100
+#dx_dy_dz: 0.001 0.001 0.001
+#time_window: 3e-9
+
+#waveform: gaussiandot 1 1e9 myWave
+#hertzian_dipole: z 0.050 0.050 0.050 myWave
+
+#rx: 0.08 0.08 0.08
+
+#geometry_objects_read: 0 0 0 full_volume_read.h5 full_volume_read_materials.txt
+#geometry_objects_write: 0.02 0.02 0.02 0.06 0.06 0.06 partial_volume
+#geometry_objects_write:  0 0 0 0.1 0.1 0.1 full_volume

reframe_tests/tests/src/geometry_object_tests/geometry_object_write.in

@@ -6,6 +6,8 @@
 #waveform: gaussiandot 1 1e9 myWave
 #hertzian_dipole: z 0.050 0.050 0.050 myWave
 
+#rx: 0.08 0.08 0.08
+
 #material: 4.9 0 1 0 myWater
 #material: 2 0 1.4 0 unusedMaterial
 #material: 3 0 2 0 boxMaterial

reframe_tests/tests/standard_tests.py

@@ -1,11 +1,17 @@
+from pathlib import Path
+
+from reframe.core.builtins import run_before
+
 from reframe_tests.tests.base_tests import GprMaxBaseTest
 from reframe_tests.tests.mixins import (
-    GeometryObjectMixin,
+    GeometryObjectsReadMixin,
+    GeometryObjectsWriteMixin,
     GeometryOnlyMixin,
     GeometryViewMixin,
     ReceiverMixin,
     SnapshotMixin,
 )
+from reframe_tests.tests.regression_checks import GeometryObjectMaterialsRegressionCheck
 
 
 class GprMaxRegressionTest(ReceiverMixin, GprMaxBaseTest):
@@ -20,9 +26,35 @@ class GprMaxGeometryViewTest(GeometryViewMixin, GeometryOnlyMixin, GprMaxBaseTes
     pass
 
 
-class GprMaxGeometryObjectTest(GeometryObjectMixin, GeometryOnlyMixin, GprMaxBaseTest):
+class GprMaxGeometryObjectsWriteTest(GeometryObjectsWriteMixin, GprMaxBaseTest):
     pass
 
 
-class GprMaxGeometryTest(GeometryObjectMixin, ReceiverMixin, GprMaxBaseTest):
-    geometry_objects = ["full_volume"]
+class GprMaxGeometryObjectsReadTest(GeometryObjectsReadMixin, GprMaxBaseTest):
+    pass
+
+
+class GprMaxGeometryObjectsReadWriteTest(
+    GeometryObjectsReadMixin, GeometryObjectsWriteMixin, GprMaxBaseTest
+):
+    @run_before("sanity")
+    def update_material_files(self):
+        checks = [
+            check
+            for check in self.regression_checks
+            if isinstance(check, GeometryObjectMaterialsRegressionCheck)
+        ]
+        for check in checks:
+            for geometry_object in self.geometry_objects_read.values():
+                material_file = Path(self.stagedir, check.output_file)
+                with open(material_file, "r") as f:
+                    lines = f.readlines()
+
+                with open(material_file, "w") as f:
+                    for line in lines:
+                        new_line = line.replace(f"{{{geometry_object}_materials}}", "")
+                        f.write(new_line)
+
+
+class GprMaxGeometryTest(GeometryObjectsWriteMixin, ReceiverMixin, GprMaxBaseTest):
+    geometry_objects_write = ["full_volume"]

reframe_tests/tests/test_geometry_objects.py

@@ -1,16 +1,20 @@
 import reframe as rfm
 from reframe.core.builtins import parameter
 
-from reframe_tests.tests.mixins import MpiMixin
-from reframe_tests.tests.standard_tests import GprMaxGeometryObjectTest
+from reframe_tests.tests.mixins import GeometryOnlyMixin, MpiMixin, ReceiverMixin
+from reframe_tests.tests.standard_tests import (
+    GprMaxGeometryObjectsReadTest,
+    GprMaxGeometryObjectsReadWriteTest,
+    GprMaxGeometryObjectsWriteTest,
+)
 
 
 @rfm.simple_test
-class TestGeometryObject(GprMaxGeometryObjectTest):
+class TestGeometryObject(ReceiverMixin, GprMaxGeometryObjectsWriteTest):
     tags = {"test", "serial", "geometry only", "geometry object"}
     sourcesdir = "src/geometry_object_tests"
     model = parameter(["geometry_object_write"])
-    geometry_objects = ["partial_volume", "full_volume"]
+    geometry_objects_write = ["partial_volume", "full_volume"]
 
 
 @rfm.simple_test
@@ -18,3 +22,59 @@ class TestGeometryObjectMPI(MpiMixin, TestGeometryObject):
     tags = {"test", "mpi", "geometry only", "geometry object"}
     mpi_layout = parameter([[2, 2, 2], [4, 4, 1]])
     test_dependency = TestGeometryObject
+
+
+@rfm.simple_test
+class TestGeometryObjectReadFullVolume(ReceiverMixin, GprMaxGeometryObjectsReadTest):
+    tags = {"test", "serial", "geometry only", "geometry object"}
+    sourcesdir = "src/geometry_object_tests"
+    model = parameter(["geometry_object_read_full_volume"])
+    geometry_objects_read = {"full_volume": "full_volume_read"}
+    test_dependency = TestGeometryObject
+
+
+@rfm.simple_test
+class TestGeometryObjectReadFullVolumeMPI(MpiMixin, TestGeometryObjectReadFullVolume):
+    tags = {"test", "mpi", "geometry only", "geometry object"}
+    mpi_layout = parameter([[2, 2, 2], [4, 4, 1]])
+    test_dependency = TestGeometryObject
+
+
+@rfm.simple_test
+class TestGeometryObjectReadWrite(GeometryOnlyMixin, GprMaxGeometryObjectsReadWriteTest):
+    tags = {"test", "serial", "geometry only", "geometry object"}
+    sourcesdir = "src/geometry_object_tests"
+    model = parameter(["geometry_object_read_write"])
+    geometry_objects_read = {
+        "full_volume": "full_volume_read",
+    }
+    geometry_objects_write = ["partial_volume", "full_volume"]
+    test_dependency = TestGeometryObject
+
+
+@rfm.simple_test
+class TestGeometryObjectReadWriteMPI(MpiMixin, TestGeometryObjectReadWrite):
+    tags = {"test", "mpi", "geometry only", "geometry object"}
+    mpi_layout = parameter([[2, 2, 2], [4, 4, 1]])
+    test_dependency = TestGeometryObject
+
+
+# TODO: This test fails in the serial implementation due to the geometry
+# object being positioned such that it overflows the grid
+# @rfm.simple_test
+class TestGeometryObjectMove(GeometryOnlyMixin, GprMaxGeometryObjectsReadWriteTest):
+    tags = {"test", "serial", "geometry only", "geometry object"}
+    sourcesdir = "src/geometry_object_tests"
+    model = parameter(["geometry_object_move"])
+    geometry_objects_read = {
+        "full_volume": "full_volume_read",
+    }
+    geometry_objects_write = ["partial_volume"]
+    test_dependency = TestGeometryObject
+
+
+@rfm.simple_test
+class TestGeometryObjectMoveMPI(MpiMixin, TestGeometryObjectMove):
+    tags = {"test", "mpi", "geometry only", "geometry object"}
+    mpi_layout = parameter([[2, 2, 2], [4, 3, 1]])
+    test_dependency = TestGeometryObject