Fixed point-to-line distance signed equation, added symmetry constraint

pymead/core/constraint_equations.py

@@ -1,25 +1,22 @@
 from collections import namedtuple
 
-import jax
 import numpy as np
-from jax import jit
-from jax import numpy as jnp
 
 
 def measure_distance(x1: float, y1: float, x2: float, y2: float):
-    return jnp.hypot(x1 - x2, y1 - y2)
+    return np.hypot(x1 - x2, y1 - y2)
 
 
 def measure_abs_angle(x1: float, y1: float, x2: float, y2: float):
-    return (jnp.arctan2(y2 - y1, x2 - x1)) % (2 * jnp.pi)
+    return (np.arctan2(y2 - y1, x2 - x1)) % (2 * np.pi)
 
 
 def measure_rel_angle3(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float):
-    return (jnp.arctan2(y1 - y2, x1 - x2) - jnp.arctan2(y3 - y2, x3 - x2)) % (2 * jnp.pi)
+    return (np.arctan2(y1 - y2, x1 - x2) - np.arctan2(y3 - y2, x3 - x2)) % (2 * np.pi)
 
 
 def measure_rel_angle4(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float, x4: float, y4: float):
-    return (jnp.arctan2(y4 - y3, x4 - x3) - jnp.arctan2(y2 - y1, x2 - x1)) % (2 * jnp.pi)
+    return (np.arctan2(y4 - y3, x4 - x3) - np.arctan2(y2 - y1, x2 - x1)) % (2 * np.pi)
 
 
 def measure_point_line_distance_signed(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float):
@@ -49,19 +46,19 @@ def measure_point_line_distance_signed(x1: float, y1: float, x2: float, y2: floa
     float
         Distance from the target point to the line
     """
-    return (x2 - x1) * (y1 - y3) - (x1 - x3) * (y2 - y1) / measure_distance(x1, y1, x2, y2)
+    return ((x2 - x1) * (y1 - y3) - (x1 - x3) * (y2 - y1)) / measure_distance(x1, y1, x2, y2)
 
 
 def measure_point_line_distance_unsigned(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float):
-    return jnp.abs(measure_point_line_distance_signed(x1, y1, x2, y2, x3, y3))
+    return np.abs(measure_point_line_distance_signed(x1, y1, x2, y2, x3, y3))
 
 
 def measure_radius_of_curvature_bezier(Lt: float, Lc: float, n: int, psi: float):
-    return jnp.abs(jnp.true_divide(Lt ** 2, Lc * (1 - 1 / n) * jnp.sin(psi)))
+    return np.abs(np.true_divide(Lt ** 2, Lc * (1 - 1 / n) * np.sin(psi)))
 
 
 def measure_curvature_bezier(Lt: float, Lc: float, n: int, psi: float):
-    return jnp.abs(jnp.true_divide(Lc * (1 - 1 / n) * jnp.sin(psi), Lt ** 2))
+    return np.abs(np.true_divide(Lc * (1 - 1 / n) * np.sin(psi), Lt ** 2))
 
 
 def measure_data_bezier_curve_joint(xy: np.ndarray, n: np.ndarray):
@@ -71,15 +68,15 @@ def measure_data_bezier_curve_joint(xy: np.ndarray, n: np.ndarray):
     theta2 = measure_abs_angle(xy[3, 0], xy[3, 1], xy[4, 0], xy[4, 1])
     psi1 = theta1 - phi1
     psi2 = theta2 - phi2
-    phi_rel = (phi1 - phi2) % (2 * jnp.pi)
+    phi_rel = (phi1 - phi2) % (2 * np.pi)
     Lt1 = measure_distance(xy[1, 0], xy[1, 1], xy[2, 0], xy[2, 1])
     Lt2 = measure_distance(xy[2, 0], xy[2, 1], xy[3, 0], xy[3, 1])
     Lc1 = measure_distance(xy[0, 0], xy[0, 1], xy[1, 0], xy[1, 1])
     Lc2 = measure_distance(xy[3, 0], xy[3, 1], xy[4, 0], xy[4, 1])
     kappa1 = measure_curvature_bezier(Lt1, Lc1, n[0], psi1)
     kappa2 = measure_curvature_bezier(Lt2, Lc2, n[1], psi2)
-    R1 = jnp.true_divide(1, kappa1)
-    R2 = jnp.true_divide(1, kappa2)
+    R1 = np.true_divide(1, kappa1)
+    R2 = np.true_divide(1, kappa2)
     n1 = n[0]
     n2 = n[1]
     field_names = ["phi1", "phi2", "theta1", "theta2", "psi1", "psi2", "phi_rel", "Lt1", "Lt2", "Lc1", "Lc2",
@@ -94,7 +91,7 @@ def measure_data_bezier_curve_joint(xy: np.ndarray, n: np.ndarray):
 def radius_of_curvature_constraint(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float, R: float, n: int):
     Lt = measure_distance(x1, y1, x2, y2)
     psi = measure_rel_angle3(x1, y1, x2, y2, x3, y3)
-    Lc = jnp.abs(jnp.true_divide(Lt ** 2, R * (1 - 1 / n) * jnp.sin(psi)))
+    Lc = np.abs(np.true_divide(Lt ** 2, R * (1 - 1 / n) * np.sin(psi)))
     return distance_constraint(x2, y2, x3, y3, Lc)
 
 
@@ -163,19 +160,19 @@ def rel_angle4_constraint(x1: float, y1: float, x2: float, y2: float, x3: float,
 
 
 def perp3_constraint(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float):
-    return measure_rel_angle3(x1, y1, x2, y2, x3, y3) - (jnp.pi / 2)
+    return measure_rel_angle3(x1, y1, x2, y2, x3, y3) - (np.pi / 2)
 
 
 def perp4_constraint(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float, x4: float, y4: float):
-    return measure_rel_angle4(x1, y1, x2, y2, x3, y3, x4, y4) - (jnp.pi / 2)
+    return measure_rel_angle4(x1, y1, x2, y2, x3, y3, x4, y4) - (np.pi / 2)
 
 
 def antiparallel3_constraint(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float):
-    return measure_rel_angle3(x1, y1, x2, y2, x3, y3) - jnp.pi
+    return measure_rel_angle3(x1, y1, x2, y2, x3, y3) - np.pi
 
 
 def antiparallel4_constraint(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float, x4: float, y4: float):
-    return measure_rel_angle4(x1, y1, x2, y2, x3, y3, x4, y4) - jnp.pi
+    return measure_rel_angle4(x1, y1, x2, y2, x3, y3, x4, y4) - np.pi
 
 
 def parallel3_constraint(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float):

pymead/core/gcs2.py

@@ -1,7 +1,9 @@
 import networkx
 import matplotlib.pyplot as plt
+import numpy as np
 
 from pymead.core.constraints import *
+from pymead.core.constraint_equations import *
 from pymead.core.point import Point
 
 
@@ -193,6 +195,9 @@ def add_constraint(self, constraint: GeoCon):
 
             self.add_edge(constraint.p2, constraint.p3, angle=constraint)
             return
+        elif isinstance(constraint, SymmetryConstraint):
+            self.solve_symmetry_constraint(constraint)
+            self.update_from_points(constraint)
 
     def remove_constraint(self, constraint: GeoCon):
         raise NotImplementedError("Constraint removal not yet implemented")
@@ -201,6 +206,7 @@ def make_point_copies(self):
         return {point: Point(x=point.x().value(), y=point.y().value()) for point in self.points.values()}
 
     def solve(self, source: GeoCon):
+        points_solved = []
         if isinstance(source, DistanceConstraint):
             edge_data_p12 = self.get_edge_data(source.p1, source.p2)
             if edge_data_p12 and edge_data_p12["distance"] is source:
@@ -217,12 +223,8 @@ def solve(self, source: GeoCon):
             for point in networkx.bfs_tree(self, source=start):
                 point.x().set_value(point.x().value() + dx)
                 point.y().set_value(point.y().value() + dy)
-        # elif isinstance(source, RelAngle3Constraint):
-        #     angle_1 = source.vertex.measure_angle(source.start_point)
-        #     dist = source.vertex.measure_distance(source.end_point)
-        #     source.end_point.x().set_value(source.vertex.x().value() + dist * np.cos(angle_1 + source.param().value()))
-        #     source.end_point.y().set_value(source.vertex.y().value() + dist * np.sin(angle_1 + source.param().value()))
-        #     starting_point = source.end_point
+                if point not in points_solved:
+                    points_solved.append(point)
         elif isinstance(source, AntiParallel3Constraint) or isinstance(source, Perp3Constraint):
             edge_data_p21 = self.get_edge_data(source.p2, source.p1)
             edge_data_p23 = self.get_edge_data(source.p2, source.p3)
@@ -250,6 +252,8 @@ def solve(self, source: GeoCon):
                 new_xy = (rotation_mat @ dx_dy + rotation_point_mat).flatten()
                 point.x().set_value(new_xy[0])
                 point.y().set_value(new_xy[1])
+                if point not in points_solved:
+                    points_solved.append(point)
 
         elif isinstance(source, RelAngle3Constraint):
             edge_data_p21 = self.get_edge_data(source.vertex, source.start_point)
@@ -278,13 +282,45 @@ def solve(self, source: GeoCon):
                 new_xy = (rotation_mat @ dx_dy + rotation_point_mat).flatten()
                 point.x().set_value(new_xy[0])
                 point.y().set_value(new_xy[1])
+                if point not in points_solved:
+                    points_solved.append(point)
+
+        elif isinstance(source, SymmetryConstraint):
+            self.solve_symmetry_constraint(source)
+
+        self.solve_symmetry_constraints(points_solved)
+
+    @staticmethod
+    def solve_symmetry_constraint(constraint: SymmetryConstraint):
+        x1, y1 = constraint.p1.x().value(), constraint.p1.y().value()
+        x2, y2 = constraint.p2.x().value(), constraint.p2.y().value()
+        x3, y3 = constraint.p3.x().value(), constraint.p3.y().value()
+        line_angle = measure_abs_angle(x1, y1, x2, y2)
+        tool_angle = measure_rel_angle3(x1, y1, x2, y2, x3, y3)
+        if tool_angle < np.pi:
+            mirror_angle = line_angle - np.pi / 2
+        elif tool_angle > np.pi:
+            mirror_angle = line_angle + np.pi / 2
+        else:
+            constraint.p4.request_move(constraint.p3.x().value(), constraint.p3.y().value(), force=True)
+            return
+
+        mirror_distance = 2 * measure_point_line_distance_unsigned(x1, y1, x2, y2, x3, y3)
+        constraint.p4.request_move(
+            x3 + mirror_distance * np.cos(mirror_angle),
+            y3 + mirror_distance * np.sin(mirror_angle), force=True
+        )
+
+    def solve_symmetry_constraints(self, points: typing.List[Point]):
+        symmetry_constraints_solved = []
+        for point in points:
+            symmetry_constraints = [geo_con for geo_con in point.geo_cons if isinstance(geo_con, SymmetryConstraint)]
+            for symmetry_constraint in symmetry_constraints:
+                if symmetry_constraint in symmetry_constraints_solved:
+                    continue
+                self.solve_symmetry_constraint(symmetry_constraint)
+                symmetry_constraints_solved.append(symmetry_constraint)
 
-        # for point in self.adj[starting_point]:
-        #     cnstr_dict = self.get_edge_data(starting_point, point)
-        #     print(f"{point = }, {starting_point = }, {cnstr_dict = }")
-        #     if len(cnstr_dict) == 1:
-        #         dist = point_copies[starting_point].measure_distance(point_copies[point])
-        #         angle = point_copies[starting_point].measure_angle(point_copies[point])
 
     def update_from_points(self, constraint: GeoCon):
         curves_to_update = []

pymead/core/point.py

@@ -83,9 +83,9 @@ def measure_distance(self, other: "Point"):
     def measure_angle(self, other: "Point"):
         return np.arctan2(other.y().value() - self.y().value(), other.x().value() - self.x().value())
 
-    def request_move(self, xp: float, yp: float):
+    def request_move(self, xp: float, yp: float, force: bool = False):
 
-        if self.gcs is None or (self.gcs is not None and len(self.geo_cons) == 0):
+        if self.gcs is None or (self.gcs is not None and len(self.geo_cons) == 0) or force:
 
             self.x().set_value(xp)
             self.y().set_value(yp)

pymead/gui/airfoil_canvas.py

@@ -444,7 +444,6 @@ def addROCurvatureConstraint(self):
         curvature_data = ROCurvatureConstraint.calculate_curvature_data(self.geo_col.selected_objects["points"][0])
         R = 0.5 * (curvature_data.R1 + curvature_data.R2)
         # R = curvature_data.R1
-        print(f"{R = }")
         R_param = self.geo_col.add_param(R, name="ROC-1", unit_type="length")
         constraint = ROCurvatureConstraint(*self.geo_col.selected_objects["points"], value=R_param)
         self.geo_col.add_constraint(constraint)

pymead/tests/core_tests/test_constraints.py

@@ -1,12 +1,7 @@
 from unittest import TestCase
 
-import numpy as np
-
-from pymead.core.gcs2 import GCS2
 from pymead.core.constraints import *
 from pymead.core.geometry_collection import GeometryCollection
-from pymead.core.point import Point
-from pymead.core.param import LengthParam, AngleParam
 
 
 class GCSTests(TestCase):
@@ -44,3 +39,14 @@ def test_rel_angle3_then_two_distance(self):
         self.assertAlmostEqual(ra_val, ra1.param().value())
         self.assertAlmostEqual(p1.measure_distance(p2), d1.param().value())
         self.assertAlmostEqual(p2.measure_distance(p3), d2.param().value())
+
+    def test_symmetry_constraint(self):
+        geo_col = GeometryCollection()
+        p1 = geo_col.add_point(0.4, -0.03)
+        p2 = geo_col.add_point(1.2, -0.03)
+        p3 = geo_col.add_point(0.5, 0.22)
+        p4 = geo_col.add_point(0.8, 0.9)
+        s1 = SymmetryConstraint(p1, p2, p3, p4)
+        geo_col.add_constraint(s1)
+        self.assertAlmostEqual(p4.x().value(), 0.5)
+        self.assertAlmostEqual(p4.y().value(), -0.28)