Feature/exclusions

ard/layout/exclusions.py

@@ -111,7 +111,7 @@ def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
                 regions=self.exclusion_regions,
             )
         )
-
+        
         outputs["exclusion_distances"] = -exclusion_distances
 
     def compute_partials(self, inputs, partials, discrete_inputs=None):

ard/utils/geometry.py

@@ -69,6 +69,10 @@ def get_nearest_polygons(
 
     return region
 
+# Pad all polygons to have the same number of vertices
+def pad_polygon(polygon, max_vertices):
+    padding = max_vertices - len(polygon)
+    return jnp.pad(polygon, ((0, padding), (0, 0)), mode="edge")
 
 def distance_multi_point_to_multi_polygon_ray_casting(
     points_x: np.ndarray[float],
@@ -97,35 +101,33 @@ def distance_multi_point_to_multi_polygon_ray_casting(
         np.ndarray: Constraint values for each turbine.
         np.ndarray (optional): Region assignments for each turbine (if `return_region` is True).
     """
-
     # Combine points_x and points_y into a single array of points
     points = jnp.stack([points_x, points_y], axis=1)
 
-    # Determine the maximum number of vertices in any polygon
-    max_vertices = max(len(polygon) for polygon in boundary_vertices)
-
-    # Pad all polygons to have the same number of vertices
-    def pad_polygon(polygon):
-        padding = max_vertices - len(polygon)
-        return jnp.pad(polygon, ((0, padding), (0, 0)), mode="edge")
-
-    padded_boundary_vertices = jnp.stack(
-        [pad_polygon(polygon) for polygon in boundary_vertices]
-    )
-
-    # Define a function to compute the distance for a single point and its assigned region
+    # Convert boundary_vertices to JAX arrays
+    boundary_vertices_jax = [jnp.asarray(poly, dtype=jnp.float32) for poly in boundary_vertices]
+
+    # Create a function for each polygon that computes distance
+    def make_distance_func(vertices):
+        def compute_for_polygon(point):
+            return distance_point_to_polygon_ray_casting(
+                point=point,
+                vertices=vertices,
+                s=s,
+                shift=tol,
+                return_distance=True,
+            )
+        return compute_for_polygon
+
+    # Create list of functions, one per polygon
+    distance_funcs = [make_distance_func(vertices) for vertices in boundary_vertices_jax]
+
+    # Define function to compute distance using jax.lax.switch
     def compute_distance(point, region_idx):
-        vertices = padded_boundary_vertices[region_idx]
-        return distance_point_to_polygon_ray_casting(
-            point=point,
-            vertices=vertices,
-            s=s,
-            shift=tol,
-            return_distance=True,
-        )
-
-    # Vectorize the computation over all points
-    distances = jax.vmap(compute_distance, in_axes=(0, 0))(points, regions)
+        return jax.lax.switch(region_idx, distance_funcs, point)
+
+    # Vectorize over all points
+    distances = jax.vmap(compute_distance)(points, regions)
 
     return distances
 
@@ -134,6 +136,61 @@ def compute_distance(point, region_idx):
     distance_multi_point_to_multi_polygon_ray_casting
 )
 
+# Define a function to process a single edge
+def process_edge(
+    edge_start: jnp.ndarray, 
+    edge_end: jnp.ndarray, 
+    point: jnp.ndarray, 
+    shift: float,
+):
+    """
+    Process a single polygon edge for ray-casting algorithm.
+
+    Determines if a vertical ray cast from the point up intersects
+    the edge, and calculates the distance from the point to the edge segment.
+
+    Args:
+        edge_start (jnp.ndarray): Start vertex of the edge ([x, y]).
+        edge_end (jnp.ndarray): End vertex of the edge ([x, y]).
+        point (jnp.ndarray): Point of interest ([x, y]).
+        shift (float): Small value added to avoid division by zero when edge is vertical.
+
+    Returns:
+        tuple: A tuple containing:
+            - is_below (bool): True if point is below the edge at the ray intersection.
+            - distance (float): Shortest distance from point to the edge segment.
+            - vertex_crossing (int): 1 if ray crosses exactly at edge_start vertex, 0 otherwise.
+    """
+
+    # Check if the x-coordinate of the point is between the x-coordinates of the edge, 
+    # including when the point is directly below a vertical edge
+    x_condition = ((edge_start[0] <= point[0]) & (point[0] < edge_end[0])) | (
+                   (edge_start[0] >= point[0]) & (point[0] > edge_end[0])) | (
+                   (jnp.isclose(edge_start[0], edge_end[0]) & jnp.isclose(point[0], edge_start[0]))
+                   )
+
+    # Calculate the y-coordinate of the edge at the x-coordinate of the point
+    y = ((edge_end[1] - edge_start[1]) / (edge_end[0] - edge_start[0] + shift)) * (
+        point[0] - edge_start[0] + shift
+    ) + edge_start[1]
+
+    # Determine if the point is below the edge
+    is_below = x_condition & (point[1] < y)
+
+    # record vertex crossings
+    vertex_crossing = jnp.where(
+        x_condition & (edge_start[0] == point[0]) & is_below,
+        1,
+        0
+    )
+
+    # Calculate the distance to the edge
+    distance = distance_point_to_lineseg_nd(point, edge_start, edge_end)
+
+    return is_below, distance, vertex_crossing
+
+# Vectorize the edge processing function
+process_edge_vec = jax.vmap(process_edge, in_axes=(0, 0, None, None))
 
 def distance_point_to_polygon_ray_casting(
     point: jnp.ndarray,
@@ -164,59 +221,45 @@ def distance_point_to_polygon_ray_casting(
     Returns:
         float: Signed distance or inside/outside status. Negative if inside, positive if outside.
     """
+
     # Ensure inputs are JAX arrays with explicit data types
     point = jnp.asarray(point, dtype=jnp.float32)
     vertices = jnp.asarray(vertices, dtype=jnp.float32)
 
     # Add the first vertex to the end to close the polygon loop
     vertices = jnp.vstack([vertices, vertices[0]])
 
-    # Define a function to process a single edge
-    def process_edge(edge_start, edge_end, point):
-        # Check if the x-coordinate of the point is between the x-coordinates of the edge
-        x_condition = ((edge_start[0] <= point[0]) & (point[0] < edge_end[0])) | (
-            (edge_start[0] >= point[0]) & (point[0] > edge_end[0])
-        )
-
-        # Calculate the y-coordinate of the edge at the x-coordinate of the point
-        y = (edge_end[1] - edge_start[1]) / (edge_end[0] - edge_start[0] + shift) * (
-            point[0] - edge_start[0]
-        ) + edge_start[1]
-
-        # Determine if the point is below the edge
-        is_below = x_condition & (point[1] < y)
-
-        # Calculate the distance to the edge
-        distance = distance_point_to_lineseg_nd(point, edge_start, edge_end)
-
-        return is_below, distance
-
-    # Vectorize the edge processing function
-    process_edge_vec = jax.vmap(process_edge, in_axes=(0, 0, None))
-
     edge_starts = vertices[:-1]
     edge_ends = vertices[1:]
-    is_below, distances = process_edge_vec(edge_starts, edge_ends, point)
+    is_below, distances, vertex_crossings = process_edge_vec(edge_starts, edge_ends, point, shift)
+
+    # Check for tangential vertex crossings
+    # Create arrays of current and previous edge differences
+    curr_edge_end_diff = edge_ends[:, 0] - point[0]
+    prev_edge_start_diff = jnp.roll(edge_starts[:, 0], 1) - point[0]  # Shift by 1 to get previous
+    # Tangential crossing occurs when vertex_crossing > 0 AND the other two connected vertices are on same side
+    tangential_vertex_crossings = jnp.where(
+        (vertex_crossings > 0) & (curr_edge_end_diff * prev_edge_start_diff > 0),
+        1,
+        0
+    )
 
-    # Count the number of intersections
-    intersection_counter = jnp.sum(is_below)
+    # Count the number of intersections, ignoring tangential vertex crossings
+    intersection_counter = jnp.sum(is_below) - jnp.sum(tangential_vertex_crossings)
 
-    # Compute the signed distance
+    # Compute the signed distance if desired
     if return_distance:
-        c = smooth_min(distances, s=s)
-        c = jax.lax.cond(
-            intersection_counter % 2 == 1,
-            lambda _: -c,
-            lambda _: c,
-            operand=None,
-        )
+        c_prime = smooth_min(distances, s=s)
     else:
-        c = jax.lax.cond(
+        c_prime = 1
+
+    c = jax.lax.cond(
             intersection_counter % 2 == 1,
-            lambda _: -1.0,
-            lambda _: 1.0,
+            lambda _: -c_prime,
+            lambda _: c_prime,
             operand=None,
         )
+
     return c
 
 

test/ard/system/geometry/test_constraints.py

@@ -116,7 +116,7 @@ def test_constraint_evaluation(self, subtests):
 
     def test_constraint_optimization(self, subtests):
         """test boundary-constrained optimization distances (yes derivatives)"""
-
+        #TODO I think this test is not formulated correctly
         # setup the working/design variables
         self.prob.model.add_design_var("spacing_target", lower=2.0, upper=13.0)
         self.prob.model.add_constraint("boundary_distances", upper=0.0)
@@ -143,8 +143,8 @@ def test_constraint_optimization(self, subtests):
             )
 
         # make sure the target spacing matches well
-        spacing_target_validation = 5.46721656  # from a run on 24 June 2025
-        area_target_validation = 10.49498327  # from a run on 24 June 2025
+        spacing_target_validation = 5.46727038  # from a run on 15 Jan 2026
+        area_target_validation = 10.4951899  # from a run on 15 Jan 2026
         with subtests.test("validation spacing matches"):
             assert np.isclose(
                 self.prob.get_val("spacing_target"), spacing_target_validation

test/ard/unit/layout/test_boundary.py

@@ -293,7 +293,8 @@ def test_single_triangle_distance(self):
 
     def test_offset_triangle_distance(self):
         """make sure boundaries not on the origin are working"""
-
+        bx = [0.0, 800.0, 0.0]
+        by = [1000.0, 1000.0, 200.0]
         # set modeling options
         modeling_options_single = {
             "windIO_plant": {
@@ -303,8 +304,8 @@ def test_offset_triangle_distance(self):
                     "boundaries": {
                         "polygons": [
                             {
-                                "x": [0.0, 800.0, 0.0],
-                                "y": [1000.0, 1000.0, 200.0],
+                                "x": bx,
+                                "y": by,
                                 # "y": [200.0, 1000.0, 1000.0],
                             },
                         ]
@@ -338,17 +339,82 @@ def test_offset_triangle_distance(self):
 
         prob_single.run_model()
 
-        expected_distances = -np.array(  # minus sign on the data from exclusion
+        expected_distances = np.array(
             [
-                -200.000000000000,  # 0
-                -424.2640687119286,  # 1
-                -707.1067811865476,  # 2
+                200.000000000000,  # 0
+                424.2640687119286,  # 1
+                707.1067811865476,  # 2
                 0.000000000000,  # 3
-                -141.4213562373095,  # 4
-                -424.2640687119286,  # 5
+                141.4213562373095,  # 4
+                424.2640687119286,  # 5
                 0.000000000000,  # 6
-                141.4213562373095,  # 7
-                -141.4213562373095,  # 8
+                -141.4213562373095,  # 7
+                141.4213562373095,  # 8
+            ]
+        )
+
+        assert np.allclose(
+            prob_single["boundary_distances"], expected_distances, atol=1e-3
+        )
+
+    def test_offset_triangle_distance_left(self):
+        """make sure boundaries not on the origin are working"""
+        bx = [800.0, 0.0, 800.0]
+        by = [1000.0, 1000.0, 200.0]
+        # set modeling options
+        modeling_options_single = {
+            "windIO_plant": {
+                "name": "unit test dummy",
+                "site": {
+                    "name": "unit test site",
+                    "boundaries": {
+                        "polygons": [
+                            {
+                                "x": bx,
+                                "y": by,
+                            },
+                        ]
+                    },
+                },
+                "wind_farm": {
+                    "turbine": {
+                        "rotor_diameter": self.D_rotor,
+                    }
+                },
+            },
+            "layout": {
+                "N_turbines": self.N_turbines,
+            },
+        }
+
+        # set up openmdao problem
+        model_single = om.Group()
+        model_single.add_subsystem(
+            "boundary",
+            boundary.FarmBoundaryDistancePolygon(
+                modeling_options=modeling_options_single,
+            ),
+            promotes=["*"],
+        )
+        prob_single = om.Problem(model_single)
+        prob_single.setup()
+
+        prob_single.set_val("x_turbines", self.x_turbines)
+        prob_single.set_val("y_turbines", self.y_turbines)
+
+        prob_single.run_model()
+
+        expected_distances = np.array(  # minus sign on the data from exclusion
+            [
+                707.1067811865476,  # 0
+                424.2640687119286,  # 1
+                200.0,  # 2
+                424.2640687119286,  # 3
+                141.4213562373095,  # 4
+                0.0,  # 5
+                141.4213562373095,  # 6
+                -141.4213562373095,  # 7
+                0.0,  # 8
             ]
         )
 
@@ -380,6 +446,12 @@ def test_multi_polygon_distance(self):
         region_assignments = np.ones(self.N_turbines, dtype=int)
         region_assignments[0:3] = 0
 
+        bx1 = boundary_vertices_0[:, 0].tolist()
+        by1 = boundary_vertices_0[:, 1].tolist()
+
+        bx2 = boundary_vertices_1[:, 0].tolist()
+        by2 = boundary_vertices_1[:, 1].tolist()
+
         # set modeling options
         modeling_options_multi = {
             "windIO_plant": {
@@ -389,12 +461,12 @@ def test_multi_polygon_distance(self):
                     "boundaries": {
                         "polygons": [
                             {
-                                "x": boundary_vertices_0[:, 0].tolist(),
-                                "y": boundary_vertices_0[:, 1].tolist(),
+                                "x": bx1,
+                                "y": by1,
                             },
                             {
-                                "x": boundary_vertices_1[:, 0].tolist(),
-                                "y": boundary_vertices_1[:, 1].tolist(),
+                                "x": bx2,
+                                "y": by2,
                             },
                         ]
                     },