XFOIL optimization re-implemented

docs/source/todo.md

@@ -18,6 +18,13 @@ Planned feature additions
 - Graphical highlighting of constraint parameters and constraints
 - Add unit selection ComboBox
 - Add renaming to parameters
+- Add metadata to save files (date, pymead version, etc.)
+- Make the "Analysis" tab focused by default after an aerodynamics analysis (possibly implement a user option to
+  override this behavior)
+- Write the XFOIL/MSES analysis code using the same `CPUBoundProcess` architecture used by optimization
+- Re-implement downsampling (method stored in Airfoil --> MEA --> Optimization)
+- Create a custom context menu for all callback plots (similar to airfoil canvas context menu) with a subset of the
+  original context menu actions
 
 Refactoring
 -----------
@@ -28,6 +35,11 @@ Bug fixes
 ---------
 - May need to eliminate the use of modulo in constraint equations to make the variable 
   space smooth
+- Fix loading bar for new version
+- Fix excessively large font of "XFOIL" hyperlink text after running XFOIL
+- XFOIL cp not working?
+- Fix blank line in Objective/Constraint setup not reverting to background color after editing and erasing
+- Remove wave/viscous drag from XFOIL drag history plots (optimization)
 
 Testing
 -------

pymead/analysis/read_aero_data.py

@@ -23,18 +23,25 @@ def read_Cp_from_file_xfoil(fname: str):
     dict
       Dictionary containing 1-D *numpy* arrays for :math:`x`, :math:`y`, and :math:`C_p`
     """
+    y_in_header = False
+    header_line_idx = None
     with open(fname, "r") as f:
-        first_line = f.readline()
+        for idx, line in enumerate(f):
+            if "Cp" not in line:
+                continue
+            header_line_idx = idx
+            if "y" in line:
+                y_in_header = True
+            break
+        else:
+            raise ValueError("Could not detect the header line in the XFOIL Cp file")
 
     # Read in the Cp data for XFOIL versions 6.93 or 6.99 (other versions are untested)
-    if "y" in first_line:  # For XFOIL 6.99
-        df = pd.read_csv(fname, skiprows=3, names=['x', 'y', 'Cp'], sep='\s+', engine='python')
-        array_ = df.to_numpy()
-        return {'x': array_[:, 0], 'y': array_[:, 1], 'Cp': array_[:, 2]}
-    else:  # For XFOIL 6.93
-        df = pd.read_csv(fname, skiprows=1, names=['x', 'Cp'], sep='\s+', engine='python')
-        array_ = df.to_numpy()
-        return {'x': array_[:, 0], 'Cp': array_[:, 1]}
+    names = ["x", "y", "Cp"] if y_in_header else ["x", "Cp"]
+    cols = {"x": 0, "y": 1, "Cp": 2} if y_in_header else {"x": 0, "Cp": 1}
+    df = pd.read_csv(fname, skiprows=header_line_idx + 1, names=names, sep='\s+', engine='python')
+    array_ = df.to_numpy()
+    return {name: array_[:, cols[name]] for name in names}
 
 
 def read_aero_data_from_xfoil(fpath: str, aero_data: dict):

pymead/core/airfoil.py

@@ -145,6 +145,14 @@ def save_coords_selig_format(self, file_name: str) -> None:
     def convert_coords_to_shapely_format(coords: np.ndarray):
         return list(map(tuple, coords))
 
+    @staticmethod
+    def create_line_string(coords_shapely_format: list):
+        return LineString(coords_shapely_format)
+
+    @staticmethod
+    def create_shapely_polygon(line_string: LineString):
+        return Polygon(line_string)
+
     def compute_area(self, airfoil_polygon: Polygon):
         """Computes the area of the airfoil as the area of a many-sided polygon enclosed by the airfoil coordinates
         using the `shapely <https://shapely.readthedocs.io/en/stable/manual.html>`_ library.
@@ -189,7 +197,7 @@ def compute_thickness(airfoil_line_string: LineString, n_lines: int = 201):
 
         Returns
         =======
-        tuple or dict
+        dict
           The list of \(x\)-values used for the thickness distribution calculation, the thickness distribution, the
           maximum value of the thickness distribution, and, if :code:`return_max_thickness_location=True`,
           the :math:`x/c`-location of the maximum thickness value.

pymead/core/constraint_equations.py

@@ -6,27 +6,22 @@
 from jax import numpy as jnp
 
 
-@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)) % (2 * jnp.pi)
 
 
-@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_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)
 
 
-@jit
 def measure_point_line_distance_signed(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float):
     """
     Measures the signed distance from a point :math:`(x_3, y_3)` to a line defined by
@@ -57,22 +52,18 @@ def measure_point_line_distance_signed(x1: float, y1: float, x2: float, y2: floa
     return (x2 - x1) * (y1 - y3) - (x1 - x3) * (y2 - y1) / measure_distance(x1, y1, x2, y2)
 
 
-@jit
 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))
 
 
-@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 measure_data_bezier_curve_joint(xy: np.ndarray, n: np.ndarray):
     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])
@@ -100,7 +91,6 @@ def measure_data_bezier_curve_joint(xy: np.ndarray, n: np.ndarray):
     return data
 
 
-@jit
 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)
@@ -109,127 +99,104 @@ def radius_of_curvature_constraint(x1: float, y1: float, x2: float, y2: float, x
     return distance_constraint(x2, y2, x3, y3, Lc)
 
 
-@jit
 def empty_constraint_weak():
     return 0.0
 
 
-@jit
 def fixed_param_constraint(p_val: float, val: float):
     return p_val - val
 
 
-@jit
 def fixed_param_constraint_weak(new_val: float, old_val: float):
     return new_val - old_val
 
 
-@jit
 def fixed_x_constraint(x: float, val: float):
     return x - val
 
 
-@jit
 def fixed_x_constraint_weak(x_new: float, x_old: float):
     return x_new - x_old
 
 
-@jit
 def fixed_y_constraint(y: float, val: float):
     return y - val
 
 
-@jit
 def fixed_y_constraint_weak(y_new: float, y_old: float):
     return y_new - y_old
 
 
-@jit
 def distance_constraint(x1: float, y1: float, x2: float, y2: float, dist: float):
     return measure_distance(x1, y1, x2, y2) - dist
 
 
-@jit
 def distance_constraint_weak(x1_new: float, y1_new: float, x2_new: float, y2_new: float,
                              x1_old: float, y1_old: float, x2_old: float, y2_old: float):
     return (measure_distance(x1_new, y1_new, x2_new, y2_new) -
             measure_distance(x1_old, y1_old, x2_old, y2_old))
 
 
-@jit
 def abs_angle_constraint(x1: float, y1: float, x2: float, y2: float, angle: float):
     return measure_abs_angle(x1, y1, x2, y2) - angle
 
 
-@jit
 def abs_angle_constraint_weak(x1_new: float, y1_new: float, x2_new: float, y2_new: float,
                               x1_old: float, y1_old: float, x2_old: float, y2_old: float):
     return (measure_abs_angle(x1_new, y1_new, x2_new, y2_new) -
             measure_abs_angle(x1_old, y1_old, x2_old, y2_old))
 
 
-@jit
 def rel_angle3_constraint(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float, angle: float):
     return measure_rel_angle3(x1, y1, x2, y2, x3, y3) - angle
 
 
-@jit
 def rel_angle3_constraint_weak(x1_new: float, y1_new: float, x2_new: float, y2_new: float, x3_new: float, y3_new: float,
                                x1_old: float, y1_old: float, x2_old: float, y2_old: float, x3_old: float,
                                y3_old: float):
     return (measure_rel_angle3(x1_new, y1_new, x2_new, y2_new, x3_new, y3_new) -
             measure_rel_angle3(x1_old, y1_old, x2_old, y2_old, x3_old, y3_old))
 
 
-@jit
 def rel_angle4_constraint(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float, x4: float, y4: float,
                           angle: float):
     return measure_rel_angle4(x1, y1, x2, y2, x3, y3, x4, y4) - angle
 
 
-@jit
 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)
 
 
-@jit
 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)
 
 
-@jit
 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
 
 
-@jit
 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
 
 
-@jit
 def parallel3_constraint(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float):
     return measure_rel_angle3(x1, y1, x2, y2, x3, y3)
 
 
-@jit
 def parallel4_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)
 
 
-@jit
 def point_on_line_constraint(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float):
     return (y1 - y2) * (x3 - x1) - (y3 - y1) * (x1 - x2)
 
 
-@jit
 def points_equidistant_from_line_constraint_unsigned(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float,
                                                      x4: float, y4: float):
     return (measure_point_line_distance_unsigned(x1, y1, x2, y2, x3, y3) -
             measure_point_line_distance_unsigned(x1, y1, x2, y2, x4, y4))
 
 
-@jit
 def points_equidistant_from_line_constraint_signed(x1: float, y1: float, x2: float, y2: float, x3: float, y3: float,
                                                    x4: float, y4: float):
     return (measure_point_line_distance_signed(x1, y1, x2, y2, x3, y3) +

pymead/core/constraint_graph.py

@@ -70,13 +70,15 @@ 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()
@@ -602,7 +604,7 @@ def equation_system(x: np.ndarray,
 
             return func_outputs
 
-        constraint.data.root_finders[method] = PymeadRootFinder(jit(equation_system), method=method)
+        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")
@@ -708,10 +710,11 @@ def update_points(self, constraint: GeoCon, new_x: np.ndarray):
 
         for airfoil in airfoils_to_update:
             airfoil.update_coords()
-            airfoil.canvas_item.generatePicture()
+            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):
+            if isinstance(node, GeoCon) and node.canvas_item is not None:
                 node.canvas_item.update()
 
 

pymead/core/geometry_collection.py

@@ -328,18 +328,22 @@ def add_pymead_obj_by_ref(self, pymead_obj: PymeadObj, assign_unique_name: bool
 
         return pymead_obj
 
-    def remove_pymead_obj(self, pymead_obj: PymeadObj):
+    def remove_pymead_obj(self, pymead_obj: PymeadObj, promotion_demotion: bool = False):
         """
         Removes a pymead object from the geometry collection.
 
         Parameters
         ==========
         pymead_obj: PymeadObj
             Pymead object to remove
+
+        promotion_demotion: bool
+            When this flag is set to ``True``, the ``ValueError`` normally raised when directly deleting a ``Param``
+            associated with a ``GeoCon`` is ignored. Default: ``False``
         """
         # Type-specific actions
         if isinstance(pymead_obj, Param):
-            if len(pymead_obj.geo_cons) != 0:
+            if len(pymead_obj.geo_cons) != 0 and not promotion_demotion:
                 raise ValueError(f"Please delete each constraint associated with this parameter ({pymead_obj.geo_cons})"
                                  f" before deleting this parameter")
 
@@ -590,7 +594,7 @@ def promote_param_to_desvar(self, param: Param or str, lower: float or None = No
 
         # Remove the parameter
         if param.point is None:
-            self.remove_pymead_obj(param)
+            self.remove_pymead_obj(param, promotion_demotion=True)
 
         return desvar
 
@@ -637,7 +641,7 @@ def demote_desvar_to_param(self, desvar: DesVar):
         param.gcs.constraint_params[self.gcs.constraint_params.index(desvar)] = param
 
         # Remove the design variable
-        self.remove_pymead_obj(desvar)
+        self.remove_pymead_obj(desvar, promotion_demotion=True)
 
         return param
 

pymead/core/mea.py

@@ -30,7 +30,8 @@ def add_airfoil(self, airfoil: Airfoil):
     def remove_airfoil(self, airfoil: Airfoil):
         self.airfoils.remove(airfoil)
 
-    def get_coords_list(self):
+    def get_coords_list(self, downsample: bool = False, ds_max_points: int = 100, ds_curve_exp: float = 2.0):
+        # TODO: implement downsampling here
         return [airfoil.get_coords_selig_format() for airfoil in self.airfoils]
 
     def write_mses_blade_file(self,
@@ -84,7 +85,7 @@ def write_to_IGES(self, file_name: str):
         """
         bez_IGES_entities = [
             [BezierIGES(np.column_stack((c.P[:, 0], np.zeros(len(c.P)), c.P[:, 1]))) for c in a.curve_list]
-            for a in self.airfoils.values()]
+            for a in self.airfoils]
         entities_flattened = list(itertools.chain.from_iterable(bez_IGES_entities))
         iges_generator = IGESGenerator(entities_flattened)
         iges_generator.generate(file_name)

pymead/core/point.py

@@ -113,7 +113,8 @@ def request_move(self, xp: float, yp: float):
 
             for airfoil in airfoils_to_update:
                 airfoil.update_coords()
-                airfoil.canvas_item.generatePicture()
+                if airfoil.canvas_item is not None:
+                    airfoil.canvas_item.generatePicture()
         # else:
         #     start_param_vec = np.array([p.value() for p in self.gcs.params])
         #