diff --git a/pymead/core/constraint_graph.py b/pymead/core/constraint_graph.py deleted file mode 100644 index e64901e0..00000000 --- a/pymead/core/constraint_graph.py +++ /dev/null @@ -1,722 +0,0 @@ -from copy import deepcopy - -import jaxopt -from jax import numpy as jnp -from jax import jit -import numpy as np -import networkx - -from pymead.core.point import Point -from pymead.core.constraints import * - - -class PymeadRootFinder(jaxopt.ScipyRootFinding): - def __init__(self, equation_system: typing.Callable, method: str = "lm"): - r""" - Parameters - ---------- - equation_system: typing.Callable - The non-linear system of equations to be solved. This system of equations represents interdependencies - (constraints) between geometric points. For example, a distance constraint with a value of 0.5 between two - points, :math:`P_1 = (x_1, y_1)` and :math:`P_2 = (x_2, y_2)` is represented by the equation - - .. math:: - \sqrt{(x_1^2 - x_2^2) + (y_1^2 - y_2^2)} = 0.5 - - method: str - Root-finding technique to use. The available methods are the subset of methods available in - ``scipy.optimize.root`` that allow for a user-specified Jacobian rather than a finite difference - approximation. Available methods: - - - ``"lm"``: Levenberg-Marquardt non-linear damped least-squares method - (`see here `_) - - ``"hybr"``: MINPACK's ``hybrd`` and ``hybrj`` routines - (`see here `_) - """ - super().__init__( - method=method, - jit=True, - has_aux=False, - optimality_fun=equation_system, - tol=1e-10, - use_jacrev=True, # Use the forward Jacobian calculation since the matrix is square - ) - - def solve(self, x0: np.ndarray, start_param_vec: np.ndarray, intermediate_param_vec: np.ndarray): - """ - Solves the compiled non-linear system of equations for this ``Point`` using Jax's wrapper for - ``scipy.optimize.root``. - - Parameters - ---------- - x0: np.ndarray - Initial guess for the solution to the system of equations - - start_param_vec: np.ndarray - The parameter vector - - intermediate_param_vec: np.ndarray - The parameter vector - - Returns - ------- - np.ndarray, dict - A two-element tuple containing the new parameter vector and information about the solution state - """ - return self.run(x0, start_param_vec, intermediate_param_vec) - - -class EquationData: - - def __init__(self): - self.root_finders = {} - self.root_finders_uncompiled = {} - self.geo_cons = [] - self.equations = [] - self.arg_idx_array = [] - self.variable_pos = [] - - def clear(self): - self.root_finders.clear() - self.root_finders_uncompiled.clear() - self.geo_cons.clear() - self.equations.clear() - self.arg_idx_array.clear() - self.variable_pos.clear() - - -class ConstraintGraph(networkx.Graph): - def __init__(self): - self.points = [] - self.constraint_params = [] - super().__init__() - - def add_point(self, point: Point): - """ - Adds a ``Point`` as a new (unconnected) node in the graph, and to the graph instance's list of points. - - Parameters - ---------- - point: Point - Point to add to the graph - - Returns - ------- - - """ - point.gcs = self - self.add_node(point) - self.points.append(point) - - def remove_point(self, point: Point): - connected_constraints = point.geo_cons.copy() - self.remove_node(point) - self.points.remove(point) - - for constraint in connected_constraints: - self.remove_constraint(constraint) - - def add_constraint(self, constraint: GeoCon, compile: bool = True, solve_and_update: bool = True): - """ - Adds the specified constraint to the graph, then analyzes the entire graph of connected points and constraints - and adds any weak constraints needed to make the problem well-constrained. The non-linear system of equations - representing the constraints is then solved to tolerance, and all points updated. - - Parameters - ---------- - constraint: GeoCon - GeoCon to add to the graph - - compile: bool - Whether to compile the non-linear system of equations after adding this constraint. Normally left as - ``True`` except for most constraints when loading an airfoil into ``pymead``. Note that if this option - is set to ``False``, the "solve and update" step will not be performed. Default: ``True`` - - solve_and_update: bool - Whether to solve the system of equations and update the points. Normally left as ``True`` except - for most constraints when loading an airfoil into ``pymead``. Default: ``True`` - - Returns - ------- - - """ - if constraint.param() is not None: - constraint.param().gcs = self - if constraint.param() not in self.constraint_params: - self.constraint_params.append(constraint.param()) - - if len(constraint.secondary_params) != 0: - for param in constraint.secondary_params: - param.gcs = self - self.constraint_params.append(param) - - self.add_node(constraint) - for child_node in constraint.child_nodes: - self.add_edge(constraint, child_node) - - # Perform depth-first search on the nodes to find if any connected node has the "data" attr set - for node in networkx.dfs_preorder_nodes(self, source=constraint): - if isinstance(node, GeoCon) and node.data is not None: - constraint.data = node.data - break - else: # If not, assign a new instance of the EquationData class to the "data" attr - constraint.data = EquationData() - - # Analyze the constraint to add the strong constraints and any weak constraints necessary to make the problem - # well-constrained - self.analyze(constraint) - - # Compile two separate root finders: one using Levenberg-Marquardt damped non-linear least squares algorithm, - # another using MINPACK's hybrd/hybrj routines - if compile: - self.compile_equation_for_entity_or_constraint(constraint, method="lm") - self.compile_equation_for_entity_or_constraint(constraint, method="hybr") - - # Solve using first the least-squares method and then MINPACK if necessary. Update the points if the solution - # falls within the tolerance specified in the PymeadRootFinder class - if compile and solve_and_update: - self.multisolve_and_update(constraint) - - def remove_constraint(self, constraint: GeoCon): - """ - Removes a constraint from the ``ConstraintGraph`` - - Parameters - ---------- - constraint: GeoCon - Geometric constraint to remove - - Returns - ------- - - """ - - # Get a copy of the constraint's neighbors - adj_points = [nbr for nbr in self.adj[constraint]].copy() - - # Remove the constraint from the constraint graph - self.remove_node(constraint) - # - # For each point that was adjacent to the constraint, - for point in adj_points: - - # Remove the constraint references stored in the Point object - point.geo_cons.remove(constraint) - - clusters = {} - - # Compile all other constraints - for point in adj_points: - - # Only compile this constraint if it is not a member of a compiled cluster - compile_point_constraint = True - for k in clusters.keys(): - if point in clusters[k]: - compile_point_constraint = False - break - if not compile_point_constraint or len(point.geo_cons) == 0: - continue - - # Create a new instance of EquationData to separate this point's first constraint's data reference from the - # disconnected clusters - point.geo_cons[0].data = EquationData() - - # Get a list of all the nodes now connecting to this constraint - descendants = networkx.algorithms.descendants(self, point) - for descendant in descendants: - if isinstance(descendant, GeoCon) and descendant is not point.geo_cons[0]: - descendant.data = point.geo_cons[0].data - - # Compile, solve, update, and add this constraint to the list of cluster parents - constraint = point.geo_cons[0] - self.analyze(constraint) - self.multicompile(constraint) - self.multisolve_and_update(constraint) - clusters[constraint] = networkx.algorithms.descendants(self, constraint) - pass - - def get_points_to_fix(self, source: GeoCon) -> typing.List[Point]: - - earliest_point = None - earliest_index = None - - points_to_fix = [] - - for node in networkx.dfs_preorder_nodes(self, source=source): - if not isinstance(node, Point): - continue - - # If the node has already set by the user to be fixed, append this point to the list and continue - if node.fixed(): - points_to_fix.append(node) - continue - - # If this point was added earlier than "earliest_point," then set this point to "earliest_point" - if earliest_point is None: - earliest_point = node - earliest_index = self.points.index(node) - else: - if self.points.index(node) > earliest_index: - continue - earliest_point = node - earliest_index = self.points.index(node) - - # If no fixed points were added by the user, we fix exactly one point (the earliest point) to satisfy - # the required degrees of freedom for a rigid 2-D body - if len(points_to_fix) == 0: - points_to_fix = [earliest_point] - - return points_to_fix - - @staticmethod - def fix_points_weak(points_to_fix: typing.List[Point]): - for point in points_to_fix: - if point.fixed(): - continue - else: - point.set_fixed_weak(True) - return points_to_fix - - @staticmethod - def free_points(points_to_free: typing.List[Point]): - for point in points_to_free: - point.set_fixed(False) - return points_to_free - - def analyze(self, constraint: GeoCon): - - constraint.data.clear() - params = self.get_params() - points = self.get_points(source_node=constraint) - for point in points: - point.set_fixed_weak(False) - points_to_fix = self.get_points_to_fix(source=constraint) - fixed_points = self.fix_points_weak(points_to_fix) - strong_constraints = self.get_strong_constraints(source_node=constraint) - strong_equations = self.get_strong_equations(strong_constraints, source_node=constraint) - abs_angle_constraints = [cnstr for cnstr in strong_constraints if isinstance(cnstr, AbsAngleConstraint)] - weak_constraints = [] - - dof = 2 * len(points) - len(strong_equations) - if len(fixed_points) == 0: - dof -= 2 - else: - dof -= 2 * len(fixed_points) - - # Only need to add an absolute angle constraint if one does not yet exist and if there is only one fixed point - if len(abs_angle_constraints) == 0 and len(fixed_points) == 1: - dof -= 1 - - # for point in points: - # if point.fixed_weak(): - # print(f"{point = } fixed weak") - # print(f"{len(points) = }, {len(strong_equations) = }, {len(fixed_points) = }, {len(abs_angle_constraints) = }") - - if dof < 0: - raise OverConstrainedError("System is over-constrained") - - points_may_need_distance = [] - points_may_need_angle = [] - a4_constraints = [] - - for node in networkx.dfs_preorder_nodes(self, source=constraint): - - if isinstance(node, Point) and not any(["d" in cnstr.kind for cnstr in self.adj[node]]): - if node not in points_may_need_distance: - points_may_need_distance.append(node) - - if isinstance(node, Point) and not any(["a" in cnstr.kind for cnstr in self.adj[node]]): - if node not in points_may_need_angle: - points_may_need_angle.append(node) - - if isinstance(node, GeoCon) and node.kind == "a4": - a4_constraints.append(node) - - for point in points_may_need_distance: - - # Weak distance constraint algorithm - - if dof == 0: - break - - # If this was the last point added, use the second-to-last point as p2. Otherwise, use the next point. - if self.points.index(point) == len(self.points) - 1: - p2 = self.points[-2] - else: - p2 = self.points[self.points.index(point) + 1] - - # Add a weak distance constraint between this point and the next candidate point, chosen somewhat - # arbitrarily for now - dist_constraint_weak = DistanceConstraintWeak(point, p2, "dcw1") - weak_constraints.append(dist_constraint_weak) - dof -= 1 - - for cnstr in a4_constraints: - - # A4 constraint algorithm - - if dof == 0: - break - - # 2. For each "a4," check that each of the subsets {{1,2},{3,4}} and {{2,3},{1,4}} is constrained w/ "a3" - # a3_constrained_points = [] - # for point in [cnstr.p1, cnstr.p2, cnstr.p3, cnstr.p4]: - # if any([c.kind == "a3" for c in point.geo_cons]): - # a3_constrained_points.append(point) - - add_a3 = True - analyze_angle_loop = True - - # 2 (cont). If less than 3 points in the a4 constraint have an a3 constraint, then we know for sure that an - # a3 constraint connecting two adjacent sides of the a4 constraint will be required - # if len(a3_constrained_points) < 3: - # analyze_angle_loop = False - - point_pairs = { - 0: [cnstr.p1, cnstr.p2], - 1: [cnstr.p3, cnstr.p4], - 2: [cnstr.p2, cnstr.p3], - 3: [cnstr.p1, cnstr.p4] - } - - common_neighbors = [] - - for idx, pair in point_pairs.items(): - common_neighbors.append( - [n for n in networkx.common_neighbors(self, pair[0], pair[1]) if n.kind == "a3"] - ) - - point_pair_n_constraints = {idx: len(neighbors) for idx, neighbors in enumerate(common_neighbors)} - - pair_combos = [[0, 2], [0, 3], [1, 2], [1, 3]] - found_combo = False - - for pair_combo in pair_combos: - if point_pair_n_constraints[pair_combo[0]] > 0 and point_pair_n_constraints[pair_combo[1]] > 0: - if len(set(common_neighbors[pair_combo[0]]).intersection(common_neighbors[pair_combo[1]])) > 0: - add_a3 = False - break - found_combo = True - - if not found_combo: - analyze_angle_loop = False - - # This next block should all go inside the next if statement - cycles = [cycle for cycle in networkx.simple_cycles(self) if all([isinstance(node, Point) or ( - isinstance(node, GeoCon) and node.kind in ["a3", "a4"]) for node in cycle])] - # for cycle in cycles: - # print(f"{cycle = }") - - # 4. Angle loop closure - if add_a3 and analyze_angle_loop: - visited_nodes = [] - - if not add_a3: - continue - - # Add the weak a3 constraint - candidate_point_lists = [[cnstr.p2, cnstr.p1, cnstr.p3], [cnstr.p4, cnstr.p2, cnstr.p1]] - check_dist_point_lists = [[cnstr.p1, cnstr.p3], [cnstr.p2, cnstr.p4]] - for candidate_point_list, check_dist_point_list in zip(candidate_point_lists, check_dist_point_lists): - common_dist_neighbors = ( - [n for n in networkx.common_neighbors( - self, check_dist_point_list[0], check_dist_point_list[1]) if "d" in n.kind] - ) - if len(common_dist_neighbors) == 0: - continue - - weak_a3_constraint = RelAngle3ConstraintWeak(*candidate_point_list, name="ra3w") - weak_constraints.append(weak_a3_constraint) - break - - for point in points_may_need_angle: - - # Weak relative angle constraint algorithm - - if dof == 0: - break - - for cnstr in point.geo_cons: - if isinstance(cnstr, DistanceConstraint) or isinstance(cnstr, DistanceConstraintWeak): - start_point = point - p2 = cnstr.p2 if cnstr.p2 is not point else cnstr.p1 - end_point = None - - for sub_cnstr in p2.geo_cons: - if sub_cnstr is cnstr: - continue - if isinstance(sub_cnstr, RelAngle3Constraint) or isinstance(sub_cnstr, RelAngle3ConstraintWeak): - end_point = sub_cnstr.p1 - elif (isinstance(sub_cnstr, Perp3Constraint) or isinstance(sub_cnstr, Parallel3Constraint) or - isinstance(sub_cnstr, AntiParallel3Constraint)): - end_point = sub_cnstr.p1 - if end_point is not None: - rel_angle3_constraint_weak = RelAngle3ConstraintWeak(start_point, p2, end_point, "ra3") - weak_constraints.append(rel_angle3_constraint_weak) - dof -= 1 - break - - for sub_cnstr in p2.geo_cons: - if sub_cnstr is cnstr: - continue - - if isinstance(sub_cnstr, DistanceConstraint) or isinstance(sub_cnstr, DistanceConstraintWeak): - end_point = sub_cnstr.p1 if sub_cnstr.p1 not in [start_point, p2] else sub_cnstr.p2 - - if end_point is not None: - rel_angle3_constraint_weak = RelAngle3ConstraintWeak(start_point, p2, end_point, "ra3") - weak_constraints.append(rel_angle3_constraint_weak) - dof -= 1 - break - - # Determine the absolute angle constraint to add to eliminate the rotational degree of freedom, if necessary - if len(abs_angle_constraints) == 0 and len(fixed_points) == 1: - start_point = fixed_points[0] - - # Sort the points in the cluster by the order they were added - candidate_points = [pt for pt in self.points if pt in points] - - # Choose the first point in the cluster by order added as the end point, unless the first point is the start - # point. In this case, choose the second point instead. - if candidate_points[0] == start_point: - end_point = candidate_points[1] - else: - end_point = candidate_points[0] - abs_angle_constraint = AbsAngleConstraintWeak(start_point, end_point, "aa1") - weak_constraints.append(abs_angle_constraint) - # print(f"{abs_angle_constraint = }, {p1 = }, {p3 = }") - - weak_equations = [] - for cnstr in weak_constraints: - weak_equations.extend(cnstr.equations) - - points_to_vary = [point for point in points if (not point.fixed() and not point.fixed_weak())] - self.add_points_as_variables(constraint.data, points_to_vary, params=params) - - for eq in strong_equations: - constraint.data.equations.append(eq) - - for cnstr in strong_constraints: - constraint.data.arg_idx_array.extend(cnstr.get_arg_idx_array(params)) - constraint.data.geo_cons.append(cnstr) - - for eq in weak_equations: - constraint.data.equations.append(eq) - - for cnstr in weak_constraints: - constraint.data.arg_idx_array.extend(cnstr.get_arg_idx_array(params)) - constraint.data.geo_cons.append(cnstr) - - pass - - @staticmethod - def verify_constraint_addition(info): - if np.any(info.fun_val > 1e-6): - raise ValueError("Could not solve the constraint problem within tolerance after constraint addition") - - def get_strong_constraints(self, source_node: Point or GeoCon): - - strong_constraints = [] - for node in networkx.dfs_preorder_nodes(self, source=source_node): - if isinstance(node, GeoCon): - strong_constraints.append(node) - - return strong_constraints - - def get_strong_equations(self, strong_constraints: typing.List[GeoCon] = None, - source_node: Point or GeoCon = None): - strong_equations = [] - - strong_constraints = self.get_strong_constraints( - source_node) if strong_constraints is None else strong_constraints - - for constraint in strong_constraints: - strong_equations.extend(constraint.equations) - - return strong_equations - - def get_points(self, source_node: Point or GeoCon): - points = [] - - for node in networkx.dfs_preorder_nodes(self, source=source_node): - if not isinstance(node, Point): - continue - - points.append(node) - - return points - - def get_params(self): - params = [] - - for point in self.points: - params.extend([point.x(), point.y()]) - - # for node in networkx.dfs_preorder_nodes(self, source=source_node): - # # if isinstance(node, Point): - # # params.extend([node.x(), node.y()]) - # if isinstance(node, GeoCon): - # if node.param() is None: - # continue - # params.append(node.param()) - - for constraint_param in self.constraint_params: - params.append(constraint_param) - - return params - - def get_param_values(self): - return np.array([p.value() for p in self.get_params()]) - - @staticmethod - def add_variable(equation_data: EquationData, variable: Param, - params: typing.List[Param]): - - var_index = params.index(variable) - - if var_index not in equation_data.variable_pos: - equation_data.variable_pos.append(var_index) - - def add_points_as_variables(self, equation_data: EquationData, points: typing.List[Point], - params: typing.List[Param]): - param_list = [] - for point in points: - param_list.extend([point.x(), point.y()]) - - for param in param_list: - self.add_variable(equation_data, param, params) - - @staticmethod - def compile_equation_for_entity_or_constraint(constraint: GeoCon, method: str = "lm"): - - def equation_system(x: np.ndarray, - start_param_vec: np.ndarray, intermediate_param_vec: np.ndarray): - final_param_vec = jnp.array([v for v in intermediate_param_vec]) - for idx, pos in enumerate(constraint.data.variable_pos): - final_param_vec = final_param_vec.at[pos].set(x[idx]) - - func_outputs = jnp.array([cnstr(*[start_param_vec[w[0]] if w[1] == 0 - else intermediate_param_vec[w[0]] if w[1] == 1 - else final_param_vec[w[0]] for w in arg_idx_array]) - for cnstr, arg_idx_array in zip(constraint.data.equations, - constraint.data.arg_idx_array) - ]) - - return func_outputs - - constraint.data.root_finders[method] = PymeadRootFinder(equation_system, method=method) - - def multicompile(self, constraint: GeoCon): - self.compile_equation_for_entity_or_constraint(constraint, method="lm") - self.compile_equation_for_entity_or_constraint(constraint, method="hybr") - - def solve(self, constraint: GeoCon, method: str): - params = self.get_params() - x0 = np.array([params[x_pos].value() for x_pos in constraint.data.variable_pos]) - v = np.array([p.value() for p in params]) - w = deepcopy(v) - x, info = constraint.data.root_finders[method].solve(x0, v, w) - return x, info - - def multisolve_and_update(self, constraint: GeoCon): - x, info = self.solve(constraint, method="lm") - - try: - self.verify_constraint_addition(info) - except ValueError: - # Update the points anyway and try to solve using the other root-finding method - self.update_points(constraint, new_x=x) - x, info = self.solve(constraint, method="hybr") - try: - self.verify_constraint_addition(info) - except ValueError: - raise ValueError("Could not converge the solution within tolerance using both root-finding methods") - - self.update_points(constraint, new_x=x) - - def compile_and_solve_first_constraint_of_all_clusters(self, constraint_list: typing.List[GeoCon]): - """ - Method used in the loading process of the GUI to compile only the first constraint (in order of addition to - the geometry collection) of each unconnected cluster in the ``ConstraintGraph``. Executing this method rather - than compiling/solving the non-linear system of equations for each constraint added during the load process - significantly decreases the load times for a ``GeometryCollection``. To be more precise, using this method - allows the load time to scale proportional to the number of unconnected clusters rather than the total number - of constraints. Compiling only once per cluster is possible during the loading process because the points - are already in the correct positions following the last equation solve prior to saving the ``.jmea`` file. - - Parameters - ---------- - constraint_list: typing.List[GeoCon] - List of constraints in the ``GeometryCollection`` in the order they were added - """ - # Compile the first constraint - constraint = constraint_list[0] - self.multicompile(constraint) - self.multisolve_and_update(constraint) - - clusters = {constraint: networkx.algorithms.descendants(self, constraint)} - - # Compile all other constraints - for constraint in constraint_list[1:]: - - # Only compile this constraint if it is not a member of a compiled cluster - compile_this_constraint = True - for k in clusters.keys(): - if constraint in clusters[k]: - compile_this_constraint = False - break - if not compile_this_constraint: - continue - - # Compile, solve, update, and add this constraint to the list of cluster parents - self.multicompile(constraint) - self.multisolve_and_update(constraint) - clusters[constraint] = networkx.algorithms.descendants(self, constraint) - - def update_points(self, constraint: GeoCon, new_x: np.ndarray): - """ - Updates the variable points and parameters in the constraint system - - Parameters - ---------- - constraint: GeoCon - Normally, this is the constraint that had its parameter updated and therefore called the ``solve()`` method. - - new_x: np.ndarray - Vector of updated parameter values - - Returns - ------- - - """ - params = self.get_params() - for idx, x in zip(constraint.data.variable_pos, new_x): - params[idx].set_value(x) - - curves_to_update = [] - for point in self.points: - if point.canvas_item is not None: - point.canvas_item.updateCanvasItem(point.x().value(), point.y().value()) - - for curve in point.curves: - if curve not in curves_to_update: - curves_to_update.append(curve) - - airfoils_to_update = [] - for curve in curves_to_update: - if curve.airfoil is not None and curve.airfoil not in airfoils_to_update: - airfoils_to_update.append(curve.airfoil) - curve.update() - - for airfoil in airfoils_to_update: - airfoil.update_coords() - if airfoil.canvas_item is not None: - airfoil.canvas_item.generatePicture() - - for node in networkx.dfs_preorder_nodes(self, source=constraint): - if isinstance(node, GeoCon) and node.canvas_item is not None: - node.canvas_item.update() - - -class OverConstrainedError(Exception): - pass diff --git a/pymead/core/constraint_solver.py b/pymead/core/constraint_solver.py deleted file mode 100644 index 4e4f0593..00000000 --- a/pymead/core/constraint_solver.py +++ /dev/null @@ -1,170 +0,0 @@ -import time - -from jax import jacfwd, jit, debug -import jax.numpy as jnp -from math import floor -import numpy as np -from scipy.optimize import root -import matplotlib.pyplot as plt - -from pymead.core.param import Param -from pymead.core.point import Point, PointSequence -from pymead.core.bezier import Bezier - - -@jit -def measure_distance(x1: float, y1: float, x2: float, y2: float): - return jnp.hypot(x1 - x2, y1 - y2) - - -@jit -def measure_abs_angle(x1: float, y1: float, x2: float, y2: float): - return jnp.arctan2(y2 - y1, x2 - x1) - - -@jit -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) - - -@jit -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))) - - -@jit -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)) - - -@jit -def f(x: np.ndarray): - f1 = (1.0 - x[0]) ** 2 + (1.0 - x[1]) ** 2 - 3.0 ** 2 - f2 = (x[0] - x[2]) ** 2 + (x[1] - x[3]) ** 2 - 4.0 ** 2 - f3 = (x[2] - 1.0) ** 2 + (x[3] - 1.0) ** 2 - 5.0 ** 2 - f4 = jnp.arctan2(x[3] - x[1], x[2] - x[0]) - np.pi / 4 - return jnp.array([f1, f2, f3, f4]) - - -@jit -def f_curvature(x: np.ndarray, x_pos: np.ndarray, xy: np.ndarray, n: np.ndarray): - - for idx, pos in enumerate(x_pos): - xy = xy.at[pos[0], pos[1]].set(x[idx]) - - debug.print(" {xy} ", xy=xy) - - phi1 = measure_abs_angle(xy[2, 0], xy[2, 1], xy[1, 0], xy[1, 1]) - phi2 = measure_abs_angle(xy[2, 0], xy[2, 1], xy[3, 0], xy[3, 1]) - theta1 = measure_abs_angle(xy[1, 0], xy[1, 1], xy[0, 0], xy[0, 1]) - theta2 = measure_abs_angle(xy[3, 0], xy[3, 1], xy[4, 0], xy[4, 1]) - psi1 = theta1 - phi1 - psi2 = theta2 - phi2 - # phi_rel = measure_rel_angle3(xy[1, 0], xy[1, 1], xy[2, 0], xy[2, 1], xy[3, 0], xy[3, 1]) - phi_rel = (phi1 - phi2) % (2 * jnp.pi) - # psi1 = measure_rel_angle3(xy[0, 0], xy[0, 1], xy[1, 0], xy[1, 1], xy[2, 0], xy[2, 1]) - # psi2 = measure_rel_angle3(xy[2, 0], xy[2, 1], xy[3, 0], xy[3, 1], xy[4, 0], xy[4, 1]) - 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) - f1 = kappa1 - kappa2 - f2 = phi_rel - jnp.pi - - debug.print(" {f1} ", f1=f1) - f_arr = jnp.array([f1, f2]) - return f_arr - - -@jit -def fp(x: np.ndarray): - jac = jacfwd(f)(x) - return jac - - -@jit -def fp_curvature(x: np.ndarray, pos: np.ndarray, xy: np.ndarray, n: np.ndarray): - return jacfwd(f_curvature)(x, pos, xy, n) - - -def main(): - x0 = np.array([4.87, -0.1, 16.221, 0.756]) - t1 = time.perf_counter() - for _ in range(3500): - res = root(f, x0=x0, jac=fp) - t2 = time.perf_counter() - print(f"Root found 3500 times in {t2 - t1} seconds. Estimated operations per second: {floor(3500 / (t2 - t1))}") - fig, ax = plt.subplots() - ax.plot([1.0, res.x[0], res.x[2]], [1.0, res.x[1], res.x[3]], ls="--", marker="o") - ax.set_aspect("equal") - plt.show() - - -def main2(): - x0 = np.array([0.25, 0.8]) - res = root(f_curvature, x0=x0, jac=fp_curvature) - print(f"{res = }") - - -def main3(): - xy_original = np.array([ - [0.0, 0.0], - [0.3, 0.5], - [0.45, 0.55], - [0.8, 0.3], - [1.1, -0.1] - ]) - original_point_seq = PointSequence.generate_from_array(xy_original) - bez1 = Bezier(point_sequence=original_point_seq[:3]) - bez2 = Bezier(point_sequence=original_point_seq[2:][::-1]) - - fig, ax = plt.subplots() - ax.plot(xy_original[:, 0], xy_original[:, 1], ls="--", marker="o", mfc="grey", mec="grey", - color="grey") - bez1_original_data = bez1.evaluate() - bez2_original_data = bez2.evaluate() - ax.plot(bez1_original_data.xy[:, 0], bez1_original_data.xy[:, 1], color="steelblue") - ax.plot(bez2_original_data.xy[:, 0], bez2_original_data.xy[:, 1], color="indianred") - - pos = np.array([[3, 1], [4, 1]]) - res = root(f_curvature, x0=np.array([xy_original[p[0], p[1]] for p in pos]), jac=fp_curvature, - args=(pos, xy_original, np.array([bez1.degree, bez2.degree]))) - print(f"{res = }") - - for idx, x in enumerate(res.x): - xy_original[pos[idx, 0], pos[idx, 1]] = x - - print(f"{xy_original = }") - - new_point_seq = PointSequence.generate_from_array(xy_original) - bez1 = Bezier(point_sequence=new_point_seq[:3]) - bez2 = Bezier(point_sequence=new_point_seq[2:]) - bez1_new_data = bez1.evaluate() - bez2_new_data = bez2.evaluate() - - ax.plot(xy_original[:, 0], xy_original[:, 1], ls="--", marker="o", mfc="magenta", mec="magenta", - color="magenta") - ax.plot(bez1_new_data.xy[:, 0], bez1_new_data.xy[:, 1], color="cornflowerblue") - ax.plot(bez2_new_data.xy[:, 0], bez2_new_data.xy[:, 1], color="hotpink") - - print(f"{bez1_new_data.R[-1] = }, {bez2_new_data.R[0] = }") - - plt.show() - - -def main4(): - - from pymead.core.constraints import GCS - p1 = Point(0.0, 0.0) - p2 = Point(0.5, 0.2) - dist = Param(0.5, "dist") - gcs = GCS(parent=p1) - gcs.add_distance_constraint(p1, p2, dist) - res = gcs.solve() - print(f"{res = }") - - -if __name__ == "__main__": - main3() diff --git a/pymead/core/gcs2.py b/pymead/core/gcs.py similarity index 100% rename from pymead/core/gcs2.py rename to pymead/core/gcs.py diff --git a/pymead/core/geometry_collection.py b/pymead/core/geometry_collection.py index a7a4231e..a6d4f4d4 100644 --- a/pymead/core/geometry_collection.py +++ b/pymead/core/geometry_collection.py @@ -7,8 +7,7 @@ from pymead.core.airfoil import Airfoil from pymead.core.bezier import Bezier from pymead.core.constraints import * -from pymead.core.constraint_graph import OverConstrainedError -from pymead.core.gcs2 import GCS2 +from pymead.core.gcs import GCS2 from pymead.core.dimensions import LengthDimension, AngleDimension, Dimension from pymead.core.mea import MEA from pymead.core.pymead_obj import DualRep, PymeadObj @@ -796,13 +795,14 @@ def add_constraint(self, constraint: GeoCon, assign_unique_name: bool = True, ** if isinstance(constraint, AntiParallel3Constraint) or isinstance(constraint, Perp3Constraint): points_solved = self.gcs.solve(constraint) self.gcs.update_canvas_items(points_solved) - except (OverConstrainedError, ValueError) as e: + except ValueError as e: self.remove_pymead_obj(constraint) self.clear_selected_objects() if self.gui_obj is not None: - self.gui_obj.showColoredMessage("Constraint cluster is over-constrained. Removing constraint...", - 4000, "#eb4034") - raise e + # self.gui_obj.showColoredMessage("Constraint cluster is over-constrained. Removing constraint...", + # 4000, "#eb4034") + self.gui_obj.disp_message_box(str(e), message_mode="error") + return return constraint def get_dict_rep(self):