Added radius of curvature constraint, made small speed-up to canvas i…

…tem update after constraint param change

pymead/core/constraint_equations.py

@@ -57,6 +57,13 @@ def measure_radius_of_curvature_bezier(Lt: float, Lc: float, n: int, psi: float)
     return np.abs(np.true_divide(Lt ** 2, Lc * (1 - 1 / n) * np.sin(psi)))
 
 
+def measure_curvature_length_bezier(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 = np.abs(np.true_divide(Lt ** 2, R * (1 - 1 / n) * np.sin(psi)))
+    return Lc
+
+
 def measure_curvature_bezier(Lt: float, Lc: float, n: int, psi: float):
     return np.abs(np.true_divide(Lc * (1 - 1 / n) * np.sin(psi), Lt ** 2))
 

pymead/core/constraints.py

@@ -588,8 +588,7 @@ def __init__(self, curve_joint: Point, value: float or LengthParam, name: str =
         self.g2_point_curve_1 = self.curve_1.point_sequence().points()[self.g2_point_index_curve_1]
         self.g2_point_curve_2 = self.curve_2.point_sequence().points()[self.g2_point_index_curve_2]
 
-        points = [self.g2_point_curve_1, self.g1_point_curve_1,
-                  self.curve_joint,
+        points = [self.curve_joint, self.g2_point_curve_1, self.g1_point_curve_1,
                   self.g1_point_curve_2, self.g2_point_curve_2]
 
         param = value if isinstance(value, Param) else LengthParam(value=value, name="ROC-1")

pymead/core/gcs2.py

@@ -1,6 +1,4 @@
 import networkx
-import matplotlib.pyplot as plt
-import numpy as np
 
 from pymead.core.constraints import *
 from pymead.core.constraint_equations import *
@@ -197,7 +195,10 @@ def add_constraint(self, constraint: GeoCon):
             return
         elif isinstance(constraint, SymmetryConstraint):
             self.solve_symmetry_constraint(constraint)
-            self.update_from_points(constraint)
+            self.update_canvas_items(constraint)
+        elif isinstance(constraint, ROCurvatureConstraint):
+            self.solve_roc_constraint(constraint)
+            self.update_canvas_items(constraint)
 
     def remove_constraint(self, constraint: GeoCon):
         raise NotImplementedError("Constraint removal not yet implemented")
@@ -287,8 +288,10 @@ def solve(self, source: GeoCon):
 
         elif isinstance(source, SymmetryConstraint):
             self.solve_symmetry_constraint(source)
+        elif isinstance(source, ROCurvatureConstraint):
+            self.solve_roc_constraint(source)
 
-        self.solve_symmetry_constraints(points_solved)
+        self.solve_other_constraints(points_solved)
 
     @staticmethod
     def solve_symmetry_constraint(constraint: SymmetryConstraint):
@@ -311,6 +314,22 @@ def solve_symmetry_constraint(constraint: SymmetryConstraint):
             y3 + mirror_distance * np.sin(mirror_angle), force=True
         )
 
+    @staticmethod
+    def solve_roc_constraint(constraint: ROCurvatureConstraint):
+
+        def solve_for_single_curve(p_g1: Point, p_g2: Point, n: int):
+            Lc = measure_curvature_length_bezier(
+                constraint.curve_joint.x().value(), constraint.curve_joint.y().value(),
+                p_g1.x().value(), p_g1.y().value(),
+                p_g2.x().value(), p_g2.y().value(), constraint.param().value(), n
+            )
+            angle = p_g1.measure_angle(p_g2)
+            p_g2.request_move(p_g1.x().value() + Lc * np.cos(angle),
+                              p_g1.y().value() + Lc * np.sin(angle), force=True)
+
+        solve_for_single_curve(constraint.g1_point_curve_1, constraint.g2_point_curve_1, constraint.curve_1.degree)
+        solve_for_single_curve(constraint.g1_point_curve_2, constraint.g2_point_curve_2, constraint.curve_2.degree)
+
     def solve_symmetry_constraints(self, points: typing.List[Point]):
         symmetry_constraints_solved = []
         for point in points:
@@ -321,10 +340,37 @@ def solve_symmetry_constraints(self, points: typing.List[Point]):
                 self.solve_symmetry_constraint(symmetry_constraint)
                 symmetry_constraints_solved.append(symmetry_constraint)
 
+    def solve_roc_constraints(self, points: typing.List[Point]):
+        roc_constraints_solved =[]
+        for point in points:
+            roc_constraints = [geo_con for geo_con in point.geo_cons if isinstance(geo_con, ROCurvatureConstraint)]
+            for roc_constraint in roc_constraints:
+                if roc_constraint in roc_constraints_solved:
+                    continue
+                self.solve_roc_constraint(roc_constraint)
+                roc_constraints_solved.append(roc_constraint)
+
+    def solve_other_constraints(self, points: typing.List[Point]):
+        self.solve_symmetry_constraints(points)
+        self.solve_roc_constraints(points)
+
+    def update_canvas_items(self, source: GeoCon or Point):
+
+        points_to_update = []
+        if isinstance(source, Point):
+            starting_points = [source]
+        elif isinstance(source, GeoCon):
+            starting_points = source.child_nodes
+        else:
+            raise ValueError("source must be of type GeoCon or of type Point")
+        for starting_point in starting_points:
+            for point in networkx.dfs_preorder_nodes(self, source=starting_point):
+                if point in points_to_update:
+                    continue
+                points_to_update.append(point)
 
-    def update_from_points(self, constraint: GeoCon):
         curves_to_update = []
-        for point in self.points.values():
+        for point in points_to_update:
             if point.canvas_item is not None:
                 point.canvas_item.updateCanvasItem(point.x().value(), point.y().value())
 
@@ -344,7 +390,7 @@ def update_from_points(self, constraint: GeoCon):
                 airfoil.canvas_item.generatePicture()
 
         constraints_to_update = []
-        for point in networkx.dfs_preorder_nodes(self, source=constraint.child_nodes[0]):
+        for point in networkx.dfs_preorder_nodes(self, source=source.child_nodes[0]):
             for geo_con in point.geo_cons:
                 if geo_con not in constraints_to_update:
                     constraints_to_update.append(geo_con)

pymead/core/geometry_collection.py

@@ -800,7 +800,7 @@ def add_constraint(self, constraint: GeoCon, assign_unique_name: bool = True, **
             self.gcs.add_constraint(constraint)
             if isinstance(constraint, AntiParallel3Constraint) or isinstance(constraint, Perp3Constraint):
                 self.gcs.solve(constraint)
-                self.gcs.update_from_points(constraint)
+                self.gcs.update_canvas_items(constraint)
         except (OverConstrainedError, ValueError) as e:
             self.remove_pymead_obj(constraint)
             self.clear_selected_objects()

pymead/core/param.py

@@ -135,7 +135,7 @@ def set_value(self, value: float or int, updated_objs: typing.List[PymeadObj] =
 
         if self.gcs is not None:
             self.gcs.solve(self.geo_cons[0])
-            self.gcs.update_from_points(self.geo_cons[0])
+            self.gcs.update_canvas_items(self.geo_cons[0])
 
         if self.tree_item is not None:
             self.tree_item.treeWidget().itemWidget(self.tree_item, 1).setValue(self.value())

pymead/tests/core_tests/test_constraints.py

@@ -15,7 +15,7 @@ def test_rel_angle3_constraint(self):
         d1 = RelAngle3Constraint(p1, p2, p3, value=3*np.pi/4)
         geo_col.add_constraint(d1)
         geo_col.gcs.solve(d1)
-        geo_col.gcs.update_from_points(d1)
+        geo_col.gcs.update_canvas_items(d1)
         a1 = p2.measure_angle(p1)
         a2 = p2.measure_angle(p3)
         ra1 = (a1 - a2) % (2 * np.pi)
@@ -32,7 +32,7 @@ def test_rel_angle3_then_two_distance(self):
         for cnstr in [ra1, d1, d2]:
             geo_col.add_constraint(cnstr)
             geo_col.gcs.solve(cnstr)
-            geo_col.gcs.update_from_points(cnstr)
+            geo_col.gcs.update_canvas_items(cnstr)
         a1 = p2.measure_angle(p1)
         a2 = p2.measure_angle(p3)
         ra_val = (a1 - a2) % (2 * np.pi)