diff --git a/.vsptime b/.vsptime new file mode 100644 index 000000000..bbfda6610 --- /dev/null +++ b/.vsptime @@ -0,0 +1 @@ +1770024222 \ No newline at end of file diff --git a/geometries/vspfiles/F-16_simplified_correct.vsp3 b/geometries/vspfiles/F-16_simplified_correct.vsp3 new file mode 100644 index 000000000..3ab81a816 --- /dev/null +++ b/geometries/vspfiles/F-16_simplified_correct.vsp3 @@ -0,0 +1,7729 @@ + + + 5 + + + OVFZHLRHIL + Vehicle + + + + + + + 0.000000000000000000e+00, 0.000000000000000000e+00, 0.000000000000000000e+00, + + + + 1.000000000000000021e-02, 5.000000000000000000e-01, 3.200000000000000067e-01, + + + + + + 1.000000000000000000e+00, 1.000000000000000000e+00, 1.000000000000000000e+00, + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 8 + + + ZMNJXJLAQX + Default + + + + + + + + + + + + + + DGWAMKODFI + Default + + + + + + + + + + + + + + BABLWABWCK + Default + + + + + + + + + + + + + + VBUBAQBJEL + Default + + + + + + + + + + + + + + NXJFDINLLV + Default + + + + + + + + + + + + + + ULVUEJKYIO + Default + + + + + + + + + + + + + + TUBCUKKVFF + Default + + + + + + + + + + + + + + LOHEESRGYH + Default + + + + + + + + + + + + + + 0 + 0 + 0 + 0 + + + + WQILYUQDOY + Mainbody + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Fuselage + 4 + 0 + NONE + + + FCUSHHANPS + + + + + + Grey Rubber + + + + GCWZCXYGWT + Default + + + + + + + + + + 0 + + + 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + + + + + + YZMMJBJUAC + Default + + + + + PZMZXIHOIK + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + XSec + + + CNREKWGWCL + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + MGJIGAOMXZ + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + XSec + + + PWHIDXPKGF + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + RCXLGWJKIB + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + XSec + + + LXXEAVJBAO + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + QPHGJHYQFO + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + XSec + + + VHVIPESJFA + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + MGQSDTUOAM + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + XSec + + + IBKAGYXTUF + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + XLRXIPTQSX + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + XSec + + + ULGWCVUMKU + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + + + + FCUSHHANPS + Main wing + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Wing + 5 + 0 + WQILYUQDOY + + + ZQJKRYZJUB + + + AISVNRZVST + + + LITPOTTKKV + + + BLUFJWYLOA + + + JMRQGETVOW + + + + + + Dark Grey Rubber + + + + DHWVTUEPBZ + Default + + + + + + + + + + 0 + + + 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + + + + + + YJFJBJQGOQ + Default + + + + + MWRRSLREKK + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSec + + 8 + 3 + 1, 5, 6, + + + + 25 + + + JBDXCQCNXT + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 11 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + UPRCKENJEJ + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSec + + 8 + 3 + 1, 5, 6, + + + + DZCDYNGLER + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 7 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + + + + ZQJKRYZJUB + Hardpoints + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Wing + 5 + 0 + FCUSHHANPS + + + + + Default + + + + LQQYFHESPO + Default + + + + + + + + + + 0 + + + 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + + + + + + MSVILOMHUO + Default + + + + + EVCLVOJSGN + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSec + + 8 + 3 + 1, 5, 6, + + + + 25 + + + WBHGPZZGUE + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 11 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + JMPCKKOPDD + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSec + + 8 + 3 + 1, 5, 6, + + + + AYYBJTRMZW + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 7 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + + + + AISVNRZVST + External Tank + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Fuselage + 4 + 0 + FCUSHHANPS + + + + + Obsidian + + + + CZJTBTTEUA + Default + + + + + + + + + + 0 + + + 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + + + + + + WCTRSSSSUV + Default + + + + + SIOTNLNPLL + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + XSec + + + RMLPDNNCHS + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + OOKHZGERYF + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + XSec + + + FAXASZNUPV + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + JJLZVPLQQM + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + XSec + + + YWSCRTZMXH + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + AGKRRNJMQD + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + XSec + + + JNXXZHWFMP + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + MYHODJQSCU + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + XSec + + + LANOGLDHOR + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + + + + LITPOTTKKV + Tail Elevon + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Wing + 5 + 0 + FCUSHHANPS + + + + + Dark Grey Rubber + + + + BHCPSHCIGA + Default + + + + + + + + + + 0 + + + 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + + + + GRJHLLAMWA + SS_CONT_0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 3 + + GDJDSPJGVF + + + + + + + + HKQEJXCAHC + Default + + + + + VZZECCQDOM + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSec + + 8 + 3 + 1, 5, 6, + + + + 25 + + + GYSRAEPVKT + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 11 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + BYIHMOFTXX + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSec + + 8 + 3 + 1, 5, 6, + + + + MGZCEUNSQE + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 7 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + + + + BLUFJWYLOA + Vetical Stablizer + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Wing + 5 + 0 + FCUSHHANPS + + + + + Dark Grey Rubber + + + + DIKYQJJCQS + Default + + + + + + + + + + 0 + + + 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + + + + BNEGVGCNTD + SS_CONT_0 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 3 + + GDJDSPJGVF + + + + + + + + GISDAJTVNZ + Default + + + + + JZLIGNEEJN + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSec + + 8 + 3 + 1, 5, 6, + + + + 25 + + + TQTWNAGKOG + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 11 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + YKQGKJFIWI + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSec + + 8 + 3 + 1, 5, 6, + + + + 25 + + + EGXCEKQUEB + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 11 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + RSUWYMZYWQ + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSec + + 8 + 3 + 1, 5, 6, + + + + DWMYKGYJMD + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 7 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + + + + JMRQGETVOW + Back-front wing connector + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Wing + 5 + 0 + FCUSHHANPS + + + + + Dark Grey Rubber + + + + IVMZEAXZYR + Default + + + + + + + + + + 0 + + + 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + + + + + + MANRLUPJOW + Default + + + + + EUOSBVBHPW + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSec + + 8 + 3 + 1, 5, 6, + + + + 25 + + + YSIVGJTROE + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 11 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + RSDJVHTDOB + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 2 + XSec + + 8 + 3 + 1, 5, 6, + + + + YJBTDTEVEM + Default + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 7 + XSecCurve + + 4 + 2 + 0, 2, + + + + + + + + + + + + Custom + 0.000000000000000000e+00, 0.000000000000000000e+00, 5.019607843137254832e-01, 1.000000000000000000e+00, + 0.000000000000000000e+00, 0.000000000000000000e+00, 5.019607843137254832e-01, 1.000000000000000000e+00, + 0.000000000000000000e+00, 0.000000000000000000e+00, 5.019607843137254832e-01, 1.000000000000000000e+00, + 0.000000000000000000e+00, 0.000000000000000000e+00, 5.019607843137254832e-01, 1.000000000000000000e+00, + 0.000000000000000000e+00 + + + wings + 5.019607843137254832e-01, 0.000000000000000000e+00, 5.019607843137254832e-01, 1.000000000000000000e+00, + 0.000000000000000000e+00, 0.000000000000000000e+00, 5.019607843137254832e-01, 1.000000000000000000e+00, + 0.000000000000000000e+00, 5.019607843137254832e-01, 5.019607843137254832e-01, 1.000000000000000000e+00, + 0.000000000000000000e+00, 0.000000000000000000e+00, 5.019607843137254832e-01, 1.000000000000000000e+00, + 0.000000000000000000e+00 + + + HS + 1.000000000000000000e+00, 0.000000000000000000e+00, 0.000000000000000000e+00, 1.000000000000000000e+00, + 5.019607843137254832e-01, 0.000000000000000000e+00, 5.019607843137254832e-01, 1.000000000000000000e+00, + 1.000000000000000000e+00, 0.000000000000000000e+00, 0.000000000000000000e+00, 1.000000000000000000e+00, + 1.000000000000000000e+00, 0.000000000000000000e+00, 0.000000000000000000e+00, 1.000000000000000000e+00, + 1.280000000000000000e+02 + + + VS + 0.000000000000000000e+00, 1.000000000000000000e+00, 0.000000000000000000e+00, 1.000000000000000000e+00, + 0.000000000000000000e+00, 0.000000000000000000e+00, 1.000000000000000000e+00, 1.000000000000000000e+00, + 5.019607843137254832e-01, 0.000000000000000000e+00, 5.019607843137254832e-01, 1.000000000000000000e+00, + 5.019607843137254832e-01, 0.000000000000000000e+00, 0.000000000000000000e+00, 1.000000000000000000e+00, + 1.280000000000000000e+02 + + + Transparent + 2.000000000000000111e-01, 2.000000000000000111e-01, 2.000000000000000111e-01, 1.000000000000000000e+00, + 8.000000000000000444e-01, 8.000000000000000444e-01, 8.000000000000000444e-01, 2.999999999999999889e-01, + 0.000000000000000000e+00, 0.000000000000000000e+00, 0.000000000000000000e+00, 1.000000000000000000e+00, + 0.000000000000000000e+00, 0.000000000000000000e+00, 0.000000000000000000e+00, 1.000000000000000000e+00, + 1.280000000000000000e+02 + + + Grey Rubber + 3.568627450980392246e-01, 3.568627450980392246e-01, 3.568627450980392246e-01, 1.000000000000000000e+00, + 1.000000000000000021e-02, 1.000000000000000021e-02, 1.000000000000000021e-02, 1.000000000000000000e+00, + 4.000000000000000222e-01, 4.000000000000000222e-01, 4.000000000000000222e-01, 1.000000000000000000e+00, + 0.000000000000000000e+00, 0.000000000000000000e+00, 0.000000000000000000e+00, 1.000000000000000000e+00, + 1.000000000000000000e+01 + + + Dark Grey Rubber + 1.019607843137254888e-01, 1.019607843137254888e-01, 1.019607843137254888e-01, 1.000000000000000000e+00, + 1.000000000000000021e-02, 1.000000000000000021e-02, 1.000000000000000021e-02, 1.000000000000000000e+00, + 4.000000000000000222e-01, 4.000000000000000222e-01, 4.000000000000000222e-01, 1.000000000000000000e+00, + 0.000000000000000000e+00, 0.000000000000000000e+00, 0.000000000000000000e+00, 1.000000000000000000e+00, + 1.000000000000000000e+01 + + + + 16 + + + + + + + + + + + + + + + + + + + + + + BISHRHCVIX + VSPAEROSettings + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + 0 + 0 + 0 + + + + + NRXSPSYWGK + + + JWGLCTITCB + CFDMeshSettings + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + UFNTJNYFCF + SurfaceIntersectSettings + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ELDNCTPUAI + CFDGridDensity + + + + + + + + + + + + + + + + + + + + + + + + + + ZQCESHPJLO + Default + + + + + + + GDJDSPJGVF + DefaultShell + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + _Al6061T6 + + + + ZYKOKSADKQ + DefaultBeam + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + _Al6061T6 + + + + + + UQSRVVPKRR + ClippingMgr + + + + + + + + + + + + + + + + + + + ZPXAYZFJPO + WaveDragSettings + + + + + + + + + + + + + + + + + + + + + + + + EOIUFUWCUB + ParasiteDragSettings + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0 + + + + + LWLNEGMDVR + AeroStructSettings + + + + + + + + 0 + + + All + Shown + Not_Shown + Set_0 + Set_1 + Set_2 + Set_3 + Set_4 + Set_5 + Set_6 + Set_7 + Set_8 + Set_9 + Set_10 + Set_11 + Set_12 + Set_13 + Set_14 + Set_15 + Set_16 + Set_17 + Set_18 + Set_19 + + + diff --git "a/src/app/01_\342\234\210\357\270\217_Geometry.py" "b/src/app/01_\342\234\210\357\270\217_Geometry.py" index 092db12ca..a5cb04776 100644 --- "a/src/app/01_\342\234\210\357\270\217_Geometry.py" +++ "b/src/app/01_\342\234\210\357\270\217_Geometry.py" @@ -12,10 +12,12 @@ # Imports import os import hashlib +import shutil import numpy as np import streamlit as st - -from ceasiompy.utils import get_wkdir +from stl import mesh +import plotly.graph_objects as go +import colorsys from cpacspy.cpacsfunctions import get_value from ceasiompy.utils.guiobjects import add_value from gmshairfoil2d.airfoil_func import get_airfoil_points @@ -39,6 +41,9 @@ render_openvsp_panel, convert_vsp3_to_cpacs, ) +from ceasiompy.stl2cpacs.func.split_stl_into_components import split_main +from ceasiompy.stl2cpacs.stl2cpacs import main as stl2cpacs_main +from ceasiompy.stl2cpacs.stl2cpacs import cpacs_component_detection from typing import Final from pathlib import Path @@ -47,7 +52,7 @@ from ceasiompy import log from constants import BLOCK_CONTAINER -from ceasiompy.utils.commonpaths import CPACS_FILES_PATH +from ceasiompy.utils.commonpaths import CPACS_FILES_PATH, get_wkdir from ceasiompy.utils.commonxpaths import ( AREA_XPATH, LENGTH_XPATH, @@ -68,28 +73,163 @@ # Methods -def _bootstrap_cli_geometry() -> None: - """Inject geometry passed from CLI into the same upload/conversion flow.""" +def _show_stl_mesh(stl_path: str, + key: str, + show_wireframe: bool | None = None + ) -> None: + """Interactive STL viewer with proper wireframe toggle refresh.""" + + # --- Cache STL loading only --- + @st.cache_data(show_spinner=False) + def load_stl(path): + return mesh.Mesh.from_file(path) + + stl_mesh = load_stl(stl_path) + + # Build unique vertices + face indices + vertices = stl_mesh.vectors.reshape(-1, 3) + verts_unique, index = np.unique(vertices, axis=0, return_inverse=True) + + x, y, z = verts_unique.T + i = index[0::3] + j = index[1::3] + k = index[2::3] + + fig = go.Figure() + + # ---- Surface mesh ---- + fig.add_trace( + go.Mesh3d( + x=x, + y=y, + z=z, + i=i, + j=j, + k=k, + color="grey", + opacity=1.0, + flatshading=True, + ) + ) - if st.session_state.get("_cli_geometry_bootstrapped"): - return + # ---- Wireframe ---- + if show_wireframe: + edge_x = [] + edge_y = [] + edge_z = [] + + for tri in stl_mesh.vectors: + for e in [(0, 1), (1, 2), (2, 0)]: + edge_x += [tri[e[0]][0], tri[e[1]][0], None] + edge_y += [tri[e[0]][1], tri[e[1]][1], None] + edge_z += [tri[e[0]][2], tri[e[1]][2], None] + + fig.add_trace( + go.Scatter3d( + x=edge_x, + y=edge_y, + z=edge_z, + mode="lines", + line=dict(color="black", width=1), + ) + ) + + fig.update_layout( + scene=dict(aspectmode="data"), + margin=dict(l=0, r=0, t=0, b=0), + ) + + # Force redraw by changing key when toggle changes + dynamic_key = f"{key}_{show_wireframe}" + + st.plotly_chart(fig, use_container_width=True, key=dynamic_key) - st.session_state["_cli_geometry_bootstrapped"] = True - cli_geometry = os.environ.get("CEASIOMPY_GEOMETRY", "").strip() - if not cli_geometry: - return - geometry_path = Path(cli_geometry).expanduser() - if not geometry_path.is_absolute(): - geometry_path = (Path.cwd() / geometry_path).resolve() +def _show_split_components_mesh( + split_dir: str | Path, + key: str, + show_wireframe: bool | None = None +) -> None: + """Display split STL components with one color per component.""" - if not geometry_path.exists(): - st.warning(f"CLI geometry path does not exist: {geometry_path}") + split_dir = Path(split_dir) + component_files = sorted(split_dir.glob("*.stl")) + if not component_files: + st.warning(f"No split component STL found in: {split_dir}") return - st.session_state["pending_default_cpacs"] = str(geometry_path) - st.session_state["uploaded_default_cpacs"] = True - st.session_state["_cli_geometry_autonext"] = True + @st.cache_data(show_spinner=False) + def load_stl(path: str): + return mesh.Mesh.from_file(path) + + n_components = len(component_files) + + def generate_distinct_colors(n: int) -> list[str]: + if n <= 0: + return [] + colors = [] + for idx in range(n): + h = idx / float(n) + r, g, b = colorsys.hsv_to_rgb(h, 0.75, 0.95) + colors.append(f"#{int(r * 255):02x}{int(g * 255):02x}{int(b * 255):02x}") + return colors + + colors = generate_distinct_colors(n_components) + + fig = go.Figure() + for idx, comp_file in enumerate(component_files): + stl_mesh = load_stl(str(comp_file)) + vertices = stl_mesh.vectors.reshape(-1, 3) + verts_unique, index = np.unique(vertices, axis=0, return_inverse=True) + + x, y, z = verts_unique.T + i = index[0::3] + j = index[1::3] + k = index[2::3] + + fig.add_trace( + go.Mesh3d( + x=x, + y=y, + z=z, + i=i, + j=j, + k=k, + color=colors[idx], + opacity=1.0, + flatshading=True, + name=comp_file.stem, + showlegend=True, + ) + ) + + if show_wireframe: + edge_x = [] + edge_y = [] + edge_z = [] + for tri in stl_mesh.vectors: + for e in [(0, 1), (1, 2), (2, 0)]: + edge_x += [tri[e[0]][0], tri[e[1]][0], None] + edge_y += [tri[e[0]][1], tri[e[1]][1], None] + edge_z += [tri[e[0]][2], tri[e[1]][2], None] + + fig.add_trace( + go.Scatter3d( + x=edge_x, + y=edge_y, + z=edge_z, + mode="lines", + line=dict(color="black", width=1), + showlegend=False, + ) + ) + + fig.update_layout( + scene=dict(aspectmode="data"), + margin=dict(l=0, r=0, t=0, b=0), + ) + dynamic_key = f"{key}_{show_wireframe}_{len(component_files)}" + st.plotly_chart(fig, use_container_width=True, key=dynamic_key) def _clean_toolspecific(cpacs: CPACS) -> CPACS: @@ -160,11 +300,13 @@ def _vector_to_str(values: list[float]) -> str: st.session_state["uploaded_default_cpacs"] = False wkdir = st.session_state.workflow.working_dir wkdir.mkdir(parents=True, exist_ok=True) - - safe_name = airfoil_name - if safe_name is None: + if airfoil_name is None: + airfoil_name = "custom" + safe_name = "".join( + char if (char.isalnum() or char in ("-", "_")) else "_" for char in airfoil_name.strip() + ) + if not safe_name: safe_name = "custom" - new_cpacs_path = Path(wkdir, f"airfoil_{safe_name}.xml") cpacs.save_cpacs(new_cpacs_path, overwrite=True) @@ -208,20 +350,6 @@ def _read_airfoil_xy(airfoil_path: Path) -> tuple[list[float], list[float]]: return x_vals, y_vals -def _to_float_or_default(value: str | int | float | None, default: float = 0.0) -> float: - """Best-effort conversion to float for Streamlit numeric widgets.""" - - if value is None: - return default - - try: - casted = float(value) - except (TypeError, ValueError): - return default - - return casted if np.isfinite(casted) else default - - def _section_generate_cpacs_airfoil() -> CPACS | None: st.markdown("#### Generate an Airfoil Profile") @@ -344,9 +472,11 @@ def _section_load_cpacs() -> CPACS | None: st.session_state["last_uploaded_digest"] = uploaded_digest st.session_state["last_uploaded_name"] = uploaded_file.name - suffix = uploaded_path.suffix.lower() - - if suffix in {".csv", ".dat", ".txt"}: + if ( + uploaded_path.suffix == ".csv" + or uploaded_path.suffix == ".dat" + or uploaded_path.suffix == ".txt" + ): try: airfoil_x, airfoil_y = _read_airfoil_xy(uploaded_path) except Exception as exc: @@ -359,7 +489,7 @@ def _section_load_cpacs() -> CPACS | None: airfoil_name=uploaded_path.stem, ) - elif suffix in {".vsp3", ".vsp"}: + elif uploaded_path.suffix == ".vsp3": should_convert = (not is_same_upload) or ( st.session_state.get("last_converted_vsp3_digest") != uploaded_digest ) @@ -421,7 +551,7 @@ def _section_load_cpacs() -> CPACS | None: st.session_state["cpacs"] = cpacs return cpacs - elif suffix == ".xml": + elif uploaded_path.suffix == ".xml": new_cpacs_path = uploaded_path st.session_state["last_converted_cpacs_path"] = str(uploaded_path) cpacs = CPACS(str(uploaded_path)) @@ -438,8 +568,309 @@ def _section_load_cpacs() -> CPACS | None: return st.session_state.get("cpacs") +def _section_stl_to_cpacs(): + + # 1. Load the stl + st.markdown("#### Load an STL file or multiple file") + st.markdown("You can upload a 3D STL model (binary or ASCII) to convert it to CPACS later.") + + # 1. Load STL files via the Streamlit uploader + uploaded_files = st.file_uploader( + label="Upload an STL geometry", + type=["stl"], + key="stl_file_uploader", + accept_multiple_files=True, + label_visibility="collapsed", + ) + + # Return early if nothing loaded + if not uploaded_files: + st.info("Please upload one or more STL files first.") + return None + + # Avoid rewriting same files, but still render preview on every rerun. + wkdir = st.session_state.workflow.working_dir + wkdir.mkdir(parents=True, exist_ok=True) + previous_uploads = st.session_state.get("last_uploaded_stl_files", {}) + current_uploads: dict[str, dict[str, str]] = {} + uploaded_paths: list[Path] = [] + had_new_upload = False + + for uploaded_file in uploaded_files: + uploaded_bytes = uploaded_file.getbuffer() + uploaded_digest = hashlib.sha256(uploaded_bytes).hexdigest() + uploaded_path = wkdir / uploaded_file.name + + previous_meta = previous_uploads.get(uploaded_file.name, {}) + is_same_upload = ( + previous_meta.get("digest") == uploaded_digest + and Path(previous_meta.get("path", "")).exists() + ) + + if is_same_upload: + uploaded_path = Path(previous_meta["path"]) + else: + with open(uploaded_path, "wb") as f: + f.write(uploaded_bytes) + had_new_upload = True + + current_uploads[uploaded_file.name] = { + "digest": uploaded_digest, + "path": str(uploaded_path), + } + uploaded_paths.append(uploaded_path) + + st.session_state["last_uploaded_stl_files"] = current_uploads + st.session_state["last_uploaded_stl_paths"] = [str(path) for path in uploaded_paths] + + selected_stl_name = st.selectbox( + "Select STL to preview", + options=[path.name for path in uploaded_paths], + key="selected_stl_for_preview", + ) + uploaded_path = next(path for path in uploaded_paths if path.name == selected_stl_name) + + # New uploads invalidate stale split results from older geometries. + if had_new_upload: + st.session_state.pop("last_split_components_dir", None) + st.session_state["show_split_tools"] = False + st.session_state["show_cpacs_conversion_tools"] = False + + # 2. visualization + wireframe_view = st.toggle("Wireframe view", value=False, key="stl_wireframe_view") + + # 3. split tools (hidden by default) + if "show_split_tools" not in st.session_state: + st.session_state["show_split_tools"] = False + previous_show_split_tools = bool(st.session_state["show_split_tools"]) + st.session_state["show_split_tools"] = st.toggle( + "Enable split tools", + value=bool(st.session_state["show_split_tools"]), + key="split_tools_toggle", + ) + split_just_enabled = st.session_state["show_split_tools"] and not previous_show_split_tools + + split_dir = st.session_state.get("last_split_components_dir") + if split_dir and Path(split_dir).exists(): + _show_split_components_mesh(split_dir, + key="stl_viewer_split", + show_wireframe=wireframe_view + ) + else: + _show_stl_mesh(str(uploaded_path), key="stl_viewer", show_wireframe=wireframe_view) + + if st.session_state["show_split_tools"]: + + if split_just_enabled: + try: + split_root = wkdir / "STL2CPACS" + if split_root.exists(): + shutil.rmtree(split_root) + + split_dir = split_main( + stl_path=str(uploaded_path), + namefile=f"{uploaded_path.stem}_", + out_dir=wkdir, + ) + st.session_state["last_split_components_dir"] = str(split_dir) + st.success("Component split completed.") + st.rerun() + except Exception as exc: + st.error(f"STL split failed: {exc}") + + # 4. Convert to CPACS using automatic component detection + st.markdown("#### Convert to CPACS") + + candidate_paths: list[Path] + if split_dir and Path(split_dir).exists(): + candidate_paths = sorted(Path(split_dir).glob("*.stl")) + else: + candidate_paths = uploaded_paths + + if not candidate_paths: + st.warning("No STL sources found in the working directory.") + return None + + current_candidates = [str(path) for path in candidate_paths] + selected_stl_files = current_candidates + try: + detected_cpacs_types = cpacs_component_detection(selected_stl_files) + except Exception as exc: + st.error(f"Failed to detect CPACS component types: {exc}") + return None + + if len(selected_stl_files) != len(detected_cpacs_types): + st.error("Detected component data is inconsistent.") + return None + + st.markdown("##### Component Settings") + settings: list[dict] = [] + for idx, (stl_path, auto_type) in enumerate(zip(selected_stl_files, detected_cpacs_types)): + normalized_type = "fuselage" if str(auto_type).lower() == "fuselage" else "wing" + st.markdown( + f"**Component {idx + 1}** - `{Path(stl_path).name}` " + f"(auto-detected: `{normalized_type}`)" + ) + + use_advanced = st.toggle( + "Advanced", + value=False, + key=f"stl_component_advanced_{idx}", + ) + + setting_dict: dict[str, float | int] = {} + if use_advanced: + if normalized_type == "wing": + adv_col1, adv_col2, adv_col3 = st.columns(3) + with adv_col1: + setting_dict["EXTREME_TOL_perc_start"] = float( + st.number_input( + "EXTREME_TOL_perc_start", + value=0.005, + min_value=0.0, + max_value=1.0, + format="%.4f", + help="Fraction of span removed at the root side start before slicing.", + key=f"stl_component_adv_extreme_tol_start_{idx}", + ) + ) + setting_dict["N_Y_SLICES"] = int( + st.number_input( + "Slices", + value=500, + min_value=2, + step=1, + help="Number of spanwise slicing planes for wing extraction.", + key=f"stl_component_adv_n_y_slices_{idx}", + ) + ) + with adv_col2: + setting_dict["EXTREME_TOL_perc_end"] = float( + st.number_input( + "EXTREME_TOL_perc_end", + value=0.005, + min_value=0.0, + max_value=1.0, + format="%.4f", + help="Fraction of span removed at the tip side end before slicing.", + key=f"stl_component_adv_extreme_tol_end_{idx}", + ) + ) + setting_dict["N_SLICE_ADDING"] = int( + st.number_input( + "N_SLICE_ADDING", + value=0, + min_value=0, + step=1, + help="Extra interpolated sections inserted where local angles change.", + key=f"stl_component_adv_n_slice_adding_{idx}", + ) + ) + with adv_col3: + setting_dict["TE_CUT"] = float( + st.number_input( + "TE_CUT", + value=0.0, + min_value=0.0, + max_value=0.5, + format="%.4f", + help=("Trailing-edge trim ratio applied on each extracted airfoil.", + " Increase if there are oscillations in the TE region."), + key=f"stl_component_adv_te_cut_{idx}", + ) + ) + setting_dict["N_BIN"] = int( + st.number_input( + "N_BIN", + value=10, + min_value=3, + step=1, + help=("Bins used to split the cloud into upper/lower surfaces.", + " Increase if you have a refined profile with sharp features,", + " or if you have oscillations on the CPACS airfoil."), + key=f"stl_component_adv_n_bin_{idx}", + ) + ) + else: + adv_col1, adv_col2 = st.columns(2) + with adv_col1: + setting_dict["EXTREME_TOL_perc_start"] = float( + st.number_input( + "EXTREME_TOL_perc_start", + value=0.005, + min_value=0.0, + max_value=1.0, + format="%.4f", + help="Fraction of fuselage length removed at the nose before slicing.", + key=f"stl_component_adv_extreme_tol_start_{idx}", + ) + ) + setting_dict["N_X_SLICES"] = int( + st.number_input( + "N_X_SLICES", + value=100, + min_value=2, + step=1, + help="Numbers of slicing planes for fuselage extraction.", + key=f"stl_component_adv_n_x_slices_{idx}", + ) + ) + with adv_col2: + setting_dict["EXTREME_TOL_perc_end"] = float( + st.number_input( + "EXTREME_TOL_perc_end", + value=0.005, + min_value=0.0, + max_value=1.0, + format="%.4f", + help="Fraction of fuselage length removed at the tail before slicing.", + key=f"stl_component_adv_extreme_tol_end_{idx}", + ) + ) + setting_dict["N_SLICE_ADDING"] = int( + st.number_input( + "N_SLICE_ADDING", + value=0, + min_value=0, + step=1, + help="Extra interpolated sections inserted at angle transitions.", + key=f"stl_component_adv_n_slice_adding_{idx}", + ) + ) + settings.append(setting_dict) + + if st.button("Convert Geometry to CPACS", key="convert_stl_to_cpacs_button", width="stretch"): + try: + cpacs_path, conversion_report = stl2cpacs_main( + stl_file=selected_stl_files, + setting=settings, + out_dir=wkdir, + ) + st.session_state["last_stl2cpacs_conversion_report"] = conversion_report + + st.session_state["last_converted_cpacs_path"] = str(cpacs_path) + close_cpacs_handles(st.session_state.get("cpacs")) + st.session_state["cpacs"] = CPACS(str(cpacs_path)) + st.success(f"CPACS generated: {cpacs_path}") + + failed_components = conversion_report.get("failed", []) + if failed_components: + st.warning( + f"{len(failed_components)} component could not be converted and were skipped." + ) + for comp in failed_components: + st.caption( + f"- `{Path(comp.get('path', '')).name}` " + f"({comp.get('type', 'unknown')}): {comp.get('error', 'Unknown error')}" + ) + except Exception as exc: + st.error(f"STL to CPACS conversion failed: {exc}") + + return st.session_state.get("cpacs") + # Functions + def section_select_cpacs() -> None: wkdir = get_wkdir() wkdir.mkdir(parents=True, exist_ok=True) @@ -448,6 +879,7 @@ def section_select_cpacs() -> None: tabs = [ "Load Geometry", "Generate Airfoil", + "STL to CPACS", ] show_openvsp = not parse_bool(os.environ.get("CEASIOMPY_CLOUD", "False")) @@ -460,8 +892,10 @@ def section_select_cpacs() -> None: _section_load_cpacs() with selected_tab[1]: _section_generate_cpacs_airfoil() + with selected_tab[2]: + _section_stl_to_cpacs() if show_openvsp: - with selected_tab[2]: + with selected_tab[3]: render_openvsp_panel() # ALWAYS use the session CPACS @@ -472,20 +906,24 @@ def section_select_cpacs() -> None: tixi = cpacs.tixi # 2D mode or else; geometry mode is expected to exist. + # STL2CPACS exports may not include CEASIOMpy geometry mode; default to 3D. + if not tixi.checkElement(GEOMETRY_MODE_XPATH): + add_value( + tixi=tixi, + xpath=GEOMETRY_MODE_XPATH, + value="3D", + ) geometry_mode = str(get_value(tixi, xpath=GEOMETRY_MODE_XPATH)) dim_mode = (geometry_mode == "2D") st.markdown("---") - ref_area: float | None = None - ref_length: float = 0.0 + ref_area = None try: if not dim_mode: - area, length = compute_aircraft_ref_values(cpacs) - ref_area = _to_float_or_default(area, 0.0) - ref_length = _to_float_or_default(length, 0.0) + ref_area, ref_length = compute_aircraft_ref_values(cpacs) if dim_mode: - ref_length = _to_float_or_default(compute_airfoil_ref_length(cpacs), 0.0) + ref_length = compute_airfoil_ref_length(cpacs) except Exception as e: ref_area, ref_length = 0.0, 0.0 @@ -502,6 +940,7 @@ def section_select_cpacs() -> None: section_3d_view( cpacs=cpacs, force_regenerate=True, + plot_key="geometry_page_cpacs_preview", ) # Once 3D view of CPACS file is done scroll down @@ -509,10 +948,10 @@ def section_select_cpacs() -> None: st.markdown("---") if tixi.checkElement(AREA_XPATH): - ref_area = _to_float_or_default(get_value(tixi, xpath=AREA_XPATH), 0.0) + ref_area = get_value(tixi, xpath=AREA_XPATH) if tixi.checkElement(LENGTH_XPATH): - ref_length = _to_float_or_default(get_value(tixi, xpath=LENGTH_XPATH), 0.0) + ref_length = get_value(tixi, xpath=LENGTH_XPATH) spec = 1 if dim_mode else 2 cols = st.columns( @@ -524,12 +963,11 @@ def section_select_cpacs() -> None: value=ref_length, min_value=0.0, ) - new_ref_area: float | None = ref_area if not dim_mode: with cols[1]: new_ref_area = st.number_input( label="Reference Area", - value=ref_area if ref_area is not None else 0.0, + value=ref_area, min_value=0.0, ) @@ -549,10 +987,8 @@ def section_select_cpacs() -> None: # 3D update if ( not dim_mode - and ref_area is not None and np.isfinite(ref_area) and ref_area > 0.0 and np.isfinite(ref_length) and ref_length > 0.0 - and new_ref_area is not None ): add_value( tixi=tixi, @@ -598,8 +1034,6 @@ def section_select_cpacs() -> None: if "workflow" not in st.session_state: st.session_state["workflow"] = Workflow() - _bootstrap_cli_geometry() - st.markdown("---") section_select_cpacs() diff --git a/src/app/streamlitutils.py b/src/app/streamlitutils.py index 6f9eaf6c5..cd966b858 100755 --- a/src/app/streamlitutils.py +++ b/src/app/streamlitutils.py @@ -408,7 +408,12 @@ def section_3d_view( try: pv_mesh = pv.read(str(vtp_file)) - surface = pv_mesh.extract_surface(algorithm="dataset_surface").triangulate().clean() + try: + surface = pv_mesh.extract_surface(algorithm="dataset_surface") + except TypeError: + # Older PyVista versions do not support the "algorithm" keyword. + surface = pv_mesh.extract_surface() + surface = surface.triangulate().clean() surface = surface.compute_normals() except Exception as exc: st.error(f"Failed to read generated preview mesh for 3D view: {exc}") @@ -467,8 +472,7 @@ def _axis_values(vmin: float, vmax: float, count: int) -> list[float]: scene=dict( aspectmode="data", xaxis=dict( - title="X", - titlefont=dict(color="black"), + title=dict(text="X", font=dict(color="black")), tickfont=dict(color="black"), range=[x_min, x_max], tickmode="array", @@ -479,8 +483,7 @@ def _axis_values(vmin: float, vmax: float, count: int) -> list[float]: backgroundcolor="white", ), yaxis=dict( - title="Y" if show_yaxis else "", - titlefont=dict(color="black"), + title=dict(text="Y" if show_yaxis else "", font=dict(color="black")), tickfont=dict(color="black"), range=[y_min, y_max] if show_yaxis else [0.0, 0.0], tickmode="array", @@ -493,8 +496,7 @@ def _axis_values(vmin: float, vmax: float, count: int) -> list[float]: backgroundcolor="white", ), zaxis=dict( - title="Z", - titlefont=dict(color="black"), + title=dict(text="Z", font=dict(color="black")), tickfont=dict(color="black"), range=[z_min, z_max], tickmode="array", diff --git a/src/ceasiompy/stl2cpacs/func/split_stl_into_components.py b/src/ceasiompy/stl2cpacs/func/split_stl_into_components.py new file mode 100644 index 000000000..391d2a85b --- /dev/null +++ b/src/ceasiompy/stl2cpacs/func/split_stl_into_components.py @@ -0,0 +1,692 @@ +"""Utilities to split a full-aircraft STL into geometric components. + +This module provides a pragmatic STL splitter for CEASIOMpy workflows: +- load an STL (ASCII or binary) +- split disconnected components using triangle connectivity +- export each detected connected part as a generic `component_i.stl +""" + +from __future__ import annotations +from dataclasses import dataclass +from pathlib import Path +import struct +from typing import Dict, List, Tuple +import numpy as np + +# Quantization tolerance used to merge numerically close vertices. +VERTEX_MERGE_TOL = 1e-6 + +# Connectivity tolerance used while grouping triangles into disconnected parts. +DEFAULT_VERTEX_TOL = VERTEX_MERGE_TOL +DEFAULT_FEATURE_ANGLE_DEG = 55.0 +SIGNIFICANT_COMPONENT_MIN_TRIS = 100 + +# ====================================================================================== +# HOW THE SPLIT WORKS +# ====================================================================================== +# 1) Read STL triangles as an array with shape (N, 3, 3) +# N triangles, each triangle has 3 vertices, each vertex has (x, y, z). +# +# 2) Build triangle-to-triangle adjacency: +# - If two triangles share at least one vertex, they are connected. +# - Vertices are quantized with `vertex_tol` to avoid floating-point noise +# causing false disconnections. +# +# 3) Run graph connected-components: +# - Each connected triangle group becomes one output component. +# +# 4) Save each group as `component_i.stl`. +# +# Result: +# - A fully connected STL => one component file. +# - A multi-part STL => one file per disconnected part. + + +@dataclass +class ComponentInfo: + """Container for one split component.""" + + name: str + component_type: str + triangles: np.ndarray + n_triangles: int + bbox_min: np.ndarray + bbox_max: np.ndarray + + +def read_ascii_stl(path: str | Path) -> np.ndarray: + """Read ASCII STL and return triangles as (N, 3, 3).""" + + tri = [] + with open(path, "r", encoding="utf-8", errors="ignore") as f: + for line in f: + if line.strip().startswith("vertex"): + _, x, y, z = line.split()[:4] + tri.append([float(x), float(y), float(z)]) + + return np.asarray(tri, dtype=float).reshape(-1, 3, 3) + + +def read_binary_stl(path: str | Path) -> np.ndarray: + """Read binary STL and return triangles as (N, 3, 3).""" + + with open(path, "rb") as f: + f.read(80) + ntri = struct.unpack(" np.ndarray: + """Auto-detect STL format and read triangles as (N, 3, 3).""" + + path = Path(path) + with open(path, "rb") as f: + start = f.read(80) + + if start[:5].lower() == b"solid": + try: + return read_ascii_stl(path) + except Exception: + return read_binary_stl(path) + + return read_binary_stl(path) + + +def write_binary_stl(path: str | Path, + triangles: np.ndarray, + solid_name: str = "component" + ) -> None: + """Write triangles to a binary STL file.""" + + tris = np.asarray(triangles, dtype=np.float32).reshape(-1, 3, 3) + with open(path, "wb") as f: + header = solid_name.encode("ascii", errors="ignore")[:80] + header = header + b" " * (80 - len(header)) + f.write(header) + f.write(struct.pack(" List[List[int]]: + """Build triangle adjacency list by shared quantized vertices.""" + + tris = np.asarray(triangles, dtype=float).reshape(-1, 3, 3) + n_tri = tris.shape[0] + if n_tri == 0: + return [] + + verts = tris.reshape(-1, 3) + qverts = np.round(verts / tol).astype(np.int64) + tri_ids = np.repeat(np.arange(n_tri, dtype=np.int64), 3) + + vert_to_tris: Dict[Tuple[int, int, int], List[int]] = {} + for tri_id, key in zip(tri_ids, map(tuple, qverts)): + vert_to_tris.setdefault(key, []).append(int(tri_id)) + + adjacency: List[set] = [set() for _ in range(n_tri)] + for tri_list in vert_to_tris.values(): + if len(tri_list) < 2: + continue + unique_tris = list(set(tri_list)) + for i in unique_tris: + adjacency[i].update(unique_tris) + adjacency[i].discard(i) + + return [list(nei) for nei in adjacency] + + +def _triangle_adjacency_from_shared_edges(triangles: np.ndarray, + tol: float = VERTEX_MERGE_TOL + ) -> List[List[int]]: + """Build triangle adjacency from shared *manifold* edges. + + Used as a fallback when vertex-connectivity yields a single component. + It prevents over-merging components that only touch at points/edges. + """ + + tris = np.asarray(triangles, dtype=float).reshape(-1, 3, 3) + n_tri = tris.shape[0] + if n_tri == 0: + return [] + + verts = tris.reshape(-1, 3) + qverts = np.round(verts / tol).astype(np.int64) + + tri_ids = np.repeat(np.arange(n_tri, dtype=np.int64), 3) + local_vid = np.tile(np.arange(3, dtype=np.int64), n_tri) + + idx_mat = np.empty((n_tri, 3), dtype=np.int64) + idx_mat[tri_ids, local_vid] = np.arange(n_tri * 3, dtype=np.int64) + + e0 = idx_mat[:, [0, 1]] + e1 = idx_mat[:, [1, 2]] + e2 = idx_mat[:, [2, 0]] + edge_idx = np.vstack([e0, e1, e2]) + edge_tri = np.repeat(np.arange(n_tri, dtype=np.int64), 3) + + va = qverts[edge_idx[:, 0]] + vb = qverts[edge_idx[:, 1]] + + swap = ( + (va[:, 0] > vb[:, 0]) + | ((va[:, 0] == vb[:, 0]) & (va[:, 1] > vb[:, 1])) + | ((va[:, 0] == vb[:, 0]) & (va[:, 1] == vb[:, 1]) & (va[:, 2] > vb[:, 2])) + ) + first = np.where(swap[:, None], vb, va) + second = np.where(swap[:, None], va, vb) + edges = np.hstack([first, second]) + + order = np.lexsort( + (edges[:, 5], edges[:, 4], edges[:, 3], edges[:, 2], edges[:, 1], edges[:, 0]) + ) + edges_sorted = edges[order] + tris_sorted = edge_tri[order] + + adjacency: List[set] = [set() for _ in range(n_tri)] + i = 0 + m = edges_sorted.shape[0] + while i < m: + j = i + 1 + while j < m and np.array_equal(edges_sorted[j], edges_sorted[i]): + j += 1 + + # Only connect through manifold edges (exactly 2 incident triangles). + if j - i == 2: + t0 = int(tris_sorted[i]) + t1 = int(tris_sorted[i + 1]) + if t0 != t1: + adjacency[t0].add(t1) + adjacency[t1].add(t0) + + i = j + + return [list(nei) for nei in adjacency] + + +def _triangle_normals(triangles: np.ndarray) -> np.ndarray: + """Compute unit triangle normals for an array of shape (N, 3, 3).""" + + tris = np.asarray(triangles, dtype=float).reshape(-1, 3, 3) + v1 = tris[:, 1] - tris[:, 0] + v2 = tris[:, 2] - tris[:, 0] + n = np.cross(v1, v2) + norm = np.linalg.norm(n, axis=1) + valid = norm > 1e-16 + out = np.zeros_like(n) + out[valid] = n[valid] / norm[valid, None] + return out + + +def _triangle_adjacency_from_smooth_shared_edges( + triangles: np.ndarray, + tol: float = VERTEX_MERGE_TOL, + max_dihedral_deg: float = DEFAULT_FEATURE_ANGLE_DEG +) -> List[List[int]]: + """Build adjacency using only manifold edges with smooth dihedral angle.""" + + tris = np.asarray(triangles, dtype=float).reshape(-1, 3, 3) + n_tri = tris.shape[0] + if n_tri == 0: + return [] + + verts = tris.reshape(-1, 3) + qverts = np.round(verts / tol).astype(np.int64) + normals = _triangle_normals(tris) + + idx_mat = np.arange(n_tri * 3, dtype=np.int64).reshape(n_tri, 3) + edge_idx = np.vstack([idx_mat[:, [0, 1]], idx_mat[:, [1, 2]], idx_mat[:, [2, 0]]]) + edge_tri = np.repeat(np.arange(n_tri, dtype=np.int64), 3) + + va = qverts[edge_idx[:, 0]] + vb = qverts[edge_idx[:, 1]] + swap = ( + (va[:, 0] > vb[:, 0]) + | ((va[:, 0] == vb[:, 0]) & (va[:, 1] > vb[:, 1])) + | ((va[:, 0] == vb[:, 0]) & (va[:, 1] == vb[:, 1]) & (va[:, 2] > vb[:, 2])) + ) + first = np.where(swap[:, None], vb, va) + second = np.where(swap[:, None], va, vb) + edges = np.hstack([first, second]) + + order = np.lexsort( + (edges[:, 5], edges[:, 4], edges[:, 3], edges[:, 2], edges[:, 1], edges[:, 0]) + ) + edges_sorted = edges[order] + tris_sorted = edge_tri[order] + + cos_thr = float(np.cos(np.deg2rad(max_dihedral_deg))) + adjacency: List[set] = [set() for _ in range(n_tri)] + i = 0 + m = edges_sorted.shape[0] + while i < m: + j = i + 1 + while j < m and np.array_equal(edges_sorted[j], edges_sorted[i]): + j += 1 + + if j - i == 2: + t0 = int(tris_sorted[i]) + t1 = int(tris_sorted[i + 1]) + if t0 != t1: + # abs(dot) makes this robust to inconsistent normal orientation. + c = abs(float(np.dot(normals[t0], normals[t1]))) + if c >= cos_thr: + adjacency[t0].add(t1) + adjacency[t1].add(t0) + + i = j + + return [list(nei) for nei in adjacency] + + +def _extract_components_from_adjacency(adjacency: List[List[int]]) -> List[np.ndarray]: + """Extract connected components from a precomputed adjacency list.""" + + n_tri = len(adjacency) + if n_tri == 0: + return [] + + visited = np.zeros(n_tri, dtype=bool) + components: List[np.ndarray] = [] + for seed in range(n_tri): + if visited[seed]: + continue + stack = [seed] + visited[seed] = True + comp_idx = [] + while stack: + t = stack.pop() + comp_idx.append(t) + for nb in adjacency[t]: + if not visited[nb]: + visited[nb] = True + stack.append(nb) + components.append(np.asarray(comp_idx, dtype=int)) + return components + + +def _connected_triangle_components(triangles: np.ndarray, + tol: float = VERTEX_MERGE_TOL + ) -> List[np.ndarray]: + """Split triangles into components with robust auto-connectivity. + + We treat triangles as nodes of a graph: + - node = one triangle + - edge = triangle connectivity relationship + """ + + tris = np.asarray(triangles, dtype=float).reshape(-1, 3, 3) + if tris.shape[0] == 0: + return [] + + # Fast path: vertex connectivity preserves legacy behavior. + vertex_adj = _triangle_adjacency_from_shared_vertices(tris, tol=tol) + vertex_components = _extract_components_from_adjacency(vertex_adj) + if len(vertex_components) > 1: + return vertex_components + + # Fallback for fully connected meshes with coincident interfaces. + edge_adj = _triangle_adjacency_from_shared_edges(tris, tol=tol) + return _extract_components_from_adjacency(edge_adj) + + +def _count_significant_components(components: List[np.ndarray], + min_triangles: int = SIGNIFICANT_COMPONENT_MIN_TRIS + ) -> int: + """Count components with enough triangles to be meaningful geometric parts.""" + + return int(sum(1 for c in components if len(c) >= min_triangles)) + + +def _feature_split_largest_component( + triangles: np.ndarray, comp_indices: List[np.ndarray], tol: float = VERTEX_MERGE_TOL +) -> List[np.ndarray]: + """Optionally split the largest component by sharp feature edges.""" + + if not comp_indices: + return comp_indices + + if len(comp_indices) > 3: + return comp_indices + + if _count_significant_components(comp_indices) >= 2: + return comp_indices + + largest_id = int(np.argmax([len(c) for c in comp_indices])) + largest = comp_indices[largest_id] + local_tris = triangles[largest] + smooth_adj = _triangle_adjacency_from_smooth_shared_edges( + local_tris, tol=tol, max_dihedral_deg=DEFAULT_FEATURE_ANGLE_DEG + ) + smooth_sub = _extract_components_from_adjacency(smooth_adj) + + if _count_significant_components(smooth_sub) < 2: + return comp_indices + + remapped = [largest[sub] for sub in smooth_sub] + for i, comp in enumerate(comp_indices): + if i != largest_id: + remapped.append(comp) + return remapped + + +def _histogram_valley_threshold(values: np.ndarray, n_bins: int = 80) -> float | None: + """Find a robust valley threshold in a 1D distribution.""" + + v = np.asarray(values, dtype=float) + if v.size < 100: + return None + + hist, edges = np.histogram(v, bins=n_bins) + if np.all(hist == 0): + return None + + smooth = np.convolve(hist, np.array([1, 2, 3, 2, 1], dtype=float), mode="same") + cdf = np.cumsum(hist) / float(np.sum(hist)) + lo = int(0.10 * n_bins) + hi = int(0.90 * n_bins) + if hi <= lo + 2: + return None + + best_i = None + best_val = None + for i in range(lo + 1, hi - 1): + left_frac = float(cdf[i]) + right_frac = 1.0 - left_frac + if left_frac < 0.15 or right_frac < 0.15: + continue + if smooth[i] <= smooth[i - 1] and smooth[i] <= smooth[i + 1]: + val = smooth[i] + if best_val is None or val < best_val: + best_val = val + best_i = i + + if best_i is None: + return None + + return float(0.5 * (edges[best_i] + edges[best_i + 1])) + + +def _induced_components_from_mask(adjacency: List[List[int]], + mask: np.ndarray + ) -> List[np.ndarray]: + """Connected components of the subgraph induced by `mask`.""" + + mask = np.asarray(mask, dtype=bool) + local_ids = np.where(mask)[0] + if local_ids.size == 0: + return [] + + g2l = -np.ones(mask.shape[0], dtype=int) + g2l[local_ids] = np.arange(local_ids.size, dtype=int) + + local_adj: List[List[int]] = [[] for _ in range(local_ids.size)] + for g in local_ids: + li = int(g2l[g]) + for nb in adjacency[int(g)]: + lnb = int(g2l[int(nb)]) + if lnb >= 0: + local_adj[li].append(lnb) + + local_components = _extract_components_from_adjacency(local_adj) + return [local_ids[c] for c in local_components] + + +def _span_split_largest_component( + triangles: np.ndarray, + comp_indices: List[np.ndarray], + tol: float = VERTEX_MERGE_TOL +) -> List[np.ndarray]: + """Split one dominant shell into inboard/outboard parts using |y| valley.""" + + if not comp_indices: + return comp_indices + if len(comp_indices) > 3: + return comp_indices + if _count_significant_components(comp_indices) >= 2: + return comp_indices + + largest_id = int(np.argmax([len(c) for c in comp_indices])) + largest = comp_indices[largest_id] + local_tris = triangles[largest] + centroids = np.mean(local_tris, axis=1) + abs_y = np.abs(centroids[:, 1]) + y_thr = _histogram_valley_threshold(abs_y) + if y_thr is None: + return comp_indices + + mask_out = abs_y > y_thr + n_out = int(np.count_nonzero(mask_out)) + n_in = int(mask_out.size - n_out) + if min(n_in, n_out) < SIGNIFICANT_COMPONENT_MIN_TRIS: + return comp_indices + + local_adj = _triangle_adjacency_from_shared_vertices(local_tris, tol=tol) + out_sub = _induced_components_from_mask(local_adj, mask_out) + in_sub = _induced_components_from_mask(local_adj, ~mask_out) + split_sub = out_sub + in_sub + if _count_significant_components(split_sub) < 2: + return comp_indices + + remapped = [largest[sub] for sub in split_sub] + for i, comp in enumerate(comp_indices): + if i != largest_id: + remapped.append(comp) + return remapped + + +def _build_generic_components(disconnected_components: List[np.ndarray], + triangles: np.ndarray + ) -> List[ComponentInfo]: + """Build generic component objects from connected triangle groups. + + This function only packages metadata and names: + component_1, component_2, ... + """ + + components: List[ComponentInfo] = [] + + for i, idxs in enumerate(disconnected_components, start=1): + # Extract triangle subset for one connected group. + comp_tris = triangles[idxs] + # Bounding box is useful for manual inspection in logs. + verts = comp_tris.reshape(-1, 3) + components.append( + ComponentInfo( + name=f"component_{i}", + component_type="component", + triangles=comp_tris, + n_triangles=comp_tris.shape[0], + bbox_min=np.min(verts, axis=0), + bbox_max=np.max(verts, axis=0), + ) + ) + + # Largest first for easier manual review. + components.sort(key=lambda c: c.n_triangles, reverse=True) + return components + + +def _symmetry_plane_to_axis_index(symmetry_plane: str) -> int: + """Map symmetry plane string to axis index normal to that plane.""" + + plane_key = symmetry_plane.strip().lower().replace(" ", "").replace("_", "-") + plane_to_axis = { + "x-z": 1, + "z-x": 1, + "xz": 1, + "zx": 1, + "x-y": 2, + "y-x": 2, + "xy": 2, + "yx": 2, + "y-z": 0, + "z-y": 0, + "yz": 0, + "zy": 0, + } + if plane_key not in plane_to_axis: + raise ValueError( + f"Unsupported symmetry_plane '{symmetry_plane}'. " + "Use one of: x-z, x-y, y-z." + ) + return int(plane_to_axis[plane_key]) + + +def find_symmetry_split_midpoint(stl_path: str | Path, symmetry_plane: str) -> float: + """Return midpoint coordinate along the normal axis of a symmetry plane. + + Examples + -------- + - symmetry_plane="x-z" -> split axis is y + - symmetry_plane="x-y" -> split axis is z + - symmetry_plane="y-z" -> split axis is x + """ + + triangles = load_stl_auto(stl_path) + if triangles.size == 0: + raise ValueError(f"No triangles found in STL: {stl_path}") + + axis_idx = _symmetry_plane_to_axis_index(symmetry_plane) + coords = triangles.reshape(-1, 3)[:, axis_idx] + axis_min = float(np.min(coords)) + axis_max = float(np.max(coords)) + return 0.5 * (axis_min + axis_max) + + +def split_stl_by_symmetry_plane( + stl_path: str | Path, + symmetry_plane: str, + output_dir: str | Path, + name: str = "component", + eps: float = 1e-10, +) -> Tuple[Path, Path]: + """Split STL into left/right parts and write two STL files. + + Triangles are assigned by centroid coordinate on the split axis: + - left: centroid <= midpoint + eps + - right: centroid >= midpoint - eps + """ + + triangles = load_stl_auto(stl_path) + if triangles.size == 0: + raise ValueError(f"No triangles found in STL: {stl_path}") + + axis_idx = _symmetry_plane_to_axis_index(symmetry_plane) + midpoint = find_symmetry_split_midpoint(stl_path, symmetry_plane) + centroids = np.mean(triangles, axis=1) + axis_values = centroids[:, axis_idx] + + left_mask = axis_values <= (midpoint + eps) + right_mask = axis_values >= (midpoint - eps) + left_tris = triangles[left_mask] + right_tris = triangles[right_mask] + + if left_tris.shape[0] == 0 or right_tris.shape[0] == 0: + raise ValueError( + "Symmetry split failed: one side is empty. " + "Check symmetry_plane or STL orientation." + ) + + out_dir = Path(output_dir) + out_dir.mkdir(parents=True, exist_ok=True) + left_path = out_dir / f"{name}_left.stl" + right_path = out_dir / f"{name}_right.stl" + write_binary_stl(left_path, left_tris, solid_name=f"{name}_left") + write_binary_stl(right_path, right_tris, solid_name=f"{name}_right") + return left_path, right_path + + +def split_aircraft_stl( + stl_path: str | Path, + output_dir: str | Path, + name: str, + vertex_tol: float = VERTEX_MERGE_TOL, +) -> List[ComponentInfo]: + """Split full-aircraft STL into generic connected components. + + Parameters + ---------- + stl_path: + Input STL path. + output_dir: + If provided, one STL per component is written into this folder. + vertex_tol: + Vertex merging tolerance used to detect triangle connectivity. + + Returns + ------- + list[ComponentInfo] + Split components (generic, no semantic classification). + """ + + # 1) Load STL triangles + triangles = load_stl_auto(stl_path) + if triangles.size == 0: + return [] + + # 2) Split by disconnected triangle connectivity. + comp_indices = _connected_triangle_components(triangles, tol=vertex_tol) + comp_indices = _feature_split_largest_component(triangles, comp_indices, tol=vertex_tol) + comp_indices = _span_split_largest_component(triangles, comp_indices, tol=vertex_tol) + + # 3) Convert connected groups to generic component metadata. + components = _build_generic_components(comp_indices, triangles) + + # 4) output files (one STL per component) + if output_dir is not None: + out_dir = Path(output_dir) + out_dir.mkdir(parents=True, exist_ok=True) + split_dir = out_dir / "STL2CPACS" + split_dir.mkdir(parents=True, exist_ok=True) + for comp in components: + out_path = split_dir / f"{name}{comp.name}.stl" + write_binary_stl(out_path, comp.triangles, solid_name=comp.name) + + return components + + +def split_main(stl_path: str | Path, namefile: str, out_dir: str | Path) -> Path: + + vertex_tol = DEFAULT_VERTEX_TOL + + if not Path(stl_path).exists(): + raise FileNotFoundError( + f"STL file not found: {stl_path}\n" + "Set INPUT_STL_PATH in splitstlgeom.py." + ) + + print("Running STL splitter...") + print(f"Input STL : {stl_path}") + split_dir = Path(out_dir) / "STL2CPACS" + print(f"Output dir: {split_dir}") + + comps = split_aircraft_stl(stl_path, output_dir=out_dir, vertex_tol=vertex_tol, name=namefile) + if not comps: + print("No triangles found.") + else: + print(f"\nWrote {len(comps)} split STL file(s) into: {split_dir}") + + return split_dir diff --git a/src/ceasiompy/stl2cpacs/func/stl2fuselage.py b/src/ceasiompy/stl2cpacs/func/stl2fuselage.py new file mode 100644 index 000000000..dbc013178 --- /dev/null +++ b/src/ceasiompy/stl2cpacs/func/stl2fuselage.py @@ -0,0 +1,1040 @@ + +# ================================================================================================= +# IMPORTS +# ================================================================================================= + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.cm as cm +from scipy.interpolate import PchipInterpolator +from ceasiompy.stl2cpacs.stl2cpacs import export_mesh, parse_cart3d_tri +from pathlib import Path + +# --------------------------- +# CONFIG +# --------------------------- + +INTERSECT_TOL = 1e-6 +SLAB_TOLS = [1e-5, 5e-5, 1e-4, 5e-4, 1e-3] +DUPLICATE_YZ_TOL = 1e-8 +# ================================================================================================= +# FUNCTIONS +# ================================================================================================= + + +def _save_debug_stl_and_slices_plot( + pts: np.ndarray, + tris: np.ndarray, + per_slice_clouds: list[np.ndarray], + ref_points: list[np.ndarray | None], + output_directory: str | Path, + name: str, + view_elev: float = 0, + view_azim: float = 0.0, + interactive: bool = False, +) -> None: + """ + Save a debug 3D figure with STL surface, first slice clouds and reference picks. + """ + + out_dir = Path(output_directory) + out_dir.mkdir(parents=True, exist_ok=True) + debug_path = out_dir / f"{name}_debug_stl_first_slices.png" + + fig = plt.figure(figsize=(11, 8)) + ax = fig.add_subplot(111, projection="3d") + + if tris.shape[0] > 0: + + tris_plot = tris + ax.plot_trisurf( + pts[:, 0], + pts[:, 1], + pts[:, 2], + triangles=tris_plot, + color="#4C78A8", + alpha=0.24, + linewidth=0.0, + shade=False, + ) + ax.plot([], [], [], color="#4C78A8", lw=4, alpha=0.85, label="STL mesh") + + cmap = cm.get_cmap("viridis", max(1, len(per_slice_clouds))) + first_slice_label = True + for i, cloud in enumerate(per_slice_clouds): + if cloud.shape[0] == 0: + continue + ax.scatter( + cloud[:, 0], + cloud[:, 1], + cloud[:, 2], + s=5, + color=cmap(i), + alpha=0.8, + label="Slice intersections" if first_slice_label else None, + ) + first_slice_label = False + + ref_valid = [pt for pt in ref_points if pt is not None] + if ref_valid: + ref_arr = np.vstack(ref_valid) + ax.scatter( + ref_arr[:, 0], + ref_arr[:, 1], + ref_arr[:, 2], + s=24, + c="red", + alpha=0.95, + label="Reference points", + ) + + # Camera setup: azim=0 gives a front view aligned with X direction. + ax.view_init(elev=view_elev, azim=view_azim) + ax.set_title("Fuselage slicing") + ax.set_xlabel("X") + ax.set_ylabel("Y") + ax.set_zlabel("Z") + ax.grid(True, alpha=0.3) + ax.legend(loc="upper right") + fig.tight_layout() + fig.savefig(debug_path, dpi=200) + if interactive: + plt.show() + plt.close(fig) + + +def resample_fuselage_cpacs( + yr, zr, + yl, zl, + y_bot, z_bot, + y_top, z_top, + n_points, + use_cosine_spacing=False, +): + """ + Regularize a fuselage cross-section using PCHIP interpolation. + + CPACS order: + bottom → right side → top → left side → bottom + """ + + yr = np.asarray(yr) + zr = np.asarray(zr) + yl = np.asarray(yl) + zl = np.asarray(zl) + + def _prepare_strictly_increasing(z_in, y_in, tol=1e-9): + """ + Sort by z (bottom -> top) and enforce strict increase with a tiny epsilon. + This keeps the original side curvature, even for complex shapes. + """ + + z_in = np.asarray(z_in, dtype=float) + y_in = np.asarray(y_in, dtype=float) + + # Stable sort preserves relative ordering for equal z values. + order = np.argsort(z_in, kind="mergesort") + z_sorted = z_in[order] + y_sorted = y_in[order] + + # Final strict-increase pass (numerical safety). + for i in range(1, z_sorted.size): + if z_sorted[i] <= z_sorted[i - 1]: + z_sorted[i] = z_sorted[i - 1] + tol + + return z_sorted, y_sorted + + # ------------------------------------------------- + # 1) Augment splines with extrema (CRITICAL) + # ------------------------------------------------- + zr_s = np.hstack([z_bot, zr, z_top]) + yr_s = np.hstack([y_bot, yr, y_top]) + + zl_s = np.hstack([z_bot, zl, z_top]) + yl_s = np.hstack([y_bot, yl, y_top]) + + # Ensure strict monotonicity for PCHIP while preserving original side curvature. + zr_s, yr_s = _prepare_strictly_increasing(zr_s, yr_s) + zl_s, yl_s = _prepare_strictly_increasing(zl_s, yl_s) + + # ------------------------------------------------- + # 2) Build splines + # ------------------------------------------------- + pchip_right = PchipInterpolator(zr_s, yr_s, extrapolate=False) + pchip_left = PchipInterpolator(zl_s, yl_s, extrapolate=False) + + # ------------------------------------------------- + # 3) Build z distribution + # ------------------------------------------------- + n_half = n_points // 2 + + if use_cosine_spacing: + beta = np.linspace(0.0, np.pi, n_half) + z_dist = z_bot + (z_top - z_bot) * 0.5 * (1 - np.cos(beta)) + else: + z_dist = np.linspace(z_bot, z_top, n_half, endpoint=True) + + # Remove extrema (added explicitly later) + z_mid = z_dist[1:-1] + + # ------------------------------------------------- + # 4) Evaluate splines + # ------------------------------------------------- + y_right = pchip_right(z_mid) + y_left = pchip_left(z_mid) + + # Enforce strict side ordering at each z: + # first branch is always +y side, second branch always -y side. + y_pos = np.maximum(y_right, y_left) + y_neg = np.minimum(y_right, y_left) + + # ------------------------------------------------- + # 5) Assemble CPACS closed loop + # ------------------------------------------------- + Airfoil = np.hstack([ + np.array([[y_bot], [z_bot]]), + np.vstack([y_pos, z_mid]), + np.array([[y_top], [z_top]]), + np.vstack([y_neg[::-1], z_mid[::-1]]), + np.array([[y_bot], [z_bot]]) + ]) + Airfoil = np.round(Airfoil, 3) + + if not hasattr(resample_fuselage_cpacs, "_debug_plot_saved"): + plt.figure() + plt.plot(yr, zr, ".", color="red", label="Input right (+y)") + plt.plot(yl, zl, ".", color="blue", label="Input left (-y)") + plt.plot(Airfoil[0, :], Airfoil[1, :], "-k", linewidth=1.2, label="Resampled profile") + plt.xlabel("y") + plt.ylabel("z") + plt.title("Fuselage Resampling") + plt.legend() + plt.axis("equal") + plt.grid() + plt.savefig("fuselage resampling") + plt.close() + resample_fuselage_cpacs._debug_plot_saved = True + + return Airfoil + + +def deduplicate_yz_points(y, z, tol=DUPLICATE_YZ_TOL): + """Remove duplicate/near-duplicate points in (y, z) with a tolerance grid.""" + + y = np.asarray(y, dtype=float) + z = np.asarray(z, dtype=float) + if y.size == 0: + return y, z + + q = np.round(np.column_stack((y, z)) / tol).astype(np.int64) + _, idx = np.unique(q, axis=0, return_index=True) + idx = np.sort(idx) + return y[idx], z[idx] + + +def find_top_bottom_near_centerline(y, z, y_tol_frac=0.08, z_band_frac=0.02): + """ + Find top and bottom points using only points close to symmetry plane (y ~= 0). + + Steps: + 1) Select points with |y| <= y_tol. + 2) Find highest and lowest z inside that subset. + 3) If multiple points are close to those extrema (z band), average them. + 4) Return indices of original points closest to these averaged top/bottom points. + """ + + y = np.asarray(y, dtype=float) + z = np.asarray(z, dtype=float) + + if y.size == 0: + return 0, 0 + + y_span = max(float(np.ptp(y)), 1e-12) + z_span = max(float(np.ptp(z)), 1e-12) + + y_tol = y_tol_frac * y_span + z_band = max(z_band_frac * z_span, 1e-12) + + # 1) points near y == 0 + center_mask = np.abs(y) <= y_tol + + # If too few points, enlarge tolerance. + if np.count_nonzero(center_mask) < 3: + center_mask = np.abs(y) <= (2.0 * y_tol) + + # If still too few, fall back to all points. + if np.count_nonzero(center_mask) < 2: + center_mask = np.ones_like(y, dtype=bool) + + zc = z[center_mask] + + # 2) raw top and bottom in center subset + z_top_raw = float(np.max(zc)) + z_bot_raw = float(np.min(zc)) + + # 3) average all points close to those extrema + top_group = center_mask & (np.abs(z - z_top_raw) <= z_band) + bot_group = center_mask & (np.abs(z - z_bot_raw) <= z_band) + + if np.count_nonzero(top_group) == 0: + top_group = center_mask & np.isclose(z, z_top_raw) + if np.count_nonzero(bot_group) == 0: + bot_group = center_mask & np.isclose(z, z_bot_raw) + + y_top = 0.0 + z_top = float(np.mean(z[top_group])) + y_bot = 0.0 + z_bot = float(np.mean(z[bot_group])) + + # 4) pick nearest original points (indices) to constrained centerline extrema + i_top = int(np.argmin((y - y_top) ** 2 + (z - z_top) ** 2)) + i_bot = int(np.argmin((y - y_bot) ** 2 + (z - z_bot) ** 2)) + + return i_top, i_bot + + +def extract_airfoil_surface_local(cloud_xyz, p0, n): + if cloud_xyz.shape[0] < 10: + return np.zeros((2, 0)), 0.0 + + # ------------------------------------------------- + # Define local in-plane basis (Y, Z) + # ------------------------------------------------- + e1 = np.array([0.0, 1.0, 0.0]) + e2 = np.array([0.0, 0.0, 1.0]) + + # ------------------------------------------------- + # Project cloud into slicing plane + # ------------------------------------------------- + local = np.array([ + [ + np.dot(p - p0, e1), + np.dot(p - p0, e2), + ] + for p in cloud_xyz + ]) + + y = local[:, 0] + z = local[:, 1] + # Left/right and top/bottom extrema detection + i_le = np.argmin(y) + i_te = np.argmax(y) + i_up, i_low = find_top_bottom_near_centerline( + y, z, y_tol_frac=0.08, z_band_frac=0.02 + ) + z_up = z[i_up] + z_low = z[i_low] + y_le = y[i_le] + y_te = y[i_te] + + width = y_te - y_le + height = z_up - z_low + + if width <= 1e-8 or height <= 1e-8: + return np.zeros((2, 0)), 0.0 + + # Normalize around section center so profile is centered at origin. + y_center = 0.5 * (y_le + y_te) + z_center = 0.5 * (z_up + z_low) + y = (y - y_center) / width + z = (z - z_center) / height + + # Enforce exact lateral normalization: min(y)=-0.5 and max(y)=+0.5 + y_min = float(np.min(y)) + y_max = float(np.max(y)) + y_span = y_max - y_min + if y_span > 1e-12: + y = (y - 0.5 * (y_max + y_min)) / y_span + + # Enforce exact extrema after normalization: + # bottom -> (y=0, z=-0.5), top -> (y=0, z=+0.5) + i_up, i_low = find_top_bottom_near_centerline( + y, z, y_tol_frac=0.08, z_band_frac=0.02 + ) + i_ymax = int(np.argmax(y)) + i_ymin = int(np.argmin(y)) + y[i_ymax] = 0.5 + y[i_ymin] = -0.5 + y[i_up] = 0.0 + y[i_low] = 0.0 + + # Remove duplicate points before splitting the surface. + y, z = deduplicate_yz_points(y, z, tol=DUPLICATE_YZ_TOL) + if y.size < 10: + return np.zeros((2, 0)), 0.0 + + # Recompute extrema after deduplication. + i_up, i_low = find_top_bottom_near_centerline( + y, z, y_tol_frac=0.08, z_band_frac=0.02 + ) + y_top_n, z_top_n = y[i_up], z[i_up] + y_bot_n, z_bot_n = y[i_low], z[i_low] + + n = 10 # number of bins for camber line + airfoil = split_fuselage_left_right_by_centerline( + y, z, y_bot_n, z_bot_n, y_top_n, z_top_n, n + ) + return airfoil, [width, height] + + +def split_fuselage_left_right_by_centerline( + y_raw, + z_raw, + y_bot, + z_bot, + y_top, + z_top, + n_bins +): + """ + Split fuselage cross-section into right (+y) and left (-y) sides + using a centerline computed from max/min y per z-bin. + + Top and bottom extrema are excluded from left/right classification. + """ + + y = np.asarray(y_raw) + z = np.asarray(z_raw) + + # Sort by z (vertical) + idx = np.argsort(z) + y = y[idx] + z = z[idx] + + # --- Detect bottom / top on symmetry plane (y ~= 0) --- + i_top, i_bot = find_top_bottom_near_centerline( + y, z, y_tol_frac=0.08, z_band_frac=0.02 + ) + z_top, y_top = z[i_top], y[i_top] + z_bot, y_bot = z[i_bot], y[i_bot] + + # --- Remove top & bottom from processing --- + mask_mid = np.ones(len(z), dtype=bool) + mask_mid[[i_bot, i_top]] = False + + y_mid = y[mask_mid] + z_mid = z[mask_mid] + + # --- Build bins only in interior --- + bins = np.linspace(z_bot, z_top, n_bins + 1) + + center_z = [] + center_y = [] + + for i in range(n_bins): + mask = (z_mid >= bins[i]) & (z_mid < bins[i + 1]) + if np.count_nonzero(mask) < 2: + continue + + y_bin = y_mid[mask] + z_bin = z_mid[mask] + + i_r = np.argmax(y_bin) + i_l = np.argmin(y_bin) + + center_y.append(0.5 * (y_bin[i_r] + y_bin[i_l])) + center_z.append(0.5 * (z_bin[i_r] + z_bin[i_l])) + + center_y = np.asarray(center_y) + center_z = np.asarray(center_z) + + # --- Interpolate centerline over interior points only --- + y_center_mid = np.interp(z_mid, center_z, center_y) + + # --- Classification (excluding extrema) --- + right_mask_mid = y_mid > y_center_mid + left_mask_mid = y_mid < y_center_mid + + # --- Rebuild full masks --- + right_mask = np.zeros(len(y), dtype=bool) + left_mask = np.zeros(len(y), dtype=bool) + + right_mask[np.where(mask_mid)[0][right_mask_mid]] = True + left_mask[np.where(mask_mid)[0][left_mask_mid]] = True + + if not hasattr(split_fuselage_left_right_by_centerline, "_debug_plot_saved"): + plt.figure() + for i, zb in enumerate(bins): + plt.hlines( + zb, + -0.6, + 0.6, + colors="black", + linestyles="--", + linewidth=1.2, + alpha=1, + label="Bin limits" if i == 0 else None, + ) + plt.plot(y[right_mask], z[right_mask], ".", color="red", label="Right (+y)") + plt.plot(y[left_mask], z[left_mask], ".", color="blue", label="Left (-y)") + plt.plot(center_y, center_z, "o", color="orange", ms=4, label="Center points") + plt.plot(y_center_mid, z_mid, "-k", label="Center line") + plt.xlabel("y") + plt.ylabel("z") + plt.title("Fuselage Points Classification") + plt.legend() + plt.axis("equal") + plt.grid() + plt.savefig("fuselage splitting") + plt.close() + split_fuselage_left_right_by_centerline._debug_plot_saved = True + + N_RESAMPLE_POINTS = 60 + return resample_fuselage_cpacs( + y[right_mask], z[right_mask], + y[left_mask], z[left_mask], + y_bot, z_bot, + y_top, z_top, + N_RESAMPLE_POINTS, + use_cosine_spacing=True, + ) + + +def intersect_triangle_with_plane_point_normal(p0, n, a, b, c, tol=INTERSECT_TOL): + da = np.dot(n, a - p0) + db = np.dot(n, b - p0) + dc = np.dot(n, c - p0) + pts = [] + + def edge_int(p1, d1, p2, d2): + if abs(d1) < tol and abs(d2) < tol: + # Both vertices lie on the plane + return [p1, p2] + if abs(d1) < tol: + # One vertex on plane + return [p1] + if abs(d2) < tol: + # One vertex above, one below. + # There is a parametric line equation P(t)= p1 + t*(p2 - p1) + return [p2] + if d1 * d2 < 0: + t = d1 / (d1 - d2) + return [p1 + t * (p2 - p1)] + # Edge does not intersect plane + return [] + + pts += edge_int(a, da, b, db) + pts += edge_int(b, db, c, dc) + pts += edge_int(c, dc, a, da) + if not pts: + return [] + uniq = [] + for p in pts: + if not any(np.linalg.norm(p - q) < 1e-10 for q in uniq): + # Sometimes the intersection produces duplicate points + uniq.append(p) + return uniq + + +def slice_mesh_rotated_YZ( + pts, + tris, + p0, + dihedral_deg, + sweep_deg, + tol=INTERSECT_TOL, + debug=False, +): + """ + Slice mesh with a plane orthogonal to local direction + defined by dihedral. + """ + + # Build local direction + a = np.deg2rad(dihedral_deg) + + Rx = np.array([ + [1, 0, 0], + [0, np.cos(a), -np.sin(a)], + [0, np.sin(a), np.cos(a)] + ]) + + RR = Rx + + e_span = RR @ np.array([1.0, 0.0, 0.0]) + e_span /= np.linalg.norm(e_span) + + # Signed distances + dverts = (pts - p0) @ e_span + dtri = dverts[tris] + + tri_min = dtri.min(axis=1) + tri_max = dtri.max(axis=1) + + hits = np.where((tri_min <= tol) & (tri_max >= -tol))[0] + if hits.size == 0: + return np.zeros((0, 3)), e_span + + # ------------------------------------------------- + # Intersections + # ------------------------------------------------- + inter = [] + for ti in hits: + i0, i1, i2 = tris[ti] + ip = intersect_triangle_with_plane_point_normal( + p0, e_span, + pts[i0], pts[i1], pts[i2], + tol=tol + ) + if ip: + inter.extend(ip) + + if not inter: + return np.zeros((0, 3)), e_span + + arr = np.vstack(inter) + + # Deduplicate + rtol = 1e-8 + key = np.round(arr / rtol).astype(np.int64) + dtype = np.dtype((np.void, key.dtype.itemsize * key.shape[1])) + _, idx = np.unique(key.view(dtype), return_index=True) + arr = arr[np.sort(idx)] + + # ------------------------------------------------- + # DEBUG PLOT + # ------------------------------------------------- + if debug: + import matplotlib.pyplot as plt + + fig = plt.figure(figsize=(9, 7)) + ax = fig.add_subplot(111, projection="3d") + + # Mesh + ax.scatter( + pts[:, 0], pts[:, 1], pts[:, 2], + s=1, alpha=0.1, color="gray", label="Mesh" + ) + + # Intersection points + ax.scatter( + arr[:, 0], arr[:, 1], arr[:, 2], + s=20, color="red", label="Slice points" + ) + + # Plane normal + L = np.linalg.norm(arr.max(axis=0) - arr.min(axis=0)) + ax.quiver( + p0[0], p0[1], p0[2], + e_span[0], e_span[1], e_span[2], + length=0.3 * L, + color="blue", + linewidth=3, + label="Plane normal" + ) + + # Plane visualization + u = np.linspace(-0.3 * L, 0.3 * L, 10) + v = np.linspace(-0.3 * L, 0.3 * L, 10) + U, V = np.meshgrid(u, v) + + # Two orthogonal vectors in plane + t1 = np.cross(e_span, [1, 0, 0]) + if np.linalg.norm(t1) < 1e-6: + t1 = np.cross(e_span, [0, 0, 1]) + t1 /= np.linalg.norm(t1) + t2 = np.cross(e_span, t1) + + Xp = p0[0] + U * t1[0] + V * t2[0] + Yp = p0[1] + U * t1[1] + V * t2[1] + Zp = p0[2] + U * t1[2] + V * t2[2] + + ax.plot_surface(Xp, Yp, Zp, alpha=0.25, color="cyan") + + ax.set_xlabel("X") + ax.set_ylabel("Y") + ax.set_zlabel("Z") + ax.set_title( + f"Rotated slice plane (dihedral={dihedral_deg:.2f} deg)" + ) + ax.legend() + ax.set_box_aspect([1, 1, 1]) + plt.tight_layout() + plt.show() + + return arr, e_span + + +def slice_mesh_at_Y(pts, tris, x_plane, tol): + """ + slicing with plane Y = y_plane + """ + + p0 = np.array([x_plane, 0.0, 0.0]) + n = np.array([1.0, 0.0, 0.0]) + dverts = (pts - p0) @ n + dtri = dverts[tris] + + tri_min = dtri.min(axis=1) + tri_max = dtri.max(axis=1) + + hits = np.where((tri_min <= tol) & (tri_max >= -tol))[0] + if hits.size == 0: + return np.zeros((0, 3)) + + inter = [] + for ti in hits: + + i0, i1, i2 = tris[ti] + ip = intersect_triangle_with_plane_point_normal( + p0, n, pts[i0], pts[i1], pts[i2], tol + ) + if ip: + inter.extend(ip) + + if not inter: + return np.zeros((0, 3)) + + arr = np.vstack(inter) + + # deduplicate + rtol = 1e-9 + key = np.round(arr / rtol).astype(np.int64) + dtype = np.dtype((np.void, key.dtype.itemsize * key.shape[1])) + _, idx = np.unique(key.view(dtype), return_index=True) + return arr[np.sort(idx)] + + +def compute_local_angles_from_ref(ref_pts): + """ + Compute sweep and dihedral from section reference points. + Sweep is defined in the XY plane. + Dihedral is defined in the XZ plane. + """ + + ref_pts = np.asarray(ref_pts) + M = ref_pts.shape[0] + if M < 2: + return np.array([]), np.array([]) + + sweep = np.zeros(M, dtype=int) + dihedral = np.zeros(M, dtype=int) + + for i in range(M - 1): + dx = ref_pts[i+1, 0] - ref_pts[i, 0] + dy = ref_pts[i+1, 1] - ref_pts[i, 1] + dz = ref_pts[i+1, 2] - ref_pts[i, 2] + + # ---- SWEEP: XY projection ---- + if abs(dx) < 1e-12 and abs(dy) < 1e-12: + sweep[i] = 0 + else: + sweep[i] = int(np.rint(np.degrees(np.arctan2(dy, dx)))) + + # ---- DIHEDRAL: XZ projection ---- + if abs(dx) < 1e-12 and abs(dz) < 1e-12: + dihedral[i] = 0 + else: + dihedral[i] = int(np.rint(np.degrees(np.arctan2(dz, dx)))) + + # copy last value + sweep[-1] = sweep[-2] + dihedral[-1] = dihedral[-2] + + return sweep, dihedral + + +def insert_slices(y_vals, sweep_deg, dihedral_deg, ref_pts, n_insert): + """ + Keep all original slices and only insert new slices at angle transitions. + """ + + y_vals = np.asarray(y_vals, dtype=float) + sweep_deg = np.asarray(sweep_deg, dtype=float) + dihedral_deg = np.asarray(dihedral_deg, dtype=float) + ref_pts = np.asarray(ref_pts, dtype=float) + + # Initialize output arrays with the first slice + y_out = [y_vals[0]] + sweep_out = [sweep_deg[0]] + dihedral_out = [dihedral_deg[0]] + ref_out = [ref_pts[0]] + is_inserted = [False] + + for i in range(len(y_vals) - 1): + # Transition if either sweep OR dihedral changes. + sweep_changed = not np.isclose(sweep_deg[i], sweep_deg[i + 1], atol=1e-9) + dihedral_changed = not np.isclose(dihedral_deg[i], dihedral_deg[i + 1], atol=1e-9) + has_transition = sweep_changed or dihedral_changed + + # Insert interpolated slices only at transitions. + if has_transition and n_insert > 0: + + for k in range(1, n_insert + 1): + t = k / (n_insert + 1) + y_new = (1 - t) * y_vals[i] + t * y_vals[i + 1] + ref_new = (1 - t) * ref_pts[i] + t * ref_pts[i + 1] + sweep_new = (1 - t) * sweep_deg[i] + t * sweep_deg[i + 1] + dihedral_new = (1 - t) * dihedral_deg[i] + t * dihedral_deg[i + 1] + y_out.append(y_new) + ref_out.append(ref_new) + sweep_out.append(sweep_new) + dihedral_out.append(dihedral_new) + is_inserted.append(True) + + # Always keep the next original slice. + y_out.append(y_vals[i + 1]) + ref_out.append(ref_pts[i + 1]) + sweep_out.append(sweep_deg[i + 1]) + dihedral_out.append(dihedral_deg[i + 1]) + is_inserted.append(False) + + return ( + np.array(y_out), + np.rint(sweep_out).astype(int), + np.rint(dihedral_out).astype(int), + np.array(ref_out), + np.array(is_inserted, dtype=bool), + ) + + +def compute_section_centers(clouds): + """Return per-slice geometric centers as Nx3 array.""" + + centers = [] + for cloud in clouds: + if cloud.shape[0] == 0: + continue + centers.append([ + np.mean(cloud[:, 0]), + 0.5 * (np.min(cloud[:, 1]) + np.max(cloud[:, 1])), + 0.5 * (np.min(cloud[:, 2]) + np.max(cloud[:, 2])), + ]) + return np.asarray(centers, dtype=float) + + +def plot_profile_diagnostics(airfoil_profiles): + """2D debug plots to identify section continuity anomalies.""" + + if not airfoil_profiles: + return + + n_sec = len(airfoil_profiles) + idx = np.arange(n_sec) + amp_y = np.array([np.max(np.abs(p[0])) for p in airfoil_profiles]) + amp_z = np.array([np.max(np.abs(p[1])) for p in airfoil_profiles]) + + fig, axes = plt.subplots(1, 2, figsize=(13, 5)) + + colors = cm.viridis(np.linspace(0, 1, n_sec)) + for i, p in enumerate(airfoil_profiles): + axes[0].plot(p[0], p[1], color=colors[i], lw=1.2, alpha=0.8, label=f"S{i}") + axes[0].plot(p[0][0], p[1][0], "o", color=colors[i], ms=2) + axes[0].set_title("Section Profiles Overlay") + axes[0].set_xlabel("y") + axes[0].set_ylabel("z") + axes[0].axis("equal") + axes[0].grid(True) + if n_sec <= 12: + axes[0].legend(fontsize=7, ncol=2) + + axes[1].plot(idx, amp_y, "-o", label="max |y|") + axes[1].plot(idx, amp_z, "-o", label="max |z|") + axes[1].set_title("Section Trend Check") + axes[1].set_xlabel("Section index") + axes[1].set_ylabel("Amplitude") + axes[1].grid(True) + axes[1].legend() + + plt.tight_layout() + plt.show() + + +# --------------------------- +# MAIN +# --------------------------- + + +def stl2fuselage_main(stl_file: str | Path, + setting: dict, + output_directory: str | Path, + name: str + ): + + tri_fname = export_mesh(output_directory, stl_file, name) + pts, tris = parse_cart3d_tri(tri_fname) + + # some initialization + airfoil_profiles = [] + Fuse_Dict = {} + per_slice_clouds = [] + bottom_points = [] + slice_x = [] + airfoil_profiles = [] + center_prev = None + per_slice_clouds_used = [] + base_idx = 0 + + # extract the setting + EXTREME_TOL_perc_start = setting['EXTREME_TOL_perc_start'] + EXTREME_TOL_perc_end = setting['EXTREME_TOL_perc_end'] + N_X_SLICES = setting['N_X_SLICES'] + N_SLICE_ADDING = setting['N_SLICE_ADDING'] + + # build X sampling positions + xmin, xmax = float(np.min(pts[:, 0])), float(np.max(pts[:, 0])) + EXTREME_TOL_start = EXTREME_TOL_perc_start * (xmax - xmin) + EXTREME_TOL_end = EXTREME_TOL_perc_end * (xmax - xmin) + x_start = xmin + EXTREME_TOL_start + x_end = xmax - EXTREME_TOL_end + x_vals_base = np.linspace(x_start, x_end, N_X_SLICES, dtype=float) + + # refine in the nose and tail regions + # Refine first and last 10% of effective fuselage length with +5 slices each. + extra_slices = 5 + refine_frac = 0.1 + x_len = x_end - x_start + nose_end = x_start + refine_frac * x_len + tail_start = x_end - refine_frac * x_len + + # Internal points only: avoid duplicating x_start/x_end. + x_nose_extra = np.linspace(x_start, nose_end, extra_slices + 2, dtype=float)[1:-1] + x_tail_extra = np.linspace(tail_start, x_end, extra_slices + 2, dtype=float)[1:-1] + + x_vals = np.unique(np.concatenate([x_vals_base, x_nose_extra, x_tail_extra])) + + # First slicing to get one reference point per slice (bottom point) + for i, x0 in enumerate(x_vals): + cloud = slice_mesh_at_Y(pts, tris, x0, INTERSECT_TOL) + + # if still empty, skip and record None + if cloud.shape[0] == 0: + per_slice_clouds.append(np.zeros((0, 3))) + bottom_points.append(None) + slice_x.append(x0) + continue + + # Bottom point: point with minimum Z in the slice + min_idx = int(np.argmin(cloud[:, 2])) + bottom_pt = cloud[min_idx].copy() + + per_slice_clouds.append(cloud) + bottom_points.append(bottom_pt) + slice_x.append(x0) + + _save_debug_stl_and_slices_plot( + pts=pts, + tris=tris, + per_slice_clouds=per_slice_clouds, + ref_points=bottom_points, + output_directory=output_directory, + name=name, + view_elev=float(setting.get("DEBUG_PLOT_ELEV", 12.0)), + view_azim=float(setting.get("DEBUG_PLOT_AZIM", 0.0)), + interactive=bool(setting.get("DEBUG_PLOT_INTERACTIVE", False)), + ) + + # build reference-point array + valid_idxs = [i for i, p in enumerate(bottom_points) if p is not None] + if len(valid_idxs) < 2: + raise RuntimeError("Too few bottom reference points found. Check mesh and N_X_SLICES.") + + per_slice_clouds_valid = [per_slice_clouds[i] for i in valid_idxs] + center_pts = compute_section_centers(per_slice_clouds_valid) + # start to build the dictionary. + Fuse_Dict["Transformation"] = { + "Name_type": "Fuselage", + "Name": str(name), + "X_Rot": [0, 0, 0], + "Symmetry": "", + "abs_system": True, + "Relative_dih": 0, + "Relative_Twist": 0, + "ParentUid": 0, + "reference_length": 0, + "idx_engine": None, + "Length": xmax - xmin + } + + # compute sweep & dihedral along section centers + sweep_deg, dihedral_deg = compute_local_angles_from_ref(center_pts) + + x_vals, sweep_deg, dihedral_deg, center_pts, is_inserted = insert_slices( + center_pts[:, 0], + sweep_deg, + dihedral_deg, + center_pts, + N_SLICE_ADDING, + ) + + for i, x0 in enumerate(x_vals): + center_ref = center_pts[i] + + if is_inserted[i]: + cloud = slice_mesh_at_Y(pts, tris, x0, INTERSECT_TOL) + if cloud.shape[0] == 0: + fb = min(max(base_idx - 1, 0), len(per_slice_clouds_valid) - 1) + cloud = per_slice_clouds_valid[fb] + else: + cloud = per_slice_clouds_valid[base_idx] + base_idx += 1 + + per_slice_clouds_used.append(cloud) + + # slice and rotate mesh + cloud_rot, n_rot = slice_mesh_rotated_YZ( + pts, + tris, + p0=center_ref, + dihedral_deg=dihedral_deg[i], + sweep_deg=sweep_deg[i], + tol=INTERSECT_TOL + ) + + airfoil_xz, Scaling = extract_airfoil_surface_local( + cloud_rot, + p0=center_ref, + n=n_rot, + ) + + # Section center in the slice plane (global coordinates) + y_center = center_ref[1] + z_center = center_ref[2] + + # Store in Fuse_Dict + if i == 0: + Fuse_Dict[f'Section{i}'] = { + 'x_scal': 1, + 'y_scal': round(Scaling[0], 2), + 'z_scal': round(Scaling[1], 2), + 'x_loc': x0, + 'y_trasl': 0.0, + 'z_trasl': 0.0, + 'x_rot': 0, + 'y_rot': 0, + 'z_rot': 0, + 'Span': 0, + 'Airfoil': 'Airfoil', + 'Airfoil_coordinates': airfoil_xz, + 'Sweep_loc': 0, + 'Sweep_angle': 0, + 'Dihedral_angle': 0 + } + Fuse_Dict["Transformation"]["X_Trasl"] = [ + x0, + y_center, + z_center, + ] + center_prev = (y_center, z_center) + + else: + Fuse_Dict[f'Section{i}'] = { + 'x_scal': 1, + 'y_scal': round(Scaling[0], 2), + 'z_scal': round(Scaling[1], 2), + 'x_loc': x0-x_start, + 'y_trasl': y_center - center_prev[0], + 'z_trasl': z_center - center_prev[1], + 'x_rot': 0, + 'y_rot': 0, + 'z_rot': 0, + 'Airfoil': 'Airfoil', + 'Airfoil_coordinates': airfoil_xz, + 'Sweep_loc': 0, + 'Sweep_angle': 0, + 'Dihedral_angle': 0 + } + + airfoil_profiles.append(airfoil_xz) + + return Fuse_Dict diff --git a/src/ceasiompy/stl2cpacs/func/stl2wing.py b/src/ceasiompy/stl2cpacs/func/stl2wing.py new file mode 100644 index 000000000..62c15ba5d --- /dev/null +++ b/src/ceasiompy/stl2cpacs/func/stl2wing.py @@ -0,0 +1,894 @@ + + +# ================================================================================================= +# IMPORTS +# ================================================================================================= + +import numpy as np +import os +import matplotlib.pyplot as plt +import struct +import matplotlib.cm as cm +from scipy.interpolate import PchipInterpolator +from ceasiompy.stl2cpacs.stl2cpacs import export_mesh, parse_cart3d_tri, load_stl_auto +from pathlib import Path + +# --------------------------- +# CONFIG +# --------------------------- + +INTERSECT_TOL = 1e-6 +SLAB_TOLS = [1e-5, 5e-5, 1e-4, 5e-4, 1e-3] +WING_AIRFOIL_ROUND_DECIMALS = 5 +WING_AIRFOIL_MAX_ROUND_DECIMALS = 8 +WING_AIRFOIL_MIN_SEG = 1e-5 +WING_CHORD_SCALE_DECIMALS = 4 +WING_ANGLE_DECIMALS = 4 + +# ================================================================================================= +# FUNCTIONS +# ================================================================================================= + + +def _save_debug_stl_and_slices_plot( + pts: np.ndarray, + tris: np.ndarray, + per_slice_clouds: list[np.ndarray], + le_points: list[np.ndarray | None], + output_directory: str | Path, + name: str, +) -> None: + """ + Save a debug 3D figure with STL points, first slice clouds, and LE picks. + """ + + out_dir = Path(output_directory) + out_dir.mkdir(parents=True, exist_ok=True) + debug_path = out_dir / f"{name}_debug_stl_first_slices.png" + + fig = plt.figure() + ax = fig.add_subplot(111, projection="3d") + + # STL as semi-transparent surface for better visibility than sparse points. + if tris.shape[0] > 0: + max_tris = 50000 + if tris.shape[0] > max_tris: + tri_idx = np.linspace(0, tris.shape[0] - 1, max_tris, dtype=int) + tris_plot = tris[tri_idx] + else: + tris_plot = tris + ax.plot_trisurf( + pts[:, 0], + pts[:, 1], + pts[:, 2], + triangles=tris_plot, + color="#4C78A8", + alpha=0.24, + linewidth=0.0, + shade=False, + ) + # Proxy handle so STL appears in legend. + ax.plot([], [], [], color="#4C78A8", lw=4, alpha=0.85, label="STL mesh") + + cmap = cm.get_cmap("viridis", max(1, len(per_slice_clouds))) + first_slice_label = True + for i, cloud in enumerate(per_slice_clouds): + if cloud.shape[0] == 0: + continue + ax.scatter( + cloud[:, 0], + cloud[:, 1], + cloud[:, 2], + s=5, + color=cmap(i), + alpha=0.8, + label="Slice intersections" if first_slice_label else None, + ) + first_slice_label = False + + ax.set_title("Wing slicing ") + ax.set_xlabel("X") + ax.set_ylabel("Y") + ax.set_zlabel("Z") + ax.grid(True, alpha=0.3) + ax.legend(loc="upper right", fontsize=18) + ax.axis('equal') + fig.savefig(debug_path, dpi=200) + plt.close(fig) + + +def _remove_consecutive_duplicate_points(poly, tol=1e-12): + """ + Remove consecutive duplicate points in a 2xN polyline while preserving closure. + """ + + if poly.shape[1] <= 2: + return poly + + keep = [0] + for i in range(1, poly.shape[1]): + if np.hypot(poly[0, i] - poly[0, keep[-1]], poly[1, i] - poly[1, keep[-1]]) > tol: + keep.append(i) + + out = poly[:, keep] + if out.shape[1] >= 2 and np.hypot(out[0, 0] - out[0, -1], out[1, 0] - out[1, -1]) > tol: + out = np.hstack([out, out[:, [0]]]) + return out + + +def _round_airfoil_safely(poly, pref_dec=4, max_dec=8, min_seg=1e-5): + """ + Round airfoil coordinates for CPACS while avoiding degenerate thin geometries. + """ + + base = _remove_consecutive_duplicate_points(np.asarray(poly, dtype=float), tol=1e-12) + if base.shape[1] < 5: + return base + + for dec in range(pref_dec, max_dec + 1): + cand = np.round(base, dec) + cand = _remove_consecutive_duplicate_points(cand, tol=10 ** (-(dec + 2))) + if cand.shape[1] < 5: + continue + seg = np.hypot(np.diff(cand[0]), np.diff(cand[1])) + if seg.size == 0: + continue + if np.min(seg) >= min_seg: + return cand + + # Fallback: keep original precision and only clean exact consecutive duplicates. + return base + + +def fix_airfoil_cpacs(x, z, tol_x): + """ + Remove duplicate / near-duplicate airfoil points. + + If two points have |x1 - x2| < tol_x → keep only one + Keep the point that is farther from the local mean + + """ + + x = np.asarray(x, dtype=float) + z = np.asarray(z, dtype=float) + if x.size == 0: + return x, z + + # Preserve raw extrema references before any cleanup. + x_te_raw = float(np.max(x)) + x_le_raw = float(np.min(x)) + z_te_raw = float(np.mean(z[np.isclose(x, x_te_raw, atol=tol_x)])) + z_le_raw = float(np.mean(z[np.isclose(x, x_le_raw, atol=tol_x)])) + + # Sort by x for robust deduplication in profile parameter space. + idx = np.argsort(x) + x = x[idx] + z = z[idx] + + x_clean = [x[0]] + z_clean = [z[0]] + for i in range(1, len(x)): + dx = abs(x[i] - x_clean[-1]) + dz = abs(z[i] - z_clean[-1]) + + # Keep TE/LE neighborhood points to avoid losing extrema. + near_te = (abs(x[i] - x_te_raw) <= tol_x) or (abs(x_clean[-1] - x_te_raw) <= tol_x) + near_le = (abs(x[i] - x_le_raw) <= tol_x) or (abs(x_clean[-1] - x_le_raw) <= tol_x) + + if dx < tol_x and dz < tol_x and not (near_te or near_le): + # Skip near-duplicate interior point. + continue + + x_clean.append(x[i]) + z_clean.append(z[i]) + + x_clean = np.asarray(x_clean, dtype=float) + z_clean = np.asarray(z_clean, dtype=float) + + # Enforce LE/TE presence explicitly. + if not np.any(np.isclose(x_clean, x_te_raw, atol=tol_x)): + x_clean = np.append(x_clean, x_te_raw) + z_clean = np.append(z_clean, z_te_raw) + if not np.any(np.isclose(x_clean, x_le_raw, atol=tol_x)): + x_clean = np.append(x_clean, x_le_raw) + z_clean = np.append(z_clean, z_le_raw) + + # Restore x-sorted order after possible appends. + idx = np.argsort(x_clean) + x_clean = x_clean[idx] + z_clean = z_clean[idx] + + return np.array(x_clean), np.array(z_clean) + + +def resample_airfoil_cpacs( + xu, zu, + xl, zl, + x_te, z_te, + x_le, z_le, + n_points, + +): + """ + Regularize an airfoil by spline-fitting upper/lower surfaces + and resampling with either uniform or cosine spacing in x. + + xu, zu : array-like + Upper surface points. + xl, zl : array-like + Lower surface points. + x_te, z_te : float + Trailing edge coordinates. + n_points : int + Total number of points for the closed airfoil polyline. + """ + + use_cosine_spacing = False + xu = np.asarray(xu) + zu = np.asarray(zu) + xl = np.asarray(xl) + zl = np.asarray(zl) + # Sort surfaces + # Upper: LE -> TE + idx_u = np.argsort(xu) + xu, zu = xu[idx_u], zu[idx_u] + xu[-1] = x_te + zu[-1] = z_te + + # Lower: TE -> LE + idx_l = np.argsort(xl)[::-1] + xl, zl = xl[idx_l], zl[idx_l] + xl[0] = x_te + zl[0] = z_te + + # Build x-distribution on LE → TE + n_half = n_points // 2 + if use_cosine_spacing: + beta = np.linspace(0.0, np.pi, n_half) + x_dist = x_le + (x_te - x_le) * 0.5 * (1 - np.cos(beta)) + else: + x_dist = np.linspace(x_le, x_te, n_half) + + x_dist = x_dist[1:-1] + # PCHIP interpolation + pchip_upper = PchipInterpolator(xu, zu, extrapolate=False) + pchip_lower = PchipInterpolator(xl[::-1], zl[::-1], extrapolate=False) + + z_upper = pchip_upper(x_dist) + z_lower = pchip_lower(x_dist) + + x_u, z_u = x_dist, z_upper + x_l, z_l = x_dist[::-1], z_lower[::-1] + # CPACS order + x_te = 1.0 + x_le = 0.0 + airfoil = np.hstack([ + np.array([[x_te], [z_te]]), + np.vstack([x_l, z_l]), + np.array([[x_le], [z_le]]), + np.vstack([x_u, z_u]), + np.array([[x_te], [z_te]]) + ]) + airfoil = _round_airfoil_safely( + airfoil, + pref_dec=WING_AIRFOIL_ROUND_DECIMALS, + max_dec=WING_AIRFOIL_MAX_ROUND_DECIMALS, + min_seg=WING_AIRFOIL_MIN_SEG, + ) + + if not hasattr(resample_airfoil_cpacs, "_debug_plot_saved"): + plt.figure() + plt.plot(airfoil[0, :np.shape(airfoil)[1]//2], airfoil[1, :np.shape(airfoil)[1]//2], + ".", color="red", + label="Input upper" + ) + plt.plot(airfoil[0, np.shape(airfoil)[1]//2:-1], + airfoil[1, np.shape(airfoil)[1]//2:-1], + ".", color="blue", + label="Input lower" + ) + plt.plot(airfoil[0, :], airfoil[1, :], "-k", linewidth=1.2, label="Resampled profile") + plt.xlabel("x/c") + plt.ylabel("z/c") + plt.title("Airfoil Resampling") + plt.legend() + plt.axis("equal") + plt.grid() + plt.savefig("airfoil resampling") + plt.close() + resample_airfoil_cpacs._debug_plot_saved = True + + return airfoil + + +def extract_airfoil_surface_local(cloud_xyz, p0, n, N_BIN, TE_CUT): + + if cloud_xyz.shape[0] < 10: + return np.zeros((2, 0)), 0.0 + + n = n / np.linalg.norm(n) + + # Local basis + ex = np.array([1.0, 0.0, 0.0]) + e1 = ex - np.dot(ex, n) * n + if np.linalg.norm(e1) < 1e-10: + return np.zeros((2, 0)), 0.0 + e1 /= np.linalg.norm(e1) + + e2 = np.cross(n, e1) + e2 /= np.linalg.norm(e2) + + # Project STL cloud + local = np.array([ + [np.dot(p - p0, e1), np.dot(p - p0, -e2)] + for p in cloud_xyz + ]) + + x = local[:, 0] + z = local[:, 1] + + # LE / TE detection + i_le = np.argmin(x) + i_te = np.argmax(x) + + x_le = x[i_le] + x_te = x[i_te] + + chord = abs(x_te - x_le) + + if chord <= 1e-8: + return np.zeros((2, 0)), 0.0 + + # Normalize ONCE + x = (x - x_le) / chord + z = z / chord + + airfoil = split_upper_lower_by_camber(x, z, N_BIN, TE_CUT) + + return airfoil, chord + + +def split_upper_lower_by_camber(x_raw, z_raw, n_bins, te_cut): + """ + upper/lower split using camber line only where reliable. + The TE zone is not split to avoid sharp-edge issues. + """ + + x, z = fix_airfoil_cpacs(x_raw, z_raw, 1e-4) + + # Sort by x + idx = np.argsort(x) + x = x[idx] + z = z[idx] + + # Detect LE / TE + i_le = np.argmin(x) + i_te = np.argmax(x) + + x_le, z_le = x[i_le], z[i_le] + x_te, z_te = x[i_te], z[i_te] + + # Define camber-valid zone + camber_mask = x < (x_te - te_cut) + x_camber = x[camber_mask] + z_camber = z[camber_mask] + + # Build coarse camber line + bins = np.linspace(x_le, x_te - te_cut, n_bins + 1) + + camber_x = [x_le] + camber_z = [z_le] + + for i in range(n_bins): + mask = (x_camber >= bins[i]) & (x_camber < bins[i + 1]) + if np.count_nonzero(mask) < 2: + continue + + x_bin = x_camber[mask] + z_bin = z_camber[mask] + # inside every bin, the camber point is done + # using the point with the maximum z and an other one wiht the minimum z + i_up = np.argmax(z_bin) + i_lo = np.argmin(z_bin) + camber_x.append(0.5 * (x_bin[i_up] + x_bin[i_lo])) + camber_z.append(0.5 * (z_bin[i_up] + z_bin[i_lo])) + + camber_x = np.array(camber_x) + camber_z = np.array(camber_z) + + # Interpolate camber on camber zone + zc = np.interp(x_camber, camber_x, camber_z) + + # Classification + upper_mask = np.zeros_like(x, dtype=bool) + lower_mask = np.zeros_like(x, dtype=bool) + + upper_mask[camber_mask] = z_camber > zc + lower_mask[camber_mask] = z_camber < zc + + # Resample + return resample_airfoil_cpacs( + x[upper_mask], z[upper_mask], + x[lower_mask], z[lower_mask], + x_te=x_te, + z_te=z_te, + x_le=x_le, + z_le=z_le, + n_points=60 + ) + + +def intersect_triangle_with_plane_point_normal(p0, n, a, b, c, tol=INTERSECT_TOL): + + da = np.dot(n, a - p0) + db = np.dot(n, b - p0) + dc = np.dot(n, c - p0) + pts = [] + + def edge_int(p1, d1, p2, d2): + if abs(d1) < tol and abs(d2) < tol: + # Both vertices lie on the plane + return [p1, p2] + if abs(d1) < tol: + # One vertex on plane + return [p1] + if abs(d2) < tol: + # One vertex above, one below. + # There is a parametric line equation P(t)= p1 + t*(p2 - p1) + return [p2] + if d1 * d2 < 0: + t = d1 / (d1 - d2) + return [p1 + t * (p2 - p1)] + # Edge does not intersect plane + return [] + pts += edge_int(a, da, b, db) + pts += edge_int(b, db, c, dc) + pts += edge_int(c, dc, a, da) + if not pts: + return [] + uniq = [] + for p in pts: + if not any(np.linalg.norm(p - q) < 1e-10 for q in uniq): + # Sometimes the intersection produces duplicate points + uniq.append(p) + return uniq + + +def slice_mesh_rotated_YZ( + pts, + tris, + p0, + dihedral_deg, + sweep_deg, + tol=INTERSECT_TOL, +): + """ + Slice mesh with a plane orthogonal to local span direction + defined by dihedral. + """ + + # Build span direction + a = np.deg2rad(dihedral_deg) + + Rx = np.array([ + [1, 0, 0], + [0, np.cos(a), -np.sin(a)], + [0, np.sin(a), np.cos(a)] + ]) + + RR = Rx + + e_span = RR @ np.array([0.0, 1.0, 0.0]) + e_span /= np.linalg.norm(e_span) + + # Signed distances + dverts = (pts - p0) @ e_span + dtri = dverts[tris] + + tri_min = dtri.min(axis=1) + tri_max = dtri.max(axis=1) + + hits = np.where((tri_min <= tol) & (tri_max >= -tol))[0] + if hits.size == 0: + return np.zeros((0, 3)), e_span + + # ------------------------------------------------- + # Intersections + # ------------------------------------------------- + inter = [] + for ti in hits: + i0, i1, i2 = tris[ti] + ip = intersect_triangle_with_plane_point_normal( + p0, e_span, + pts[i0], pts[i1], pts[i2], + tol=tol + ) + if ip: + inter.extend(ip) + + if not inter: + return np.zeros((0, 3)), e_span + + arr = np.vstack(inter) + + # Deduplicate + rtol = 1e-8 + key = np.round(arr / rtol).astype(np.int64) + dtype = np.dtype((np.void, key.dtype.itemsize * key.shape[1])) + _, idx = np.unique(key.view(dtype), return_index=True) + arr = arr[np.sort(idx)] + + return arr, e_span + + +def slice_mesh_at_Y(pts, tris, y_plane, tol): + """ + Slicing with plane Y = y_plane + """ + + p0 = np.array([0.0, y_plane, 0.0]) + n = np.array([0.0, 1.0, 0.0]) + dverts = (pts - p0) @ n + dtri = dverts[tris] + tri_min = dtri.min(axis=1) + tri_max = dtri.max(axis=1) + + hits = np.where((tri_min <= tol) & (tri_max >= -tol))[0] + if hits.size == 0: + return np.zeros((0, 3)) + + inter = [] + for ti in hits: + + i0, i1, i2 = tris[ti] + ip = intersect_triangle_with_plane_point_normal( + p0, n, pts[i0], pts[i1], pts[i2], tol + ) + if ip: + inter.extend(ip) + + if not inter: + return np.zeros((0, 3)) + + arr = np.vstack(inter) + + # deduplicate + rtol = 1e-9 + key = np.round(arr / rtol).astype(np.int64) + dtype = np.dtype((np.void, key.dtype.itemsize * key.shape[1])) + _, idx = np.unique(key.view(dtype), return_index=True) + return arr[np.sort(idx)] + + +def compute_local_angles_from_le(le_pts): + """ + Compute sweep and dihedral from LE points. + Sweep is defined in the XY plane. + Dihedral is defined in the YZ plane. + """ + + le_pts = np.asarray(le_pts) + M = le_pts.shape[0] + if M < 2: + return np.array([]), np.array([]) + + sweep = np.zeros(M, dtype=int) + dihedral = np.zeros(M, dtype=int) + + for i in range(M - 1): + dx = le_pts[i+1, 0] - le_pts[i, 0] + dy = le_pts[i+1, 1] - le_pts[i, 1] + dz = le_pts[i+1, 2] - le_pts[i, 2] + + # ---- SWEEP: XY projection ---- + if abs(dy) < 1e-12: + sweep[i] = 0 + else: + sweep[i] = np.round((np.rint(np.degrees(np.arctan(dx / np.sqrt(dy**2 + dz**2))))), + WING_ANGLE_DECIMALS + ) + # ---- DIHEDRAL: YZ projection ---- + if abs(dy) < 1e-12: + dihedral[i] = 0 + else: + dihedral[i] = np.round((np.rint(np.degrees(np.arctan(dz / dy)))), + WING_ANGLE_DECIMALS + ) + + # copy last value + sweep[-1] = sweep[-2] + dihedral[-1] = dihedral[-2] + + return sweep, dihedral + + +def filter_and_insert(y_vals, sweep_deg, dihedral_deg, le_pts, n_insert): + """ + Refine wing sections for CPACS generation. + + Behavior: + - Keep first slice. + - In constant-angle regions, keep only boundary slices (filter interior). + - When angles change between i and i+1, insert n_insert interpolated slices. + - Keep the slice at the middle to be sure that also symmetrical wings are well treated. + - Keep last slice. + """ + + y_vals = np.asarray(y_vals, dtype=float) + sweep_deg = np.asarray(sweep_deg, dtype=float) + dihedral_deg = np.asarray(dihedral_deg, dtype=float) + le_pts = np.asarray(le_pts, dtype=float) + + # Initialize output arrays with the first slice + y_out = [y_vals[0]] + sweep_out = [sweep_deg[0]] + dihedral_out = [dihedral_deg[0]] + le_out = [le_pts[0]] + + for i in range(len(y_vals) - 1): + boll_add_slice = ( + np.isclose(sweep_deg[i], sweep_deg[i + 1], atol=0.1, rtol=0.0) and + np.isclose(dihedral_deg[i], dihedral_deg[i + 1], atol=0.1, rtol=0.0) and + (i < len(y_vals) // 2 - 1 or i > len(y_vals) // 2 + 1) + ) + + if not boll_add_slice and n_insert > 0: + # Interpolate transition slices with linear angle evolution. + for k in range(1, n_insert + 1): + t = k / (n_insert + 1) + y_new = (1 - t) * y_vals[i] + t * y_vals[i + 1] + le_new = (1 - t) * le_pts[i] + t * le_pts[i + 1] + sweep_new = (1 - t) * sweep_deg[i] + t * sweep_deg[i + 1] + dihedral_new = (1 - t) * dihedral_deg[i] + t * dihedral_deg[i + 1] + y_out.append(y_new) + le_out.append(le_new) + sweep_out.append(sweep_new) + dihedral_out.append(dihedral_new) + + # Keep i+1 when: + # - entering or leaving a transition region + # - always for final slice + is_last_pair = (i == len(y_vals) - 2) + next_changes = False if is_last_pair else not ( + np.isclose(sweep_deg[i + 1], sweep_deg[i + 2]) and + np.isclose(dihedral_deg[i + 1], dihedral_deg[i + 2]) + ) + keep_next = (not boll_add_slice) or next_changes or is_last_pair + + if keep_next: + y_out.append(y_vals[i + 1]) + le_out.append(le_pts[i + 1]) + sweep_out.append(sweep_deg[i + 1]) + dihedral_out.append(dihedral_deg[i + 1]) + + return ( + np.array(y_out), + np.array(sweep_out, dtype=float), + np.array(dihedral_out, dtype=float), + np.array(le_out), + ) + + +def rotate_vertical_tail(stl_file): + """Rotate STL 90 deg around X while keeping the lowest point fixed.""" + + stl_path = Path(stl_file) + if not stl_path.exists(): + raise FileNotFoundError(f"STL file not found: {stl_path}") + + tris = load_stl_auto(str(stl_path)).astype(float) + if tris.size == 0: + return str(stl_path) + + pts = tris.reshape(-1, 3) + # Use the lowest point (minimum global Z) as rotation pivot + # so the geometry does not drift in space after rotation. + pivot = pts[int(np.argmin(pts[:, 2]))].copy() + + rot_x_neg_90 = np.array( + [ + [1.0, 0.0, 0.0], + [0.0, 0.0, 1.0], + [0.0, -1.0, 0.0], + ], + dtype=float, + ) + + rotated_pts = (pts - pivot) @ rot_x_neg_90.T + pivot + rotated_tris = rotated_pts.reshape(tris.shape) + + out_path = stl_path.with_name(f"{stl_path.stem}_rotated_vertical_tail.stl") + with open(out_path, "wb") as f: + header = b"CEASIOMpy rotated vertical tail" + f.write(header.ljust(80, b"\0")) + f.write(struct.pack(" 0.0: + n = n / n_norm + else: + n = np.array([0.0, 0.0, 0.0], dtype=float) + + f.write(struct.pack(" list: + """ + Automatically detects if STL represents a wing, vertical tail, or fuselage + using PCA. + + Criteria: + - Wing: maximum variance approximately aligned with global Y axis + - Vertical tail: maximum variance approximately aligned with global Z axis + - Fuselage: maximum variance approximately aligned with global X axis + + Returns: + List[str]: ["wing"], ["wing_vertical_tail"], or ["fuselage"] + """ + + cpacs_component = [] + + if isinstance(stl_file, (str, Path)): + stl_file = [stl_file] + + for path in stl_file: + + # --------------------------------------------------------- + # 1) Load STL and extract unique vertices + # --------------------------------------------------------- + tris = load_stl_auto(path) + points = tris.reshape(-1, 3) + points = np.unique(points, axis=0) + + # --------------------------------------------------------- + # 2) Center the geometry + # --------------------------------------------------------- + # PCA must be applied on centered data. + # We subtract the centroid so that variance describes + # only shape distribution and not absolute position. + centroid = np.mean(points, axis=0) + centered = points - centroid + + # --------------------------------------------------------- + # 3) Principal Component Analysis (PCA) + # --------------------------------------------------------- + # Covariance matrix describes how geometry spreads in space. + # Eigenvectors give principal directions of variance. + # Eigenvalues give amount of variance in each direction. + cov = np.cov(centered.T) + eigvals, eigvecs = np.linalg.eigh(cov) + + # Sort eigenvalues descending so eigvecs[:,0] corresponds + # to direction of maximum geometric development. + order = np.argsort(eigvals)[::-1] + eigvals = eigvals[order] + eigvecs = eigvecs[:, order] + + # --------------------------------------------------------- + # 4) Principal direction (maximum variance direction) + # --------------------------------------------------------- + # This vector tells us along which axis the STL is + # developing the most (longest geometric dimension). + principal_dir = eigvecs[:, 0] + + # Normalize to make it a unit vector. + # IMPORTANT: + # For a unit vector v = [vx, vy, vz], + # each component is a DIRECTION COSINE: + # + # vx = cos(angle with global X axis) + # vy = cos(angle with global Y axis) + # vz = cos(angle with global Z axis) + # + # This comes from the dot product definition: + # v · ex = |v||ex| cos(theta_x) + # Since both are unit vectors, v · ex = cos(theta_x). + # + # Therefore the components of a normalized vector + # directly measure angular alignment with axes. + principal_dir = principal_dir / np.linalg.norm(principal_dir) + + # --------------------------------------------------------- + # 5) Axis alignment (direction cosines) + # --------------------------------------------------------- + # We take absolute value because orientation sign + # does not matter: + # +X and -X are both fuselage directions. + # + # Example: + # align_x > 0.75 -> angle < acos(0.75) ≈ 41° + # align_y > 0.6 -> angle < acos(0.6) ≈ 53° + # + # So we are effectively testing whether the principal + # direction lies inside a cone around a global axis. + align_x = abs(principal_dir[0]) + align_y = abs(principal_dir[1]) + align_z = abs(principal_dir[2]) + + # --------------------------------------------------------- + # 6) Variance ratio (elongation check) + # --------------------------------------------------------- + # For a fuselage, geometry is strongly elongated in X. + # So first eigenvalue should be significantly larger + # than the second. + # + # If ratio is close to 1 → geometry is more isotropic. + var_ratio = eigvals[0] / eigvals[1] + + # --------------------------------------------------------- + # 7) Classification logic + # --------------------------------------------------------- + + # Fuselage: + # - Strong elongation + # - Principal direction mostly aligned with X + if align_x > 0.75 and var_ratio > 1.2: + cpacs_component.append("fuselage") + + # Vertical tail: + # - Dominant direction along Z + elif align_z > 0.6: + cpacs_component.append("wing_vertical_tail") + + # Wing: + # - Dominant direction along Y + # - Works even with sweep or winglets + elif align_y > 0.6: + cpacs_component.append("wing") + + # Fallback: choose axis with strongest alignment + else: + if align_x > align_y and align_x > align_z: + cpacs_component.append("fuselage") + elif align_z > align_y: + cpacs_component.append("wing_vertical_tail") + else: + cpacs_component.append("wing") + + return cpacs_component + + +def load_stl_auto(path): + with open(path, "rb") as f: + start = f.read(80) + if start[:5].lower() == b"solid": + try: + return read_ascii_stl(path) + except ValueError: + return read_binary_stl(path) + return read_binary_stl(path) + + +def write_cart3d_tri(filename, triangles): + """ + Saves triangles to Cart3D .tri format + """ + + verts = triangles.reshape(-1, 3) + uniq, inverse = np.unique(verts, axis=0, return_inverse=True) + tri_idx = inverse.reshape(-1, 3) + 1 # Cart3D uses 1-based indexing + + with open(filename, "w") as f: + f.write(f"{uniq.shape[0]} {tri_idx.shape[0]}\n") # first line + for v in uniq: + f.write(f"{v[0]:.9g} {v[1]:.9g} {v[2]:.9g}\n") # vertices + for t in tri_idx: + f.write(f"{t[0]} {t[1]} {t[2]}\n") # triangle + + return filename + + +def read_ascii_stl(path): + """Reads ASCII STL and returns Nx3x3 triangle array""" + + tri = [] + with open(path, "r") as f: + for line in f: + if line.strip().startswith("vertex"): + _, x, y, z = line.split() + tri.append([float(x), float(y), float(z)]) + tri = np.array(tri).reshape(-1, 3, 3) + return tri + + +def read_binary_stl(path): + """Reads binary STL and returns Nx3x3 triangle array""" + with open(path, "rb") as f: + # Binary STL has a mandatory 80-byte header; consume it to align stream. + _ = f.read(80) + ntri = struct.unpack(" tuple[Path, dict[str, list[dict[str, str]]]]: + """Convert STL components to one CPACS file and return output XML path + report.""" + + # Local imports although it shows errors + from ceasiompy.stl2cpacs.func.stl2wing import stl2wing_main + from ceasiompy.stl2cpacs.func.stl2fuselage import stl2fuselage_main + + cpacs_component = cpacs_component_detection(stl_file=stl_file) + + if not (len(stl_file) == len(cpacs_component) == len(setting)): + raise ValueError( + "stl_file, cpacs_component and setting must have the same length." + ) + + output_dir = Path(out_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + dict_exportcpacs = {} + conversion_report: dict[str, list[dict[str, str]]] = { + "converted": [], + "failed": [], + } + for idx, (stl_path, item, comp_setting) in enumerate( + zip(stl_file, cpacs_component, setting), start=1 + ): + item_norm = str(item).strip().lower() + + comp_name = Path(stl_path).stem + try: + if item_norm == "wing" or item_norm == "wing_vertical_tail": + effective_setting = { + "EXTREME_TOL_perc_start": 0.02, + "EXTREME_TOL_perc_end": 0.02, + "N_Y_SLICES": 50, + "N_SLICE_ADDING": 0, + "TE_CUT": 0.0, + "N_BIN": 10, + "vertical_tail": True if item_norm == "wing_vertical_tail" else False, + } + effective_setting.update(comp_setting) + dict_exportcpacs[f"component_{idx}"] = stl2wing_main( + stl_path, effective_setting, output_dir, comp_name + ) + conversion_report["converted"].append( + {"name": comp_name, "path": str(stl_path), "type": item_norm} + ) + elif item_norm == "fuselage": + effective_setting = { + "EXTREME_TOL_perc_start": 0.02, + "EXTREME_TOL_perc_end": 0.02, + "N_X_SLICES": 50, + "N_SLICE_ADDING": 0, + } + effective_setting.update(comp_setting) + dict_exportcpacs[f"component_{idx}"] = stl2fuselage_main( + stl_path, effective_setting, output_dir, comp_name + ) + conversion_report["converted"].append( + {"name": comp_name, "path": str(stl_path), "type": item_norm} + ) + else: + conversion_report["failed"].append( + { + "name": comp_name, + "path": str(stl_path), + "type": item_norm, + "error": f"Unsupported detected type: {item_norm}", + } + ) + except Exception as exc: + conversion_report["failed"].append( + { + "name": comp_name, + "path": str(stl_path), + "type": item_norm, + "error": str(exc), + } + ) + + if not dict_exportcpacs: + raise RuntimeError( + "No component could be converted to CPACS. " + "Check detected types and STL quality." + ) + exporter = Export_CPACS(dict_exportcpacs, "stl2cpacs_file", output_dir) + return exporter.run(), conversion_report diff --git a/src/ceasiompy/utils/exportcpacs.py b/src/ceasiompy/utils/exportcpacs.py new file mode 100644 index 000000000..ec8de2bd7 --- /dev/null +++ b/src/ceasiompy/utils/exportcpacs.py @@ -0,0 +1,1077 @@ + +""" +CEASIOMpy: Conceptual Aircraft Design Software + +Developed by CFS ENGINEERING, 1015 Lausanne, Switzerland + +The geometry is built in OpenVSP, saved as a .vsp3 file, and then selected in the GUI. +It is subsequently processed by this module to generate a CPACS file. + +| Author: Nicolo Perasso +| Creation: 23/12/2025 +""" + +# Imports +import re +import defusedxml +import numpy as np +import xml.etree.ElementTree as ET # nosec B405 + +from pathlib import Path +from defusedxml import ElementTree as DefusedElementTree + +from ceasiompy import log + + +defuse = getattr(defusedxml, "defuse_stdlib", None) +if defuse is None: # pragma: no cover + raise ImportError("defusedxml does not support defuse_stdlib in this version.") +defuse() + + +# Functions +def safe_parse_xml(*args, **kwargs): + return DefusedElementTree.parse(*args, **kwargs) + + +def safe_fromstring_xml(text: str): + return DefusedElementTree.fromstring(text) + + +def make(doc, name, parent=None, text=None, **attrs): + """ + Create an XML element, optionally add attributes, text value, + and append to a parent element. + """ + node = doc.createElement(name) + + # set attributes + for k, v in attrs.items(): + node.setAttribute(k, str(v)) + + # optional text + if text is not None: + node.appendChild(doc.createTextNode(str(text))) + + # optional parent + if parent is not None: + parent.appendChild(node) + + return node + + +def add_text(doc, name, parent, value): + make(doc, name, parent, text=value) + + +def segment(doc, Parent, Name_wing, Name_Element, Name_Element_before, idx): + segment = doc.createElement('segment') + segment.setAttribute('uID', f'{Name_wing}Seg{idx}') + Parent.appendChild(segment) + name = doc.createElement('name') + name.appendChild(doc.createTextNode(f'{Name_wing}Seg{idx}')) + segment.appendChild(name) + FromElement = doc.createElement('fromElementUID') + FromElement.appendChild(doc.createTextNode(Name_Element_before)) + segment.appendChild(FromElement) + toElement = doc.createElement('toElementUID') + toElement.appendChild(doc.createTextNode(Name_Element)) + segment.appendChild(toElement) + + +def initialization(doc, data, name_file): + cpacs = doc.createElement('cpacs') + cpacs.setAttribute("xmlns:xsi", "http://www.w3.org/2001/XMLSchema-instance") + cpacs.setAttribute("xsi_noNamespaceSchemaLocation", "CPACS_21_Schema.xsd") + doc.appendChild(cpacs) + + # --- HEADER --- + header = make(doc, 'header', cpacs) + make(doc, 'name', header, name_file) + make(doc, 'version', header, '0.2') + make(doc, 'cpacsVersion', header, '3.0') + + versionInfos = make(doc, 'versionInfos', header) + versionInfo = make(doc, 'versionInfo', versionInfos, version="version1") + make(doc, 'description', versionInfo, 'Geometry from openVSP') + make(doc, 'creator', versionInfo, 'cpacs2to3') + make(doc, 'version', versionInfo, '0.3') + make(doc, 'cpacsVersion', versionInfo, '3.0') + make(doc, 'timestamp', versionInfo, '2025-10-30T13:40:44') + + # --- VEHICLES --- + vehicles = make(doc, 'vehicles', cpacs) + + keys = list(data.keys()) + for i, k in enumerate(keys): + + if data[k]["Transformation"]["idx_engine"] is not None: + + transformation = data[k]["Transformation"] + + # CASE 1 : Normal nacelle + prev_trans = data[keys[i - 1]]["Transformation"] + if prev_trans["Name_type"] != "Duct": + NameEngine = f"{transformation['Name']}_Engine{transformation['idx_engine']}" + + # ENGINES + engines = make(doc, 'engines', vehicles) + + # ENGINE + engine = make(doc, 'engine', engines, uID=NameEngine) + make(doc, 'name', engine, NameEngine) + + # NACELLE + nacelle = make(doc, 'nacelle', engine, uID=f'{NameEngine}_Nacelle') + fanCowl = make(doc, 'fanCowl', nacelle, uID=f'{NameEngine}_Nacelle_fanCowl') + + sections = make(doc, 'sections', fanCowl) + section = make(doc, 'section', sections, uID=f'{NameEngine}_Nacelle_fanCowl_Sec1') + make( + doc, + 'profileUID', section, + f"{NameEngine}_Nacelle_fanCowl_{data[k]['Airfoil']}" + ) + + Transformation( + doc, + section, + f"{NameEngine}_Nacelle_fanCowl_Sec1", + [ + data[k]["Transformation"]["chord"], + 1, + data[k]["Transformation"]["chord"] / 2, + ], + [0, 0, 0], + [0, 0, 0.5], + ) + + rotationCurve = make( + doc, + 'rotationCurve', + fanCowl, + uID=f'{NameEngine}_Nacelle_fanCowl_RotCurve', + ) + make(doc, 'referenceSectionUID', rotationCurve, + f'{NameEngine}_Nacelle_fanCowl_Sec1') + make(doc, 'curveProfileUID', rotationCurve, f'{NameEngine}_fanCowlRotationCurve') + make(doc, 'startZeta', rotationCurve, '-0.28') + make(doc, 'endZeta', rotationCurve, '-0.25') + make(doc, 'startZetaBlending', rotationCurve, '-0.3') + make(doc, 'endZetaBlending', rotationCurve, '-0.23') + + # CASE 2 : Duct + Duct + elif ( + prev_trans["Name_type"] == "Duct" + and data[k]["Transformation"]["Name_type"] == "Duct" + ): + NameEngine = f"{transformation['Name']}_Engine{transformation['idx_engine']}" + + coreCowl = make(doc, 'coreCowl', nacelle, uID=f'{NameEngine}_coreCowl') + + sections = make(doc, 'sections', coreCowl) + section = make(doc, 'section', sections, uID=f'{NameEngine}_coreCowl_Sec1') + make(doc, 'profileUID', + section, f"{NameEngine}_Nacelle_coreCowl_{data[k]['Airfoil']}") + + Transformation( + doc, + section, + f"{NameEngine}_coreCowl_Sec1", + [ + data[k]["Transformation"]["chord"], + 1, + data[k]["Transformation"]["chord"] / 2, + ], + [0, 0, 0], + [ + 0, + 0, + (0.5 * data[k]["Transformation"]["height"]) + / data[keys[i - 1]]["Transformation"]["height"], + ], + ) + + rotationCurve = make(doc, 'rotationCurve', + coreCowl, uID=f'{NameEngine}_coreCowl_RotCurve') + make(doc, 'referenceSectionUID', rotationCurve, f'{NameEngine}_coreCowl_Sec1') + make(doc, 'curveProfileUID', rotationCurve, f'{NameEngine}_coreCowlRotationCurve') + make(doc, 'startZeta', rotationCurve, '-0.28') + make(doc, 'endZeta', rotationCurve, '-0.25') + make(doc, 'startZetaBlending', rotationCurve, '-0.3') + make(doc, 'endZetaBlending', rotationCurve, '-0.23') + + # CASE 3 : Duct + Duct + Pod + elif data[k]['Transformation']['Name_type'] == 'Pod': + NameEngine = f"{transformation['Name']}_Engine{transformation['idx_engine']}" + + centerCowl = make(doc, 'centerCowl', nacelle, uID=f'{NameEngine}_centerCowl') + make(doc, 'curveUID', centerCowl, f'{NameEngine}_centerCowlRotationCurve') + + offset_val = ( + data[keys[i - 2]]["Transformation"]["X_Trasl"][0] + - data[keys[i - 2]]["Transformation"]["chord"] / 2 + - data[k]["Transformation"]["X_Trasl"][0] + - data[keys[i - 2]]["Transformation"]["chord"] / 2 + ) + + make(doc, 'xOffset', centerCowl, offset_val) + + # --- AIRCRAFT BODY --- + aircraft = make(doc, 'aircraft', vehicles) + model = make(doc, 'model', aircraft, uID='CPACS_test') + make(doc, 'name', model, 'CPACS_test') + + # --- REFERENCE DATA --- + reference = make(doc, 'reference', model) + make(doc, 'area', reference, '1') + make(doc, 'length', reference, '1') + + point = make(doc, 'point', reference, uID='Point') + make(doc, 'x', point, '0') + make(doc, 'y', point, '0') + make(doc, 'z', point, '0') + + return doc, model, vehicles + + +def Transformation(doc, Parent, Name, X_scal, X_Rot, X_Trasl, abs_or_rel=None): + transformation = make(doc, 'transformation', Parent, uID=f'{Name}_Tr') + + txt = ['x', 'y', 'z'] + x_Transf = np.column_stack((X_scal, X_Rot, X_Trasl)) + + labels = ['scaling', 'rotation', 'translation'] + labels_uID = ['Scal', 'Rot', 'Transl'] + + for idx, label in enumerate(labels): + + child = make(doc, label, transformation, uID=f'{Name}_Tr_{labels_uID[idx]}') + + # Only translation may contain refType attribute + if idx == 2 and abs_or_rel is not None: + child.setAttribute('refType', 'absGlobal' if abs_or_rel else 'absLocal') + + # x, y, z coordinates + for i, ax in enumerate(txt): + make(doc, ax, child, x_Transf[i, idx]) + + +def Wing_section(doc, Parent, Name_wing, Section_key, Sections_parameters): + + # Extract numeric index at the end of the section key + match = re.search(r'\d+$', Section_key) + number = int(match.group()) if match else 0 + + # --- Section UID --- + uID_section_root = f'{Name_wing}Sec{number}' + + # Create section element + section_root = make(doc, 'section', Parent, uID=uID_section_root) + + # Section name + make(doc, 'name', section_root, uID_section_root) + + # --- Transformation parameters --- + x_Scal = [ + Sections_parameters[Section_key]["x_scal"], + Sections_parameters[Section_key]["y_scal"], + Sections_parameters[Section_key]["z_scal"], + ] + x_Rot = np.zeros(3) + x_Trasl = [ + Sections_parameters[Section_key]["x_trasl"], + 0, + 0, + ] + + # --- Transformation --- + Transformation( + doc, + section_root, + uID_section_root, + x_Scal, + x_Rot, + x_Trasl, + Sections_parameters["Transformation"]["abs_system"], + ) + # --- Element definition for the section --- + Element( + doc, + section_root, + uID_section_root, + Section_key, + Sections_parameters + ) + + +def Fuse_section(doc, Parent, Name, Section_key, Sections_parameters): + + # Extract numeric index at the end of the section key + match = re.search(r'\d+$', Section_key) + number = int(match.group()) if match else 0 + + # --- Section UID --- + uID_section = f'{Name}Sec{number}' + + # Create section element + section = make(doc, 'section', Parent, uID=uID_section) + + # Section name + make(doc, 'name', section, uID_section) + + # --- Transformation parameters --- + X_Scal = [ + Sections_parameters[Section_key]["x_scal"], + Sections_parameters[Section_key]["y_scal"], + Sections_parameters[Section_key]["z_scal"], + ] + + X_Rot = [ + Sections_parameters[Section_key]["x_rot"], + Sections_parameters[Section_key]["y_rot"], + Sections_parameters[Section_key]["z_rot"], + ] + + X_Trasl = [ + 0, + Sections_parameters[Section_key]["y_trasl"], + Sections_parameters[Section_key]["z_trasl"], + ] + + # --- Transformation --- + Transformation( + doc, + section, + uID_section, + X_Scal, + X_Rot, + X_Trasl, + Sections_parameters["Transformation"]["abs_system"], + ) + + # --- Element definition of the fuselage section --- + Fuse_Element( + doc, + section, + uID_section, + Section_key, + Sections_parameters + ) + + +def Fuse_Element(doc, Parent, Name_section, Section_key, Section_parameters): + + # Element UID based on the section name + Name_element = f'{Name_section}Elem' + + # Create container + Elements = make(doc, 'elements', Parent) + + # Create with uID + Element = make(doc, 'element', Elements, uID=Name_element) + + # Element name tag + make(doc, 'name', Element, Name_element) + + # Refers to fuselage airfoil profile + profile_uid_str = f"{Name_section}_{Section_parameters[Section_key]['Airfoil']}" + make(doc, 'profileUID', Element, profile_uid_str) + + # Identity transformation for the element (external function) + x_Scal = np.ones(3) + x_Rot = np.zeros(3) + x_Trasl = np.zeros(3) + + Transformation( + doc, + Element, + Name_element, + x_Scal, + x_Rot, + x_Trasl, + Section_parameters['Transformation']['abs_system'] + ) + + +def Fuse_positioning(doc, Parent, Name, Section_key, Section_parameters, length_before_perc): + # Extract only the trailing section index (safe when fuselage name contains digits). + match = re.search(r'(.*?)(\d+)$', Name) + if match: + prefix = match.group(1) + number = int(match.group(2)) + else: + prefix = Name + number = 0 + # Create element with uID + positioning = make(doc, 'positioning', Parent, uID=f'{Name}GenPos') + + # Name tag + make(doc, 'name', positioning, f'{Name}GenPos') + + # Fixed sweep and dihedral values for fuselage + make(doc, 'sweepAngle', positioning, '90') + make(doc, 'dihedralAngle', positioning, '0') + + # First fuselage section (no length) + if number == 0: + make(doc, 'length', positioning, '0') + make(doc, 'toSectionUID', positioning, Name) + length_before_perc.append(0) + else: + # Section longitudinal spacing + section_length = float( + Section_parameters[Section_key]['x_loc'] + ) - float(length_before_perc[-1]) + length_before_perc.append(float(Section_parameters[Section_key]['x_loc'])) + + make(doc, 'length', positioning, str(section_length)) + + # Connect from previous to current fuselage section + Name_prev = f"{prefix}{number - 1}" + make(doc, 'fromSectionUID', positioning, Name_prev) + make(doc, 'toSectionUID', positioning, Name) + + +def Fuse_Profile(doc, Parent, Section_parameters, Section_key, uid): + + # Create fuselage profile element with uID + fuselageProfile = make(doc, 'fuselageProfile', Parent, uID=str(uid)) + + # Name and description + make(doc, 'name', fuselageProfile, str(uid)) + make(doc, 'description', fuselageProfile, 'Profile from OpenVSP') + + # Coordinates container + pointList = make(doc, 'pointList', fuselageProfile) + + txt = ['x', 'y', 'z'] + y_values = Section_parameters[Section_key]['Airfoil_coordinates'][0] + z_values = Section_parameters[Section_key]['Airfoil_coordinates'][1] + x_values = np.zeros(np.shape(y_values)) + + # CPACS vector formatting + x_cpacs = np.column_stack((x_values, y_values, z_values)) + + for tag, elem in enumerate(txt): + values_str = ' '.join(str(x_cpacs[:, tag]).strip().split()) \ + .replace('[ ', '').replace( + '[', '').replace(' ]', '').replace(']', '').replace(' ', ';') + child = make(doc, elem, pointList, values_str) + child.setAttribute('mapType', 'vector') + + +def Element(doc, Parent, Name_section, Section_key, Section_parameters): + + Name_element = f'{Name_section}Elem' + + # Create elements container + Elements = make(doc, 'elements', Parent) + + # Create element with uID + Element = make(doc, 'element', Elements, uID=Name_element) + + # Name tag + make(doc, 'name', Element, Name_element) + + # Airfoil UID reference + airfoil_str = f"{Name_section}_{Section_parameters[Section_key]['Airfoil']}" + make(doc, 'airfoilUID', Element, airfoil_str) + + # Identity transformation (external function) + x_Scal = np.ones(3) + x_Rot = np.zeros(3) + x_Trasl = np.zeros(3) + + Transformation( + doc, + Element, + Name_element, + x_Scal, + x_Rot, + x_Trasl, + Section_parameters['Transformation']['abs_system'] + ) + + +def Wing_positioning(doc, Parent, Name, Section_key, Sections_parameters, dih_list): + # Create positioning element with uID + positioning = make(doc, 'positioning', Parent, uID=f'{Name}GenPos') + make(doc, 'name', positioning, f'{Name}GenPos') + + # First wing section + import re + if Name.endswith("Sec0"): + make(doc, 'length', positioning, '0') + make(doc, 'sweepAngle', positioning, '0') + make(doc, 'dihedralAngle', positioning, '0') + make(doc, 'toSectionUID', positioning, Name) + + else: + length_val = Sections_parameters[Section_key]['Span'] / ( + np.cos(np.deg2rad(Sections_parameters[Section_key]['Sweep_angle'])) + ) + make(doc, 'length', positioning, str(length_val)) + + sweep_angle = Sections_parameters[Section_key]['Sweep_angle'] + dihedral_angle = Sections_parameters[Section_key].get('Dihedral_angle', 0.0) + + make(doc, 'sweepAngle', positioning, str(sweep_angle)) + make(doc, 'dihedralAngle', positioning, str(float(dihedral_angle) - dih_list[-1])) + + match = re.search(r'(.*?)(\d+)$', Name) + prev_name = f"{match.group(1)}{int(match.group(2)) - 1}" + + make(doc, 'fromSectionUID', positioning, prev_name) + make(doc, 'toSectionUID', positioning, Name) + # Update dihedral history list + dih_val = Sections_parameters[Section_key].get('Dihedral_angle', 0.0) + if Sections_parameters["Transformation"]["Relative_dih"]: + dih_list.append(dih_val) + else: + dih_list.append(0) + + +def wingAirfoil(doc, Parent, Section_key, Section_parameters, uid): + + # Create wing airfoil element with uID + wingAirfoil = make(doc, 'wingAirfoil', Parent, uID=str(uid)) + make(doc, 'name', wingAirfoil, str(uid)) + + # List of coordinate vectors + pointList = make(doc, 'pointList', wingAirfoil) + + txt = ['x', 'y', 'z'] + x_values = Section_parameters[Section_key]['Airfoil_coordinates'][0] + y_values = np.zeros(np.shape(x_values)) + z_values = Section_parameters[Section_key]['Airfoil_coordinates'][1] + + x_cpacs = np.column_stack((x_values, y_values, z_values)) + + for tag, elem in enumerate(txt): + values_str = ' '.join(str(x_cpacs[:, tag]).strip().split()) \ + .replace('[ ', '').replace('[', '').replace( + ' ]', '').replace(']', '').replace(' ', ';') + child = make(doc, elem, pointList, values_str) + child.setAttribute('mapType', 'vector') + + +def Wing_to_CPACS( + WingData, + doc, + Parent_wing, + Parent_prof, + name_file, + output_dir: Path | None = None, +): + + # ---- keys of the dictionary( number of sections, trasformation of the main wing...) ----# + keys = list(WingData.keys()) + Name_wing = WingData[keys[0]]['Name'] + + # + wings = make(doc, 'wings', Parent_wing) + + # + wing = make(doc, 'wing', wings, uID=Name_wing) + make(doc, 'name', wing, Name_wing) + + if WingData[keys[0]]['Symmetry'] != '0': + wing.setAttribute( + 'symmetry', WingData[keys[0]]['Symmetry']) + + # name - description - transformation # + Transformation(doc, wing, Name_wing, np.ones( + 3).T, np.array(WingData[keys[0]]['X_Rot']), + WingData[keys[0]]['X_Trasl'], WingData[keys[0]]['abs_system']) + + # + Name_section = {} + Name_element = {} + Name_airfoil = {} + sections = make(doc, 'sections', wing) + + for Section_key in keys[1:]: + #
f'{Name_wing}Sec{Section_key[-1]}' + Wing_section( + doc, sections, Name_wing, Section_key, WingData) + + match = re.search(r'\d+$', Section_key) + number = int(match.group()) if match else 0 + + Name_section[Section_key] = { + 'Name': f'{Name_wing}Sec{number}' + } + Name_element[Section_key] = { + 'Name': f"{Name_wing}Sec{number}Elem" + } + Name_airfoil[Section_key] = { + 'Name': f"{Name_wing}Sec{number}_{WingData[Section_key]['Airfoil']}" + } + + # + positionings = make(doc, 'positionings', wing) + dih_list = [] + for Section_key in keys[1:]: + Wing_positioning( + doc, positionings, Name_section[Section_key]['Name'], Section_key, WingData, dih_list) + + # + segments = make(doc, 'segments', wing) + Name_element_before = None + for Section_key in keys[1:]: + + match = re.search(r'\d+$', Section_key) + number = int(match.group()) if match else 0 + + if Section_key != 'Section0' and Name_element_before is not None: + segment( + doc, segments, Name_wing, + Name_element[Section_key]['Name'], Name_element_before, number + ) + Name_element_before = Name_element[Section_key]['Name'] + Name_element_before = Name_element[Section_key]['Name'] + + # + profiles = make(doc, 'profiles', Parent_prof) + wingAirfoils = make(doc, 'wingAirfoils', profiles) + + for Section_key in keys[1:]: + wingAirfoil( + doc, wingAirfoils, Section_key, WingData, Name_airfoil[Section_key]['Name']) + + return Save_CPACS_file(doc, name_file, output_dir) + + +def Fuselage_to_CPACS( + FuseData, + doc, + Parent_Fuse, + Parent_prof, + name_file, + output_dir: Path | None = None, +): + + # ---- keys of the dictionary( number of sections, trasformation of the fuselage...) ----# + keys = list(FuseData.keys()) + Fuse_name = FuseData[keys[0]]['Name'] + # + fuselages = make(doc, 'fuselages', Parent_Fuse) + + # + fuselage = make(doc, 'fuselage', fuselages, uID=Fuse_name) + make(doc, 'name', fuselage, Fuse_name) + + if FuseData[keys[0]]['Symmetry'] != '0': + fuselage.setAttribute( + 'symmetry', FuseData[keys[0]]['Symmetry']) + + # name - description - transformation # + Transformation(doc, fuselage, Fuse_name, np.ones( + 3).T, np.array(FuseData[keys[0]]['X_Rot']), + FuseData[keys[0]]['X_Trasl'], FuseData[keys[0]]['abs_system']) + + # + Name_section = {} + Name_element = {} + Name_airfoil = {} + sections = make(doc, 'sections', fuselage) + + for Section_key in keys[1:]: + #
+ Fuse_section( + doc, sections, Fuse_name, Section_key, FuseData) + + match = re.search(r'\d+$', Section_key) + + number = int(match.group()) if match else 0 + + Name_section[Section_key] = { + 'Name': f"{Fuse_name}Sec{number}" + } + Name_element[Section_key] = { + 'Name': f"{Fuse_name}Sec{number}Elem" + } + Name_airfoil[Section_key] = { + 'Name': f"{Fuse_name}Sec{number}_{FuseData[Section_key]['Airfoil']}" + } + + # + positionings = make(doc, 'positionings', fuselage) + + Store_perc_length = [] + for Section_key in keys[1:]: + Fuse_positioning( + doc, positionings, + Name_section[Section_key]['Name'], Section_key, FuseData, Store_perc_length) + + # + segments = make(doc, 'segments', fuselage) + + Name_element_before = None + for Section_key in keys[1:]: + + match = re.search(r'\d+$', Section_key) + number = int(match.group()) if match else 0 + + if Section_key != "Section0" and Name_element_before is not None: + segment( + doc, + segments, + Fuse_name, + Name_element[Section_key]["Name"], + Name_element_before, + number, + ) + + Name_element_before = Name_element[Section_key]["Name"] + + # + profiles = make(doc, 'profiles', Parent_prof) + fuselageProfiles = make(doc, 'fuselageProfiles', profiles) + + for Section_key in keys[1:]: + Fuse_Profile( + doc, fuselageProfiles, FuseData, Section_key, Name_airfoil[Section_key]['Name']) + return Save_CPACS_file(doc, name_file, output_dir) + + +def Engine_profile(doc, vehicle, data, name, i): + # Create container + profiles = make(doc, 'profiles', vehicle) + + # + curveProfiles = make(doc, 'curveProfiles', profiles) + if len(i) == 1: + curveProfile = make(doc, 'curveProfile', curveProfiles, uID=f"{name}_fanCowlRotationCurve") + pointList = make(doc, 'pointList', curveProfile) + + make( + doc, + 'x', + pointList, + '0;0.05;0.1;0.15;0.2;0.3;0.4;0.5;0.6;0.7;0.8;0.85;0.9;0.95;1' + ) + + make( + doc, + "y", + pointList, + ( + "0;-0.006;-0.0085;-0.0102;-0.0113;" + "-0.0125;-0.0132;-0.0135;-0.013;" + "-0.0118;-0.0098;-0.0082;" + "-0.0062;-0.0035;0" + ), + ) + + elif len(i) == 2: + curveProfile = make(doc, 'curveProfile', + curveProfiles, uID=f"{name}_coreCowlRotationCurve") + pointList = make(doc, 'pointList', curveProfile) + + make( + doc, + 'x', + pointList, + '0.27108433734939763;0.2720883534136546' + ) + + make( + doc, + 'y', + pointList, + '-0.0252;-0.0252' + ) + + elif len(i) == 3: + curveProfile = make(doc, 'curveProfile', + curveProfiles, uID=f"{name}_centerCowlRotationCurve") + pointList = make(doc, 'pointList', curveProfile) + cp = " ".join(str(data["Transformation"]["curveProfile"][0]).strip().split()) + cp = cp.replace("[ ", "").replace( + "[", "").replace(" ]", "").replace("]", "").replace(" ", ";") + make(doc, "x", pointList, cp) + cp = " ".join(str(data["Transformation"]["curveProfile"][1]).strip().split()) + cp = cp.replace("[ ", "").replace( + "[", "").replace(" ]", "").replace("]", "").replace(" ", ";") + make(doc, "y", pointList, cp) + + # + nacelleProfiles = make(doc, 'nacelleProfiles', profiles) + if len(i) == 1: + nacelleProfile = make(doc, 'nacelleProfile', + nacelleProfiles, uID=f"{name}_Nacelle_fanCowl_{data['Airfoil']}") + make(doc, 'name', nacelleProfile, data['Airfoil']) + pointList = make(doc, 'pointList', nacelleProfile) + coord_x = " ".join(str(data["Airfoil_coordinates"][0]).strip().split()) + coord_x = coord_x.replace("[ ", "").replace("[", "").replace( + " ]", "").replace("]", "").replace(" ", ";") + make(doc, "x", pointList, coord_x) + coord_y = " ".join(str(data["Airfoil_coordinates"][1]).strip().split()) + coord_y = coord_y.replace("[ ", "").replace("[", "").replace( + " ]", "").replace("]", "").replace(" ", ";") + make(doc, "y", pointList, coord_y) + elif len(i) == 2: + nacelleProfile = make(doc, 'nacelleProfile', + nacelleProfiles, uID=f"{name}_Nacelle_coreCowl_{data['Airfoil']}") + make(doc, 'name', nacelleProfile, data['Airfoil']) + pointList = make(doc, 'pointList', nacelleProfile) + coord_x = " ".join(str(data["Airfoil_coordinates"][0]).strip().split()) + coord_x = coord_x.replace("[ ", "").replace( + "[", "").replace(" ]", "").replace("]", "").replace(" ", ";") + make(doc, "x", pointList, coord_x) + coord_y = " ".join(str(data["Airfoil_coordinates"][1]).strip().split()) + coord_y = coord_y.replace("[ ", "").replace("[", "").replace( + " ]", "").replace("]", "").replace(" ", ";") + make(doc, "y", pointList, coord_y) + + +def Engine_to_CPACS( + EngineData, + doc, + Parent_engine, + Parent_prof, + i, + name_file, + output_dir: Path | None = None, +): + # ---- keys of the dictionary( number of sections, trasformation of the fuselage...) ----# + NameEngine = f"""{ + EngineData['Transformation']['Name'] + }_Engine{EngineData['Transformation']['idx_engine']}""" + + if len(i) == 1: + # + engines = make(doc, 'engines', Parent_engine) + + # + engine = make(doc, 'engine', engines, uID=f"{NameEngine}_Positioning") + if EngineData['Transformation']['Symmetry'] != '0': + engine.setAttribute('symmetry', EngineData['Transformation']['Symmetry']) + + # + make(doc, 'parentUID', engine, EngineData['Transformation']['Child_to_Parent']) + + # + make(doc, 'engineUID', engine, NameEngine) + + # Transformation for engine + Transformation(doc, engine, NameEngine, + [1, EngineData['Transformation']['width'], + EngineData['Transformation']['height']], + EngineData['Transformation']['X_Rot'], + EngineData['Transformation']['X_Trasl']) + + # + Engine_profile(doc, Parent_prof, EngineData, NameEngine, i) + + return Save_CPACS_file(doc, name_file, output_dir) + + +def merge_elements(document, parent_tag, target_tag): + + # unify more element inside the arent node , only if exist. + parent_nodes = document.getElementsByTagName(parent_tag) + if not parent_nodes: + return document + parent = parent_nodes[0] + elements = list(parent.getElementsByTagName(target_tag)) + if not elements: + return document + first = elements[0] + if len(elements) > 1: + for extra in elements[1:]: + for c in list(extra.childNodes): + extra.removeChild(c) + first.appendChild(c) + extra.parentNode.removeChild(extra) + return document + + +def Save_CPACS_file(Document, name_file, output_dir: Path | None = None) -> Path: + + merge_elements(Document, 'model', 'wings') + merge_elements(Document, 'model', 'fuselages') + merge_elements(Document, 'vehicles', 'profiles') + merge_elements(Document, 'profiles', 'fuselageProfiles') + merge_elements(Document, 'profiles', 'wingAirfoils') + merge_elements(Document, 'profiles', 'curveProfiles') + merge_elements(Document, 'profiles', 'nacelleProfiles') + + xml_str = Document.toprettyxml(indent=" ") + + # The file will be saved either in the provided output directory or, by default, + # in the same folder as this module. + if output_dir is None: + module_dir = Path(__file__).parent + output_dir = module_dir.parent + else: + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + + output_path = output_dir / f'{name_file}.xml' + + with open(output_path, 'w') as xml_file: + xml_file.write(xml_str) + + return output_path + + +# Classes +class _ETElement: + def __init__(self, element, parent: "_ETElement | None" = None): + self._element = element + self._parent = parent + + @property + def parentNode(self) -> "_ETElement | None": + return self._parent + + @property + def childNodes(self) -> list["_ETElement"]: + return [_ETElement(child, parent=self) for child in list(self._element)] + + def setAttribute(self, key: str, value: str) -> None: + self._element.set(key, value) + + def appendChild(self, node): + if isinstance(node, str): + text_value = str(node) + if self._element.text is None: + self._element.text = text_value + else: + self._element.text += text_value + return node + if isinstance(node, _ETElement): + self._element.append(node._element) + node._parent = self + return node + raise TypeError(f"Unsupported node type for appendChild: {type(node)!r}") + + def removeChild(self, node: "_ETElement"): + if not isinstance(node, _ETElement): + raise TypeError(f"Unsupported node type for removeChild: {type(node)!r}") + self._element.remove(node._element) + node._parent = None + return node + + def getElementsByTagName(self, tag: str) -> list["_ETElement"]: + matches: list[_ETElement] = [] + + def walk(element, parent_wrapper: _ETElement) -> None: + for child in list(element): + child_wrapper = _ETElement(child, parent=parent_wrapper) + if child.tag == tag: + matches.append(child_wrapper) + walk(child, child_wrapper) + + walk(self._element, self) + return matches + + +class _ETDocument: + def __init__(self): + self._root: _ETElement | None = None + + def createElement(self, name: str) -> _ETElement: + return _ETElement(ET.Element(name)) + + def createTextNode(self, value: str) -> str: + return str(value) + + def appendChild(self, node: _ETElement) -> _ETElement: + if self._root is not None: + raise ValueError("Document already has a root element") + self._root = node + return node + + def getElementsByTagName(self, tag: str) -> list[_ETElement]: + if self._root is None: + return [] + if self._root._element.tag == tag: + return [self._root, *self._root.getElementsByTagName(tag)] + return self._root.getElementsByTagName(tag) + + def toprettyxml(self, indent: str = " ") -> str: + if self._root is None: + return "" + ET.indent(self._root._element, space=indent) + return ET.tostring(self._root._element, encoding="unicode", xml_declaration=True) + + +class Export_CPACS: + def __init__(self, Data, name_file, output_dir: Path | None = None): + self.Data = Data + self.name_file = name_file + self.output_dir = output_dir + self.output_path: Path | None = None + + def run(self) -> Path: + # Create the document + Doc = _ETDocument() + + # find the keys + keys = list(self.Data.keys()) + + # dummy variable for the engine + dummy_idx_engine = [] + + # CPACS's initialization + Doc, model, vehicles = initialization(Doc, self.Data, self.name_file) + + # Loop to connect the components inside the CPACS + for item in keys: + log.info(f"Component {item} {self.Data[f'{item}']['Transformation']['Name_type']} ") + if self.Data[f'{item}']['Transformation']['Name_type'] == 'Wing': + self.output_path = Wing_to_CPACS( + self.Data[f'{item}'], + Doc, + model, + vehicles, + self.name_file, + self.output_dir, + ) + if self.Data[f'{item}']['Transformation']['Name_type'] == 'Fuselage': + self.output_path = Fuselage_to_CPACS( + self.Data[f'{item}'], + Doc, + model, + vehicles, + self.name_file, + self.output_dir, + ) + if ( + self.Data[f'{item}']['Transformation']['Name_type'] == 'Pod' + and self.Data[f'{item}']['Transformation']['idx_engine'] is None + ): + self.output_path = Fuselage_to_CPACS( + self.Data[f'{item}'], + Doc, + model, + vehicles, + self.name_file, + self.output_dir, + ) + if ( + self.Data[f'{item}']['Transformation']['Name_type'] == 'Duct' + or self.Data[f'{item}']['Transformation']['idx_engine'] is not None + ): + if len(dummy_idx_engine) < 3: + dummy_idx_engine.append( + self.Data[f'{item}']['Transformation']['idx_engine'] + ) + self.output_path = Engine_to_CPACS( + self.Data[f'{item}'], + Doc, + model, + vehicles, + dummy_idx_engine, + self.name_file, + self.output_dir, + ) + + if self.output_path is None: + raise RuntimeError("Failed to export CPACS file from OpenVSP geometry.") + + return self.output_path diff --git a/src/ceasiompy/vsp2cpacs/func/wing.py b/src/ceasiompy/vsp2cpacs/func/wing.py index d004531b4..1e41ba651 100644 --- a/src/ceasiompy/vsp2cpacs/func/wing.py +++ b/src/ceasiompy/vsp2cpacs/func/wing.py @@ -857,10 +857,15 @@ def get_coord_from_file(xsec_id, n): Reads upper and lower surface points from the selected x-section, merges them into a closed profile, and returns the coordinates along with the section name and chord scaling. +<<<<<<< HEAD:src/ceasiompy/VSP2CPACS/func/wing.py + """ + geom_parm = get_params_by_name(xsec_id,['Chord','ThickChord','BaseThickChord']) +======= """ geom_parm = get_params_by_name(xsec_id,['Chord','ThickChord','BaseThickChord']) # Name +>>>>>>> origin/main:src/ceasiompy/vsp2cpacs/func/wing.py Name = vsp.GetXSecCurveAlias(xsec_id).replace(' ','_') # import the points @@ -869,9 +874,14 @@ def get_coord_from_file(xsec_id, n): upper_coords = np.array([[p.x(), p.y(), p.z()] for p in upper_pts]) lower_coords = np.array([[p.x(), p.y(), p.z()] for p in lower_pts]) +<<<<<<< HEAD:src/ceasiompy/VSP2CPACS/func/wing.py + x_u, y_u = upper_coords[:, 0], upper_coords[:, 1] + x_l, y_l = lower_coords[:, 0], lower_coords[:, 1] +======= x_u, y_u = upper_coords[:,0], upper_coords[:,1] x_l, y_l = lower_coords[:,0], lower_coords[:,1] +>>>>>>> origin/main:src/ceasiompy/vsp2cpacs/func/wing.py x = np.concatenate((x_l[::-1], x_u), axis=0) y = np.concatenate((y_l[::-1], y_u), axis=0) * (geom_parm['ThickChord']/geom_parm['BaseThickChord']) @@ -1181,10 +1191,15 @@ def Get_coordinates_profile(idx, *args, **kwargs): the trailing edge. """ - + breakpoint() Airfoil_name_type = vsp.GetXSecShape(idx) +<<<<<<< HEAD:src/ceasiompy/VSP2CPACS/func/wing.py + print(Airfoil_name_type) + func = profile_mapping[Airfoil_name_type] +======= func = profile_mapping()[Airfoil_name_type] print(f'Working with {idx}') +>>>>>>> origin/main:src/ceasiompy/vsp2cpacs/func/wing.py try: return func(idx, *args, **kwargs) except TypeError: @@ -1196,8 +1211,12 @@ def Get_coordinates_profile(idx, *args, **kwargs): def get_profile_section(Component,xsec_id, idx, Twist_val, Twist_loc, Rel, Twist_list): # Tess_W control how many points you need to define the shape of the profile +<<<<<<< HEAD:src/ceasiompy/VSP2CPACS/func/wing.py + Tess_W = int(vsp.GetParmVal(Component, "Tess_W", "Shape")) +======= Tess_W = int(vsp.GetParmVal(Component,'Tess_W','Shape')) +>>>>>>> origin/main:src/ceasiompy/vsp2cpacs/func/wing.py # get profile x, y, Airfoil_name,Scaling,shift = Get_coordinates_profile(xsec_id,Tess_W) @@ -1234,6 +1253,12 @@ def get_profile_section(Component,xsec_id, idx, Twist_val, Twist_loc, Rel, Twist Coord_rot = np.dot(RotationMatrix, Coord_shift) + Origin_shift x = Coord_rot[0, :] y = Coord_rot[1, :] + + # LE duplicates from twist part. + zero_idx = np.where(np.isclose(x, 0.0, atol=1e-12))[0] + if len(zero_idx) > 1: + x = np.delete(x, zero_idx[1:]) + y = np.delete(y, zero_idx[1:]) zero_idx = np.where(np.isclose(x, 0.0, atol=1e-12))[0]