Update reading geometry objects for parallel build

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/src/geometry_object_tests/geometry_object_move.in

@@ -8,5 +8,5 @@
 
 #rx: 0.08 0.08 0.08
 
-#geometry_objects_read: 0.05 0.05 0.05 partial_volume_read.h5 partial_volume_read_materials.txt
+#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_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/standard_tests.py

@@ -1,3 +1,7 @@
+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 (
     GeometryObjectsReadMixin,
@@ -7,6 +11,7 @@ from reframe_tests.tests.mixins import (
     ReceiverMixin,
     SnapshotMixin,
 )
+from reframe_tests.tests.regression_checks import GeometryObjectMaterialsRegressionCheck
 
 
 class GprMaxRegressionTest(ReceiverMixin, GprMaxBaseTest):
@@ -32,7 +37,23 @@ class GprMaxGeometryObjectsReadTest(GeometryObjectsReadMixin, GprMaxBaseTest):
 class GprMaxGeometryObjectsReadWriteTest(
     GeometryObjectsReadMixin, GeometryObjectsWriteMixin, GprMaxBaseTest
 ):
-    pass
+    @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):

reframe_tests/tests/test_geometry_objects.py

@@ -30,7 +30,6 @@ class TestGeometryObjectReadFullVolume(ReceiverMixin, GprMaxGeometryObjectsReadT
     sourcesdir = "src/geometry_object_tests"
     model = parameter(["geometry_object_read_full_volume"])
     geometry_objects_read = {"full_volume": "full_volume_read"}
-    geometry_objects_write = ["full_volume"]
     test_dependency = TestGeometryObject
 
 
@@ -45,9 +44,11 @@ class TestGeometryObjectReadFullVolumeMPI(MpiMixin, TestGeometryObjectReadFullVo
 class TestGeometryObjectReadWrite(GeometryOnlyMixin, GprMaxGeometryObjectsReadWriteTest):
     tags = {"test", "serial", "geometry only", "geometry object"}
     sourcesdir = "src/geometry_object_tests"
-    model = parameter(["geometry_object_read_write_partial_volume"])
-    geometry_objects_read = {"partial_volume": "partial_volume_read"}
-    geometry_objects_write = ["partial_volume"]
+    model = parameter(["geometry_object_read_write"])
+    geometry_objects_read = {
+        "full_volume": "full_volume_read",
+    }
+    geometry_objects_write = ["partial_volume", "full_volume"]
     test_dependency = TestGeometryObject
 
 
@@ -56,3 +57,24 @@ 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