From 6529f35f8658389539774157cfd44bde11a1d279 Mon Sep 17 00:00:00 2001 From: mlauer154 Date: Thu, 1 Feb 2024 13:07:30 -0600 Subject: [PATCH] Fixed point-to-line distance signed equation, added symmetry constraint --- pymead/core/constraint_equations.py | 35 ++++++------ pymead/core/gcs2.py | 60 ++++++++++++++++----- pymead/core/point.py | 4 +- pymead/gui/airfoil_canvas.py | 1 - pymead/tests/core_tests/test_constraints.py | 16 ++++-- 5 files changed, 77 insertions(+), 39 deletions(-) diff --git a/pymead/core/constraint_equations.py b/pymead/core/constraint_equations.py index 26ed8605..66957494 100644 --- a/pymead/core/constraint_equations.py +++ b/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): diff --git a/pymead/core/gcs2.py b/pymead/core/gcs2.py index 28be0b76..26a9d803 100644 --- a/pymead/core/gcs2.py +++ b/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 = [] diff --git a/pymead/core/point.py b/pymead/core/point.py index 5a0f5417..c0dd1330 100644 --- a/pymead/core/point.py +++ b/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) diff --git a/pymead/gui/airfoil_canvas.py b/pymead/gui/airfoil_canvas.py index c4f7ee2c..c3e33169 100644 --- a/pymead/gui/airfoil_canvas.py +++ b/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) diff --git a/pymead/tests/core_tests/test_constraints.py b/pymead/tests/core_tests/test_constraints.py index 9c6db869..a3e65148 100644 --- a/pymead/tests/core_tests/test_constraints.py +++ b/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)