diff --git a/src/ceasiompy/utils/cpacs_preprocessing.py b/src/ceasiompy/utils/cpacs_preprocessing.py new file mode 100644 index 000000000..3c264a56e --- /dev/null +++ b/src/ceasiompy/utils/cpacs_preprocessing.py @@ -0,0 +1,276 @@ +from pathlib import Path +import xml.etree.ElementTree as ET +import math + +from cpacspy.cpacspy import CPACS + +def _strip_ns(tag: str) -> str: + """Strip XML namespace from a tag name.""" + return tag.split("}", 1)[-1] if "}" in tag else tag + + +def _find_child(parent: ET.Element, name: str): + """Find direct child by local (namespace-free) tag name.""" + for child in parent: + if _strip_ns(child.tag) == name: + return child + return None + + +def _find_path(root: ET.Element, parts: list[str]): + """Follow a direct-child path, e.g. ['vehicles', 'profiles', 'wingAirfoils'].""" + cur = root + for p in parts: + cur = _find_child(cur, p) + if cur is None: + return None + return cur + + +def _parse_vector(text: str | None) -> list[float]: + """Parse CPACS mapType='vector' text into list[float].""" + if not text: + return [] + out = [] + for t in text.replace("\n", "").split(";"): + t = t.strip() + if t: + out.append(float(t)) + return out + + +def _to_vector(values: list[float]) -> str: + """Convert list[float] back to CPACS vector string format.""" + return ";".join(f"{v:.8g}" for v in values) + ";" + + +def _interp_at_x(x0, y0, z0, x1, y1, z1, x_target): + if math.isclose(x1, x0, rel_tol=0.0, abs_tol=1e-15): + return x_target, y0, z0 + t = (x_target - x0) / (x1 - x0) + return x_target, y0 + t * (y1 - y0), z0 + t * (z1 - z0) + + +def _clip_lower_te_to_le(x_vals, y_vals, z_vals, x_cut): + # Lower branch follows CPACS order TE -> LE. + # Remove aft part (x > x_cut), then add one interpolated cut point at x_cut. + i_in = next((i for i, x in enumerate(x_vals) if x <= x_cut), None) + if i_in is None: + return [x_cut], [y_vals[-1]], [z_vals[-1]] + if i_in == 0: + return x_vals[:], y_vals[:], z_vals[:] + + xi, yi, zi = _interp_at_x( + x_vals[i_in - 1], + y_vals[i_in - 1], + z_vals[i_in - 1], + x_vals[i_in], + y_vals[i_in], + z_vals[i_in], + x_cut, + ) + return [xi] + x_vals[i_in:], [yi] + y_vals[i_in:], [zi] + z_vals[i_in:] + + +def _clip_upper_le_to_te(x_vals, y_vals, z_vals, x_cut): + # Upper branch follows CPACS order LE -> TE. + # Remove aft part (x > x_cut), then add one interpolated cut point at x_cut. + i_out = next((i for i, x in enumerate(x_vals) if x > x_cut), None) + if i_out is None: + return x_vals[:], y_vals[:], z_vals[:] + if i_out == 0: + return [x_cut], [y_vals[0]], [z_vals[0]] + + xi, yi, zi = _interp_at_x( + x_vals[i_out - 1], + y_vals[i_out - 1], + z_vals[i_out - 1], + x_vals[i_out], + y_vals[i_out], + z_vals[i_out], + x_cut, + ) + return x_vals[:i_out] + [xi], y_vals[:i_out] + [yi], z_vals[:i_out] + [zi] + + +def cut_trailing_edge_airfoils(airfoils, te_cut): + """ + Cut trailing edge at x_cut = x_te - te_cut for each airfoil and rebuild + each profile as: + TE -> lower side -> LE -> upper side -> TE. + + Keep distinct lower/upper TE cut points and close the polyline by + repeating TE-lower at the end. + Coordinates are then renormalized so x_cut -> 1 for each profile. + """ + tol = 1e-12 + + def clip_lower(points, x_cut): + """Clip lower branch ordered TE->LE to x <= x_cut.""" + i_in = next((i for i, (xp, _, _) in enumerate(points) if xp <= x_cut), None) + if i_in is None: + _, yp, zp = points[-1] + return [(x_cut, yp, zp)] + if i_in == 0: + return points[:] + cut = _interp_at_x(*points[i_in - 1], *points[i_in], x_cut) + return [cut] + points[i_in:] + + def clip_upper(points, x_cut): + """Clip upper branch ordered LE->TE to x <= x_cut.""" + i_out = next((i for i, (xp, _, _) in enumerate(points) if xp > x_cut), None) + if i_out is None: + return points[:] + if i_out == 0: + _, yp, zp = points[0] + return [(x_cut, yp, zp)] + cut = _interp_at_x(*points[i_out - 1], *points[i_out], x_cut) + return points[:i_out] + [cut] + + out = [] + for airfoil in airfoils: + points = list(zip(airfoil["x"], airfoil["y"], airfoil["z"])) + if len(points) < 2: + out.append(dict(airfoil)) + continue + + x_te = points[0][0] + x_cut = x_te - te_cut + if x_cut <= 0.0: + raise ValueError( + f"Invalid te_cut={te_cut} for airfoil {airfoil.get('uid')}: " + f"x_te={x_te}, x_cut={x_cut}." + ) + + i_le = min(range(len(points)), key=lambda i: points[i][0]) + lower = clip_lower(points[: i_le + 1], x_cut) # TE -> LE + upper = clip_upper(points[i_le:], x_cut) # LE -> TE + + same_le = ( + math.isclose(lower[-1][0], upper[0][0], rel_tol=0.0, abs_tol=tol) + and math.isclose(lower[-1][1], upper[0][1], rel_tol=0.0, abs_tol=tol) + and math.isclose(lower[-1][2], upper[0][2], rel_tol=0.0, abs_tol=tol) + ) + merged = lower + (upper[1:] if same_le else upper) + + # Close polyline at the trailing edge: + # ... -> TE_upper -> TE_lower + te_lower = lower[0] + merged.append(te_lower) + + scale = 1.0 + new_x = [xp * scale for xp, _, _ in merged] + new_y = [yp * scale for _, yp, _ in merged] + new_z = [zp * scale for _, _, zp in merged] + import matplotlib.pyplot as plt + plt.figure() + for i in range(len(new_x) - 1): + plt.plot(new_x[i : i + 2], new_z[i : i + 2], "-k", linewidth=1.0) + plt.plot(new_x[0], new_z[0], "go", markersize=6, label="start") + plt.plot(new_x[-1], new_z[-1], "ro", markersize=6, label="end") + plt.title(f"Airfoil {airfoil.get('uid')} after TE cut") + plt.xlabel('x-coordinate') + plt.ylabel('z-coordinate') + plt.grid() + plt.axis('equal') + plt.legend() + plt.show() + + new_airfoil = dict(airfoil) + new_airfoil["x"] = new_x + new_airfoil["y"] = new_y + new_airfoil["z"] = new_z + new_airfoil["points"] = list(zip(new_x, new_y, new_z)) + out.append(new_airfoil) + + return out + + +def load_wing_airfoils(root: ET.Element): + """Load wing airfoils from CPACS profile section into Python dictionaries.""" + airfoils_node = _find_path(root, ["vehicles", "profiles", "wingAirfoils"]) + if airfoils_node is None: + return [] + + data = [] + for wa in airfoils_node: + if _strip_ns(wa.tag) != "wingAirfoil": + continue + + uid = wa.attrib.get("uID") + point_list = _find_child(wa, "pointList") + if point_list is None: + continue + + x_node = _find_child(point_list, "x") + y_node = _find_child(point_list, "y") + z_node = _find_child(point_list, "z") + if x_node is None or y_node is None or z_node is None: + continue + + x = _parse_vector(x_node.text) + y = _parse_vector(y_node.text) + z = _parse_vector(z_node.text) + + data.append( + { + "uid": uid, + "x": x, + "y": y, + "z": z, + "points": list(zip(x, y, z)), + "x_node": x_node, + "y_node": y_node, + "z_node": z_node, + } + ) + return data + + +def write_airfoils_to_cpacs(tree: ET.ElementTree, airfoils_cut, cpacs_out: Path): + """ + Update / in the XML tree and save a new CPACS file. + Matching is done by wingAirfoil uID. + """ + root = tree.getroot() + airfoils_now = load_wing_airfoils(root) + by_uid = {af["uid"]: af for af in airfoils_cut} + updated = 0 + for af in airfoils_now: + uid = af["uid"] + if uid not in by_uid: + continue + new_airfoil = by_uid[uid] + af["x_node"].text = _to_vector(new_airfoil["x"]) + af["y_node"].text = _to_vector(new_airfoil["y"]) + af["z_node"].text = _to_vector(new_airfoil["z"]) + updated += 1 + + tree.write(cpacs_out, encoding="utf-8", xml_declaration=True) + return updated + + +def truncate_trailing_edge(cpacs: CPACS, te_cut: float = 0.4) -> CPACS: + tree = ET.parse(cpacs) + root = tree.getroot() + + # import the airfoils + airfoils = load_wing_airfoils(root) + print(f"Loaded {len(airfoils)} airfoils") + + # te_cut -> fraction of chord removed at TE (x_cut = 1 - te_cut) + airfoils_cut = cut_trailing_edge_airfoils(airfoils, te_cut) + + + + + # save as a new CPACS + cpacs_out = cpacs.with_name(f"{cpacs.stem}_tecut_straight_New.xml") + n_updated = write_airfoils_to_cpacs(tree, airfoils_cut, cpacs_out) + print(f"Updated {n_updated} airfoils in CPACS") + print(f"Saved new CPACS: {cpacs_out}") + return cpacs_out + + +if __name__ == "__main__": + truncate_trailing_edge(Path(__file__).with_name("00_ToolInput(1).xml"))