From e9bc7d3978b2151903431d20bcfdd223f99ee291 Mon Sep 17 00:00:00 2001 From: Nathan Mannall Date: Fri, 21 Feb 2025 18:29:11 +0000 Subject: [PATCH] Fix bug when objects start in a positive halo Objects with a thickness (triangles and circular sectors) can no longer extend beyond the grid in the dimension they are aligned to. Building these objects where the lower extent is aligned with a positive halo, no longer corrupts memory due to out of bounds memory access --- gprMax/user_inputs.py | 106 ++++++++++++------ .../cmds_geometry/cylindrical_sector.py | 4 +- gprMax/user_objects/cmds_geometry/triangle.py | 34 +++--- .../cylindrical_sector_z_outside_boundary.in | 2 +- 4 files changed, 89 insertions(+), 57 deletions(-) diff --git a/gprMax/user_inputs.py b/gprMax/user_inputs.py index 1f712bbf..6a550890 100644 --- a/gprMax/user_inputs.py +++ b/gprMax/user_inputs.py @@ -140,6 +140,9 @@ class UserInput(Generic[GridType]): """ return self.discretise_point(point) * self.grid.dl + def discrete_to_continuous(self, point: npt.NDArray[np.int32]) -> npt.NDArray[np.float64]: + return point * self.grid.dl + class MainGridUserInput(UserInput[GridType]): """Handles (x, y, z) points supplied by the user in the main grid.""" @@ -199,30 +202,29 @@ class MainGridUserInput(UserInput[GridType]): p2: Tuple[float, float, float], p3: Tuple[float, float, float], cmd_str: str, - ) -> Tuple[bool, npt.NDArray[np.int32], npt.NDArray[np.int32], npt.NDArray[np.int32]]: - p1_within_grid, p1_checked = self.check_point(p1, cmd_str, name="vertex_1") - p2_within_grid, p2_checked = self.check_point(p2, cmd_str, name="vertex_2") - p3_within_grid, p3_checked = self.check_point(p3, cmd_str, name="vertex_3") + ) -> Tuple[npt.NDArray[np.int32], npt.NDArray[np.int32], npt.NDArray[np.int32]]: + # We only care if the point are in the global grid (an error + # will be thrown if that is not the case). + _, p1_checked = self.check_point(p1, cmd_str, name="vertex_1") + _, p2_checked = self.check_point(p2, cmd_str, name="vertex_2") + _, p3_checked = self.check_point(p3, cmd_str, name="vertex_3") - return ( - p1_within_grid and p2_within_grid and p3_within_grid, - p1_checked, - p2_checked, - p3_checked, - ) + return p1_checked, p2_checked, p3_checked def check_thickness( - self, dimension: str, start: float, thickness: float + self, + dimension: str, + lower_extent: float, + thickness: float, + cmd_str: str, ) -> Tuple[bool, float, float]: """Check the thickness of an object in a specified dimension. Args: dimension: Dimension to check the thickness value for. This must have value x, y, or z. - start: Start coordinate of the object in the specified - dimension. This must be in the local grid coordinate - system - i.e. previously translated using the - round_to_grid function. + lower_extent: Lower extent of the object in the specified + dimension. thickness: Thickness of the object. Raises: @@ -231,35 +233,71 @@ class MainGridUserInput(UserInput[GridType]): Returns: within_grid: True if part of the object is within the current grid. False otherwise. - start: Start value limited to the bounds of the grid. - thickness: Thickness value such that start + thickness is - within the bounds of the grid. + lower_extent: Lower extent limited to the bounds of the + grid. + thickness: Thickness value such that lower_extent + + thickness is within the bounds of the grid. """ - # TODO: This should throw an error if the bounds are outside the - # global grid. + if thickness < 0: + raise ValueError(f"'{cmd_str}' requires a non negative thickness") + + if lower_extent < 0: + raise ValueError( + f"'{cmd_str}' lower extent should be non negative in the {dimension} dimension" + ) + + upper_extent = lower_extent + thickness + if dimension == "x": - grid_size = self.grid.nx * self.grid.dx + lower_point = self.discretise_point((lower_extent, 0, 0)) + upper_point = self.discretise_point((upper_extent, 0, 0)) + thickness_point = self.discretise_static_point((thickness, 0, 0)) + index = 0 elif dimension == "y": - grid_size = self.grid.ny * self.grid.dy + lower_point = self.discretise_point((0, lower_extent, 0)) + upper_point = self.discretise_point((0, upper_extent, 0)) + thickness_point = self.discretise_static_point((0, thickness, 0)) + index = 1 elif dimension == "z": - grid_size = self.grid.nz * self.grid.dz + lower_point = self.discretise_point((0, 0, lower_extent)) + upper_point = self.discretise_point((0, 0, upper_extent)) + thickness_point = self.discretise_static_point((0, 0, thickness)) + index = 2 else: raise ValueError("Dimension should have value x, y, or z") - end = start + thickness + try: + self.grid.within_bounds(upper_point) + except ValueError: + raise ValueError( + f"'{cmd_str}' extends beyond the size of the model in the {dimension} dimension" + ) - if start > grid_size: - return False, start, thickness - elif end < 0: - return False, start, thickness + # Work with discretised (int) values as reduces imprecision due + # to floating point calculations + size = self.grid.size[index] + lower_extent = lower_point[index] + upper_extent = upper_point[index] + thickness = thickness_point[index] - if start < 0: - thickness += start - start = 0 - if end > grid_size: - thickness -= end - grid_size + # These should only trigger for MPIGrids. + # TODO: Can this be structured so these checks happen in the + # MPIGridUserInput object? + if lower_extent < 0: + thickness += lower_extent + lower_extent = 0 + if upper_extent > size: + thickness -= upper_extent - size - return True, start, thickness + dl = self.grid.dl[index] + + return ( + lower_extent <= size + and upper_extent >= 0 + and not (upper_extent > size and thickness <= 0), + lower_extent * dl, + thickness * dl, + ) class MPIUserInput(MainGridUserInput[MPIGrid]): diff --git a/gprMax/user_objects/cmds_geometry/cylindrical_sector.py b/gprMax/user_objects/cmds_geometry/cylindrical_sector.py index 4562739e..97da02ca 100644 --- a/gprMax/user_objects/cmds_geometry/cylindrical_sector.py +++ b/gprMax/user_objects/cmds_geometry/cylindrical_sector.py @@ -90,7 +90,9 @@ class CylindricalSector(GeometryUserObject): elif normal == "z": ctr1, ctr2, level = uip.round_to_grid((ctr1, ctr2, extent1)) - sector_within_grid, level, thickness = uip.check_thickness(normal, level, thickness) + sector_within_grid, level, thickness = uip.check_thickness( + normal, extent1, thickness, self.__str__() + ) # Exit early if none of the cylindrical sector is in this grid # as there is nothing else to do. diff --git a/gprMax/user_objects/cmds_geometry/triangle.py b/gprMax/user_objects/cmds_geometry/triangle.py index c2d28abe..27d55edc 100644 --- a/gprMax/user_objects/cmds_geometry/triangle.py +++ b/gprMax/user_objects/cmds_geometry/triangle.py @@ -96,37 +96,35 @@ class Triangle(RotatableMixin, GeometryUserObject): logger.exception(f"{self.__str__()} no materials have been specified") raise - if thickness < 0: - logger.exception(f"{self.__str__()} requires a positive value for thickness") - raise ValueError - uip = self._create_uip(grid) # Check whether points are valid against grid - uip.check_tri_points(up1, up2, up3, self.__str__()) + dp1, dp2, dp3 = uip.check_tri_points(up1, up2, up3, self.__str__()) # Convert points to metres - x1, y1, z1 = uip.round_to_grid(up1) - x2, y2, z2 = uip.round_to_grid(up2) - x3, y3, z3 = uip.round_to_grid(up3) + x1, y1, z1 = uip.discrete_to_continuous(dp1) + x2, y2, z2 = uip.discrete_to_continuous(dp2) + x3, y3, z3 = uip.discrete_to_continuous(dp3) # Check for valid orientations # yz-plane triangle if x1 == x2 == x3: normal = "x" - start = x1 + lower_extent = up1[0] # xz-plane triangle elif y1 == y2 == y3: normal = "y" - start = y1 + lower_extent = up1[1] # xy-plane triangle elif z1 == z2 == z3: normal = "z" - start = z1 + lower_extent = up1[2] else: logger.exception(f"{self.__str__()} the triangle is not specified correctly") raise ValueError - triangle_within_grid, start, thickness = uip.check_thickness(normal, start, thickness) + triangle_within_grid, lower_extent, thickness = uip.check_thickness( + normal, lower_extent, thickness, self.__str__() + ) # Exit early if none of the triangle is in this grid as there is # nothing else to do. @@ -136,19 +134,13 @@ class Triangle(RotatableMixin, GeometryUserObject): # Update start bound of the triangle # yz-plane triangle if normal == "x": - x1 = start - x2 = start - x3 = start + x1 = x2 = x3 = lower_extent # xz-plane triangle elif normal == "y": - y1 = start - y2 = start - y3 = start + y1 = y2 = y3 = lower_extent # xy-plane triangle elif normal == "z": - z1 = start - z2 = start - z3 = start + z1 = z2 = z3 = lower_extent # Look up requested materials in existing list of material instances materials = [y for x in materialsrequested for y in grid.materials if y.ID == x] diff --git a/reframe_tests/tests/src/geometry_tests/cylindrical_sector_geometry/cylindrical_sector_z_outside_boundary.in b/reframe_tests/tests/src/geometry_tests/cylindrical_sector_geometry/cylindrical_sector_z_outside_boundary.in index 3950a615..2a609408 100644 --- a/reframe_tests/tests/src/geometry_tests/cylindrical_sector_geometry/cylindrical_sector_z_outside_boundary.in +++ b/reframe_tests/tests/src/geometry_tests/cylindrical_sector_geometry/cylindrical_sector_z_outside_boundary.in @@ -8,6 +8,6 @@ #rx: 0.080 0.080 0.080 #material: 8 0 1 0 half_space -#cylindrical_sector: z 0.07 0.06 0.04 0.12 0.05 30 240 half_space +#cylindrical_sector: z 0.07 0.06 0.04 0.1 0.05 30 240 half_space #geometry_objects_write: 0 0 0 0.1 0.1 0.1 full_volume