From cc9f5eec323d749ff34af8ceb355054b196ecc01 Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Tue, 7 Mar 2023 15:19:35 +0100
Subject: [PATCH] Add compressibility and unsteady sdpm
---
sdpm/_src/sdpmw.i | 40 +-
sdpm/models/lann.geo | 1730 ++++++++++++++++++++++++++++++++++++++
sdpm/models/oneraM6.geo | 10 +-
sdpm/src/sdpm.h | 6 +
sdpm/src/sdpmBody.cpp | 108 ++-
sdpm/src/sdpmBody.h | 54 +-
sdpm/src/sdpmBuilder.cpp | 39 +-
sdpm/src/sdpmBuilder.h | 1 +
sdpm/src/sdpmElement.cpp | 16 +-
sdpm/src/sdpmElement.h | 27 +-
sdpm/src/sdpmElement.inl | 59 ++
sdpm/src/sdpmLine2.cpp | 16 +-
sdpm/src/sdpmLine2.h | 2 +-
sdpm/src/sdpmMotion.cpp | 43 +
sdpm/src/sdpmMotion.h | 48 ++
sdpm/src/sdpmProblem.cpp | 60 +-
sdpm/src/sdpmProblem.h | 41 +-
sdpm/src/sdpmQuad4.cpp | 24 +-
sdpm/src/sdpmQuad4.h | 2 +-
sdpm/src/sdpmSolver.cpp | 340 ++++++--
sdpm/src/sdpmSolver.h | 39 +-
sdpm/src/sdpmWake.cpp | 176 ++--
sdpm/src/sdpmWake.h | 15 +-
sdpm/src/sdpm_extras.h | 30 +-
sdpm/tests/lann.py | 86 ++
sdpm/tests/m6.py | 9 +-
sdpm/tests/n12.py | 8 +-
sdpm/tests/n12tail.py | 30 +-
28 files changed, 2715 insertions(+), 344 deletions(-)
create mode 100644 sdpm/models/lann.geo
create mode 100644 sdpm/src/sdpmElement.inl
create mode 100644 sdpm/src/sdpmMotion.cpp
create mode 100644 sdpm/src/sdpmMotion.h
create mode 100644 sdpm/tests/lann.py
diff --git a/sdpm/_src/sdpmw.i b/sdpm/_src/sdpmw.i
index 8a384a2..b7342f7 100644
--- a/sdpm/_src/sdpmw.i
+++ b/sdpm/_src/sdpmw.i
@@ -51,6 +51,7 @@ threads="1"
%include "std_string.i"
%include "std_vector.i"
%include "std_map.i"
+%include "std_complex.i"
// EXCEPTION
// from: http://swig.10945.n7.nabble.com/Trapping-Swig-DirectorException-td6013.html
@@ -86,6 +87,15 @@ threads="1"
}
}
+// SDPM templates
+%template(std_vector_double) std::vector<double>;
+%template(std_vector_complex) std::vector<std::complex<double>>;
+%template(std_vector_vector_complex) std::vector<std::vector<std::complex<double>>>;
+%template(std_vector_Nodep) std::vector<sdpm::Node *>;
+%template(std_vector_Elementp) std::vector<sdpm::Element *>;
+%template(std_vector_Tagp) std::vector<sdpm::Tag *>;
+%template(std_map_string_Tagp) std::map<std::string, sdpm::Tag *>;
+
// SDPM
%include "sdpm.h"
%include "sdpmObject.h"
@@ -110,12 +120,32 @@ threads="1"
%include "sdpmCppBuf2Py.h"
%warnfilter(+401);
-// SDPM template
-%template(std_vector_Nodep) std::vector<sdpm::Node *>;
-%template(std_vector_Elementp) std::vector<sdpm::Element *>;
-%template(std_vector_Tagp) std::vector<sdpm::Tag *>;
-%template(std_map_string_Tagp) std::map<std::string, sdpm::Tag *>;
+// Interface to SDPM types
+class sdpmVector3d {};
+
+%extend sdpmVector3d {
+ double getitem(int i)
+ {
+ return (*self)(i);
+ }
+ void setitem(int i, double a)
+ {
+ (*self)(i) = a;
+ }
+ std::string __str__()
+ {
+ std::ostringstream out; out << *self; return out.str();
+ }
+
+ %pythoncode {
+def __getitem__(self, idx):
+ return self.getitem(idx)
+def __setitem__(self, idx, v):
+ return self.setitem(idx, v)
+ }
+}
+// Interface to SDPM classes
namespace sdpm {
// So that print(obj) use << *obj
diff --git a/sdpm/models/lann.geo b/sdpm/models/lann.geo
new file mode 100644
index 0000000..e61924e
--- /dev/null
+++ b/sdpm/models/lann.geo
@@ -0,0 +1,1730 @@
+/******************************************/
+/* Gmsh geometry for lann */
+/* Generated by geoGen.py */
+/* Adrien Crovato */
+/* ULiege, 2018-2019 */
+/******************************************/
+
+// --- Wing geometry ---
+// Number of spanwise stations: 8
+// Spanwise stations normalized coordinate: 0.000000 0.200000 0.325000 0.475000 0.650000 0.825000 0.950000 1.000000
+// Chord lengths: 0.360600 0.317650 0.290710 0.258060 0.220290 0.182350 0.155340 0.144450
+// Half-wing area: 0.252692
+// Half-wing span: 1.000000
+// Full-wing aspect ratio: 7.914766
+
+// --- Options ---
+DefineConstant[ nC = { 50, Name "Number of panel along chord" }
+ bC = { 0.05, Name "Progression along chord" }
+ rS = { 1, Name "Refinement factor on span" } ];
+
+
+// --- Wing points ---
+// -- Airfoil 0
+Point(1) = {0.360600,0.000000,-0.008200};
+Point(2) = {0.360513,0.000000,-0.008196};
+Point(3) = {0.360247,0.000000,-0.008176};
+Point(4) = {0.359802,0.000000,-0.008142};
+Point(5) = {0.359180,0.000000,-0.008095};
+Point(6) = {0.358380,0.000000,-0.008036};
+Point(7) = {0.357405,0.000000,-0.007968};
+Point(8) = {0.356257,0.000000,-0.007892};
+Point(9) = {0.354935,0.000000,-0.007809};
+Point(10) = {0.353441,0.000000,-0.007721};
+Point(11) = {0.351775,0.000000,-0.007597};
+Point(12) = {0.349941,0.000000,-0.007478};
+Point(13) = {0.347939,0.000000,-0.007369};
+Point(14) = {0.345771,0.000000,-0.007274};
+Point(15) = {0.343440,0.000000,-0.007201};
+Point(16) = {0.340948,0.000000,-0.007154};
+Point(17) = {0.338298,0.000000,-0.007140};
+Point(18) = {0.335491,0.000000,-0.007162};
+Point(19) = {0.332532,0.000000,-0.007224};
+Point(20) = {0.329423,0.000000,-0.007330};
+Point(21) = {0.326166,0.000000,-0.007484};
+Point(22) = {0.322765,0.000000,-0.007690};
+Point(23) = {0.319223,0.000000,-0.007946};
+Point(24) = {0.315545,0.000000,-0.008253};
+Point(25) = {0.311733,0.000000,-0.008608};
+Point(26) = {0.307791,0.000000,-0.009011};
+Point(27) = {0.303724,0.000000,-0.009461};
+Point(28) = {0.299534,0.000000,-0.009954};
+Point(29) = {0.295227,0.000000,-0.010489};
+Point(30) = {0.290807,0.000000,-0.011062};
+Point(31) = {0.286278,0.000000,-0.011670};
+Point(32) = {0.281644,0.000000,-0.012310};
+Point(33) = {0.276909,0.000000,-0.012979};
+Point(34) = {0.272080,0.000000,-0.013672};
+Point(35) = {0.267160,0.000000,-0.014386};
+Point(36) = {0.262154,0.000000,-0.015117};
+Point(37) = {0.257068,0.000000,-0.015861};
+Point(38) = {0.251906,0.000000,-0.016612};
+Point(39) = {0.246673,0.000000,-0.017364};
+Point(40) = {0.241374,0.000000,-0.018112};
+Point(41) = {0.236016,0.000000,-0.018848};
+Point(42) = {0.230602,0.000000,-0.019567};
+Point(43) = {0.225139,0.000000,-0.020261};
+Point(44) = {0.219631,0.000000,-0.020924};
+Point(45) = {0.214085,0.000000,-0.021548};
+Point(46) = {0.208505,0.000000,-0.022125};
+Point(47) = {0.202898,0.000000,-0.022649};
+Point(48) = {0.197268,0.000000,-0.023116};
+Point(49) = {0.191621,0.000000,-0.023526};
+Point(50) = {0.185963,0.000000,-0.023874};
+Point(51) = {0.180300,0.000000,-0.024160};
+Point(52) = {0.174637,0.000000,-0.024381};
+Point(53) = {0.168979,0.000000,-0.024540};
+Point(54) = {0.163332,0.000000,-0.024641};
+Point(55) = {0.157702,0.000000,-0.024685};
+Point(56) = {0.152095,0.000000,-0.024675};
+Point(57) = {0.146515,0.000000,-0.024615};
+Point(58) = {0.140969,0.000000,-0.024504};
+Point(59) = {0.135461,0.000000,-0.024342};
+Point(60) = {0.129998,0.000000,-0.024128};
+Point(61) = {0.124584,0.000000,-0.023861};
+Point(62) = {0.119226,0.000000,-0.023540};
+Point(63) = {0.113927,0.000000,-0.023167};
+Point(64) = {0.108694,0.000000,-0.022741};
+Point(65) = {0.103532,0.000000,-0.022264};
+Point(66) = {0.098446,0.000000,-0.021736};
+Point(67) = {0.093440,0.000000,-0.021157};
+Point(68) = {0.088520,0.000000,-0.020531};
+Point(69) = {0.083691,0.000000,-0.019860};
+Point(70) = {0.078956,0.000000,-0.019146};
+Point(71) = {0.074322,0.000000,-0.018394};
+Point(72) = {0.069793,0.000000,-0.017606};
+Point(73) = {0.065372,0.000000,-0.016787};
+Point(74) = {0.061066,0.000000,-0.015943};
+Point(75) = {0.056876,0.000000,-0.015080};
+Point(76) = {0.052809,0.000000,-0.014202};
+Point(77) = {0.048867,0.000000,-0.013314};
+Point(78) = {0.045055,0.000000,-0.012423};
+Point(79) = {0.041377,0.000000,-0.011530};
+Point(80) = {0.037835,0.000000,-0.010640};
+Point(81) = {0.034434,0.000000,-0.009757};
+Point(82) = {0.031177,0.000000,-0.008884};
+Point(83) = {0.028068,0.000000,-0.008022};
+Point(84) = {0.025109,0.000000,-0.007169};
+Point(85) = {0.022302,0.000000,-0.006325};
+Point(86) = {0.019652,0.000000,-0.005491};
+Point(87) = {0.017160,0.000000,-0.004665};
+Point(88) = {0.014829,0.000000,-0.003851};
+Point(89) = {0.012661,0.000000,-0.003050};
+Point(90) = {0.010659,0.000000,-0.002265};
+Point(91) = {0.008825,0.000000,-0.001498};
+Point(92) = {0.007159,0.000000,-0.000749};
+Point(93) = {0.005665,0.000000,-0.000001};
+Point(94) = {0.004342,0.000000,0.000764};
+Point(95) = {0.003194,0.000000,0.001565};
+Point(96) = {0.002220,0.000000,0.002420};
+Point(97) = {0.001422,0.000000,0.003344};
+Point(98) = {0.000801,0.000000,0.004327};
+Point(99) = {0.000356,0.000000,0.005353};
+Point(100) = {0.000089,0.000000,0.006406};
+Point(101) = {0.000000,0.000000,0.007472};
+Point(102) = {0.000089,0.000000,0.008535};
+Point(103) = {0.000356,0.000000,0.009582};
+Point(104) = {0.000801,0.000000,0.010599};
+Point(105) = {0.001422,0.000000,0.011576};
+Point(106) = {0.002220,0.000000,0.012497};
+Point(107) = {0.003194,0.000000,0.013354};
+Point(108) = {0.004342,0.000000,0.014150};
+Point(109) = {0.005665,0.000000,0.014892};
+Point(110) = {0.007159,0.000000,0.015585};
+Point(111) = {0.008825,0.000000,0.016236};
+Point(112) = {0.010659,0.000000,0.016849};
+Point(113) = {0.012661,0.000000,0.017424};
+Point(114) = {0.014829,0.000000,0.017956};
+Point(115) = {0.017160,0.000000,0.018442};
+Point(116) = {0.019652,0.000000,0.018878};
+Point(117) = {0.022302,0.000000,0.019263};
+Point(118) = {0.025109,0.000000,0.019597};
+Point(119) = {0.028068,0.000000,0.019887};
+Point(120) = {0.031177,0.000000,0.020135};
+Point(121) = {0.034434,0.000000,0.020345};
+Point(122) = {0.037835,0.000000,0.020521};
+Point(123) = {0.041377,0.000000,0.020667};
+Point(124) = {0.045055,0.000000,0.020785};
+Point(125) = {0.048867,0.000000,0.020877};
+Point(126) = {0.052809,0.000000,0.020948};
+Point(127) = {0.056876,0.000000,0.020999};
+Point(128) = {0.061066,0.000000,0.021032};
+Point(129) = {0.065372,0.000000,0.021045};
+Point(130) = {0.069793,0.000000,0.021038};
+Point(131) = {0.074322,0.000000,0.021013};
+Point(132) = {0.078956,0.000000,0.020968};
+Point(133) = {0.083691,0.000000,0.020903};
+Point(134) = {0.088520,0.000000,0.020816};
+Point(135) = {0.093440,0.000000,0.020708};
+Point(136) = {0.098446,0.000000,0.020578};
+Point(137) = {0.103532,0.000000,0.020425};
+Point(138) = {0.108694,0.000000,0.020249};
+Point(139) = {0.113927,0.000000,0.020053};
+Point(140) = {0.119226,0.000000,0.019836};
+Point(141) = {0.124584,0.000000,0.019600};
+Point(142) = {0.129998,0.000000,0.019347};
+Point(143) = {0.135461,0.000000,0.019074};
+Point(144) = {0.140969,0.000000,0.018781};
+Point(145) = {0.146515,0.000000,0.018466};
+Point(146) = {0.152095,0.000000,0.018129};
+Point(147) = {0.157702,0.000000,0.017769};
+Point(148) = {0.163332,0.000000,0.017384};
+Point(149) = {0.168979,0.000000,0.016974};
+Point(150) = {0.174637,0.000000,0.016537};
+Point(151) = {0.180300,0.000000,0.016074};
+Point(152) = {0.185963,0.000000,0.015582};
+Point(153) = {0.191621,0.000000,0.015065};
+Point(154) = {0.197268,0.000000,0.014524};
+Point(155) = {0.202898,0.000000,0.013962};
+Point(156) = {0.208505,0.000000,0.013380};
+Point(157) = {0.214085,0.000000,0.012782};
+Point(158) = {0.219631,0.000000,0.012168};
+Point(159) = {0.225139,0.000000,0.011539};
+Point(160) = {0.230602,0.000000,0.010896};
+Point(161) = {0.236016,0.000000,0.010240};
+Point(162) = {0.241374,0.000000,0.009571};
+Point(163) = {0.246673,0.000000,0.008893};
+Point(164) = {0.251906,0.000000,0.008207};
+Point(165) = {0.257068,0.000000,0.007517};
+Point(166) = {0.262154,0.000000,0.006824};
+Point(167) = {0.267160,0.000000,0.006131};
+Point(168) = {0.272080,0.000000,0.005440};
+Point(169) = {0.276909,0.000000,0.004752};
+Point(170) = {0.281644,0.000000,0.004068};
+Point(171) = {0.286278,0.000000,0.003389};
+Point(172) = {0.290807,0.000000,0.002717};
+Point(173) = {0.295227,0.000000,0.002055};
+Point(174) = {0.299534,0.000000,0.001404};
+Point(175) = {0.303724,0.000000,0.000767};
+Point(176) = {0.307791,0.000000,0.000145};
+Point(177) = {0.311733,0.000000,-0.000459};
+Point(178) = {0.315545,0.000000,-0.001045};
+Point(179) = {0.319223,0.000000,-0.001611};
+Point(180) = {0.322765,0.000000,-0.002156};
+Point(181) = {0.326166,0.000000,-0.002681};
+Point(182) = {0.329423,0.000000,-0.003184};
+Point(183) = {0.332532,0.000000,-0.003664};
+Point(184) = {0.335491,0.000000,-0.004121};
+Point(185) = {0.338298,0.000000,-0.004554};
+Point(186) = {0.340948,0.000000,-0.004963};
+Point(187) = {0.343440,0.000000,-0.005348};
+Point(188) = {0.345771,0.000000,-0.005706};
+Point(189) = {0.347939,0.000000,-0.006040};
+Point(190) = {0.349941,0.000000,-0.006347};
+Point(191) = {0.351775,0.000000,-0.006627};
+Point(192) = {0.353441,0.000000,-0.006881};
+Point(193) = {0.354935,0.000000,-0.007144};
+Point(194) = {0.356257,0.000000,-0.007382};
+Point(195) = {0.357405,0.000000,-0.007593};
+Point(196) = {0.358380,0.000000,-0.007776};
+Point(197) = {0.359180,0.000000,-0.007928};
+Point(198) = {0.359802,0.000000,-0.008049};
+Point(199) = {0.360247,0.000000,-0.008135};
+Point(200) = {0.360513,0.000000,-0.008186};
+// -- Airfoil 1
+Point(501) = {0.421732,0.200000,-0.005789};
+Point(502) = {0.421656,0.200000,-0.005794};
+Point(503) = {0.421422,0.200000,-0.005780};
+Point(504) = {0.421030,0.200000,-0.005751};
+Point(505) = {0.420481,0.200000,-0.005708};
+Point(506) = {0.419777,0.200000,-0.005654};
+Point(507) = {0.418918,0.200000,-0.005590};
+Point(508) = {0.417906,0.200000,-0.005519};
+Point(509) = {0.416742,0.200000,-0.005442};
+Point(510) = {0.415426,0.200000,-0.005362};
+Point(511) = {0.413959,0.200000,-0.005274};
+Point(512) = {0.412343,0.200000,-0.005196};
+Point(513) = {0.410579,0.200000,-0.005130};
+Point(514) = {0.408670,0.200000,-0.005077};
+Point(515) = {0.406616,0.200000,-0.005039};
+Point(516) = {0.404421,0.200000,-0.005016};
+Point(517) = {0.402087,0.200000,-0.005011};
+Point(518) = {0.399614,0.200000,-0.005031};
+Point(519) = {0.397007,0.200000,-0.005081};
+Point(520) = {0.394268,0.200000,-0.005168};
+Point(521) = {0.391399,0.200000,-0.005299};
+Point(522) = {0.388404,0.200000,-0.005479};
+Point(523) = {0.385284,0.200000,-0.005707};
+Point(524) = {0.382044,0.200000,-0.005980};
+Point(525) = {0.378686,0.200000,-0.006296};
+Point(526) = {0.375213,0.200000,-0.006654};
+Point(527) = {0.371630,0.200000,-0.007050};
+Point(528) = {0.367940,0.200000,-0.007484};
+Point(529) = {0.364146,0.200000,-0.007956};
+Point(530) = {0.360252,0.200000,-0.008466};
+Point(531) = {0.356262,0.200000,-0.009012};
+Point(532) = {0.352180,0.200000,-0.009595};
+Point(533) = {0.348010,0.200000,-0.010208};
+Point(534) = {0.343756,0.200000,-0.010848};
+Point(535) = {0.339422,0.200000,-0.011508};
+Point(536) = {0.335012,0.200000,-0.012183};
+Point(537) = {0.330532,0.200000,-0.012868};
+Point(538) = {0.325984,0.200000,-0.013558};
+Point(539) = {0.321375,0.200000,-0.014249};
+Point(540) = {0.316707,0.200000,-0.014936};
+Point(541) = {0.311987,0.200000,-0.015614};
+Point(542) = {0.307218,0.200000,-0.016281};
+Point(543) = {0.302406,0.200000,-0.016927};
+Point(544) = {0.297554,0.200000,-0.017545};
+Point(545) = {0.292668,0.200000,-0.018130};
+Point(546) = {0.287753,0.200000,-0.018671};
+Point(547) = {0.282814,0.200000,-0.019165};
+Point(548) = {0.277854,0.200000,-0.019609};
+Point(549) = {0.272880,0.200000,-0.020000};
+Point(550) = {0.267896,0.200000,-0.020339};
+Point(551) = {0.262907,0.200000,-0.020623};
+Point(552) = {0.257919,0.200000,-0.020852};
+Point(553) = {0.252935,0.200000,-0.021028};
+Point(554) = {0.247961,0.200000,-0.021151};
+Point(555) = {0.243001,0.200000,-0.021224};
+Point(556) = {0.238062,0.200000,-0.021248};
+Point(557) = {0.233146,0.200000,-0.021224};
+Point(558) = {0.228261,0.200000,-0.021154};
+Point(559) = {0.223409,0.200000,-0.021038};
+Point(560) = {0.218597,0.200000,-0.020878};
+Point(561) = {0.213828,0.200000,-0.020673};
+Point(562) = {0.209107,0.200000,-0.020426};
+Point(563) = {0.204440,0.200000,-0.020136};
+Point(564) = {0.199830,0.200000,-0.019802};
+Point(565) = {0.195283,0.200000,-0.019425};
+Point(566) = {0.190802,0.200000,-0.019004};
+Point(567) = {0.186393,0.200000,-0.018540};
+Point(568) = {0.182059,0.200000,-0.018035};
+Point(569) = {0.177805,0.200000,-0.017491};
+Point(570) = {0.173634,0.200000,-0.016910};
+Point(571) = {0.169552,0.200000,-0.016296};
+Point(572) = {0.165563,0.200000,-0.015651};
+Point(573) = {0.161668,0.200000,-0.014979};
+Point(574) = {0.157875,0.200000,-0.014287};
+Point(575) = {0.154184,0.200000,-0.013579};
+Point(576) = {0.150601,0.200000,-0.012859};
+Point(577) = {0.147129,0.200000,-0.012133};
+Point(578) = {0.143771,0.200000,-0.011402};
+Point(579) = {0.140531,0.200000,-0.010670};
+Point(580) = {0.137411,0.200000,-0.009938};
+Point(581) = {0.134415,0.200000,-0.009208};
+Point(582) = {0.131546,0.200000,-0.008482};
+Point(583) = {0.128807,0.200000,-0.007761};
+Point(584) = {0.126200,0.200000,-0.007045};
+Point(585) = {0.123728,0.200000,-0.006334};
+Point(586) = {0.121393,0.200000,-0.005629};
+Point(587) = {0.119198,0.200000,-0.004931};
+Point(588) = {0.117145,0.200000,-0.004240};
+Point(589) = {0.115236,0.200000,-0.003557};
+Point(590) = {0.113472,0.200000,-0.002883};
+Point(591) = {0.111856,0.200000,-0.002219};
+Point(592) = {0.110389,0.200000,-0.001564};
+Point(593) = {0.109072,0.200000,-0.000908};
+Point(594) = {0.107907,0.200000,-0.000236};
+Point(595) = {0.106896,0.200000,0.000462};
+Point(596) = {0.106038,0.200000,0.001200};
+Point(597) = {0.105335,0.200000,0.001986};
+Point(598) = {0.104788,0.200000,0.002812};
+Point(599) = {0.104396,0.200000,0.003664};
+Point(600) = {0.104161,0.200000,0.004530};
+Point(601) = {0.104082,0.200000,0.005397};
+Point(602) = {0.104161,0.200000,0.006254};
+Point(603) = {0.104396,0.200000,0.007099};
+Point(604) = {0.104788,0.200000,0.007933};
+Point(605) = {0.105335,0.200000,0.008754};
+Point(606) = {0.106038,0.200000,0.009563};
+Point(607) = {0.106896,0.200000,0.010359};
+Point(608) = {0.107907,0.200000,0.011133};
+Point(609) = {0.109072,0.200000,0.011877};
+Point(610) = {0.110389,0.200000,0.012581};
+Point(611) = {0.111856,0.200000,0.013236};
+Point(612) = {0.113472,0.200000,0.013835};
+Point(613) = {0.115236,0.200000,0.014380};
+Point(614) = {0.117145,0.200000,0.014875};
+Point(615) = {0.119198,0.200000,0.015322};
+Point(616) = {0.121393,0.200000,0.015728};
+Point(617) = {0.123728,0.200000,0.016095};
+Point(618) = {0.126200,0.200000,0.016425};
+Point(619) = {0.128807,0.200000,0.016722};
+Point(620) = {0.131546,0.200000,0.016987};
+Point(621) = {0.134415,0.200000,0.017222};
+Point(622) = {0.137411,0.200000,0.017431};
+Point(623) = {0.140531,0.200000,0.017614};
+Point(624) = {0.143771,0.200000,0.017773};
+Point(625) = {0.147129,0.200000,0.017908};
+Point(626) = {0.150601,0.200000,0.018023};
+Point(627) = {0.154184,0.200000,0.018116};
+Point(628) = {0.157875,0.200000,0.018191};
+Point(629) = {0.161668,0.200000,0.018247};
+Point(630) = {0.165563,0.200000,0.018286};
+Point(631) = {0.169552,0.200000,0.018309};
+Point(632) = {0.173634,0.200000,0.018317};
+Point(633) = {0.177805,0.200000,0.018309};
+Point(634) = {0.182059,0.200000,0.018283};
+Point(635) = {0.186393,0.200000,0.018238};
+Point(636) = {0.190802,0.200000,0.018173};
+Point(637) = {0.195283,0.200000,0.018086};
+Point(638) = {0.199830,0.200000,0.017978};
+Point(639) = {0.204440,0.200000,0.017850};
+Point(640) = {0.209107,0.200000,0.017703};
+Point(641) = {0.213828,0.200000,0.017539};
+Point(642) = {0.218597,0.200000,0.017357};
+Point(643) = {0.223409,0.200000,0.017158};
+Point(644) = {0.228261,0.200000,0.016940};
+Point(645) = {0.233146,0.200000,0.016704};
+Point(646) = {0.238062,0.200000,0.016448};
+Point(647) = {0.243001,0.200000,0.016172};
+Point(648) = {0.247961,0.200000,0.015875};
+Point(649) = {0.252935,0.200000,0.015558};
+Point(650) = {0.257919,0.200000,0.015219};
+Point(651) = {0.262907,0.200000,0.014858};
+Point(652) = {0.267896,0.200000,0.014475};
+Point(653) = {0.272880,0.200000,0.014072};
+Point(654) = {0.277854,0.200000,0.013647};
+Point(655) = {0.282814,0.200000,0.013203};
+Point(656) = {0.287753,0.200000,0.012740};
+Point(657) = {0.292668,0.200000,0.012258};
+Point(658) = {0.297554,0.200000,0.011760};
+Point(659) = {0.302406,0.200000,0.011244};
+Point(660) = {0.307218,0.200000,0.010713};
+Point(661) = {0.311987,0.200000,0.010168};
+Point(662) = {0.316707,0.200000,0.009609};
+Point(663) = {0.321375,0.200000,0.009038};
+Point(664) = {0.325984,0.200000,0.008459};
+Point(665) = {0.330532,0.200000,0.007875};
+Point(666) = {0.335012,0.200000,0.007287};
+Point(667) = {0.339422,0.200000,0.006698};
+Point(668) = {0.343756,0.200000,0.006109};
+Point(669) = {0.348010,0.200000,0.005520};
+Point(670) = {0.352180,0.200000,0.004934};
+Point(671) = {0.356262,0.200000,0.004349};
+Point(672) = {0.360252,0.200000,0.003767};
+Point(673) = {0.364146,0.200000,0.003190};
+Point(674) = {0.367940,0.200000,0.002621};
+Point(675) = {0.371630,0.200000,0.002060};
+Point(676) = {0.375213,0.200000,0.001510};
+Point(677) = {0.378686,0.200000,0.000974};
+Point(678) = {0.382044,0.200000,0.000451};
+Point(679) = {0.385284,0.200000,-0.000057};
+Point(680) = {0.388404,0.200000,-0.000549};
+Point(681) = {0.391399,0.200000,-0.001024};
+Point(682) = {0.394268,0.200000,-0.001482};
+Point(683) = {0.397007,0.200000,-0.001921};
+Point(684) = {0.399614,0.200000,-0.002337};
+Point(685) = {0.402087,0.200000,-0.002730};
+Point(686) = {0.404421,0.200000,-0.003097};
+Point(687) = {0.406616,0.200000,-0.003437};
+Point(688) = {0.408670,0.200000,-0.003750};
+Point(689) = {0.410579,0.200000,-0.004037};
+Point(690) = {0.412343,0.200000,-0.004299};
+Point(691) = {0.413959,0.200000,-0.004538};
+Point(692) = {0.415426,0.200000,-0.004753};
+Point(693) = {0.416742,0.200000,-0.004960};
+Point(694) = {0.417906,0.200000,-0.005149};
+Point(695) = {0.418918,0.200000,-0.005318};
+Point(696) = {0.419777,0.200000,-0.005465};
+Point(697) = {0.420481,0.200000,-0.005587};
+Point(698) = {0.421030,0.200000,-0.005683};
+Point(699) = {0.421422,0.200000,-0.005750};
+Point(700) = {0.421656,0.200000,-0.005786};
+// -- Airfoil 2
+Point(1001) = {0.459844,0.325000,-0.004603};
+Point(1002) = {0.459774,0.325000,-0.004604};
+Point(1003) = {0.459560,0.325000,-0.004590};
+Point(1004) = {0.459201,0.325000,-0.004562};
+Point(1005) = {0.458699,0.325000,-0.004524};
+Point(1006) = {0.458054,0.325000,-0.004477};
+Point(1007) = {0.457268,0.325000,-0.004423};
+Point(1008) = {0.456342,0.325000,-0.004363};
+Point(1009) = {0.455276,0.325000,-0.004300};
+Point(1010) = {0.454072,0.325000,-0.004237};
+Point(1011) = {0.452730,0.325000,-0.004154};
+Point(1012) = {0.451251,0.325000,-0.004080};
+Point(1013) = {0.449636,0.325000,-0.004017};
+Point(1014) = {0.447889,0.325000,-0.003968};
+Point(1015) = {0.446010,0.325000,-0.003933};
+Point(1016) = {0.444001,0.325000,-0.003914};
+Point(1017) = {0.441864,0.325000,-0.003915};
+Point(1018) = {0.439602,0.325000,-0.003938};
+Point(1019) = {0.437216,0.325000,-0.003988};
+Point(1020) = {0.434709,0.325000,-0.004070};
+Point(1021) = {0.432083,0.325000,-0.004188};
+Point(1022) = {0.429342,0.325000,-0.004346};
+Point(1023) = {0.426487,0.325000,-0.004544};
+Point(1024) = {0.423521,0.325000,-0.004782};
+Point(1025) = {0.420448,0.325000,-0.005061};
+Point(1026) = {0.417270,0.325000,-0.005381};
+Point(1027) = {0.413991,0.325000,-0.005741};
+Point(1028) = {0.410614,0.325000,-0.006141};
+Point(1029) = {0.407141,0.325000,-0.006578};
+Point(1030) = {0.403578,0.325000,-0.007052};
+Point(1031) = {0.399926,0.325000,-0.007560};
+Point(1032) = {0.396191,0.325000,-0.008102};
+Point(1033) = {0.392374,0.325000,-0.008672};
+Point(1034) = {0.388480,0.325000,-0.009267};
+Point(1035) = {0.384514,0.325000,-0.009882};
+Point(1036) = {0.380479,0.325000,-0.010511};
+Point(1037) = {0.376378,0.325000,-0.011152};
+Point(1038) = {0.372216,0.325000,-0.011798};
+Point(1039) = {0.367997,0.325000,-0.012447};
+Point(1040) = {0.363726,0.325000,-0.013094};
+Point(1041) = {0.359406,0.325000,-0.013735};
+Point(1042) = {0.355041,0.325000,-0.014364};
+Point(1043) = {0.350637,0.325000,-0.014977};
+Point(1044) = {0.346197,0.325000,-0.015565};
+Point(1045) = {0.341726,0.325000,-0.016122};
+Point(1046) = {0.337227,0.325000,-0.016641};
+Point(1047) = {0.332707,0.325000,-0.017117};
+Point(1048) = {0.328168,0.325000,-0.017547};
+Point(1049) = {0.323616,0.325000,-0.017930};
+Point(1050) = {0.319054,0.325000,-0.018265};
+Point(1051) = {0.314489,0.325000,-0.018550};
+Point(1052) = {0.309923,0.325000,-0.018784};
+Point(1053) = {0.305362,0.325000,-0.018969};
+Point(1054) = {0.300810,0.325000,-0.019105};
+Point(1055) = {0.296271,0.325000,-0.019195};
+Point(1056) = {0.291750,0.325000,-0.019239};
+Point(1057) = {0.287252,0.325000,-0.019239};
+Point(1058) = {0.282781,0.325000,-0.019196};
+Point(1059) = {0.278340,0.325000,-0.019110};
+Point(1060) = {0.273936,0.325000,-0.018984};
+Point(1061) = {0.269572,0.325000,-0.018818};
+Point(1062) = {0.265252,0.325000,-0.018612};
+Point(1063) = {0.260980,0.325000,-0.018369};
+Point(1064) = {0.256761,0.325000,-0.018086};
+Point(1065) = {0.252600,0.325000,-0.017765};
+Point(1066) = {0.248499,0.325000,-0.017406};
+Point(1067) = {0.244463,0.325000,-0.017008};
+Point(1068) = {0.240497,0.325000,-0.016575};
+Point(1069) = {0.236604,0.325000,-0.016108};
+Point(1070) = {0.232787,0.325000,-0.015609};
+Point(1071) = {0.229051,0.325000,-0.015081};
+Point(1072) = {0.225400,0.325000,-0.014526};
+Point(1073) = {0.221836,0.325000,-0.013949};
+Point(1074) = {0.218364,0.325000,-0.013352};
+Point(1075) = {0.214987,0.325000,-0.012740};
+Point(1076) = {0.211707,0.325000,-0.012117};
+Point(1077) = {0.208530,0.325000,-0.011486};
+Point(1078) = {0.205457,0.325000,-0.010850};
+Point(1079) = {0.202491,0.325000,-0.010211};
+Point(1080) = {0.199636,0.325000,-0.009572};
+Point(1081) = {0.196894,0.325000,-0.008934};
+Point(1082) = {0.194269,0.325000,-0.008300};
+Point(1083) = {0.191762,0.325000,-0.007671};
+Point(1084) = {0.189376,0.325000,-0.007048};
+Point(1085) = {0.187113,0.325000,-0.006434};
+Point(1086) = {0.184977,0.325000,-0.005828};
+Point(1087) = {0.182968,0.325000,-0.005232};
+Point(1088) = {0.181089,0.325000,-0.004640};
+Point(1089) = {0.179341,0.325000,-0.004048};
+Point(1090) = {0.177727,0.325000,-0.003449};
+Point(1091) = {0.176248,0.325000,-0.002838};
+Point(1092) = {0.174905,0.325000,-0.002212};
+Point(1093) = {0.173701,0.325000,-0.001579};
+Point(1094) = {0.172635,0.325000,-0.000948};
+Point(1095) = {0.171709,0.325000,-0.000330};
+Point(1096) = {0.170923,0.325000,0.000266};
+Point(1097) = {0.170280,0.325000,0.000836};
+Point(1098) = {0.169779,0.325000,0.001404};
+Point(1099) = {0.169421,0.325000,0.002000};
+Point(1100) = {0.169206,0.325000,0.002653};
+Point(1101) = {0.169134,0.325000,0.003395};
+Point(1102) = {0.169206,0.325000,0.004244};
+Point(1103) = {0.169421,0.325000,0.005168};
+Point(1104) = {0.169779,0.325000,0.006124};
+Point(1105) = {0.170280,0.325000,0.007069};
+Point(1106) = {0.170923,0.325000,0.007959};
+Point(1107) = {0.171709,0.325000,0.008761};
+Point(1108) = {0.172635,0.325000,0.009481};
+Point(1109) = {0.173701,0.325000,0.010135};
+Point(1110) = {0.174905,0.325000,0.010738};
+Point(1111) = {0.176248,0.325000,0.011306};
+Point(1112) = {0.177727,0.325000,0.011853};
+Point(1113) = {0.179341,0.325000,0.012374};
+Point(1114) = {0.181089,0.325000,0.012866};
+Point(1115) = {0.182968,0.325000,0.013323};
+Point(1116) = {0.184977,0.325000,0.013738};
+Point(1117) = {0.187113,0.325000,0.014110};
+Point(1118) = {0.189376,0.325000,0.014440};
+Point(1119) = {0.191762,0.325000,0.014735};
+Point(1120) = {0.194269,0.325000,0.015000};
+Point(1121) = {0.196894,0.325000,0.015239};
+Point(1122) = {0.199636,0.325000,0.015458};
+Point(1123) = {0.202491,0.325000,0.015658};
+Point(1124) = {0.205457,0.325000,0.015840};
+Point(1125) = {0.208530,0.325000,0.016003};
+Point(1126) = {0.211707,0.325000,0.016150};
+Point(1127) = {0.214987,0.325000,0.016279};
+Point(1128) = {0.218364,0.325000,0.016391};
+Point(1129) = {0.221836,0.325000,0.016487};
+Point(1130) = {0.225400,0.325000,0.016566};
+Point(1131) = {0.229051,0.325000,0.016627};
+Point(1132) = {0.232787,0.325000,0.016670};
+Point(1133) = {0.236604,0.325000,0.016695};
+Point(1134) = {0.240497,0.325000,0.016701};
+Point(1135) = {0.244463,0.325000,0.016688};
+Point(1136) = {0.248499,0.325000,0.016655};
+Point(1137) = {0.252600,0.325000,0.016601};
+Point(1138) = {0.256761,0.325000,0.016528};
+Point(1139) = {0.260980,0.325000,0.016437};
+Point(1140) = {0.265252,0.325000,0.016329};
+Point(1141) = {0.269572,0.325000,0.016207};
+Point(1142) = {0.273936,0.325000,0.016070};
+Point(1143) = {0.278340,0.325000,0.015918};
+Point(1144) = {0.282781,0.325000,0.015751};
+Point(1145) = {0.287252,0.325000,0.015567};
+Point(1146) = {0.291750,0.325000,0.015365};
+Point(1147) = {0.296271,0.325000,0.015144};
+Point(1148) = {0.300810,0.325000,0.014905};
+Point(1149) = {0.305362,0.325000,0.014645};
+Point(1150) = {0.309923,0.325000,0.014366};
+Point(1151) = {0.314489,0.325000,0.014066};
+Point(1152) = {0.319054,0.325000,0.013746};
+Point(1153) = {0.323616,0.325000,0.013405};
+Point(1154) = {0.328168,0.325000,0.013044};
+Point(1155) = {0.332707,0.325000,0.012666};
+Point(1156) = {0.337227,0.325000,0.012269};
+Point(1157) = {0.341726,0.325000,0.011856};
+Point(1158) = {0.346197,0.325000,0.011426};
+Point(1159) = {0.350637,0.325000,0.010981};
+Point(1160) = {0.355041,0.325000,0.010521};
+Point(1161) = {0.359406,0.325000,0.010047};
+Point(1162) = {0.363726,0.325000,0.009560};
+Point(1163) = {0.367997,0.325000,0.009061};
+Point(1164) = {0.372216,0.325000,0.008554};
+Point(1165) = {0.376378,0.325000,0.008038};
+Point(1166) = {0.380479,0.325000,0.007518};
+Point(1167) = {0.384514,0.325000,0.006993};
+Point(1168) = {0.388480,0.325000,0.006465};
+Point(1169) = {0.392374,0.325000,0.005936};
+Point(1170) = {0.396191,0.325000,0.005405};
+Point(1171) = {0.399926,0.325000,0.004873};
+Point(1172) = {0.403578,0.325000,0.004341};
+Point(1173) = {0.407141,0.325000,0.003812};
+Point(1174) = {0.410614,0.325000,0.003288};
+Point(1175) = {0.413991,0.325000,0.002771};
+Point(1176) = {0.417270,0.325000,0.002264};
+Point(1177) = {0.420448,0.325000,0.001769};
+Point(1178) = {0.423521,0.325000,0.001286};
+Point(1179) = {0.426487,0.325000,0.000817};
+Point(1180) = {0.429342,0.325000,0.000360};
+Point(1181) = {0.432083,0.325000,-0.000084};
+Point(1182) = {0.434709,0.325000,-0.000514};
+Point(1183) = {0.437216,0.325000,-0.000928};
+Point(1184) = {0.439602,0.325000,-0.001322};
+Point(1185) = {0.441864,0.325000,-0.001693};
+Point(1186) = {0.444001,0.325000,-0.002038};
+Point(1187) = {0.446010,0.325000,-0.002354};
+Point(1188) = {0.447889,0.325000,-0.002644};
+Point(1189) = {0.449636,0.325000,-0.002908};
+Point(1190) = {0.451251,0.325000,-0.003149};
+Point(1191) = {0.452730,0.325000,-0.003368};
+Point(1192) = {0.454072,0.325000,-0.003567};
+Point(1193) = {0.455276,0.325000,-0.003771};
+Point(1194) = {0.456342,0.325000,-0.003957};
+Point(1195) = {0.457268,0.325000,-0.004124};
+Point(1196) = {0.458054,0.325000,-0.004269};
+Point(1197) = {0.458699,0.325000,-0.004391};
+Point(1198) = {0.459201,0.325000,-0.004488};
+Point(1199) = {0.459560,0.325000,-0.004557};
+Point(1200) = {0.459774,0.325000,-0.004596};
+// -- Airfoil 3
+Point(1501) = {0.505256,0.475000,-0.003178};
+Point(1502) = {0.505194,0.475000,-0.003180};
+Point(1503) = {0.505003,0.475000,-0.003169};
+Point(1504) = {0.504685,0.475000,-0.003147};
+Point(1505) = {0.504239,0.475000,-0.003114};
+Point(1506) = {0.503667,0.475000,-0.003073};
+Point(1507) = {0.502969,0.475000,-0.003025};
+Point(1508) = {0.502147,0.475000,-0.002972};
+Point(1509) = {0.501201,0.475000,-0.002916};
+Point(1510) = {0.500132,0.475000,-0.002858};
+Point(1511) = {0.498940,0.475000,-0.002777};
+Point(1512) = {0.497628,0.475000,-0.002701};
+Point(1513) = {0.496195,0.475000,-0.002634};
+Point(1514) = {0.494643,0.475000,-0.002578};
+Point(1515) = {0.492975,0.475000,-0.002534};
+Point(1516) = {0.491192,0.475000,-0.002504};
+Point(1517) = {0.489295,0.475000,-0.002492};
+Point(1518) = {0.487287,0.475000,-0.002501};
+Point(1519) = {0.485169,0.475000,-0.002535};
+Point(1520) = {0.482944,0.475000,-0.002599};
+Point(1521) = {0.480613,0.475000,-0.002697};
+Point(1522) = {0.478179,0.475000,-0.002833};
+Point(1523) = {0.475645,0.475000,-0.003007};
+Point(1524) = {0.473012,0.475000,-0.003219};
+Point(1525) = {0.470284,0.475000,-0.003468};
+Point(1526) = {0.467463,0.475000,-0.003753};
+Point(1527) = {0.464553,0.475000,-0.004076};
+Point(1528) = {0.461554,0.475000,-0.004435};
+Point(1529) = {0.458472,0.475000,-0.004830};
+Point(1530) = {0.455309,0.475000,-0.005261};
+Point(1531) = {0.452068,0.475000,-0.005728};
+Point(1532) = {0.448751,0.475000,-0.006231};
+Point(1533) = {0.445363,0.475000,-0.006763};
+Point(1534) = {0.441907,0.475000,-0.007318};
+Point(1535) = {0.438386,0.475000,-0.007888};
+Point(1536) = {0.434804,0.475000,-0.008469};
+Point(1537) = {0.431164,0.475000,-0.009054};
+Point(1538) = {0.427470,0.475000,-0.009641};
+Point(1539) = {0.423725,0.475000,-0.010229};
+Point(1540) = {0.419933,0.475000,-0.010817};
+Point(1541) = {0.416098,0.475000,-0.011404};
+Point(1542) = {0.412224,0.475000,-0.011987};
+Point(1543) = {0.408314,0.475000,-0.012561};
+Point(1544) = {0.404372,0.475000,-0.013118};
+Point(1545) = {0.400403,0.475000,-0.013650};
+Point(1546) = {0.396410,0.475000,-0.014149};
+Point(1547) = {0.392397,0.475000,-0.014608};
+Point(1548) = {0.388368,0.475000,-0.015026};
+Point(1549) = {0.384327,0.475000,-0.015402};
+Point(1550) = {0.380278,0.475000,-0.015734};
+Point(1551) = {0.376226,0.475000,-0.016021};
+Point(1552) = {0.372173,0.475000,-0.016262};
+Point(1553) = {0.368124,0.475000,-0.016459};
+Point(1554) = {0.364083,0.475000,-0.016612};
+Point(1555) = {0.360054,0.475000,-0.016723};
+Point(1556) = {0.356041,0.475000,-0.016793};
+Point(1557) = {0.352048,0.475000,-0.016823};
+Point(1558) = {0.348079,0.475000,-0.016815};
+Point(1559) = {0.344137,0.475000,-0.016770};
+Point(1560) = {0.340228,0.475000,-0.016688};
+Point(1561) = {0.336353,0.475000,-0.016573};
+Point(1562) = {0.332518,0.475000,-0.016424};
+Point(1563) = {0.328727,0.475000,-0.016242};
+Point(1564) = {0.324982,0.475000,-0.016027};
+Point(1565) = {0.321287,0.475000,-0.015780};
+Point(1566) = {0.317647,0.475000,-0.015501};
+Point(1567) = {0.314065,0.475000,-0.015189};
+Point(1568) = {0.310544,0.475000,-0.014847};
+Point(1569) = {0.307088,0.475000,-0.014477};
+Point(1570) = {0.303700,0.475000,-0.014079};
+Point(1571) = {0.300384,0.475000,-0.013658};
+Point(1572) = {0.297142,0.475000,-0.013213};
+Point(1573) = {0.293979,0.475000,-0.012749};
+Point(1574) = {0.290897,0.475000,-0.012268};
+Point(1575) = {0.287899,0.475000,-0.011774};
+Point(1576) = {0.284988,0.475000,-0.011269};
+Point(1577) = {0.282167,0.475000,-0.010757};
+Point(1578) = {0.279439,0.475000,-0.010241};
+Point(1579) = {0.276806,0.475000,-0.009722};
+Point(1580) = {0.274272,0.475000,-0.009203};
+Point(1581) = {0.271838,0.475000,-0.008688};
+Point(1582) = {0.269507,0.475000,-0.008177};
+Point(1583) = {0.267282,0.475000,-0.007671};
+Point(1584) = {0.265164,0.475000,-0.007166};
+Point(1585) = {0.263156,0.475000,-0.006662};
+Point(1586) = {0.261259,0.475000,-0.006157};
+Point(1587) = {0.259476,0.475000,-0.005649};
+Point(1588) = {0.257808,0.475000,-0.005141};
+Point(1589) = {0.256257,0.475000,-0.004635};
+Point(1590) = {0.254824,0.475000,-0.004133};
+Point(1591) = {0.253511,0.475000,-0.003639};
+Point(1592) = {0.252319,0.475000,-0.003151};
+Point(1593) = {0.251249,0.475000,-0.002659};
+Point(1594) = {0.250303,0.475000,-0.002145};
+Point(1595) = {0.249481,0.475000,-0.001597};
+Point(1596) = {0.248784,0.475000,-0.000998};
+Point(1597) = {0.248213,0.475000,-0.000338};
+Point(1598) = {0.247768,0.475000,0.000368};
+Point(1599) = {0.247450,0.475000,0.001103};
+Point(1600) = {0.247259,0.475000,0.001847};
+Point(1601) = {0.247196,0.475000,0.002583};
+Point(1602) = {0.247259,0.475000,0.003295};
+Point(1603) = {0.247450,0.475000,0.003985};
+Point(1604) = {0.247768,0.475000,0.004656};
+Point(1605) = {0.248213,0.475000,0.005312};
+Point(1606) = {0.248784,0.475000,0.005959};
+Point(1607) = {0.249481,0.475000,0.006598};
+Point(1608) = {0.250303,0.475000,0.007225};
+Point(1609) = {0.251249,0.475000,0.007834};
+Point(1610) = {0.252319,0.475000,0.008418};
+Point(1611) = {0.253511,0.475000,0.008972};
+Point(1612) = {0.254824,0.475000,0.009490};
+Point(1613) = {0.256257,0.475000,0.009972};
+Point(1614) = {0.257808,0.475000,0.010421};
+Point(1615) = {0.259476,0.475000,0.010839};
+Point(1616) = {0.261259,0.475000,0.011226};
+Point(1617) = {0.263156,0.475000,0.011586};
+Point(1618) = {0.265164,0.475000,0.011920};
+Point(1619) = {0.267282,0.475000,0.012229};
+Point(1620) = {0.269507,0.475000,0.012515};
+Point(1621) = {0.271838,0.475000,0.012779};
+Point(1622) = {0.274272,0.475000,0.013023};
+Point(1623) = {0.276806,0.475000,0.013248};
+Point(1624) = {0.279439,0.475000,0.013455};
+Point(1625) = {0.282167,0.475000,0.013645};
+Point(1626) = {0.284988,0.475000,0.013820};
+Point(1627) = {0.287899,0.475000,0.013979};
+Point(1628) = {0.290897,0.475000,0.014122};
+Point(1629) = {0.293979,0.475000,0.014249};
+Point(1630) = {0.297142,0.475000,0.014358};
+Point(1631) = {0.300384,0.475000,0.014450};
+Point(1632) = {0.303700,0.475000,0.014523};
+Point(1633) = {0.307088,0.475000,0.014580};
+Point(1634) = {0.310544,0.475000,0.014621};
+Point(1635) = {0.314065,0.475000,0.014649};
+Point(1636) = {0.317647,0.475000,0.014665};
+Point(1637) = {0.321287,0.475000,0.014670};
+Point(1638) = {0.324982,0.475000,0.014663};
+Point(1639) = {0.328727,0.475000,0.014644};
+Point(1640) = {0.332518,0.475000,0.014609};
+Point(1641) = {0.336353,0.475000,0.014558};
+Point(1642) = {0.340228,0.475000,0.014489};
+Point(1643) = {0.344137,0.475000,0.014402};
+Point(1644) = {0.348079,0.475000,0.014297};
+Point(1645) = {0.352048,0.475000,0.014174};
+Point(1646) = {0.356041,0.475000,0.014035};
+Point(1647) = {0.360054,0.475000,0.013878};
+Point(1648) = {0.364083,0.475000,0.013704};
+Point(1649) = {0.368124,0.475000,0.013512};
+Point(1650) = {0.372173,0.475000,0.013302};
+Point(1651) = {0.376226,0.475000,0.013073};
+Point(1652) = {0.380278,0.475000,0.012826};
+Point(1653) = {0.384327,0.475000,0.012561};
+Point(1654) = {0.388368,0.475000,0.012280};
+Point(1655) = {0.392397,0.475000,0.011983};
+Point(1656) = {0.396410,0.475000,0.011673};
+Point(1657) = {0.400403,0.475000,0.011348};
+Point(1658) = {0.404372,0.475000,0.011011};
+Point(1659) = {0.408314,0.475000,0.010659};
+Point(1660) = {0.412224,0.475000,0.010293};
+Point(1661) = {0.416098,0.475000,0.009911};
+Point(1662) = {0.419933,0.475000,0.009514};
+Point(1663) = {0.423725,0.475000,0.009104};
+Point(1664) = {0.427470,0.475000,0.008682};
+Point(1665) = {0.431164,0.475000,0.008250};
+Point(1666) = {0.434804,0.475000,0.007811};
+Point(1667) = {0.438386,0.475000,0.007366};
+Point(1668) = {0.441907,0.475000,0.006916};
+Point(1669) = {0.445363,0.475000,0.006460};
+Point(1670) = {0.448751,0.475000,0.005999};
+Point(1671) = {0.452068,0.475000,0.005532};
+Point(1672) = {0.455309,0.475000,0.005060};
+Point(1673) = {0.458472,0.475000,0.004586};
+Point(1674) = {0.461554,0.475000,0.004114};
+Point(1675) = {0.464553,0.475000,0.003646};
+Point(1676) = {0.467463,0.475000,0.003187};
+Point(1677) = {0.470284,0.475000,0.002738};
+Point(1678) = {0.473012,0.475000,0.002302};
+Point(1679) = {0.475645,0.475000,0.001878};
+Point(1680) = {0.478179,0.475000,0.001466};
+Point(1681) = {0.480613,0.475000,0.001065};
+Point(1682) = {0.482944,0.475000,0.000677};
+Point(1683) = {0.485169,0.475000,0.000303};
+Point(1684) = {0.487287,0.475000,-0.000054};
+Point(1685) = {0.489295,0.475000,-0.000393};
+Point(1686) = {0.491192,0.475000,-0.000712};
+Point(1687) = {0.492975,0.475000,-0.001008};
+Point(1688) = {0.494643,0.475000,-0.001283};
+Point(1689) = {0.496195,0.475000,-0.001536};
+Point(1690) = {0.497628,0.475000,-0.001769};
+Point(1691) = {0.498940,0.475000,-0.001982};
+Point(1692) = {0.500132,0.475000,-0.002175};
+Point(1693) = {0.501201,0.475000,-0.002375};
+Point(1694) = {0.502147,0.475000,-0.002558};
+Point(1695) = {0.502969,0.475000,-0.002720};
+Point(1696) = {0.503667,0.475000,-0.002861};
+Point(1697) = {0.504239,0.475000,-0.002978};
+Point(1698) = {0.504685,0.475000,-0.003070};
+Point(1699) = {0.505003,0.475000,-0.003136};
+Point(1700) = {0.505194,0.475000,-0.003172};
+// -- Airfoil 4
+Point(2001) = {0.558558,0.650000,-0.001008};
+Point(2002) = {0.558505,0.650000,-0.001006};
+Point(2003) = {0.558342,0.650000,-0.000993};
+Point(2004) = {0.558070,0.650000,-0.000969};
+Point(2005) = {0.557690,0.650000,-0.000937};
+Point(2006) = {0.557202,0.650000,-0.000896};
+Point(2007) = {0.556606,0.650000,-0.000850};
+Point(2008) = {0.555904,0.650000,-0.000798};
+Point(2009) = {0.555097,0.650000,-0.000742};
+Point(2010) = {0.554184,0.650000,-0.000684};
+Point(2011) = {0.553167,0.650000,-0.000602};
+Point(2012) = {0.552046,0.650000,-0.000526};
+Point(2013) = {0.550823,0.650000,-0.000459};
+Point(2014) = {0.549499,0.650000,-0.000401};
+Point(2015) = {0.548075,0.650000,-0.000356};
+Point(2016) = {0.546553,0.650000,-0.000327};
+Point(2017) = {0.544933,0.650000,-0.000315};
+Point(2018) = {0.543219,0.650000,-0.000324};
+Point(2019) = {0.541411,0.650000,-0.000356};
+Point(2020) = {0.539511,0.650000,-0.000415};
+Point(2021) = {0.537522,0.650000,-0.000505};
+Point(2022) = {0.535444,0.650000,-0.000627};
+Point(2023) = {0.533281,0.650000,-0.000781};
+Point(2024) = {0.531034,0.650000,-0.000967};
+Point(2025) = {0.528705,0.650000,-0.001185};
+Point(2026) = {0.526297,0.650000,-0.001434};
+Point(2027) = {0.523812,0.650000,-0.001713};
+Point(2028) = {0.521253,0.650000,-0.002023};
+Point(2029) = {0.518622,0.650000,-0.002365};
+Point(2030) = {0.515921,0.650000,-0.002740};
+Point(2031) = {0.513154,0.650000,-0.003149};
+Point(2032) = {0.510323,0.650000,-0.003592};
+Point(2033) = {0.507431,0.650000,-0.004064};
+Point(2034) = {0.504481,0.650000,-0.004558};
+Point(2035) = {0.501475,0.650000,-0.005067};
+Point(2036) = {0.498417,0.650000,-0.005587};
+Point(2037) = {0.495310,0.650000,-0.006111};
+Point(2038) = {0.492157,0.650000,-0.006639};
+Point(2039) = {0.488960,0.650000,-0.007170};
+Point(2040) = {0.485723,0.650000,-0.007706};
+Point(2041) = {0.482449,0.650000,-0.008246};
+Point(2042) = {0.479142,0.650000,-0.008789};
+Point(2043) = {0.475805,0.650000,-0.009328};
+Point(2044) = {0.472440,0.650000,-0.009856};
+Point(2045) = {0.469052,0.650000,-0.010365};
+Point(2046) = {0.465643,0.650000,-0.010846};
+Point(2047) = {0.462218,0.650000,-0.011293};
+Point(2048) = {0.458778,0.650000,-0.011704};
+Point(2049) = {0.455329,0.650000,-0.012077};
+Point(2050) = {0.451872,0.650000,-0.012412};
+Point(2051) = {0.448413,0.650000,-0.012708};
+Point(2052) = {0.444953,0.650000,-0.012963};
+Point(2053) = {0.441497,0.650000,-0.013178};
+Point(2054) = {0.438047,0.650000,-0.013356};
+Point(2055) = {0.434608,0.650000,-0.013495};
+Point(2056) = {0.431182,0.650000,-0.013599};
+Point(2057) = {0.427773,0.650000,-0.013668};
+Point(2058) = {0.424385,0.650000,-0.013703};
+Point(2059) = {0.421021,0.650000,-0.013707};
+Point(2060) = {0.417683,0.650000,-0.013681};
+Point(2061) = {0.414376,0.650000,-0.013628};
+Point(2062) = {0.411102,0.650000,-0.013548};
+Point(2063) = {0.407866,0.650000,-0.013443};
+Point(2064) = {0.404669,0.650000,-0.013312};
+Point(2065) = {0.401515,0.650000,-0.013155};
+Point(2066) = {0.398408,0.650000,-0.012973};
+Point(2067) = {0.395350,0.650000,-0.012765};
+Point(2068) = {0.392344,0.650000,-0.012533};
+Point(2069) = {0.389394,0.650000,-0.012279};
+Point(2070) = {0.386502,0.650000,-0.012002};
+Point(2071) = {0.383671,0.650000,-0.011705};
+Point(2072) = {0.380904,0.650000,-0.011390};
+Point(2073) = {0.378204,0.650000,-0.011057};
+Point(2074) = {0.375573,0.650000,-0.010711};
+Point(2075) = {0.373013,0.650000,-0.010353};
+Point(2076) = {0.370528,0.650000,-0.009985};
+Point(2077) = {0.368120,0.650000,-0.009610};
+Point(2078) = {0.365792,0.650000,-0.009230};
+Point(2079) = {0.363545,0.650000,-0.008849};
+Point(2080) = {0.361381,0.650000,-0.008469};
+Point(2081) = {0.359304,0.650000,-0.008092};
+Point(2082) = {0.357314,0.650000,-0.007721};
+Point(2083) = {0.355414,0.650000,-0.007352};
+Point(2084) = {0.353606,0.650000,-0.006982};
+Point(2085) = {0.351892,0.650000,-0.006609};
+Point(2086) = {0.350273,0.650000,-0.006228};
+Point(2087) = {0.348751,0.650000,-0.005838};
+Point(2088) = {0.347327,0.650000,-0.005439};
+Point(2089) = {0.346002,0.650000,-0.005031};
+Point(2090) = {0.344779,0.650000,-0.004617};
+Point(2091) = {0.343659,0.650000,-0.004197};
+Point(2092) = {0.342641,0.650000,-0.003773};
+Point(2093) = {0.341728,0.650000,-0.003343};
+Point(2094) = {0.340920,0.650000,-0.002904};
+Point(2095) = {0.340219,0.650000,-0.002454};
+Point(2096) = {0.339624,0.650000,-0.001991};
+Point(2097) = {0.339136,0.650000,-0.001512};
+Point(2098) = {0.338757,0.650000,-0.001015};
+Point(2099) = {0.338485,0.650000,-0.000499};
+Point(2100) = {0.338322,0.650000,0.000039};
+Point(2101) = {0.338268,0.650000,0.000601};
+Point(2102) = {0.338322,0.650000,0.001187};
+Point(2103) = {0.338485,0.650000,0.001791};
+Point(2104) = {0.338757,0.650000,0.002406};
+Point(2105) = {0.339136,0.650000,0.003023};
+Point(2106) = {0.339624,0.650000,0.003634};
+Point(2107) = {0.340219,0.650000,0.004234};
+Point(2108) = {0.340920,0.650000,0.004817};
+Point(2109) = {0.341728,0.650000,0.005380};
+Point(2110) = {0.342641,0.650000,0.005919};
+Point(2111) = {0.343659,0.650000,0.006431};
+Point(2112) = {0.344779,0.650000,0.006913};
+Point(2113) = {0.346002,0.650000,0.007366};
+Point(2114) = {0.347327,0.650000,0.007792};
+Point(2115) = {0.348751,0.650000,0.008192};
+Point(2116) = {0.350273,0.650000,0.008569};
+Point(2117) = {0.351892,0.650000,0.008924};
+Point(2118) = {0.353606,0.650000,0.009259};
+Point(2119) = {0.355414,0.650000,0.009576};
+Point(2120) = {0.357314,0.650000,0.009876};
+Point(2121) = {0.359304,0.650000,0.010161};
+Point(2122) = {0.361381,0.650000,0.010433};
+Point(2123) = {0.363545,0.650000,0.010692};
+Point(2124) = {0.365792,0.650000,0.010936};
+Point(2125) = {0.368120,0.650000,0.011165};
+Point(2126) = {0.370528,0.650000,0.011379};
+Point(2127) = {0.373013,0.650000,0.011577};
+Point(2128) = {0.375573,0.650000,0.011760};
+Point(2129) = {0.378204,0.650000,0.011927};
+Point(2130) = {0.380904,0.650000,0.012081};
+Point(2131) = {0.383671,0.650000,0.012222};
+Point(2132) = {0.386502,0.650000,0.012350};
+Point(2133) = {0.389394,0.650000,0.012465};
+Point(2134) = {0.392344,0.650000,0.012566};
+Point(2135) = {0.395350,0.650000,0.012653};
+Point(2136) = {0.398408,0.650000,0.012725};
+Point(2137) = {0.401515,0.650000,0.012783};
+Point(2138) = {0.404669,0.650000,0.012826};
+Point(2139) = {0.407866,0.650000,0.012856};
+Point(2140) = {0.411102,0.650000,0.012873};
+Point(2141) = {0.414376,0.650000,0.012878};
+Point(2142) = {0.417683,0.650000,0.012871};
+Point(2143) = {0.421021,0.650000,0.012852};
+Point(2144) = {0.424385,0.650000,0.012821};
+Point(2145) = {0.427773,0.650000,0.012775};
+Point(2146) = {0.431182,0.650000,0.012715};
+Point(2147) = {0.434608,0.650000,0.012639};
+Point(2148) = {0.438047,0.650000,0.012547};
+Point(2149) = {0.441497,0.650000,0.012439};
+Point(2150) = {0.444953,0.650000,0.012315};
+Point(2151) = {0.448413,0.650000,0.012175};
+Point(2152) = {0.451872,0.650000,0.012018};
+Point(2153) = {0.455329,0.650000,0.011846};
+Point(2154) = {0.458778,0.650000,0.011659};
+Point(2155) = {0.462218,0.650000,0.011457};
+Point(2156) = {0.465643,0.650000,0.011243};
+Point(2157) = {0.469052,0.650000,0.011015};
+Point(2158) = {0.472440,0.650000,0.010775};
+Point(2159) = {0.475805,0.650000,0.010524};
+Point(2160) = {0.479142,0.650000,0.010260};
+Point(2161) = {0.482449,0.650000,0.009985};
+Point(2162) = {0.485723,0.650000,0.009699};
+Point(2163) = {0.488960,0.650000,0.009400};
+Point(2164) = {0.492157,0.650000,0.009087};
+Point(2165) = {0.495310,0.650000,0.008758};
+Point(2166) = {0.498417,0.650000,0.008413};
+Point(2167) = {0.501475,0.650000,0.008050};
+Point(2168) = {0.504481,0.650000,0.007673};
+Point(2169) = {0.507431,0.650000,0.007284};
+Point(2170) = {0.510323,0.650000,0.006888};
+Point(2171) = {0.513154,0.650000,0.006487};
+Point(2172) = {0.515921,0.650000,0.006084};
+Point(2173) = {0.518622,0.650000,0.005682};
+Point(2174) = {0.521253,0.650000,0.005281};
+Point(2175) = {0.523812,0.650000,0.004882};
+Point(2176) = {0.526297,0.650000,0.004487};
+Point(2177) = {0.528705,0.650000,0.004097};
+Point(2178) = {0.531034,0.650000,0.003713};
+Point(2179) = {0.533281,0.650000,0.003341};
+Point(2180) = {0.535444,0.650000,0.002981};
+Point(2181) = {0.537522,0.650000,0.002637};
+Point(2182) = {0.539511,0.650000,0.002311};
+Point(2183) = {0.541411,0.650000,0.002003};
+Point(2184) = {0.543219,0.650000,0.001710};
+Point(2185) = {0.544933,0.650000,0.001432};
+Point(2186) = {0.546553,0.650000,0.001168};
+Point(2187) = {0.548075,0.650000,0.000916};
+Point(2188) = {0.549499,0.650000,0.000678};
+Point(2189) = {0.550823,0.650000,0.000455};
+Point(2190) = {0.552046,0.650000,0.000248};
+Point(2191) = {0.553167,0.650000,0.000059};
+Point(2192) = {0.554184,0.650000,-0.000111};
+Point(2193) = {0.555097,0.650000,-0.000289};
+Point(2194) = {0.555904,0.650000,-0.000450};
+Point(2195) = {0.556606,0.650000,-0.000594};
+Point(2196) = {0.557202,0.650000,-0.000719};
+Point(2197) = {0.557690,0.650000,-0.000823};
+Point(2198) = {0.558070,0.650000,-0.000905};
+Point(2199) = {0.558342,0.650000,-0.000965};
+Point(2200) = {0.558505,0.650000,-0.000999};
+// -- Airfoil 5
+Point(2501) = {0.611690,0.825000,0.001248};
+Point(2502) = {0.611646,0.825000,0.001246};
+Point(2503) = {0.611511,0.825000,0.001251};
+Point(2504) = {0.611286,0.825000,0.001262};
+Point(2505) = {0.610971,0.825000,0.001279};
+Point(2506) = {0.610567,0.825000,0.001300};
+Point(2507) = {0.610074,0.825000,0.001324};
+Point(2508) = {0.609493,0.825000,0.001350};
+Point(2509) = {0.608825,0.825000,0.001376};
+Point(2510) = {0.608069,0.825000,0.001402};
+Point(2511) = {0.607227,0.825000,0.001459};
+Point(2512) = {0.606300,0.825000,0.001512};
+Point(2513) = {0.605287,0.825000,0.001558};
+Point(2514) = {0.604191,0.825000,0.001596};
+Point(2515) = {0.603012,0.825000,0.001625};
+Point(2516) = {0.601752,0.825000,0.001643};
+Point(2517) = {0.600412,0.825000,0.001648};
+Point(2518) = {0.598993,0.825000,0.001638};
+Point(2519) = {0.597496,0.825000,0.001611};
+Point(2520) = {0.595924,0.825000,0.001563};
+Point(2521) = {0.594277,0.825000,0.001491};
+Point(2522) = {0.592557,0.825000,0.001395};
+Point(2523) = {0.590766,0.825000,0.001272};
+Point(2524) = {0.588906,0.825000,0.001122};
+Point(2525) = {0.586978,0.825000,0.000942};
+Point(2526) = {0.584985,0.825000,0.000732};
+Point(2527) = {0.582928,0.825000,0.000492};
+Point(2528) = {0.580810,0.825000,0.000221};
+Point(2529) = {0.578632,0.825000,-0.000078};
+Point(2530) = {0.576396,0.825000,-0.000405};
+Point(2531) = {0.574106,0.825000,-0.000757};
+Point(2532) = {0.571763,0.825000,-0.001133};
+Point(2533) = {0.569369,0.825000,-0.001531};
+Point(2534) = {0.566926,0.825000,-0.001950};
+Point(2535) = {0.564439,0.825000,-0.002386};
+Point(2536) = {0.561907,0.825000,-0.002837};
+Point(2537) = {0.559335,0.825000,-0.003302};
+Point(2538) = {0.556725,0.825000,-0.003778};
+Point(2539) = {0.554078,0.825000,-0.004263};
+Point(2540) = {0.551399,0.825000,-0.004755};
+Point(2541) = {0.548689,0.825000,-0.005253};
+Point(2542) = {0.545952,0.825000,-0.005753};
+Point(2543) = {0.543189,0.825000,-0.006250};
+Point(2544) = {0.540404,0.825000,-0.006738};
+Point(2545) = {0.537599,0.825000,-0.007211};
+Point(2546) = {0.534778,0.825000,-0.007661};
+Point(2547) = {0.531942,0.825000,-0.008084};
+Point(2548) = {0.529095,0.825000,-0.008478};
+Point(2549) = {0.526240,0.825000,-0.008841};
+Point(2550) = {0.523379,0.825000,-0.009172};
+Point(2551) = {0.520515,0.825000,-0.009469};
+Point(2552) = {0.517651,0.825000,-0.009732};
+Point(2553) = {0.514790,0.825000,-0.009960};
+Point(2554) = {0.511934,0.825000,-0.010156};
+Point(2555) = {0.509087,0.825000,-0.010319};
+Point(2556) = {0.506252,0.825000,-0.010451};
+Point(2557) = {0.503430,0.825000,-0.010552};
+Point(2558) = {0.500626,0.825000,-0.010626};
+Point(2559) = {0.497840,0.825000,-0.010673};
+Point(2560) = {0.495078,0.825000,-0.010697};
+Point(2561) = {0.492340,0.825000,-0.010701};
+Point(2562) = {0.489630,0.825000,-0.010686};
+Point(2563) = {0.486951,0.825000,-0.010652};
+Point(2564) = {0.484305,0.825000,-0.010600};
+Point(2565) = {0.481694,0.825000,-0.010530};
+Point(2566) = {0.479122,0.825000,-0.010441};
+Point(2567) = {0.476591,0.825000,-0.010335};
+Point(2568) = {0.474103,0.825000,-0.010211};
+Point(2569) = {0.471661,0.825000,-0.010071};
+Point(2570) = {0.469267,0.825000,-0.009916};
+Point(2571) = {0.466923,0.825000,-0.009746};
+Point(2572) = {0.464633,0.825000,-0.009563};
+Point(2573) = {0.462398,0.825000,-0.009368};
+Point(2574) = {0.460220,0.825000,-0.009161};
+Point(2575) = {0.458101,0.825000,-0.008946};
+Point(2576) = {0.456044,0.825000,-0.008723};
+Point(2577) = {0.454051,0.825000,-0.008494};
+Point(2578) = {0.452123,0.825000,-0.008259};
+Point(2579) = {0.450263,0.825000,-0.008021};
+Point(2580) = {0.448472,0.825000,-0.007779};
+Point(2581) = {0.446753,0.825000,-0.007536};
+Point(2582) = {0.445106,0.825000,-0.007291};
+Point(2583) = {0.443533,0.825000,-0.007046};
+Point(2584) = {0.442037,0.825000,-0.006799};
+Point(2585) = {0.440618,0.825000,-0.006551};
+Point(2586) = {0.439277,0.825000,-0.006301};
+Point(2587) = {0.438017,0.825000,-0.006048};
+Point(2588) = {0.436838,0.825000,-0.005790};
+Point(2589) = {0.435742,0.825000,-0.005524};
+Point(2590) = {0.434730,0.825000,-0.005245};
+Point(2591) = {0.433802,0.825000,-0.004952};
+Point(2592) = {0.432960,0.825000,-0.004641};
+Point(2593) = {0.432204,0.825000,-0.004312};
+Point(2594) = {0.431536,0.825000,-0.003967};
+Point(2595) = {0.430955,0.825000,-0.003605};
+Point(2596) = {0.430462,0.825000,-0.003228};
+Point(2597) = {0.430059,0.825000,-0.002836};
+Point(2598) = {0.429745,0.825000,-0.002426};
+Point(2599) = {0.429520,0.825000,-0.001997};
+Point(2600) = {0.429385,0.825000,-0.001544};
+Point(2601) = {0.429340,0.825000,-0.001067};
+Point(2602) = {0.429385,0.825000,-0.000563};
+Point(2603) = {0.429520,0.825000,-0.000040};
+Point(2604) = {0.429745,0.825000,0.000495};
+Point(2605) = {0.430059,0.825000,0.001032};
+Point(2606) = {0.430462,0.825000,0.001563};
+Point(2607) = {0.430955,0.825000,0.002082};
+Point(2608) = {0.431536,0.825000,0.002585};
+Point(2609) = {0.432204,0.825000,0.003072};
+Point(2610) = {0.432960,0.825000,0.003543};
+Point(2611) = {0.433802,0.825000,0.003995};
+Point(2612) = {0.434730,0.825000,0.004429};
+Point(2613) = {0.435742,0.825000,0.004844};
+Point(2614) = {0.436838,0.825000,0.005241};
+Point(2615) = {0.438017,0.825000,0.005621};
+Point(2616) = {0.439277,0.825000,0.005982};
+Point(2617) = {0.440618,0.825000,0.006326};
+Point(2618) = {0.442037,0.825000,0.006655};
+Point(2619) = {0.443533,0.825000,0.006970};
+Point(2620) = {0.445106,0.825000,0.007273};
+Point(2621) = {0.446753,0.825000,0.007566};
+Point(2622) = {0.448472,0.825000,0.007850};
+Point(2623) = {0.450263,0.825000,0.008125};
+Point(2624) = {0.452123,0.825000,0.008391};
+Point(2625) = {0.454051,0.825000,0.008645};
+Point(2626) = {0.456044,0.825000,0.008889};
+Point(2627) = {0.458101,0.825000,0.009121};
+Point(2628) = {0.460220,0.825000,0.009341};
+Point(2629) = {0.462398,0.825000,0.009550};
+Point(2630) = {0.464633,0.825000,0.009747};
+Point(2631) = {0.466923,0.825000,0.009932};
+Point(2632) = {0.469267,0.825000,0.010106};
+Point(2633) = {0.471661,0.825000,0.010268};
+Point(2634) = {0.474103,0.825000,0.010419};
+Point(2635) = {0.476591,0.825000,0.010558};
+Point(2636) = {0.479122,0.825000,0.010685};
+Point(2637) = {0.481694,0.825000,0.010801};
+Point(2638) = {0.484305,0.825000,0.010905};
+Point(2639) = {0.486951,0.825000,0.010998};
+Point(2640) = {0.489630,0.825000,0.011078};
+Point(2641) = {0.492340,0.825000,0.011147};
+Point(2642) = {0.495078,0.825000,0.011204};
+Point(2643) = {0.497840,0.825000,0.011251};
+Point(2644) = {0.500626,0.825000,0.011286};
+Point(2645) = {0.503430,0.825000,0.011311};
+Point(2646) = {0.506252,0.825000,0.011326};
+Point(2647) = {0.509087,0.825000,0.011331};
+Point(2648) = {0.511934,0.825000,0.011325};
+Point(2649) = {0.514790,0.825000,0.011308};
+Point(2650) = {0.517651,0.825000,0.011277};
+Point(2651) = {0.520515,0.825000,0.011232};
+Point(2652) = {0.523379,0.825000,0.011173};
+Point(2653) = {0.526240,0.825000,0.011099};
+Point(2654) = {0.529095,0.825000,0.011011};
+Point(2655) = {0.531942,0.825000,0.010911};
+Point(2656) = {0.534778,0.825000,0.010800};
+Point(2657) = {0.537599,0.825000,0.010679};
+Point(2658) = {0.540404,0.825000,0.010547};
+Point(2659) = {0.543189,0.825000,0.010405};
+Point(2660) = {0.545952,0.825000,0.010251};
+Point(2661) = {0.548689,0.825000,0.010086};
+Point(2662) = {0.551399,0.825000,0.009909};
+Point(2663) = {0.554078,0.825000,0.009720};
+Point(2664) = {0.556725,0.825000,0.009517};
+Point(2665) = {0.559335,0.825000,0.009299};
+Point(2666) = {0.561907,0.825000,0.009065};
+Point(2667) = {0.564439,0.825000,0.008814};
+Point(2668) = {0.566926,0.825000,0.008548};
+Point(2669) = {0.569369,0.825000,0.008267};
+Point(2670) = {0.571763,0.825000,0.007973};
+Point(2671) = {0.574106,0.825000,0.007667};
+Point(2672) = {0.576396,0.825000,0.007350};
+Point(2673) = {0.578632,0.825000,0.007026};
+Point(2674) = {0.580810,0.825000,0.006698};
+Point(2675) = {0.582928,0.825000,0.006369};
+Point(2676) = {0.584985,0.825000,0.006043};
+Point(2677) = {0.586978,0.825000,0.005722};
+Point(2678) = {0.588906,0.825000,0.005407};
+Point(2679) = {0.590766,0.825000,0.005100};
+Point(2680) = {0.592557,0.825000,0.004801};
+Point(2681) = {0.594277,0.825000,0.004510};
+Point(2682) = {0.595924,0.825000,0.004228};
+Point(2683) = {0.597496,0.825000,0.003955};
+Point(2684) = {0.598993,0.825000,0.003694};
+Point(2685) = {0.600412,0.825000,0.003444};
+Point(2686) = {0.601752,0.825000,0.003207};
+Point(2687) = {0.603012,0.825000,0.002983};
+Point(2688) = {0.604191,0.825000,0.002774};
+Point(2689) = {0.605287,0.825000,0.002578};
+Point(2690) = {0.606300,0.825000,0.002398};
+Point(2691) = {0.607227,0.825000,0.002233};
+Point(2692) = {0.608069,0.825000,0.002083};
+Point(2693) = {0.608825,0.825000,0.001915};
+Point(2694) = {0.609493,0.825000,0.001763};
+Point(2695) = {0.610074,0.825000,0.001628};
+Point(2696) = {0.610567,0.825000,0.001511};
+Point(2697) = {0.610971,0.825000,0.001414};
+Point(2698) = {0.611286,0.825000,0.001338};
+Point(2699) = {0.611511,0.825000,0.001284};
+Point(2700) = {0.611646,0.825000,0.001254};
+// -- Airfoil 6
+Point(3001) = {0.649731,0.950000,0.002532};
+Point(3002) = {0.649694,0.950000,0.002533};
+Point(3003) = {0.649579,0.950000,0.002537};
+Point(3004) = {0.649388,0.950000,0.002544};
+Point(3005) = {0.649119,0.950000,0.002553};
+Point(3006) = {0.648775,0.950000,0.002563};
+Point(3007) = {0.648355,0.950000,0.002574};
+Point(3008) = {0.647860,0.950000,0.002585};
+Point(3009) = {0.647291,0.950000,0.002594};
+Point(3010) = {0.646647,0.950000,0.002601};
+Point(3011) = {0.645930,0.950000,0.002649};
+Point(3012) = {0.645139,0.950000,0.002694};
+Point(3013) = {0.644277,0.950000,0.002734};
+Point(3014) = {0.643343,0.950000,0.002768};
+Point(3015) = {0.642339,0.950000,0.002795};
+Point(3016) = {0.641266,0.950000,0.002813};
+Point(3017) = {0.640124,0.950000,0.002821};
+Point(3018) = {0.638915,0.950000,0.002815};
+Point(3019) = {0.637640,0.950000,0.002795};
+Point(3020) = {0.636300,0.950000,0.002756};
+Point(3021) = {0.634897,0.950000,0.002698};
+Point(3022) = {0.633432,0.950000,0.002618};
+Point(3023) = {0.631907,0.950000,0.002515};
+Point(3024) = {0.630322,0.950000,0.002390};
+Point(3025) = {0.628680,0.950000,0.002242};
+Point(3026) = {0.626982,0.950000,0.002071};
+Point(3027) = {0.625230,0.950000,0.001877};
+Point(3028) = {0.623425,0.950000,0.001657};
+Point(3029) = {0.621570,0.950000,0.001409};
+Point(3030) = {0.619666,0.950000,0.001128};
+Point(3031) = {0.617715,0.950000,0.000812};
+Point(3032) = {0.615718,0.950000,0.000460};
+Point(3033) = {0.613679,0.950000,0.000078};
+Point(3034) = {0.611598,0.950000,-0.000328};
+Point(3035) = {0.609479,0.950000,-0.000751};
+Point(3036) = {0.607323,0.950000,-0.001183};
+Point(3037) = {0.605132,0.950000,-0.001618};
+Point(3038) = {0.602908,0.950000,-0.002057};
+Point(3039) = {0.600653,0.950000,-0.002499};
+Point(3040) = {0.598371,0.950000,-0.002946};
+Point(3041) = {0.596062,0.950000,-0.003399};
+Point(3042) = {0.593730,0.950000,-0.003857};
+Point(3043) = {0.591377,0.950000,-0.004316};
+Point(3044) = {0.589004,0.950000,-0.004772};
+Point(3045) = {0.586615,0.950000,-0.005218};
+Point(3046) = {0.584211,0.950000,-0.005650};
+Point(3047) = {0.581796,0.950000,-0.006062};
+Point(3048) = {0.579371,0.950000,-0.006450};
+Point(3049) = {0.576938,0.950000,-0.006813};
+Point(3050) = {0.574501,0.950000,-0.007146};
+Point(3051) = {0.572061,0.950000,-0.007447};
+Point(3052) = {0.569622,0.950000,-0.007714};
+Point(3053) = {0.567184,0.950000,-0.007947};
+Point(3054) = {0.564752,0.950000,-0.008151};
+Point(3055) = {0.562326,0.950000,-0.008327};
+Point(3056) = {0.559911,0.950000,-0.008478};
+Point(3057) = {0.557507,0.950000,-0.008607};
+Point(3058) = {0.555118,0.950000,-0.008715};
+Point(3059) = {0.552745,0.950000,-0.008802};
+Point(3060) = {0.550392,0.950000,-0.008870};
+Point(3061) = {0.548060,0.950000,-0.008919};
+Point(3062) = {0.545751,0.950000,-0.008950};
+Point(3063) = {0.543469,0.950000,-0.008965};
+Point(3064) = {0.541215,0.950000,-0.008964};
+Point(3065) = {0.538991,0.950000,-0.008950};
+Point(3066) = {0.536800,0.950000,-0.008924};
+Point(3067) = {0.534643,0.950000,-0.008887};
+Point(3068) = {0.532524,0.950000,-0.008840};
+Point(3069) = {0.530444,0.950000,-0.008782};
+Point(3070) = {0.528404,0.950000,-0.008712};
+Point(3071) = {0.526408,0.950000,-0.008631};
+Point(3072) = {0.524457,0.950000,-0.008539};
+Point(3073) = {0.522552,0.950000,-0.008435};
+Point(3074) = {0.520697,0.950000,-0.008323};
+Point(3075) = {0.518893,0.950000,-0.008203};
+Point(3076) = {0.517140,0.950000,-0.008077};
+Point(3077) = {0.515442,0.950000,-0.007945};
+Point(3078) = {0.513800,0.950000,-0.007810};
+Point(3079) = {0.512216,0.950000,-0.007673};
+Point(3080) = {0.510690,0.950000,-0.007536};
+Point(3081) = {0.509225,0.950000,-0.007400};
+Point(3082) = {0.507822,0.950000,-0.007265};
+Point(3083) = {0.506482,0.950000,-0.007128};
+Point(3084) = {0.505208,0.950000,-0.006984};
+Point(3085) = {0.503999,0.950000,-0.006826};
+Point(3086) = {0.502857,0.950000,-0.006651};
+Point(3087) = {0.501783,0.950000,-0.006455};
+Point(3088) = {0.500779,0.950000,-0.006238};
+Point(3089) = {0.499845,0.950000,-0.006006};
+Point(3090) = {0.498983,0.950000,-0.005761};
+Point(3091) = {0.498193,0.950000,-0.005506};
+Point(3092) = {0.497475,0.950000,-0.005246};
+Point(3093) = {0.496831,0.950000,-0.004976};
+Point(3094) = {0.496262,0.950000,-0.004697};
+Point(3095) = {0.495767,0.950000,-0.004406};
+Point(3096) = {0.495347,0.950000,-0.004101};
+Point(3097) = {0.495004,0.950000,-0.003780};
+Point(3098) = {0.494736,0.950000,-0.003444};
+Point(3099) = {0.494545,0.950000,-0.003094};
+Point(3100) = {0.494430,0.950000,-0.002730};
+Point(3101) = {0.494391,0.950000,-0.002353};
+Point(3102) = {0.494430,0.950000,-0.001965};
+Point(3103) = {0.494545,0.950000,-0.001565};
+Point(3104) = {0.494736,0.950000,-0.001153};
+Point(3105) = {0.495004,0.950000,-0.000729};
+Point(3106) = {0.495347,0.950000,-0.000293};
+Point(3107) = {0.495767,0.950000,0.000155};
+Point(3108) = {0.496262,0.950000,0.000609};
+Point(3109) = {0.496831,0.950000,0.001062};
+Point(3110) = {0.497475,0.950000,0.001509};
+Point(3111) = {0.498193,0.950000,0.001942};
+Point(3112) = {0.498983,0.950000,0.002357};
+Point(3113) = {0.499845,0.950000,0.002754};
+Point(3114) = {0.500779,0.950000,0.003135};
+Point(3115) = {0.501783,0.950000,0.003503};
+Point(3116) = {0.502857,0.950000,0.003859};
+Point(3117) = {0.503999,0.950000,0.004205};
+Point(3118) = {0.505208,0.950000,0.004542};
+Point(3119) = {0.506482,0.950000,0.004869};
+Point(3120) = {0.507822,0.950000,0.005187};
+Point(3121) = {0.509225,0.950000,0.005496};
+Point(3122) = {0.510690,0.950000,0.005796};
+Point(3123) = {0.512216,0.950000,0.006087};
+Point(3124) = {0.513800,0.950000,0.006369};
+Point(3125) = {0.515442,0.950000,0.006642};
+Point(3126) = {0.517140,0.950000,0.006905};
+Point(3127) = {0.518893,0.950000,0.007159};
+Point(3128) = {0.520697,0.950000,0.007404};
+Point(3129) = {0.522552,0.950000,0.007639};
+Point(3130) = {0.524457,0.950000,0.007865};
+Point(3131) = {0.526408,0.950000,0.008082};
+Point(3132) = {0.528404,0.950000,0.008289};
+Point(3133) = {0.530444,0.950000,0.008486};
+Point(3134) = {0.532524,0.950000,0.008673};
+Point(3135) = {0.534643,0.950000,0.008850};
+Point(3136) = {0.536800,0.950000,0.009017};
+Point(3137) = {0.538991,0.950000,0.009173};
+Point(3138) = {0.541215,0.950000,0.009318};
+Point(3139) = {0.543469,0.950000,0.009453};
+Point(3140) = {0.545751,0.950000,0.009579};
+Point(3141) = {0.548060,0.950000,0.009695};
+Point(3142) = {0.550392,0.950000,0.009803};
+Point(3143) = {0.552745,0.950000,0.009901};
+Point(3144) = {0.555118,0.950000,0.009990};
+Point(3145) = {0.557507,0.950000,0.010069};
+Point(3146) = {0.559911,0.950000,0.010137};
+Point(3147) = {0.562326,0.950000,0.010195};
+Point(3148) = {0.564752,0.950000,0.010242};
+Point(3149) = {0.567184,0.950000,0.010280};
+Point(3150) = {0.569622,0.950000,0.010307};
+Point(3151) = {0.572061,0.950000,0.010325};
+Point(3152) = {0.574501,0.950000,0.010335};
+Point(3153) = {0.576938,0.950000,0.010334};
+Point(3154) = {0.579371,0.950000,0.010324};
+Point(3155) = {0.581796,0.950000,0.010303};
+Point(3156) = {0.584211,0.950000,0.010271};
+Point(3157) = {0.586615,0.950000,0.010227};
+Point(3158) = {0.589004,0.950000,0.010171};
+Point(3159) = {0.591377,0.950000,0.010103};
+Point(3160) = {0.593730,0.950000,0.010023};
+Point(3161) = {0.596062,0.950000,0.009930};
+Point(3162) = {0.598371,0.950000,0.009825};
+Point(3163) = {0.600653,0.950000,0.009706};
+Point(3164) = {0.602908,0.950000,0.009573};
+Point(3165) = {0.605132,0.950000,0.009425};
+Point(3166) = {0.607323,0.950000,0.009262};
+Point(3167) = {0.609479,0.950000,0.009083};
+Point(3168) = {0.611598,0.950000,0.008887};
+Point(3169) = {0.613679,0.950000,0.008674};
+Point(3170) = {0.615718,0.950000,0.008445};
+Point(3171) = {0.617715,0.950000,0.008199};
+Point(3172) = {0.619666,0.950000,0.007937};
+Point(3173) = {0.621570,0.950000,0.007662};
+Point(3174) = {0.623425,0.950000,0.007380};
+Point(3175) = {0.625230,0.950000,0.007094};
+Point(3176) = {0.626982,0.950000,0.006809};
+Point(3177) = {0.628680,0.950000,0.006529};
+Point(3178) = {0.630322,0.950000,0.006254};
+Point(3179) = {0.631907,0.950000,0.005985};
+Point(3180) = {0.633432,0.950000,0.005723};
+Point(3181) = {0.634897,0.950000,0.005468};
+Point(3182) = {0.636300,0.950000,0.005221};
+Point(3183) = {0.637640,0.950000,0.004981};
+Point(3184) = {0.638915,0.950000,0.004751};
+Point(3185) = {0.640124,0.950000,0.004531};
+Point(3186) = {0.641266,0.950000,0.004322};
+Point(3187) = {0.642339,0.950000,0.004126};
+Point(3188) = {0.643343,0.950000,0.003941};
+Point(3189) = {0.644277,0.950000,0.003769};
+Point(3190) = {0.645139,0.950000,0.003610};
+Point(3191) = {0.645930,0.950000,0.003464};
+Point(3192) = {0.646647,0.950000,0.003332};
+Point(3193) = {0.647291,0.950000,0.003172};
+Point(3194) = {0.647860,0.950000,0.003028};
+Point(3195) = {0.648355,0.950000,0.002900};
+Point(3196) = {0.648775,0.950000,0.002790};
+Point(3197) = {0.649119,0.950000,0.002698};
+Point(3198) = {0.649388,0.950000,0.002625};
+Point(3199) = {0.649579,0.950000,0.002573};
+Point(3200) = {0.649694,0.950000,0.002542};
+// -- Airfoil 7
+Point(3501) = {0.664862,1.000000,0.002901};
+Point(3502) = {0.664827,1.000000,0.002902};
+Point(3503) = {0.664721,1.000000,0.002906};
+Point(3504) = {0.664542,1.000000,0.002910};
+Point(3505) = {0.664293,1.000000,0.002915};
+Point(3506) = {0.663973,1.000000,0.002920};
+Point(3507) = {0.663582,1.000000,0.002925};
+Point(3508) = {0.663122,1.000000,0.002930};
+Point(3509) = {0.662592,1.000000,0.002933};
+Point(3510) = {0.661994,1.000000,0.002936};
+Point(3511) = {0.661327,1.000000,0.002986};
+Point(3512) = {0.660592,1.000000,0.003039};
+Point(3513) = {0.659790,1.000000,0.003091};
+Point(3514) = {0.658922,1.000000,0.003138};
+Point(3515) = {0.657988,1.000000,0.003179};
+Point(3516) = {0.656990,1.000000,0.003209};
+Point(3517) = {0.655928,1.000000,0.003226};
+Point(3518) = {0.654804,1.000000,0.003227};
+Point(3519) = {0.653618,1.000000,0.003212};
+Point(3520) = {0.652373,1.000000,0.003178};
+Point(3521) = {0.651068,1.000000,0.003125};
+Point(3522) = {0.649706,1.000000,0.003050};
+Point(3523) = {0.648287,1.000000,0.002954};
+Point(3524) = {0.646813,1.000000,0.002833};
+Point(3525) = {0.645286,1.000000,0.002689};
+Point(3526) = {0.643708,1.000000,0.002518};
+Point(3527) = {0.642078,1.000000,0.002321};
+Point(3528) = {0.640400,1.000000,0.002100};
+Point(3529) = {0.638675,1.000000,0.001854};
+Point(3530) = {0.636904,1.000000,0.001586};
+Point(3531) = {0.635090,1.000000,0.001298};
+Point(3532) = {0.633233,1.000000,0.000990};
+Point(3533) = {0.631337,1.000000,0.000662};
+Point(3534) = {0.629402,1.000000,0.000317};
+Point(3535) = {0.627431,1.000000,-0.000047};
+Point(3536) = {0.625426,1.000000,-0.000429};
+Point(3537) = {0.623389,1.000000,-0.000828};
+Point(3538) = {0.621321,1.000000,-0.001241};
+Point(3539) = {0.619225,1.000000,-0.001667};
+Point(3540) = {0.617102,1.000000,-0.002103};
+Point(3541) = {0.614955,1.000000,-0.002546};
+Point(3542) = {0.612787,1.000000,-0.002994};
+Point(3543) = {0.610598,1.000000,-0.003442};
+Point(3544) = {0.608392,1.000000,-0.003885};
+Point(3545) = {0.606170,1.000000,-0.004320};
+Point(3546) = {0.603935,1.000000,-0.004740};
+Point(3547) = {0.601689,1.000000,-0.005142};
+Point(3548) = {0.599434,1.000000,-0.005523};
+Point(3549) = {0.597172,1.000000,-0.005881};
+Point(3550) = {0.594905,1.000000,-0.006213};
+Point(3551) = {0.592637,1.000000,-0.006517};
+Point(3552) = {0.590368,1.000000,-0.006792};
+Point(3553) = {0.588102,1.000000,-0.007038};
+Point(3554) = {0.585840,1.000000,-0.007257};
+Point(3555) = {0.583585,1.000000,-0.007451};
+Point(3556) = {0.581338,1.000000,-0.007622};
+Point(3557) = {0.579103,1.000000,-0.007771};
+Point(3558) = {0.576881,1.000000,-0.007899};
+Point(3559) = {0.574675,1.000000,-0.008007};
+Point(3560) = {0.572487,1.000000,-0.008096};
+Point(3561) = {0.570318,1.000000,-0.008167};
+Point(3562) = {0.568171,1.000000,-0.008221};
+Point(3563) = {0.566049,1.000000,-0.008259};
+Point(3564) = {0.563953,1.000000,-0.008283};
+Point(3565) = {0.561885,1.000000,-0.008296};
+Point(3566) = {0.559847,1.000000,-0.008298};
+Point(3567) = {0.557842,1.000000,-0.008292};
+Point(3568) = {0.555871,1.000000,-0.008276};
+Point(3569) = {0.553937,1.000000,-0.008251};
+Point(3570) = {0.552040,1.000000,-0.008217};
+Point(3571) = {0.550184,1.000000,-0.008172};
+Point(3572) = {0.548370,1.000000,-0.008116};
+Point(3573) = {0.546599,1.000000,-0.008051};
+Point(3574) = {0.544874,1.000000,-0.007977};
+Point(3575) = {0.543195,1.000000,-0.007895};
+Point(3576) = {0.541566,1.000000,-0.007807};
+Point(3577) = {0.539987,1.000000,-0.007714};
+Point(3578) = {0.538460,1.000000,-0.007616};
+Point(3579) = {0.536987,1.000000,-0.007514};
+Point(3580) = {0.535568,1.000000,-0.007408};
+Point(3581) = {0.534206,1.000000,-0.007298};
+Point(3582) = {0.532901,1.000000,-0.007185};
+Point(3583) = {0.531655,1.000000,-0.007068};
+Point(3584) = {0.530470,1.000000,-0.006947};
+Point(3585) = {0.529346,1.000000,-0.006819};
+Point(3586) = {0.528284,1.000000,-0.006686};
+Point(3587) = {0.527286,1.000000,-0.006545};
+Point(3588) = {0.526352,1.000000,-0.006396};
+Point(3589) = {0.525484,1.000000,-0.006235};
+Point(3590) = {0.524682,1.000000,-0.006061};
+Point(3591) = {0.523947,1.000000,-0.005872};
+Point(3592) = {0.523280,1.000000,-0.005667};
+Point(3593) = {0.522681,1.000000,-0.005445};
+Point(3594) = {0.522151,1.000000,-0.005209};
+Point(3595) = {0.521691,1.000000,-0.004958};
+Point(3596) = {0.521301,1.000000,-0.004693};
+Point(3597) = {0.520981,1.000000,-0.004416};
+Point(3598) = {0.520732,1.000000,-0.004124};
+Point(3599) = {0.520554,1.000000,-0.003813};
+Point(3600) = {0.520447,1.000000,-0.003481};
+Point(3601) = {0.520412,1.000000,-0.003124};
+Point(3602) = {0.520447,1.000000,-0.002742};
+Point(3603) = {0.520554,1.000000,-0.002338};
+Point(3604) = {0.520732,1.000000,-0.001916};
+Point(3605) = {0.520981,1.000000,-0.001482};
+Point(3606) = {0.521301,1.000000,-0.001041};
+Point(3607) = {0.521691,1.000000,-0.000599};
+Point(3608) = {0.522151,1.000000,-0.000159};
+Point(3609) = {0.522681,1.000000,0.000273};
+Point(3610) = {0.523280,1.000000,0.000691};
+Point(3611) = {0.523947,1.000000,0.001092};
+Point(3612) = {0.524682,1.000000,0.001471};
+Point(3613) = {0.525484,1.000000,0.001833};
+Point(3614) = {0.526352,1.000000,0.002185};
+Point(3615) = {0.527286,1.000000,0.002533};
+Point(3616) = {0.528284,1.000000,0.002884};
+Point(3617) = {0.529346,1.000000,0.003241};
+Point(3618) = {0.530470,1.000000,0.003601};
+Point(3619) = {0.531655,1.000000,0.003956};
+Point(3620) = {0.532901,1.000000,0.004301};
+Point(3621) = {0.534206,1.000000,0.004628};
+Point(3622) = {0.535568,1.000000,0.004934};
+Point(3623) = {0.536987,1.000000,0.005220};
+Point(3624) = {0.538460,1.000000,0.005493};
+Point(3625) = {0.539987,1.000000,0.005756};
+Point(3626) = {0.541566,1.000000,0.006015};
+Point(3627) = {0.543195,1.000000,0.006273};
+Point(3628) = {0.544874,1.000000,0.006530};
+Point(3629) = {0.546599,1.000000,0.006782};
+Point(3630) = {0.548370,1.000000,0.007027};
+Point(3631) = {0.550184,1.000000,0.007263};
+Point(3632) = {0.552040,1.000000,0.007487};
+Point(3633) = {0.553937,1.000000,0.007700};
+Point(3634) = {0.555871,1.000000,0.007903};
+Point(3635) = {0.557842,1.000000,0.008094};
+Point(3636) = {0.559847,1.000000,0.008275};
+Point(3637) = {0.561885,1.000000,0.008446};
+Point(3638) = {0.563953,1.000000,0.008606};
+Point(3639) = {0.566049,1.000000,0.008758};
+Point(3640) = {0.568171,1.000000,0.008899};
+Point(3641) = {0.570318,1.000000,0.009032};
+Point(3642) = {0.572487,1.000000,0.009156};
+Point(3643) = {0.574675,1.000000,0.009271};
+Point(3644) = {0.576881,1.000000,0.009376};
+Point(3645) = {0.579103,1.000000,0.009472};
+Point(3646) = {0.581338,1.000000,0.009559};
+Point(3647) = {0.583585,1.000000,0.009636};
+Point(3648) = {0.585840,1.000000,0.009704};
+Point(3649) = {0.588102,1.000000,0.009763};
+Point(3650) = {0.590368,1.000000,0.009812};
+Point(3651) = {0.592637,1.000000,0.009851};
+Point(3652) = {0.594905,1.000000,0.009881};
+Point(3653) = {0.597172,1.000000,0.009902};
+Point(3654) = {0.599434,1.000000,0.009914};
+Point(3655) = {0.601689,1.000000,0.009916};
+Point(3656) = {0.603935,1.000000,0.009910};
+Point(3657) = {0.606170,1.000000,0.009894};
+Point(3658) = {0.608392,1.000000,0.009868};
+Point(3659) = {0.610598,1.000000,0.009831};
+Point(3660) = {0.612787,1.000000,0.009783};
+Point(3661) = {0.614955,1.000000,0.009721};
+Point(3662) = {0.617102,1.000000,0.009646};
+Point(3663) = {0.619225,1.000000,0.009556};
+Point(3664) = {0.621321,1.000000,0.009451};
+Point(3665) = {0.623389,1.000000,0.009331};
+Point(3666) = {0.625426,1.000000,0.009194};
+Point(3667) = {0.627431,1.000000,0.009041};
+Point(3668) = {0.629402,1.000000,0.008871};
+Point(3669) = {0.631337,1.000000,0.008684};
+Point(3670) = {0.633233,1.000000,0.008481};
+Point(3671) = {0.635090,1.000000,0.008262};
+Point(3672) = {0.636904,1.000000,0.008027};
+Point(3673) = {0.638675,1.000000,0.007780};
+Point(3674) = {0.640400,1.000000,0.007523};
+Point(3675) = {0.642078,1.000000,0.007259};
+Point(3676) = {0.643708,1.000000,0.006993};
+Point(3677) = {0.645286,1.000000,0.006727};
+Point(3678) = {0.646813,1.000000,0.006463};
+Point(3679) = {0.648287,1.000000,0.006203};
+Point(3680) = {0.649706,1.000000,0.005948};
+Point(3681) = {0.651068,1.000000,0.005700};
+Point(3682) = {0.652373,1.000000,0.005461};
+Point(3683) = {0.653618,1.000000,0.005231};
+Point(3684) = {0.654804,1.000000,0.005011};
+Point(3685) = {0.655928,1.000000,0.004803};
+Point(3686) = {0.656990,1.000000,0.004606};
+Point(3687) = {0.657988,1.000000,0.004422};
+Point(3688) = {0.658922,1.000000,0.004250};
+Point(3689) = {0.659790,1.000000,0.004090};
+Point(3690) = {0.660592,1.000000,0.003942};
+Point(3691) = {0.661327,1.000000,0.003805};
+Point(3692) = {0.661994,1.000000,0.003680};
+Point(3693) = {0.662592,1.000000,0.003522};
+Point(3694) = {0.663122,1.000000,0.003381};
+Point(3695) = {0.663582,1.000000,0.003257};
+Point(3696) = {0.663973,1.000000,0.003151};
+Point(3697) = {0.664293,1.000000,0.003062};
+Point(3698) = {0.664542,1.000000,0.002993};
+Point(3699) = {0.664721,1.000000,0.002942};
+Point(3700) = {0.664827,1.000000,0.002912};
+
+// --- Wing lines ---
+// -- Airfoil 0
+Spline(1) = {1:101};
+Spline(2) = {101:200, 1};
+// -- Airfoil 1
+Spline(11) = {501:601};
+Spline(12) = {601:700, 501};
+// -- Airfoil 2
+Spline(21) = {1001:1101};
+Spline(22) = {1101:1200, 1001};
+// -- Airfoil 3
+Spline(31) = {1501:1601};
+Spline(32) = {1601:1700, 1501};
+// -- Airfoil 4
+Spline(41) = {2001:2101};
+Spline(42) = {2101:2200, 2001};
+// -- Airfoil 5
+Spline(51) = {2501:2601};
+Spline(52) = {2601:2700, 2501};
+// -- Airfoil 6
+Spline(61) = {3001:3101};
+Spline(62) = {3101:3200, 3001};
+// -- Airfoil 7
+Spline(71) = {3501:3601};
+Spline(72) = {3601:3700, 3501};
+// -- Planform 0
+Line(101) = {1,501};
+Line(102) = {101,601};
+// -- Planform 1
+Line(111) = {501,1001};
+Line(112) = {601,1101};
+// -- Planform 2
+Line(121) = {1001,1501};
+Line(122) = {1101,1601};
+// -- Planform 3
+Line(131) = {1501,2001};
+Line(132) = {1601,2101};
+// -- Planform 4
+Line(141) = {2001,2501};
+Line(142) = {2101,2601};
+// -- Planform 5
+Line(151) = {2501,3001};
+Line(152) = {2601,3101};
+// -- Planform 6
+Line(161) = {3001,3501};
+Line(162) = {3101,3601};
+
+// --- Wing line loops and surfaces ---
+// -- Planform 0
+Line Loop(1) = {1,102,-11,-101};
+Line Loop(2) = {2,101,-12,-102};
+Surface(1) = {1};
+Surface(2) = {2};
+// -- Planform 1
+Line Loop(11) = {11,112,-21,-111};
+Line Loop(12) = {12,111,-22,-112};
+Surface(11) = {11};
+Surface(12) = {12};
+// -- Planform 2
+Line Loop(21) = {21,122,-31,-121};
+Line Loop(22) = {22,121,-32,-122};
+Surface(21) = {21};
+Surface(22) = {22};
+// -- Planform 3
+Line Loop(31) = {31,132,-41,-131};
+Line Loop(32) = {32,131,-42,-132};
+Surface(31) = {31};
+Surface(32) = {32};
+// -- Planform 4
+Line Loop(41) = {41,142,-51,-141};
+Line Loop(42) = {42,141,-52,-142};
+Surface(41) = {41};
+Surface(42) = {42};
+// -- Planform 5
+Line Loop(51) = {51,152,-61,-151};
+Line Loop(52) = {52,151,-62,-152};
+Surface(51) = {51};
+Surface(52) = {52};
+// -- Planform 6
+Line Loop(61) = {61,162,-71,-161};
+Line Loop(62) = {62,161,-72,-162};
+Surface(61) = {61};
+Surface(62) = {62};
+
+// --- Mesh ---
+Transfinite Line{1, 2, 11, 12, 21, 22, 31, 32, 41, 42, 51, 52, 61, 62, 71, 72} = nC + 1 Using Bump bC;
+Transfinite Line{101, 102} = 5 * rS + 1;
+Transfinite Line{111, 112} = 3 * rS + 1;
+Transfinite Line{121, 122} = 4 * rS + 1;
+Transfinite Line{131, 132} = 4 * rS + 1;
+Transfinite Line{141, 142} = 4 * rS + 1;
+Transfinite Line{151, 152} = 3 * rS + 1;
+Transfinite Line{161, 162} = 1 * rS + 1;
+Transfinite Surface{1, 2, 11, 12, 21, 22, 31, 32, 41, 42, 51, 52, 61, 62};
+Recombine Surface{1, 2, 11, 12, 21, 22, 31, 32, 41, 42, 51, 52, 61, 62};
+
+// --- Wing physical groups ---
+Physical Line("te") = {101, 111, 121, 131, 141, 151, 161};
+Physical Surface("wing") = {2, 12, 22, 32, 42, 52, 62, 1, 11, 21, 31, 41, 51, 61};
diff --git a/sdpm/models/oneraM6.geo b/sdpm/models/oneraM6.geo
index 418d940..3821d39 100644
--- a/sdpm/models/oneraM6.geo
+++ b/sdpm/models/oneraM6.geo
@@ -2,8 +2,8 @@
// Mesh parameters
// domain and mesh
-DefineConstant[ wNc = { 50, Name "Number of panel along chord" }
- wNs = { 10, Name "Number of panel along span" }
+DefineConstant[ nC = { 50, Name "Number of panel along chord" }
+ nS = { 10, Name "Number of panel along span" }
bC = { 0.05, Name "Progression along chord" }
pS = { 1.0, Name "Progression along span" } ];
@@ -315,13 +315,13 @@ Surface(11) = {-11};
Surface(12) = {-12};
//// MESH
-Transfinite Line{1,2,3,4} = wNc + 1 Using Bump bC;
-Transfinite Line{61,62} = wNs + 1 Using Progression pS;
+Transfinite Line{1,2,3,4} = nC + 1 Using Bump bC;
+Transfinite Line{61,62} = nS + 1 Using Progression pS;
Transfinite Surface{11,12};
Recombine Surface{11,12};
//// PHYSICAL GROUPS
// Trailing edge
-Physical Line("wTe") = {61};
+Physical Line("te") = {61};
// Wing:
Physical Surface("wing") = {11,12};
diff --git a/sdpm/src/sdpm.h b/sdpm/src/sdpm.h
index 4fe47a5..c580d78 100644
--- a/sdpm/src/sdpm.h
+++ b/sdpm/src/sdpm.h
@@ -32,12 +32,17 @@
#endif
#include <Eigen/Dense>
+#include <complex>
using sdpmDouble = double;
+using sdpmComplex = std::complex<double>;
using sdpmMatrixXd = Eigen::Matrix<sdpmDouble, Eigen::Dynamic, Eigen::Dynamic>;
using sdpmVectorXd = Eigen::Matrix<sdpmDouble, Eigen::Dynamic, 1>;
using sdpmMatrix3d = Eigen::Matrix<sdpmDouble, 3, 3>;
using sdpmVector3d = Eigen::Matrix<sdpmDouble, 3, 1>;
+using sdpmMatrixXcd = Eigen::Matrix<sdpmComplex, Eigen::Dynamic, Eigen::Dynamic>;
+using sdpmVectorXcd = Eigen::Matrix<sdpmComplex, Eigen::Dynamic, 1>;
+using sdpmVector3cd = Eigen::Matrix<sdpmComplex, 3, 1>;
namespace sdpm
{
@@ -60,6 +65,7 @@ class GmshExport;
class Problem;
class Body;
class Wake;
+class Motion;
class Edge;
class Builder;
class Solver;
diff --git a/sdpm/src/sdpmBody.cpp b/sdpm/src/sdpmBody.cpp
index 9f9e496..1cae6e3 100644
--- a/sdpm/src/sdpmBody.cpp
+++ b/sdpm/src/sdpmBody.cpp
@@ -21,10 +21,19 @@
#include "sdpmQuad4.h"
#include "sdpmTag.h"
#include "sdpmGmshExport.h"
+#include <unordered_set>
#include <iostream>
using namespace sdpm;
-Body::Body(Mesh &msh, std::string const &name, std::string const &teName, double xF) : Group(msh, name), cl(0), cd(0), cs(0), cm(0)
+Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
+ double shedLgth, size_t nShedDiv,
+ size_t nM, size_t nF,
+ double ux) : Group(msh, name),
+ _cl(0), _cd(0), _cs(0), _cm(0),
+ _cl1(nM, std::vector<sdpmComplex>(nF, sdpmComplex(0., 0.))),
+ _cd1(nM, std::vector<sdpmComplex>(nF, sdpmComplex(0., 0.))),
+ _cs1(nM, std::vector<sdpmComplex>(nF, sdpmComplex(0., 0.))),
+ _cm1(nM, std::vector<sdpmComplex>(nF, sdpmComplex(0., 0.)))
{
// Create wake
std::cout << "Creating wake... " << std::flush;
@@ -43,30 +52,34 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName, double
auto it = tags.find(wkName);
if (it == tags.end())
{
+ // sanity check
+ if (shedLgth <= 0.)
+ {
+ std::stringstream err;
+ err << "sdpm::Body: zero or negative characteristic chord length given for " << *_tag << "!\n";
+ throw std::runtime_error(err.str());
+ }
// create tag and add it to the mesh
Tag *tagp = new Tag(tags.size() + 1, wkName, 2);
_msh.addTag(wkName, tagp);
- // translate TE nodes (along x-coordinate)
- std::vector<Node *> wkNodes;
- for (auto n : teNodes)
+ // create wake nodes and elements
+ std::vector<Node *> wkNodes = teNodes;
+ for (size_t i = 0; i < 10 * nShedDiv; ++i)
{
- sdpmVector3d pos = n->getCoords();
- // sanity check
- if (xF <= pos(0))
+ // translate TE nodes along x-coordinate
+ sdpmVector3d delta((i + 1) * shedLgth / nShedDiv, 0., 0.);
+ for (auto n : teNodes)
{
- std::stringstream err;
- err << "sdpm::Body: zero or negative length wake requested for " << *_tag << "!\n";
- throw std::runtime_error(err.str());
+ Node *nN = new Node(_msh.getNodes().back()->getId() + 1, n->getCoords() + delta);
+ wkNodes.push_back(nN);
+ _msh.addNode(nN);
+ }
+ // create elements
+ for (size_t j = 0; j < teNodes.size() - 1; ++j)
+ {
+ std::vector<Node *> qnodes = {wkNodes[i * teNodes.size() + j], wkNodes[(i + 1) * teNodes.size() + j], wkNodes[(i + 1) * teNodes.size() + j + 1], wkNodes[i * teNodes.size() + j + 1]};
+ _msh.addElement(new Quad4(_msh.getElements().back()->getId() + 1, tagp, qnodes));
}
- Node *nN = new Node(_msh.getNodes().back()->getId() + 1, sdpmVector3d(xF, pos(1), pos(2)));
- _msh.addNode(nN);
- wkNodes.push_back(nN);
- }
- // create wake elements
- for (size_t i = 0; i < wkNodes.size() - 1; ++i)
- {
- std::vector<Node *> qnodes = {teNodes[i], wkNodes[i], wkNodes[i + 1], teNodes[i + 1]};
- _msh.addElement(new Quad4(_msh.getElements().back()->getId() + 1, tagp, qnodes));
}
// duplicate TE nodes
std::map<Node *, Node *> teMap;
@@ -77,17 +90,19 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName, double
teMap[n] = nN;
}
// create wake
- _wake = new Wake(_msh, wkName, _tag);
+ _wake = new Wake(_msh, wkName, _tag, teNodes.size() - 1, ux);
// swap nodes
- for (auto p : _wake->getMap())
+ std::unordered_set<Element *> lwEl; // get set of unique wing TE elements on the pressure side
+ for (auto p : _wake->getElMap())
+ lwEl.insert(p.second.second);
+ for (auto e : lwEl)
{
- Element *lwEl = p.second.second;
- std::vector<Node *> lwNods = lwEl->getNodes();
+ std::vector<Node *> lwNods = e->getNodes();
for (size_t i = 0; i < lwNods.size(); ++i)
{
try
{
- lwEl->setNode(i, teMap.at(lwNods[i]));
+ e->setNode(i, teMap.at(lwNods[i]));
}
catch (const std::out_of_range &)
{
@@ -111,10 +126,21 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName, double
for (auto ten : teNodes)
for (auto wn : _nodes)
if (ten->getCoords() == wn->getCoords() && ten->getId() != wn->getId())
+ {
iTeMap[wn] = ten;
+ break;
+ }
// create wake and print
- _wake = new Wake(_msh, wkName, _tag, iTeMap);
+ _wake = new Wake(_msh, wkName, _tag, teNodes.size() - 1, ux, iTeMap);
std::cout << *_wake << "already existing, nothing done." << std::endl;
+ // check if parameters match existing wake geometry
+ std::vector<Element *> wEls = _wake->getElements();
+ if (std::abs(wEls[0]->getNodes()[1]->getCoords()(0) - wEls[0]->getNodes()[0]->getCoords()(0) - shedLgth / nShedDiv) > 1e-12 || wEls.size() != 10 * nShedDiv * (teNodes.size() - 1))
+ {
+ std::stringstream err;
+ err << "sdpm::Body: parameters \'shedLgth\' and \'nShedDiv\' do not match existing wake geometry!\n";
+ throw std::runtime_error(err.str());
+ }
}
// Associate each node to its adjacent elements
@@ -122,16 +148,44 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName, double
for (auto n : e->getNodes())
_neMap[n].push_back(e);
- // Size load vectors
- nLoads.resize(_nodes.size());
+ // Size vectors
+ _nLoads.resize(_nodes.size());
+ _nLoads1.resize(nM, std::vector<std::vector<sdpmVector3cd>>(nF, std::vector<sdpmVector3cd>(_nodes.size(), sdpmVector3cd::Zero())));
}
Body::~Body()
{
delete _wake;
+ for (auto m : _motions)
+ delete m;
+ _motions.clear();
std::cout << "~sdpm::Body()\n";
}
+/**
+ * @brief Add motion
+ */
+void Body::addMotion(double ux, double amp, double xr)
+{
+ _motions.push_back(new Motion(ux, amp, xr));
+}
+
+/**
+ * @brief Compute the motion-induced velocity on the body surface
+ */
+void Body::computeVelocity(size_t m)
+{
+ std::vector<Element *> elems = getElements();
+ std::vector<sdpmVector3cd> ucS(elems.size()), ucH(elems.size());
+ for (size_t i = 0; i < elems.size(); ++i)
+ {
+ ucS[i] = _motions[m]->computeSteady(*elems[i]);
+ ucH[i] = _motions[m]->computeHarmonic(*elems[i]);
+ }
+ _iVelocityS.push_back(ucS);
+ _iVelocityH.push_back(ucH);
+}
+
/**
* @brief Return a vector of unique nodes belonging to the tag
*/
diff --git a/sdpm/src/sdpmBody.h b/sdpm/src/sdpmBody.h
index c3fa2c8..ae7e924 100644
--- a/sdpm/src/sdpmBody.h
+++ b/sdpm/src/sdpmBody.h
@@ -20,6 +20,7 @@
#include "sdpm.h"
#include "sdpmGroup.h"
#include "sdpmWake.h"
+#include "sdpmMotion.h"
#include <map>
namespace sdpm
@@ -29,20 +30,31 @@ namespace sdpm
* @brief Manage a lifting body immersed in the fluid
* @authors Adrien Crovato
*/
-class SDPM_API Body : public sdpm::Group
+class SDPM_API Body : public Group
{
+ // Geometry
std::vector<Node *> _nodes; ///< nodes of the surface
std::map<Node *, std::vector<Element *>> _neMap; ///< map between nodes and adjacent elements
Wake *_wake; ///< wake attached to the lifting body
+ // Motion
+ std::vector<Motion *> _motions; ///< body motions
+ std::vector<std::vector<sdpmVector3cd>> _iVelocityS; ///< steady motion induced velocity
+ std::vector<std::vector<sdpmVector3cd>> _iVelocityH; ///< harmonic motion induced velocity
+ // Nodal aerodynamic load (normalized by dynamic pressure)
+ std::vector<sdpmVector3d> _nLoads; ///< steady nodal aerodynamic load
+ std::vector<std::vector<std::vector<sdpmVector3cd>>> _nLoads1; ///< unsteady first harmonic nodal aerodynamic load
+ // Coefficients
+ sdpmDouble _cl; ///< steady lift coefficient
+ sdpmDouble _cd; ///< steady drag coefficient
+ sdpmDouble _cs; ///< steady side force coefficient
+ sdpmDouble _cm; ///< steady pitch moment coefficient (positive nose-up)
+ std::vector<std::vector<sdpmComplex>> _cl1; ///< unsteady first harmonic lift coefficient
+ std::vector<std::vector<sdpmComplex>> _cd1; ///< unsteady first harmonic drag coefficient
+ std::vector<std::vector<sdpmComplex>> _cs1; ///< unsteady first harmonic side force coefficient
+ std::vector<std::vector<sdpmComplex>> _cm1; ///< unsteady first harmonic pitch moment coefficient (positive nose-up)
public:
- sdpmDouble cl; ///< lift coefficient
- sdpmDouble cd; ///< drag coefficient
- sdpmDouble cs; ///< sideforce coefficient
- sdpmDouble cm; ///< pitch moment coefficient (positive nose-up)
- std::vector<sdpmVector3d> nLoads; ///< nodal aerodynamic load (normalized by dynamic pressure)
-
- Body(Mesh &msh, std::string const &name, std::string const &teName, double xF);
+ Body(Mesh &msh, std::string const &name, std::string const &teName, double shedLgth, size_t nShedDiv = 1, size_t nM = 0, size_t nF = 0, double ux = 1.);
~Body();
private:
@@ -50,12 +62,38 @@ private:
public:
std::vector<Node *> const &getNodes() const { return _nodes; }
+ std::vector<sdpmVector3d> const &getLoads() const { return _nLoads; }
+ std::vector<sdpmVector3cd> const &getUnsteadyLoads(size_t imd, size_t ifq) const { return _nLoads1[imd][ifq]; }
+ sdpmDouble getLiftCoeff() const { return _cl; }
+ sdpmDouble getDragCoeff() const { return _cd; }
+ sdpmDouble getSideforceCoeff() const { return _cs; }
+ sdpmDouble getMomentCoeff() const { return _cm; }
+ sdpmComplex getUnsteadyLiftCoeff(size_t imd, size_t ifq) const { return _cl1[imd][ifq]; }
+ sdpmComplex getUnsteadyDragCoeff(size_t imd, size_t ifq) const { return _cd1[imd][ifq]; }
+ sdpmComplex getUnsteadySideforceCoeff(size_t imd, size_t ifq) const { return _cs1[imd][ifq]; }
+ sdpmComplex getUnsteadyMomentCoeff(size_t imd, size_t ifq) const { return _cm1[imd][ifq]; }
+ void addMotion(double ux, double amp, double xr);
#ifndef SWIG
std::map<Node *, std::vector<Element *>> const &getMap() const
{
return _neMap;
}
Wake const &getWake() const { return *_wake; }
+ size_t getNumMotions() const { return _motions.size(); }
+ sdpmDouble getAmplitude(size_t imd) const { return _motions[imd]->getAmplitude(); }
+ std::vector<sdpmVector3cd> const &getSteadyVelocity(size_t imd) const { return _iVelocityS[imd]; };
+ std::vector<sdpmVector3cd> const &getHarmonicVelocity(size_t imd) const { return _iVelocityH[imd]; };
+ void setLoads(std::vector<sdpmVector3d> const &nl) { _nLoads = nl; }
+ void setUnsteadyLoads(size_t imd, size_t ifq, std::vector<sdpmVector3cd> const &nl) { _nLoads1[imd][ifq] = nl; }
+ void setLiftCoeff(sdpmDouble cl) { _cl = cl; }
+ void setDragCoeff(sdpmDouble cd) { _cd = cd; }
+ void setSideforceCoeff(sdpmDouble cs) { _cs = cs; }
+ void setMomentCoeff(sdpmDouble cm) { _cm = cm; }
+ void setUnsteadyLiftCoeff(size_t imd, size_t ifq, sdpmComplex cl) { _cl1[imd][ifq] = cl; }
+ void setUnsteadyDragCoeff(size_t imd, size_t ifq, sdpmComplex cd) { _cd1[imd][ifq] = cd; }
+ void setUnsteadySideforceCoeff(size_t imd, size_t ifq, sdpmComplex cs) { _cs1[imd][ifq] = cs; }
+ void setUnsteadyMomentCoeff(size_t imd, size_t ifq, sdpmComplex cm) { _cm1[imd][ifq] = cm; }
+ void computeVelocity(size_t imd);
virtual void write(std::ostream &out) const override;
#endif
};
diff --git a/sdpm/src/sdpmBuilder.cpp b/sdpm/src/sdpmBuilder.cpp
index 367ee52..d720fcc 100644
--- a/sdpm/src/sdpmBuilder.cpp
+++ b/sdpm/src/sdpmBuilder.cpp
@@ -26,10 +26,10 @@ using namespace sdpm;
sdpmMatrix3d Builder::computeTransfoMat(bool symy, Element const *ej)
{
// Compute tangents and normal
- sdpmVector3d u = (ej->getNodes()[0]->getCoords() - ej->getCg()).normalized();
- sdpmVector3d v_ = ej->getNodes()[1]->getCoords() - ej->getCg();
+ sdpmVector3d u = (ej->getCompressibleCoords()[0] - ej->getCompressibleCg()).normalized();
+ sdpmVector3d v_ = ej->getCompressibleCoords()[1] - ej->getCompressibleCg();
sdpmVector3d v = (v_ - (u.dot(v_)) * u).normalized();
- sdpmVector3d n = ej->getNormal();
+ sdpmVector3d n = ej->getCompressibleNormal();
// Assemble into transformation matrix
sdpmMatrix3d g2l;
@@ -49,32 +49,43 @@ sdpmMatrix3d Builder::computeTransfoMat(bool symy, Element const *ej)
*/
sdpmMatrixXd Builder::gatherCoords(bool symy, Element const *ej)
{
- // Get nodes
- std::vector<Node *> nods = ej->getNodes();
- size_t nV = nods.size();
+ // Get compressible coordinates
+ std::vector<sdpmVector3d> ccoords = ej->getCompressibleCoords();
+ size_t nV = ccoords.size();
// Get coordinates
Eigen::MatrixXd coords(3, nV);
for (size_t j = 0; j < nV; ++j)
- coords.col(j) = nods[nV - 1 - j]->getCoords();
+ coords.col(j) = ccoords[nV - 1 - j];
if (symy)
coords.row(1) = -coords.row(1);
return coords;
}
+/**
+ * @brief Compute vector joining body panel ej to body panel ei
+ */
+sdpmVector3d Builder::computeDistVec(bool symy, Element const *ei, Element const *ej)
+{
+ sdpmVector3d cgSrc = ej->getCompressibleCg();
+ sdpmVector3d cgTgt = ei->getCompressibleCg();
+ sdpmVector3d dist = cgTgt - cgSrc;
+ if (symy)
+ dist(1) = cgTgt(1) + cgSrc(1);
+ return dist;
+}
+
/**
* @brief Compute doublet AIC and source AIC from body panel ej to body panel ei
+ *
+ * Katz and Plotkin, Low-speed aerodynamics, 2001
+ * Hess and Smith, Calculation of potential flow about arbitrary bodies, 1967
*/
std::vector<sdpmDouble> Builder::computeAIC(bool symy, Element const *ei, Element const *ej, Eigen::MatrixXd const &sCoords, sdpmMatrix3d const &g2l)
{
// Compute distance between panels and target panel coordinates in local axes
- sdpmVector3d cgSrc = ej->getCg();
- sdpmVector3d cgTgt = ei->getCg();
- sdpmVector3d distance = cgTgt - cgSrc;
- if (symy)
- distance(1) = cgTgt(1) + cgSrc(1);
- distance = g2l * distance;
- sdpmVector3d tCoords = g2l * ei->getCg();
+ sdpmVector3d distance = g2l * computeDistVec(symy, ei, ej);
+ sdpmVector3d tCoords = g2l * ei->getCompressibleCg();
// Compute coefficients
int nV = static_cast<int>(ej->getNodes().size());
diff --git a/sdpm/src/sdpmBuilder.h b/sdpm/src/sdpmBuilder.h
index d569343..a816fd0 100644
--- a/sdpm/src/sdpmBuilder.h
+++ b/sdpm/src/sdpmBuilder.h
@@ -33,6 +33,7 @@ class SDPM_API Builder
public:
static sdpmMatrix3d computeTransfoMat(bool symy, Element const *ej);
static sdpmMatrixXd gatherCoords(bool symy, Element const *ej);
+ static sdpmVector3d computeDistVec(bool symy, Element const *ei, Element const *ej);
static std::vector<sdpmDouble> computeAIC(bool symy, Element const *ei, Element const *ej, sdpmMatrixXd const &trsfCoordSrc, sdpmMatrix3d const &g2l);
};
diff --git a/sdpm/src/sdpmElement.cpp b/sdpm/src/sdpmElement.cpp
index 53b092f..06bd015 100644
--- a/sdpm/src/sdpmElement.cpp
+++ b/sdpm/src/sdpmElement.cpp
@@ -38,20 +38,24 @@ SDPM_API std::ostream &operator<<(std::ostream &out, ElType const &obj)
}
} // namespace sdpm
-Element::Element(size_t n, Tag *tag, std::vector<Node *> &nodes) : sdpmObject(), _id(n), _tag(tag), _nodes(nodes)
+Element::Element(size_t n, Tag *tag, std::vector<Node *> &nodes) : sdpmObject(), _id(n), _tag(tag), _nodes(nodes), _toDate(false)
{
+ _coords.resize(_nodes.size());
_tag->addElements(this);
}
/**
* @brief Update precomputed values
*/
-void Element::update()
+void Element::update(sdpmDouble beta)
{
- // Center of gravity
- _cg = Eigen::Vector3d::Zero();
- for (auto n : _nodes)
- _cg += n->getCoords();
+ // Scale nodes and compute center of gravity
+ _cg = sdpmVector3d::Zero();
+ for (size_t i = 0; i < _nodes.size(); ++i)
+ {
+ _coords[i] = _nodes[i]->getCoords().cwiseQuotient(sdpmVector3d(beta, 1., 1.));
+ _cg += _coords[i];
+ }
_cg /= static_cast<double>(_nodes.size());
}
diff --git a/sdpm/src/sdpmElement.h b/sdpm/src/sdpmElement.h
index 938aad7..a42c410 100644
--- a/sdpm/src/sdpmElement.h
+++ b/sdpm/src/sdpmElement.h
@@ -51,10 +51,15 @@ protected:
size_t _id; ///< element number
Tag *_tag; ///< physical tag to which element belongs
std::vector<Node *> _nodes; ///< nodes of the element
+ bool _toDate; ///< whether the element values are up-to-date
// Geometry
- sdpmVector3d _cg; ///< element center of gravity
- sdpmDouble _vol; ///< element volume
- sdpmVector3d _nrm; ///< unit normal vector
+ std::vector<sdpmVector3d> _coords; ///< scaled (compressible) coordinates
+ sdpmVector3d _cg0; ///< element center of gravity
+ sdpmDouble _srf0; ///< element surface area
+ sdpmVector3d _nrm0; ///< unit normal vector
+ sdpmVector3d _cg; ///< compressible element center of gravity
+ sdpmDouble _srf; ///< compressible element surface area
+ sdpmVector3d _nrm; ///< compressible unit normal vector
public:
#ifndef SWIG
@@ -69,7 +74,7 @@ public:
#ifndef SWIG
// Management
- virtual void update();
+ virtual void update(sdpmDouble beta);
// Setters
void setNode(size_t i, Node *n) { _nodes[i] = n; }
// Getters
@@ -77,15 +82,23 @@ public:
size_t const &getIdRef() const { return _id; }
Tag const *getTag() const { return _tag; }
std::vector<Node *> const &getNodes() const { return _nodes; }
- sdpmVector3d const &getCg() const { return _cg; }
- sdpmDouble getVolume() const { return _vol; }
- sdpmVector3d const &getNormal() const { return _nrm; }
+ sdpmVector3d const &getCg() const { return _cg0; }
+ sdpmDouble getSurfaceArea() const { return _srf0; }
+ sdpmVector3d const &getNormal() const { return _nrm0; }
+ inline std::vector<sdpmVector3d> const &getCompressibleCoords() const;
+ inline sdpmVector3d const &getCompressibleCg() const;
+ inline sdpmDouble getCompressibleSurfaceArea() const;
+ inline sdpmVector3d const &getCompressibleNormal() const;
// Utilities
virtual sdpmVector3d computeGradient(std::vector<sdpmDouble> const &u) const;
virtual void write(std::ostream &out) const override;
#endif
};
+#ifndef SWIG
+#include "sdpmElement.inl"
+#endif
+
} // namespace sdpm
#endif // SDPMELEMENT_H
diff --git a/sdpm/src/sdpmElement.inl b/sdpm/src/sdpmElement.inl
new file mode 100644
index 0000000..33ef73b
--- /dev/null
+++ b/sdpm/src/sdpmElement.inl
@@ -0,0 +1,59 @@
+/*
+ * Copyright 2023 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/**
+ * @brief Return the compressible coordinates
+ */
+inline std::vector<sdpmVector3d> const &Element::getCompressibleCoords() const
+{
+ if (_toDate)
+ return _coords;
+ else
+ throw std::runtime_error("Element::getCompressibleCoords is invalid because Element is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the center of gravity (compressible coordinates)
+ */
+inline sdpmVector3d const &Element::getCompressibleCg() const
+{
+ if (_toDate)
+ return _cg;
+ else
+ throw std::runtime_error("Element::getCompressibleCg is invalid because Element is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the element volume (compressible coordinates)
+ */
+inline sdpmDouble Element::getCompressibleSurfaceArea() const
+{
+ if (_toDate)
+ return _srf;
+ else
+ throw std::runtime_error("Element::getCompressibleSurfaceArea is invalid because Element is not up-to-date!\n");
+}
+
+/**
+ * @brief Return the unit normal vector (compressible coordinates)
+ */
+inline sdpmVector3d const &Element::getCompressibleNormal() const
+{
+ if (_toDate)
+ return _nrm;
+ else
+ throw std::runtime_error("Element::getCompressibleNormal is invalid because Element is not up-to-date!\n");
+}
\ No newline at end of file
diff --git a/sdpm/src/sdpmLine2.cpp b/sdpm/src/sdpmLine2.cpp
index 21f7b6b..efa473a 100644
--- a/sdpm/src/sdpmLine2.cpp
+++ b/sdpm/src/sdpmLine2.cpp
@@ -20,15 +20,21 @@ using namespace sdpm;
Line2::Line2(size_t n, Tag *tag, std::vector<Node *> &nodes) : Element(n, tag, nodes)
{
- update();
+ // Compute unscaled values
+ update(1.);
+ _cg0 = _cg;
+ _srf0 = _srf;
+ _nrm0 = _nrm;
+ _toDate = false;
}
/**
* @brief Update precomputed values and gradients
*/
-void Line2::update()
+void Line2::update(double beta)
{
- Element::update(); // cg
- _vol = (_nodes[1]->getCoords() - _nodes[0]->getCoords()).norm(); // volume
- _nrm = sdpmVector3d::Zero(); // unit normal
+ Element::update(beta); // cg
+ _srf = (_coords[1] - _coords[0]).norm(); // volume
+ _nrm = sdpmVector3d::Zero(); // unit normal
+ _toDate = true;
}
diff --git a/sdpm/src/sdpmLine2.h b/sdpm/src/sdpmLine2.h
index 8eea05e..1022666 100644
--- a/sdpm/src/sdpmLine2.h
+++ b/sdpm/src/sdpmLine2.h
@@ -35,7 +35,7 @@ public:
virtual ElType getType() const override { return ElType::LINE2; }
#ifndef SWIG
- virtual void update() override;
+ virtual void update(double beta) override;
#endif
};
diff --git a/sdpm/src/sdpmMotion.cpp b/sdpm/src/sdpmMotion.cpp
new file mode 100644
index 0000000..1339b92
--- /dev/null
+++ b/sdpm/src/sdpmMotion.cpp
@@ -0,0 +1,43 @@
+/*
+ * Copyright 2023 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#include "sdpm.h"
+#include "sdpmMotion.h"
+#include "sdpmElement.h"
+
+using namespace sdpm;
+
+Motion::Motion(sdpmDouble xvel, sdpmDouble amp, sdpmDouble xref) : sdpmObject(),
+ _xVel(xvel), _amp(amp), _xRef(xref)
+{
+}
+
+/**
+ * @brief Compute the steady part of the motion induced velocity on an element of the body
+ */
+sdpmVector3cd Motion::computeSteady(Element const &e)
+{
+ return sdpmVector3cd(0., 0., _xVel * _amp);
+}
+
+/**
+ * @brief Compute the harmonic part of the motion induced velocity on an element of the body
+ */
+sdpmVector3cd Motion::computeHarmonic(Element const &e)
+{
+ sdpmVector3d cg = e.getCompressibleCg();
+ return sdpmVector3cd(-_amp * cg(2), 0., _amp * (cg(0) - _xRef));
+}
diff --git a/sdpm/src/sdpmMotion.h b/sdpm/src/sdpmMotion.h
new file mode 100644
index 0000000..3973d43
--- /dev/null
+++ b/sdpm/src/sdpmMotion.h
@@ -0,0 +1,48 @@
+/*
+ * Copyright 2023 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef SDPMMOTION_H
+#define SDPMMOTION_H
+
+#include "sdpm.h"
+#include "sdpmObject.h"
+
+namespace sdpm
+{
+
+/**
+ * @brief Hold the parameters defning the motion of a body and compute the induced velocity
+ * @authors Adrien Crovato
+ * @todo support several types of motion and combinations
+ */
+class SDPM_API Motion : public sdpmObject
+{
+ sdpmDouble _xVel;
+ sdpmDouble _amp;
+ sdpmDouble _xRef;
+
+public:
+ Motion(sdpmDouble xvel, sdpmDouble amp, sdpmDouble xref);
+ ~Motion(){};
+
+ sdpmDouble getAmplitude() const { return _amp; }
+ sdpmVector3cd computeSteady(Element const &e);
+ sdpmVector3cd computeHarmonic(Element const &e);
+};
+
+} // namespace sdpm
+
+#endif // SDPMMOTION_H
diff --git a/sdpm/src/sdpmProblem.cpp b/sdpm/src/sdpmProblem.cpp
index ad141a3..50fcb8f 100644
--- a/sdpm/src/sdpmProblem.cpp
+++ b/sdpm/src/sdpmProblem.cpp
@@ -21,14 +21,13 @@
#include <iostream>
using namespace sdpm;
-Problem::Problem(Mesh &msh, double aoa, double aos, double mch,
- double sref, double cref, double xref, double yref, double zref,
- bool ysym) : _msh(msh), _mach(mch), _sRef(sref), _cRef(cref), _ySym(ysym)
+Problem::Problem(Mesh &msh, size_t nmods, std::vector<double> const &fs) : _msh(msh),
+ _aoa(0.), _aos(0.), _mach(0.), _beta(1.0), _sound(0.), _uinf(sdpmVector3d::Zero()),
+ _dDir(sdpmVector3d::Zero()), _sDir(sdpmVector3d::Zero()), _lDir(sdpmVector3d::Zero()),
+ _sRef(0.), _cRef(0.), _xRef(sdpmVector3d::Zero()), _ySym(false),
+ _nModes(nmods), _omega(fs),
+ _hasAero(false), _hasGeo(false)
{
- updateFreestream(aoa, aos);
- _xRef(0) = xref;
- _xRef(1) = yref;
- _xRef(2) = zref;
}
Problem::~Problem()
@@ -45,18 +44,38 @@ void Problem::addBody(Body &b)
}
/**
- * @brief Update the freestream values
+ * @brief Update the freestream aerodynamic values
*/
-void Problem::updateFreestream(double aoa, double aos)
+void Problem::updateFreestream(double aoa, double aos, double mch, double ui)
{
- // Angles and velocity
+ // Angles
_aoa = aoa;
_aos = aos;
- _uinf = sdpmVector3d(cos(_aoa) * cos(_aos), sin(_aos), sin(_aoa) * cos(_aos));
- // Directions
+ // Compressiblity
+ _mach = mch;
+ _beta = std::sqrt(1 - _mach * _mach);
+ _sound = ui / _mach;
+ // Velocity and directions
+ _uinf = ui * sdpmVector3d(cos(_aoa) * cos(_aos), sin(_aos), sin(_aoa) * cos(_aos));
_dDir = sdpmVector3d(cos(_aoa) * cos(_aos), sin(_aos), sin(_aoa) * cos(_aos));
_sDir = sdpmVector3d(-cos(_aoa) * sin(_aos), cos(_aos), -sin(_aoa) * sin(_aos));
- _lDir = sdpmVector3d(-sin(_aoa), 0, cos(_aoa));
+ _lDir = sdpmVector3d(-sin(_aoa), 0., cos(_aoa));
+ // Flag
+ _hasAero = true;
+}
+
+/**
+ * @brief Set the reference geometric values
+ */
+void Problem::setGeometry(double sref, double cref, double xref, double yref, double zref, bool ysym)
+{
+ // Reference values
+ _sRef = sref;
+ _cRef = cref;
+ _xRef = sdpmVector3d(xref, yref, zref);
+ _ySym = ysym;
+ // Flag
+ _hasGeo = true;
}
/**
@@ -64,13 +83,12 @@ void Problem::updateFreestream(double aoa, double aos)
*/
void Problem::updateElements()
{
- // Update surface Jacobian
for (auto b : _bodies)
{
for (auto e : b->getElements())
- e->update();
+ e->update(_beta);
for (auto e : b->getWake().getElements())
- e->update();
+ e->update(_beta);
}
}
@@ -80,10 +98,15 @@ void Problem::updateElements()
void Problem::check() const
{
// Sanity checks
+ if (!_hasAero)
+ throw std::runtime_error("sdpm::Problem: aerodynamic parameters have not been provided!\n");
+ if (!_hasGeo)
+ throw std::runtime_error("sdpm::Problem: geometric parameters have not been provided!\n");
if (_bodies.empty())
throw std::runtime_error("sdpm::Problem: no lifting bodies provided!\n");
- if (_mach != 0)
- throw std::runtime_error("sdpm::Problem: compressible computation feature not implemented yet!\n");
+ for (auto b : _bodies)
+ if (b->getNumMotions() != _nModes)
+ throw std::runtime_error("sdpm::Problem: Body and Problem do not have the same number of modes!\n");
// Surface elements
for (auto b : _bodies)
@@ -102,6 +125,7 @@ void Problem::write(std::ostream &out) const
out << "sdpm::Problem definition"
<< "\n\tAngle of attack: " << _aoa
<< "\n\tAngle of sideslip: " << _aos
+ << "\n\tMach number: " << _mach
<< "\n\tReference area: " << _sRef
<< "\n\tReference length: " << _cRef
<< "\n\tReference point: [" << _xRef.transpose() << "]"
diff --git a/sdpm/src/sdpmProblem.h b/sdpm/src/sdpmProblem.h
index 4a2dac8..bee2d6f 100644
--- a/sdpm/src/sdpmProblem.h
+++ b/sdpm/src/sdpmProblem.h
@@ -35,22 +35,29 @@ class SDPM_API Problem : public sdpmObject
Mesh &_msh; ///< mesh data structure
std::vector<Body *> _bodies; ///< bodies in fluid
// Aerodynamic
- sdpmDouble _aoa; ///< Angle of attack
- sdpmDouble _aos; ///< Angle of sideslip
- double _mach; ///< Mach number
+ sdpmDouble _aoa; ///< Angle of attack
+ sdpmDouble _aos; ///< Angle of sideslip
+ sdpmDouble _mach; ///< Mach number
+ sdpmDouble _beta; ///< Compressibility factor
+ sdpmDouble _sound; ///< Speed of sound
+ sdpmVector3d _uinf; ///< Freestream velocity
+ sdpmVector3d _dDir; ///< Direction of drag force
+ sdpmVector3d _sDir; ///< Direction of side force
+ sdpmVector3d _lDir; ///< Direction of lift force
// Geometry
- double _sRef; ///< Reference surface area
- double _cRef; ///< Reference chord length
- Eigen::Vector3d _xRef; ///< Reference center point for moment computation
- bool _ySym; ///< True if mesh is symmetric with respect to y-plane
- // Freestream
- sdpmVector3d _uinf; ///< freestream velocity
- sdpmVector3d _dDir; ///< direction of drag force
- sdpmVector3d _sDir; ///< direction of side force
- sdpmVector3d _lDir; ///< direction of lift force
+ sdpmDouble _sRef; ///< Reference surface area
+ sdpmDouble _cRef; ///< Reference chord length
+ sdpmVector3d _xRef; ///< Reference center point for moment computation
+ bool _ySym; ///< True if mesh is symmetric with respect to y-plane
+ // Unsteady
+ size_t _nModes; ///< Number of modal motions
+ std::vector<double> _omega; ///< Circular frequencies
+ // Check flags
+ bool _hasAero; ///< True if aerodynamic parameters have been updated
+ bool _hasGeo; ///< True if geometric parameters have been updated
public:
- Problem(Mesh &msh, double aoa, double aos, double mch, double sref, double cref, double xref, double yref, double zref, bool ysym);
+ Problem(Mesh &msh, size_t nmods = 0, std::vector<double> const &fs = {});
~Problem();
#ifndef SWIG
@@ -62,17 +69,23 @@ public:
std::vector<Element *> const &getElements() const { return _msh.getElements(); }
sdpmDouble getAoa() const { return _aoa; }
sdpmDouble getAos() const { return _aos; }
+ sdpmDouble getMach() const { return _mach; }
+ sdpmDouble getCompressibility() const { return _beta; }
+ sdpmDouble getSpeedOfSound() const { return _sound; }
sdpmDouble getRefArea() const { return _sRef; }
sdpmDouble getRefLgth() const { return _cRef; }
sdpmVector3d const &getRefCtr() const { return _xRef; }
sdpmDouble getSym() const { return _ySym; }
+ std::vector<double> const &getFreqs() const { return _omega; }
+ size_t getNumModes() const { return _nModes; }
sdpmVector3d const &getFrstrmVel() const { return _uinf; }
sdpmVector3d const &getDragDir() const { return _dDir; }
sdpmVector3d const &getSideDir() const { return _sDir; }
sdpmVector3d const &getLiftDir() const { return _lDir; }
#endif
void addBody(Body &b);
- void updateFreestream(double aoa, double aos);
+ void updateFreestream(double aoa, double aos, double mch, double ui);
+ void setGeometry(double sref, double cref, double xref, double yref, double zref, bool ysym);
#ifndef SWIG
void check() const;
diff --git a/sdpm/src/sdpmQuad4.cpp b/sdpm/src/sdpmQuad4.cpp
index 8281830..3425915 100644
--- a/sdpm/src/sdpmQuad4.cpp
+++ b/sdpm/src/sdpmQuad4.cpp
@@ -20,19 +20,25 @@ using namespace sdpm;
Quad4::Quad4(size_t n, Tag *tag, std::vector<Node *> &nodes) : Element(n, tag, nodes)
{
- update();
+ // Compute unscaled values
+ update(1.);
+ _cg0 = _cg;
+ _srf0 = _srf;
+ _nrm0 = _nrm;
+ _toDate = false;
}
/**
* @brief Update precomputed values
*/
-void Quad4::update()
+void Quad4::update(double beta)
{
- Element::update(); // cg
- _tgt0 = _nodes[2]->getCoords() - _nodes[0]->getCoords(); // unit tangent (n0->n2)
- _tgt1 = _nodes[3]->getCoords() - _nodes[1]->getCoords(); // unit tangent (n1->n4)
- _vol = 0.5 * _tgt0.cross(_tgt1).norm(); // volume
- _nrm = _tgt0.cross(_tgt1).normalized(); // unit normal
+ Element::update(beta); // cg
+ _tgt0 = _coords[2] - _coords[0]; // unit tangent (n1->n3)
+ _tgt1 = _coords[3] - _coords[1]; // unit tangent (n2->n4)
+ _srf = 0.5 * _tgt0.cross(_tgt1).norm(); // volume
+ _nrm = _tgt0.cross(_tgt1).normalized(); // unit normal
+ _toDate = true;
}
/**
@@ -57,14 +63,14 @@ sdpmVector3d Quad4::computeGradient(std::vector<sdpmDouble> const &u) const
sdpmMatrix3d jac = sdpmMatrix3d::Zero(), ijac = sdpmMatrix3d::Zero();
// 1) compute first 2 rows: dXi, dEta (same as 2D, but with 3D coordinates)
for (size_t i = 0; i < _nodes.size(); ++i)
- jac.block(0, 0, 2, 3) += dsf.col(i) * _nodes[i]->getCoords().transpose();
+ jac.block(0, 0, 2, 3) += dsf.col(i) * _coords[i].transpose();
// 2) compute thrid rows: dZeta = dXi X dEta
jac.row(2) = jac.row(0).cross(jac.row(1));
// 3) inverse Jacobian
ijac = jac.inverse();
// Get solution at element nodes
- Eigen::VectorXd ue(_nodes.size());
+ sdpmVectorXd ue(_nodes.size());
for (size_t i = 0; i < _nodes.size(); ++i)
ue(i) = u[_nodes[i]->getId() - 1];
diff --git a/sdpm/src/sdpmQuad4.h b/sdpm/src/sdpmQuad4.h
index 4d9a948..98ee436 100644
--- a/sdpm/src/sdpmQuad4.h
+++ b/sdpm/src/sdpmQuad4.h
@@ -38,7 +38,7 @@ public:
virtual ElType getType() const override { return ElType::QUAD4; }
#ifndef SWIG
- virtual void update() override;
+ virtual void update(double beta) override;
virtual sdpmVector3d computeGradient(std::vector<sdpmDouble> const &u) const override;
#endif
};
diff --git a/sdpm/src/sdpmSolver.cpp b/sdpm/src/sdpmSolver.cpp
index ca8e0bd..e830ff1 100644
--- a/sdpm/src/sdpmSolver.cpp
+++ b/sdpm/src/sdpmSolver.cpp
@@ -21,6 +21,7 @@
#include "sdpmMesh.h"
#include "sdpmGmshExport.h"
#include "sdpmResults.h"
+#include "sdpm_extras.h"
#include <iostream>
#include <iomanip>
#include <fstream>
@@ -37,27 +38,41 @@ Solver::Solver(Problem &pbl) : _pbl(pbl), _cl(0), _cd(0), _cs(0), _cm(0)
std::cout << "*******************************************************************************" << std::endl;
// Check the problem
- pbl.check();
+ _pbl.check();
- // Count panel, set up global->local id map and resize AIC matrices
+ // Count panel and set up global->local id map
_nP = 0;
int i = 0;
- for (auto body : pbl.getBodies())
+ for (auto body : _pbl.getBodies())
{
_nP += body->getElements().size();
for (auto e : body->getElements())
_rows[e] = i++;
}
- _A = sdpmMatrixXd::Zero(_nP, _nP);
- _B = sdpmMatrixXd::Zero(_nP, _nP);
- // Set up variables
+ // Get sizes
+ size_t nm = _pbl.getNumModes();
+ size_t nf = _pbl.getFreqs().size();
size_t ne = _pbl.getElements().size();
+
+ // Resize AIC matrics
+ _A0 = sdpmMatrixXd(_nP, _nP);
+ _B0 = sdpmMatrixXd(_nP, _nP);
+ _A.resize(nf, sdpmMatrixXcd::Zero(_nP, _nP));
+ _B.resize(nf, sdpmMatrixXcd::Zero(_nP, _nP));
+
+ // Set up variables
_phi.resize(ne, 0.);
_tau.resize(ne, 0.);
_mu.resize(ne, 0.);
_u.resize(ne, sdpmVector3d::Zero());
_cp.resize(ne, 0.);
+ _u1.resize(nm, std::vector<std::vector<sdpmVector3cd>>(nf, std::vector<sdpmVector3cd>(ne, sdpmVector3cd::Zero())));
+ _cp1.resize(nm, std::vector<std::vector<sdpmComplex>>(nf, std::vector<sdpmComplex>(ne, sdpmComplex(0., 0.))));
+ _cl1.resize(nm, std::vector<sdpmComplex>(nf, sdpmComplex(0., 0.)));
+ _cd1.resize(nm, std::vector<sdpmComplex>(nf, sdpmComplex(0., 0.)));
+ _cs1.resize(nm, std::vector<sdpmComplex>(nf, sdpmComplex(0., 0.)));
+ _cm1.resize(nm, std::vector<sdpmComplex>(nf, sdpmComplex(0., 0.)));
// Display configuration
std::cout << "--- Physical groups ---\n";
@@ -67,12 +82,16 @@ Solver::Solver(Problem &pbl) : _pbl(pbl), _cl(0), _cd(0), _cs(0), _cm(0)
std::cout << "--- Freestream conditions ---\n"
<< std::setprecision(5)
<< "Angle of attack: " << _pbl.getAoa() * 180 / 3.14159 << " deg\n"
- << "Angle of sideslip: " << _pbl.getAos() * 180 / 3.14159 << " deg\n";
+ << "Angle of sideslip: " << _pbl.getAos() * 180 / 3.14159 << " deg\n"
+ << "Mach number: " << _pbl.getMach() << "\n"
+ << "Velocity: " << _pbl.getFrstrmVel().norm() << "\n";
std::cout << "--- Reference data ---\n"
- << "Reference area: " << _pbl.getRefArea() << "\n"
- << "Reference length: " << _pbl.getRefLgth() << "\n"
- << "Reference origin: (" << _pbl.getRefCtr().transpose() << ")"
- << "\n"
+ << "Surface area: " << _pbl.getRefArea() << "\n"
+ << "Chord length: " << _pbl.getRefLgth() << "\n"
+ << "Origin: (" << _pbl.getRefCtr().transpose() << ")\n";
+ std::cout << "--- Unsteady parameters ---\n"
+ << "Number of frequencies: " << _pbl.getFreqs().size() << "\n"
+ << "Number of motions: " << _pbl.getNumModes() << "\n"
<< std::endl;
}
@@ -86,34 +105,51 @@ Solver::~Solver()
*/
void Solver::update()
{
- // Update elements
+ // Update elements and get values
_pbl.updateElements();
+ std::vector<double> omega = _pbl.getFreqs();
+ sdpmDouble mach = _pbl.getMach();
+ sdpmDouble sound = _pbl.getSpeedOfSound();
+ sdpmDouble beta = _pbl.getCompressibility();
+ sdpmComplex img(0, 1);
// Compute AIC
std::cout << "Computing AIC matrices... " << std::flush;
std::vector<Body *> bodies = _pbl.getBodies();
- _A.setZero();
- _B.setZero();
+ _A0.setZero();
+ _B0.setZero();
+ std::fill(_A.begin(), _A.end(), sdpmMatrixXcd::Zero(_nP, _nP));
+ std::fill(_B.begin(), _B.end(), sdpmMatrixXcd::Zero(_nP, _nP));
for (auto jbody : bodies)
{
// body
for (auto ej : jbody->getElements())
{
+ sdpmDouble nx = ej->getCompressibleNormal()(0); /// x-component of normal vector
sdpmMatrix3d g2l = Builder::computeTransfoMat(false, ej); // global to local transformation matrix
sdpmMatrixXd sCoords = g2l * Builder::gatherCoords(false, ej); // coordinates of source panel in local axes
for (auto ibody : bodies)
{
for (auto ei : ibody->getElements())
{
- std::vector<sdpmDouble> mutau = Builder::computeAIC(false, ei, ej, sCoords, g2l);
- _A(_rows.at(ei), _rows.at(ej)) = mutau[0];
- _B(_rows.at(ei), _rows.at(ej)) = mutau[1];
+ std::vector<sdpmDouble> mutau = Builder::computeAIC(false, ei, ej, sCoords, g2l); // steady AIC
+ _A0(_rows.at(ei), _rows.at(ej)) = mutau[0];
+ _B0(_rows.at(ei), _rows.at(ej)) = mutau[1];
+ for (size_t fk = 0; fk < omega.size(); ++fk)
+ {
+ sdpmDouble om = omega[fk] / sound / beta;
+ sdpmVector3d r = Builder::computeDistVec(false, ei, ej);
+ sdpmComplex exp = std::exp(-img * om * (r.norm() - mach * r(0)));
+ _A[fk](_rows.at(ei), _rows.at(ej)) = (mutau[0] * (1. + img * om * r.norm()) - mutau[1] * img * om * mach * nx) * exp;
+ _B[fk](_rows.at(ei), _rows.at(ej)) = mutau[1] * exp;
+ }
}
}
}
// wake
Wake const &wake = jbody->getWake();
- std::map<Element *, std::pair<Element *, Element *>> wwMap = wake.getMap();
+ std::map<Element *, std::pair<Element *, Element *>> wwMap = wake.getElMap();
+ std::map<Element *, double> lagMap = wake.getLagMap();
for (auto ew : wake.getElements())
{
sdpmMatrix3d g2l = Builder::computeTransfoMat(false, ew); // global to local transformation matrix
@@ -122,9 +158,18 @@ void Solver::update()
{
for (auto ei : ibody->getElements())
{
- std::vector<sdpmDouble> mutau = Builder::computeAIC(false, ei, ew, sCoords, g2l);
- _A(_rows.at(ei), _rows.at(wwMap.at(ew).first)) += mutau[0];
- _A(_rows.at(ei), _rows.at(wwMap.at(ew).second)) -= mutau[0];
+ std::vector<sdpmDouble> mutau = Builder::computeAIC(false, ei, ew, sCoords, g2l); // steady AIC
+ _A0(_rows.at(ei), _rows.at(wwMap.at(ew).first)) += mutau[0];
+ _A0(_rows.at(ei), _rows.at(wwMap.at(ew).second)) -= mutau[0];
+ for (size_t fk = 0; fk < omega.size(); ++fk)
+ {
+ sdpmDouble om = omega[fk] / sound / beta;
+ sdpmVector3d r = Builder::computeDistVec(false, ei, ew);
+ sdpmComplex exp = std::exp(-img * om * (r.norm() - mach * r(0)));
+ sdpmComplex waktrm = mutau[0] * (1. + img * om * r.norm()) * exp * std::exp(-img * omega[fk] * lagMap.at(ew));
+ _A[fk](_rows.at(ei), _rows.at(wwMap.at(ew).first)) += waktrm;
+ _A[fk](_rows.at(ei), _rows.at(wwMap.at(ew).second)) -= waktrm;
+ }
}
}
}
@@ -137,21 +182,31 @@ void Solver::update()
// body
for (auto ej : jbody->getElements())
{
+ sdpmDouble nx = ej->getCompressibleNormal()(0); /// x-component of normal vector
sdpmMatrix3d g2l = Builder::computeTransfoMat(true, ej); // global to local transformation matrix
sdpmMatrixXd sCoords = g2l * Builder::gatherCoords(true, ej); // coordinates of source panel in local axes
for (auto ibody : bodies)
{
for (auto ei : ibody->getElements())
{
- std::vector<sdpmDouble> mutau = Builder::computeAIC(true, ei, ej, sCoords, g2l);
- _A(_rows.at(ei), _rows.at(ej)) -= mutau[0];
- _B(_rows.at(ei), _rows.at(ej)) -= mutau[1];
+ std::vector<sdpmDouble> mutau = Builder::computeAIC(true, ei, ej, sCoords, g2l); // steady AIC
+ _A0(_rows.at(ei), _rows.at(ej)) -= mutau[0];
+ _B0(_rows.at(ei), _rows.at(ej)) -= mutau[1];
+ for (size_t fk = 0; fk < omega.size(); ++fk)
+ {
+ sdpmDouble om = omega[fk] / sound / beta;
+ sdpmVector3d r = Builder::computeDistVec(true, ei, ej);
+ sdpmComplex exp = std::exp(-img * om * (r.norm() - mach * r(0)));
+ _A[fk](_rows.at(ei), _rows.at(ej)) -= (mutau[0] * (1. + img * om * r.norm()) - mutau[1] * img * om * mach * nx) * exp;
+ _B[fk](_rows.at(ei), _rows.at(ej)) -= mutau[1] * exp;
+ }
}
}
}
// wake
Wake const &wake = jbody->getWake();
- std::map<Element *, std::pair<Element *, Element *>> wwMap = wake.getMap();
+ std::map<Element *, std::pair<Element *, Element *>> wwMap = wake.getElMap();
+ std::map<Element *, double> lagMap = wake.getLagMap();
for (auto ew : wake.getElements())
{
sdpmMatrix3d g2l = Builder::computeTransfoMat(true, ew); // global to local transformation matrix
@@ -160,9 +215,18 @@ void Solver::update()
{
for (auto ei : ibody->getElements())
{
- std::vector<sdpmDouble> mutau = Builder::computeAIC(true, ei, ew, sCoords, g2l);
- _A(_rows.at(ei), _rows.at(wwMap.at(ew).first)) -= mutau[0];
- _A(_rows.at(ei), _rows.at(wwMap.at(ew).second)) += mutau[0];
+ std::vector<sdpmDouble> mutau = Builder::computeAIC(true, ei, ew, sCoords, g2l); // steady AIC
+ _A0(_rows.at(ei), _rows.at(wwMap.at(ew).first)) -= mutau[0];
+ _A0(_rows.at(ei), _rows.at(wwMap.at(ew).second)) += mutau[0];
+ for (size_t fk = 0; fk < omega.size(); ++fk)
+ {
+ sdpmDouble om = omega[fk] / sound / beta;
+ sdpmVector3d r = Builder::computeDistVec(true, ei, ew);
+ sdpmComplex exp = std::exp(-img * om * (r.norm() - mach * r(0)));
+ sdpmComplex waktrm = mutau[0] * (1. + img * om * r.norm()) * exp * std::exp(-img * omega[fk] * lagMap.at(ew));
+ _A[fk](_rows.at(ei), _rows.at(wwMap.at(ew).first)) -= waktrm;
+ _A[fk](_rows.at(ei), _rows.at(wwMap.at(ew).second)) += waktrm;
+ }
}
}
}
@@ -172,7 +236,7 @@ void Solver::update()
}
/**
- * @brief Solve the potential equation
+ * @brief Solve the steady potential equation
*/
void Solver::run()
{
@@ -182,21 +246,61 @@ void Solver::run()
sdpmVectorXd etau(_nP);
for (auto body : _pbl.getBodies())
for (auto e : body->getElements())
- etau(_rows[e]) = ui.dot(e->getNormal());
+ etau(_rows[e]) = ui.dot(e->getCompressibleNormal().cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.)));
std::cout << "done." << std::endl;
std::cout << "solving for doublets... " << std::flush;
sdpmVectorXd emu(_nP);
- emu = Eigen::PartialPivLU<sdpmMatrixXd>(_A).solve(_B * etau);
+ emu = Eigen::PartialPivLU<sdpmMatrixXd>(_A0).solve(_B0 * etau);
std::cout << "done." << std::endl;
std::cout << "computing flow and loads on bodies... " << std::flush;
- computeFlow(etau, emu);
- computeLoad();
+ post(etau, emu);
std::cout << "done." << std::endl
<< std::endl;
}
+/**
+ * @brief Solve the unsteady potential equation
+ * @todo remove output
+ */
+void Solver::runUnsteady()
+{
+ std::cout << "Running unsteady SDPM solver... " << std::endl;
+ for (size_t imd = 0; imd < _pbl.getNumModes(); ++imd)
+ {
+ // Compute motion-induced velocity
+ for (auto body : _pbl.getBodies())
+ body->computeVelocity(imd);
+ // Solve
+ for (size_t ifq = 0; ifq < _pbl.getFreqs().size(); ++ifq)
+ {
+ std::cout << "Mode #" << imd << ", frequency #" << ifq << std::endl;
+ std::cout << "computing sources... " << std::flush;
+ sdpmVectorXcd etau(_nP);
+ for (auto body : _pbl.getBodies())
+ {
+ std::vector<sdpmVector3cd> uc0 = body->getSteadyVelocity(imd);
+ std::vector<sdpmVector3cd> uc1 = body->getHarmonicVelocity(imd);
+ std::vector<Element *> elems = body->getElements();
+ for (size_t i = 0; i < elems.size(); ++i)
+ etau(_rows[elems[i]]) = (uc0[i] + uc1[i] * sdpmComplex(0., 1.) * _pbl.getFreqs()[ifq]).conjugate().dot(elems[i]->getCompressibleNormal().cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.)));
+ }
+ std::cout << "done." << std::endl;
+
+ std::cout << "solving for doublets... " << std::flush;
+ sdpmVectorXcd emu(_nP);
+ emu = Eigen::PartialPivLU<sdpmMatrixXcd>(_A[ifq]).solve(_B[ifq] * etau);
+ std::cout << "done." << std::endl;
+
+ std::cout << "computing flow and loads on bodies... " << std::flush;
+ postUnsteady(imd, ifq, etau, emu);
+ std::cout << "done." << std::endl;
+ }
+ }
+ std::cout << std::endl;
+}
+
/**
* @brief Write the results
*/
@@ -221,10 +325,10 @@ void Solver::save(GmshExport const &writer, std::string const &suffix)
{
outfile << "# Body - " << b->getName() << " (" << b->getElements().size() << " elements)" << std::endl;
outfile << std::fixed << std::setprecision(5)
- << "Cl = " << b->cl << "\n"
- << "Cd = " << b->cd << "\n"
- << "Cs = " << b->cs << "\n"
- << "Cm = " << b->cm << std::endl
+ << "Cl = " << b->getLiftCoeff() << "\n"
+ << "Cd = " << b->getDragCoeff() << "\n"
+ << "Cs = " << b->getSideforceCoeff() << "\n"
+ << "Cm = " << b->getMomentCoeff() << std::endl
<< std::endl;
}
outfile.close();
@@ -232,12 +336,12 @@ void Solver::save(GmshExport const &writer, std::string const &suffix)
}
/**
- * @brief Compute flow solution on bodies
+ * @brief Compute steady flow functionals and loads on bodies
*/
-void Solver::computeFlow(sdpmVectorXd const &etau, sdpmVectorXd const &emu)
+void Solver::post(sdpmVectorXd const &etau, sdpmVectorXd const &emu)
{
// Store source and doublet strength and compute perturbation potential on elements
- sdpmVectorXd ephi = -_A * emu - _B * etau;
+ sdpmVectorXd ephi = -_A0 * emu - _B0 * etau;
std::vector<Body *> bodies = _pbl.getBodies();
for (auto body : bodies)
{
@@ -255,35 +359,27 @@ void Solver::computeFlow(sdpmVectorXd const &etau, sdpmVectorXd const &emu)
std::vector<sdpmDouble> muNodes(_pbl.getNodes().size(), 0);
for (auto b : bodies)
{
- // associate the value to its element
- std::map<Element *, sdpmDouble> mapMu;
- for (auto e : b->getElements())
- mapMu[e] = _mu[e->getId() - 1];
- // average at nodes
std::map<Node *, std::vector<Element *>> neMap = b->getMap();
for (auto n : b->getNodes())
for (auto e : neMap.at(n))
- muNodes[n->getId() - 1] += mapMu.at(e) / neMap.at(n).size();
+ muNodes[n->getId() - 1] += _mu[e->getId() - 1] / neMap.at(n).size();
}
// Compute surface flow functionals
sdpmVector3d ui = _pbl.getFrstrmVel(); // freestream velocity
+ sdpmDouble mi = _pbl.getMach(); // freestream speed of sound
for (auto body : bodies)
{
for (auto e : body->getElements())
{
- sdpmVector3d u = ui - _tau[e->getId() - 1] * e->getNormal() - e->computeGradient(muNodes);
- _u[e->getId() - 1] = u; // velocity
- _cp[e->getId() - 1] = 1 - u.dot(u); // pressure coefficient
+ sdpmVector3d u = ui - (_tau[e->getId() - 1] * e->getCompressibleNormal() + e->computeGradient(muNodes)).cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.));
+ sdpmDouble cpl = 1 - u.dot(u) / ui.dot(ui);
+ _u[e->getId() - 1] = u; // velocity
+ _cp[e->getId() - 1] = cpl + 0.5 * 0.5 * cpl * cpl * mi * mi; // nonlinear pressure coefficient
}
}
-}
-/**
- * @brief Compute aerodynamic load coefficients on bodies
- */
-void Solver::computeLoad()
-{
+ // Compute aerodynamic loads
sdpmDouble sRef = _pbl.getRefArea(); // reference surface area
sdpmDouble cRef = _pbl.getRefLgth(); // reference length
sdpmVector3d xRef = _pbl.getRefCtr(); // reference center
@@ -296,36 +392,140 @@ void Solver::computeLoad()
_cm = 0;
for (auto b : _pbl.getBodies())
{
- // Reset
- std::fill(b->nLoads.begin(), b->nLoads.end(), sdpmVector3d::Zero());
-
// Compute pressure loads coefficient (normalized by freestream dynamic pressure) on each element and average at node
std::vector<Node *> nods = b->getNodes();
std::map<Node *, std::vector<Element *>> neMap = b->getMap();
+ std::vector<sdpmVector3d> nLoads(nods.size(), sdpmVector3d::Zero());
for (size_t i = 0; i < nods.size(); ++i)
for (auto e : neMap.at(nods[i]))
- b->nLoads[i] += -_cp[e->getId() - 1] * e->getVolume() * e->getNormal() / e->getNodes().size();
+ nLoads[i] += -_cp[e->getId() - 1] * e->getSurfaceArea() * e->getNormal() / e->getNodes().size();
+ b->setLoads(nLoads);
// Compute integrated aerodynamic load coefficients (normalized by freestream dynamic pressure and reference area)
- // force coefficients along x and vertical (y in 2D, z in 3D) directions and pitching moment (along -z in 2D, y in 3D)
- b->cm = 0;
+ // force coefficients along x, y and z directions and pitching moment around y axis
+ sdpmDouble bcm = 0;
sdpmVector3d cf(0, 0, 0);
for (size_t i = 0; i < nods.size(); ++i)
{
sdpmVector3d l = nods[i]->getCoords() - xRef; // lever arm
- sdpmVector3d cfi = b->nLoads[i]; // normalized force
+ sdpmVector3d cfi = nLoads[i]; // normalized force
cf += cfi;
- b->cm += (l(2) * cfi(0) - l(0) * cfi(2)) / (sRef * cRef); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
+ bcm += (l(2) * cfi(0) - l(0) * cfi(2)) / (sRef * cRef); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
}
+ b->setMomentCoeff(bcm);
// rotate to flow direction
- b->cd = cf.dot(dragDir) / sRef;
- b->cs = cf.dot(sideDir) / sRef;
- b->cl = cf.dot(liftDir) / sRef;
+ b->setDragCoeff(cf.dot(dragDir) / sRef);
+ b->setSideforceCoeff(cf.dot(sideDir) / sRef);
+ b->setLiftCoeff(cf.dot(liftDir) / sRef);
+ // compute total
+ _cl += b->getLiftCoeff();
+ _cd += b->getDragCoeff();
+ _cs += b->getSideforceCoeff();
+ _cm += bcm;
+ }
+}
+
+/**
+ * @brief Compute unsteady flow functionals and loads on bodies
+ */
+void Solver::postUnsteady(size_t imd, size_t ifq, sdpmVectorXcd const &etau, sdpmVectorXcd const &emu)
+{
+ // Transfer doublet at nodes
+ std::vector<sdpmDouble> muNodesRe(_pbl.getNodes().size(), 0), muNodesIm(_pbl.getNodes().size(), 0);
+ for (auto b : _pbl.getBodies())
+ {
+ std::map<Node *, std::vector<Element *>> neMap = b->getMap();
+ for (auto n : b->getNodes())
+ {
+ for (auto e : neMap.at(n))
+ {
+ muNodesRe[n->getId() - 1] += emu(_rows.at(e)).real() / neMap.at(n).size();
+ muNodesIm[n->getId() - 1] += emu(_rows.at(e)).imag() / neMap.at(n).size();
+ }
+ }
+ }
+
+ // Compute flow and loads
+ sdpmVector3d ui = _pbl.getFrstrmVel();
+ double omega = _pbl.getFreqs()[ifq];
+ sdpmDouble sRef = _pbl.getRefArea();
+ sdpmDouble cRef = _pbl.getRefLgth();
+ sdpmVector3d xRef = _pbl.getRefCtr();
+ double ca = cos(_pbl.getAos());
+ double sa = sin(_pbl.getAos());
+ _cl1[imd][ifq] = sdpmComplex(0., 0.);
+ _cd1[imd][ifq] = sdpmComplex(0., 0.);
+ _cs1[imd][ifq] = sdpmComplex(0., 0.);
+ _cm1[imd][ifq] = sdpmComplex(0., 0.);
+ for (auto b : _pbl.getBodies())
+ {
+ // compute velocity
+ std::vector<sdpmVector3cd> uc0 = b->getSteadyVelocity(imd);
+ std::vector<sdpmVector3cd> uc1 = b->getHarmonicVelocity(imd);
+ std::vector<Element *> elems = b->getElements();
+ for (size_t i = 0; i < elems.size(); ++i)
+ _u1[imd][ifq][elems[i]->getId() - 1] = uc0[i] + uc1[i] * sdpmComplex(0., 1.) * omega - (etau[_rows[elems[i]]] * elems[i]->getCompressibleNormal() + (elems[i]->computeGradient(muNodesRe) + sdpmComplex(0, 1) * elems[i]->computeGradient(muNodesIm))).cwiseQuotient(sdpmVector3d(_pbl.getCompressibility(), 1., 1.));
+ // compute pressure
+ std::vector<std::vector<sdpmComplex>> cp012(elems.size(), std::vector<sdpmComplex>(3));
+ for (size_t i = 0; i < elems.size(); ++i)
+ {
+ // si = phi_t + 1/2 conv(u_k,u_k)_k - 1/2 u_inf^2
+ sdpmVectorXcd si(5);
+ si << 0., std::conj(-sdpmComplex(0, 1) * omega * emu(_rows[elems[i]])), -0.5 * ui.dot(ui), -sdpmComplex(0, 1) * omega * emu(_rows[elems[i]]), 0.;
+ for (size_t j = 0; j < 3; ++j)
+ {
+ sdpmVector3cd utrm(_u1[imd][ifq][elems[i]->getId() - 1].conjugate()(j), _u[elems[i]->getId() - 1](j), _u1[imd][ifq][elems[i]->getId() - 1](j));
+ si += 0.5 * convolve(utrm, utrm);
+ }
+ // ecp = -2 / u_inf^2 * si + 1/u_inf^2/a_inf^2 * conv(si, si)
+ sdpmVectorXcd ecp(9);
+ ecp << 0., 0., -2 * si / ui.dot(ui), 0., 0.;
+ ecp += convolve(si, si) / (_pbl.getSpeedOfSound() * _pbl.getSpeedOfSound() * ui.dot(ui));
+ for (size_t u = 0; u < 3; ++u)
+ cp012[i][u] = ecp(4 + u);
+ _cp1[imd][ifq][elems[i]->getId() - 1] = cp012[i][1];
+ }
+
+ // compute loads on nodes
+ std::vector<Node *> nods = b->getNodes();
+ std::map<Node *, std::vector<Element *>> neMap = b->getMap();
+ std::vector<std::vector<sdpmVector3cd>> nLoads012(nods.size(), std::vector<sdpmVector3cd>(3, sdpmVector3cd::Zero()));
+ for (size_t i = 0; i < nods.size(); ++i)
+ for (auto e : neMap.at(nods[i]))
+ for (size_t u = 0; u < 3; ++u)
+ nLoads012[i][u] += -cp012[_rows.at(e)][u] * e->getSurfaceArea() * e->getNormal() / e->getNodes().size();
+ std::vector<sdpmVector3cd> nLoads1(nods.size());
+ for (size_t i = 0; i < nods.size(); ++i)
+ nLoads1[i] = nLoads012[i][1];
+ b->setUnsteadyLoads(imd, ifq, nLoads1);
+
+ // compute moment coefficient and total load on body
+ sdpmComplex bcm1(0., 0.);
+ std::vector<sdpmVector3cd> cf(3, sdpmVector3cd::Zero());
+ for (size_t i = 0; i < nods.size(); ++i)
+ {
+ sdpmVector3d l = nods[i]->getCoords() - xRef;
+ for (size_t u = 0; u < 3; ++u)
+ cf[u] += nLoads012[i][u];
+ bcm1 += (l(2) * nLoads012[i][1](0) - l(0) * nLoads012[i][1](2)) / (sRef * cRef);
+ }
+ b->setUnsteadyMomentCoeff(imd, ifq, bcm1);
+
+ // compute load coefficients, only valid for small AoAs
+ sdpmVector3cd alpha(b->getAmplitude(imd), _pbl.getAoa(), b->getAmplitude(imd));
+ sdpmVectorXcd cfx(5);
+ cfx << cf[2].conjugate()(0), cf[1].conjugate()(0), cf[0](0), cf[1](0), cf[2](0);
+ sdpmVectorXcd cfz(5);
+ cfz << cf[2].conjugate()(2), cf[1].conjugate()(2), cf[0](2), cf[1](2), cf[2](2);
+ b->setUnsteadyDragCoeff(imd, ifq, (cf[1](0) * ca + cf[1](1) * sa + convolve(cfz, alpha)(4) * ca) / sRef);
+ b->setUnsteadySideforceCoeff(imd, ifq, (-cf[1](0) * sa + cf[1](1) * ca - convolve(cfz, alpha)(4) * sa) / sRef);
+ b->setUnsteadyLiftCoeff(imd, ifq, (-convolve(cfx, alpha)(4) + cf[1](2)) / sRef);
+
// compute total
- _cl += b->cl;
- _cd += b->cd;
- _cs += b->cs;
- _cm += b->cm;
+ _cd1[imd][ifq] += b->getUnsteadyDragCoeff(imd, ifq);
+ _cs1[imd][ifq] += b->getUnsteadySideforceCoeff(imd, ifq);
+ _cl1[imd][ifq] += b->getUnsteadyLiftCoeff(imd, ifq);
+ _cm1[imd][ifq] += b->getUnsteadyMomentCoeff(imd, ifq);
}
}
diff --git a/sdpm/src/sdpmSolver.h b/sdpm/src/sdpmSolver.h
index 16e0ebc..a412044 100644
--- a/sdpm/src/sdpmSolver.h
+++ b/sdpm/src/sdpmSolver.h
@@ -36,19 +36,27 @@ class SDPM_API Solver : public sdpmObject
// AICs
size_t _nP; ///< number of body panels
std::map<Element *, int> _rows; ///< map linking an element to a matrix cell
- sdpmMatrixXd _A; ///< body-body doublet matrix
- sdpmMatrixXd _B; ///< body-body source matrix
+ sdpmMatrixXd _A0; ///< zero-frequency body-body doublet matrix
+ sdpmMatrixXd _B0; ///< zero-frequency body-body source matrix
+ std::vector<sdpmMatrixXcd> _A; ///< unsteady body-body doublet matrix
+ std::vector<sdpmMatrixXcd> _B; ///< unsteady body-body source matrix
// Solution
- std::vector<sdpmDouble> _phi; ///< perturbation potential
- std::vector<sdpmDouble> _tau; ///< source strength
- std::vector<sdpmDouble> _mu; ///< doublet strength
- std::vector<sdpmVector3d> _u; ///< velocity
- std::vector<sdpmDouble> _cp; ///< pressure coefficient
+ std::vector<sdpmDouble> _phi; ///< perturbation potential
+ std::vector<sdpmDouble> _tau; ///< source strength
+ std::vector<sdpmDouble> _mu; ///< doublet strength
+ std::vector<sdpmVector3d> _u; ///< steady velocity
+ std::vector<sdpmDouble> _cp; ///< steady pressure coefficient
+ std::vector<std::vector<std::vector<sdpmVector3cd>>> _u1; ///< unsteady first harmonic velocity
+ std::vector<std::vector<std::vector<sdpmComplex>>> _cp1; ///< unsteady first harmonic pressure coefficient
// Aerodynamic loads
- sdpmDouble _cl; ///< lift coefficient
- sdpmDouble _cd; ///< drag coefficient
- sdpmDouble _cs; ///< sideforce coefficient
- sdpmDouble _cm; ///< pitch moment coefficient (positive nose-up)
+ sdpmDouble _cl; ///< steady lift coefficient
+ sdpmDouble _cd; ///< steady drag coefficient
+ sdpmDouble _cs; ///< steady sideforce coefficient
+ sdpmDouble _cm; ///< steady pitch moment coefficient (positive nose-up)
+ std::vector<std::vector<sdpmComplex>> _cl1; ///< unsteady first harmonic lift coefficient
+ std::vector<std::vector<sdpmComplex>> _cd1; ///< unsteady first harmonic drag coefficient
+ std::vector<std::vector<sdpmComplex>> _cs1; ///< unsteady first harmonic side force coefficient
+ std::vector<std::vector<sdpmComplex>> _cm1; ///< unsteady first harmonic pitch moment coefficient (positive nose-up)
public:
Solver(Problem &pbl);
@@ -58,9 +66,14 @@ public:
sdpmDouble getDragCoeff() const { return _cd; }
sdpmDouble getSideforceCoeff() const { return _cs; }
sdpmDouble getMomentCoeff() const { return _cm; }
+ sdpmComplex getUnsteadyLiftCoeff(size_t m, size_t k) const { return _cl1[m][k]; }
+ sdpmComplex getUnsteadyDragCoeff(size_t m, size_t k) const { return _cd1[m][k]; }
+ sdpmComplex getUnsteadySideforceCoeff(size_t m, size_t k) const { return _cs1[m][k]; }
+ sdpmComplex getUnsteadyMomentCoeff(size_t m, size_t k) const { return _cm1[m][k]; }
void update();
void run();
+ void runUnsteady();
void save(GmshExport const &writer, std::string const &suffix = "");
#ifndef SWIG
@@ -68,8 +81,8 @@ public:
#endif
private:
- void computeFlow(sdpmVectorXd const &etau, sdpmVectorXd const &emu);
- void computeLoad();
+ void post(sdpmVectorXd const &etau, sdpmVectorXd const &emu);
+ void postUnsteady(size_t m, size_t k, sdpmVectorXcd const &etau, sdpmVectorXcd const &emu);
};
} // namespace sdpm
diff --git a/sdpm/src/sdpmWake.cpp b/sdpm/src/sdpmWake.cpp
index 27a78aa..cfc9358 100644
--- a/sdpm/src/sdpmWake.cpp
+++ b/sdpm/src/sdpmWake.cpp
@@ -22,122 +22,97 @@
#include <unordered_set>
using namespace sdpm;
-#include <iostream>
-
-/**
- * @brief Create the wake in the case where the elements on the pressure side still have their original (suction side) trailing edge nodes
- */
-Wake::Wake(Mesh &msh, std::string const &name, Tag *const &wingTag) : Group(msh, name)
+Wake::Wake(Mesh &msh, std::string const &name, Tag *const &wingTag, size_t nTe, double ux, std::map<Node *, Node *> const &iTeMap) : Group(msh, name)
{
// Create TE edges
- createEdges();
+ std::vector<Element *> elems = _tag->getElements();
+ std::unordered_set<Edge *, EquEdge, EquEdge> teEdges;
+ for (size_t i = 0; i < nTe; ++i)
+ {
+ std::vector<Node *> nods = elems[i]->getNodes();
+ Edge *ed = new Edge({nods[0], nods[3]});
+ ed->el0 = elems[i];
+ teEdges.insert(ed);
+ }
- // Find wing elements having an edge on the TE
- for (auto e : wingTag->getElements())
+ // Find wing elements having an edge on the TE...
+ // ... the elements on the pressure side still have their original (suction side) trailing edge nodes
+ if (iTeMap.empty())
{
- std::vector<Node *> nods = e->getNodes();
- size_t idx = 0;
- for (size_t i = 0; i < nods.size(); ++i)
+ for (auto e : wingTag->getElements())
{
- // build wing edge and check if wing edge is a TE edge
- std::vector<Node *> eN(2);
- for (size_t j = 0; j < 2; ++j)
- eN[j] = nods[(idx + j) % nods.size()];
- Edge ed(eN);
- auto it = _teEdges.find(&ed);
- if (it != _teEdges.end())
+ std::vector<Node *> nods = e->getNodes();
+ size_t idx = 0;
+ for (size_t i = 0; i < nods.size(); ++i)
{
- // check normal direction
- if (e->getNormal().dot((*it)->el0->getNormal()) > 0)
- (*it)->el1 = e; // suction side
- else
- (*it)->el2 = e; // pressure side
- break;
+ // build wing edge and check if wing edge is a TE edge
+ std::vector<Node *> eN(2);
+ for (size_t j = 0; j < 2; ++j)
+ eN[j] = nods[(idx + j) % nods.size()];
+ Edge ed(eN);
+ auto it = teEdges.find(&ed);
+ if (it != teEdges.end())
+ {
+ // check normal direction
+ if (e->getNormal().dot((*it)->el0->getNormal()) > 0)
+ (*it)->el1 = e; // suction side
+ else
+ (*it)->el2 = e; // pressure side
+ break;
+ }
+ idx++;
}
- idx++;
}
}
-
- // Map wake to wing elements
- createMap(wingTag);
-}
-
-/**
- * @brief Create the wake in the case where the elements on the pressure side already have their correct (pressure side) trailing edge nodes
- */
-Wake::Wake(Mesh &msh, std::string const &name, Tag *const &wingTag, std::map<Node *, Node *> const &iTeMap) : Group(msh, name)
-{
- // Create TE edges
- createEdges();
-
- // Find wing elements having an edge on the TE
- for (auto e : wingTag->getElements())
+ // ... the elements on the pressure side already have their correct (pressure side) trailing edge nodes
+ else
{
- std::vector<Node *> nods = e->getNodes();
- size_t idx = 0;
- for (size_t i = 0; i < nods.size(); ++i, ++idx)
+ for (auto e : wingTag->getElements())
{
- // build wing edge and check if wing edge is a TE edge
- std::vector<Node *> eN(2);
- for (size_t j = 0; j < 2; ++j)
- eN[j] = nods[(idx + j) % nods.size()];
- Edge ed(eN);
- auto it = _teEdges.find(&ed);
- // try with actual nodes
- if (it != _teEdges.end())
- {
- (*it)->el1 = e; // suction side
- break;
- }
- // otherwise try again, but by replacing the nodes by their corresponding upper nodes
- else
+ std::vector<Node *> nods = e->getNodes();
+ size_t idx = 0;
+ for (size_t i = 0; i < nods.size(); ++i, ++idx)
{
// build wing edge and check if wing edge is a TE edge
std::vector<Node *> eN(2);
- try
- {
- for (size_t j = 0; j < 2; ++j)
- eN[j] = iTeMap.at(nods[(idx + j) % nods.size()]);
- }
- catch (const std::out_of_range &)
- {
- continue;
- }
+ for (size_t j = 0; j < 2; ++j)
+ eN[j] = nods[(idx + j) % nods.size()];
Edge ed(eN);
- auto it = _teEdges.find(&ed);
- if (it != _teEdges.end())
+ auto it = teEdges.find(&ed);
+ // try with actual nodes
+ if (it != teEdges.end())
{
- (*it)->el2 = e; // pressure side
+ (*it)->el1 = e; // suction side
break;
}
+ // otherwise try again, but by replacing the nodes by their corresponding upper nodes
+ else
+ {
+ // build wing edge and check if wing edge is a TE edge
+ std::vector<Node *> eN(2);
+ try
+ {
+ for (size_t j = 0; j < 2; ++j)
+ eN[j] = iTeMap.at(nods[(idx + j) % nods.size()]);
+ }
+ catch (const std::out_of_range &)
+ {
+ continue;
+ }
+ Edge ed(eN);
+ auto it = teEdges.find(&ed);
+ if (it != teEdges.end())
+ {
+ (*it)->el2 = e; // pressure side
+ break;
+ }
+ }
}
}
}
- // Map wake to wing elements
- createMap(wingTag);
-}
-
-/**
- * @brief Create set of TE edges
- */
-void Wake::createEdges()
-{
- for (auto e : _tag->getElements())
- {
- std::vector<Node *> nods = e->getNodes();
- Edge *ed = new Edge({nods[0], nods[3]});
- ed->el0 = e;
- _teEdges.insert(ed);
- }
-}
-
-/**
- * @brief Create a map between the wake and wing elements
- */
-void Wake::createMap(Tag const *wingTag)
-{
- for (auto ed : _teEdges)
+ // Map wake elements on the TE to wing elements
+ for (auto ed : teEdges)
{
if (ed->el1 == NULL || ed->el2 == NULL)
{
@@ -148,6 +123,19 @@ void Wake::createMap(Tag const *wingTag)
_wwMap[ed->el0] = std::pair<Element *, Element *>(ed->el1, ed->el2);
delete ed;
}
+ // Map the rest of the wake and fill time lag map
+ sdpmDouble dt = (elems[0]->getNodes()[1]->getCoords()(0) - elems[0]->getNodes()[0]->getCoords()(0)) / ux;
+ size_t nw = elems.size() / nTe; // number of shed wake rows
+ for (size_t j = 0; j < nTe; ++j)
+ {
+ _lagMap[elems[j]] = dt;
+ std::pair<Element *, Element *> p = _wwMap.at(elems[j]);
+ for (size_t i = 1; i < nw; ++i)
+ {
+ _wwMap[elems[i * nTe + j]] = p;
+ _lagMap[elems[i * nTe + j]] = (i + 1) * dt;
+ }
+ }
}
void Wake::write(std::ostream &out) const
diff --git a/sdpm/src/sdpmWake.h b/sdpm/src/sdpmWake.h
index 338892b..61966bc 100644
--- a/sdpm/src/sdpmWake.h
+++ b/sdpm/src/sdpmWake.h
@@ -20,8 +20,6 @@
#include "sdpm.h"
#include "sdpmGroup.h"
#include "sdpmTag.h"
-#include "sdpmEdge.h"
-#include <unordered_set>
#include <map>
namespace sdpm
@@ -33,24 +31,19 @@ namespace sdpm
*/
class SDPM_API Wake : public Group
{
- std::unordered_set<Edge *, EquEdge, EquEdge> _teEdges; ///< trailing edge edges
std::map<Element *, std::pair<Element *, Element *>> _wwMap; ///< wing-wake map
+ std::map<Element *, double> _lagMap; ///< time lag map
public:
- Wake(Mesh &msh, std::string const &name, Tag *const &wingTag);
- Wake(Mesh &msh, std::string const &name, Tag *const &wingTag, std::map<Node *, Node *> const &iTeMap);
+ Wake(Mesh &msh, std::string const &name, Tag *const &wingTag, size_t nTe, double ux, std::map<Node *, Node *> const &iTeMap = {});
~Wake() {}
-private:
- void createEdges();
- void createMap(Tag const *tag);
-
-public:
#ifndef SWIG
- std::map<Element *, std::pair<Element *, Element *>> const &getMap() const
+ std::map<Element *, std::pair<Element *, Element *>> const &getElMap() const
{
return _wwMap;
}
+ std::map<Element *, double> const &getLagMap() const { return _lagMap; }
virtual void write(std::ostream &out) const override;
#endif
};
diff --git a/sdpm/src/sdpm_extras.h b/sdpm/src/sdpm_extras.h
index 981800a..ef7b770 100644
--- a/sdpm/src/sdpm_extras.h
+++ b/sdpm/src/sdpm_extras.h
@@ -18,12 +18,8 @@
#define SDPM_EXTRAS_H
#include "sdpm.h"
-
-#include <iostream>
-#include <fstream>
+#include <ostream>
#include <vector>
-#include <string>
-#include <Eigen/Sparse>
namespace std
{
@@ -48,19 +44,21 @@ ostream &operator<<(ostream &out, vector<T> const &v)
namespace sdpm
{
/**
- * @brief save CSR matrix to a file
- *
- * matlab usage: load K.txt; K=spconvert(K)
+ * @brief convolve two vectors
*/
-SDPM_API void tomatlab(std::string const &fname, Eigen::SparseMatrix<double> const &A)
+SDPM_API sdpmVectorXcd convolve(sdpmVectorXcd const &u, sdpmVectorXcd const &v)
{
- std::ofstream f(fname.c_str());
-
- for (int k = 0; k < A.outerSize(); ++k)
- for (Eigen::SparseMatrix<double>::InnerIterator it(A, k); it; ++it)
- f << it.row() + 1 << ' ' << it.col() + 1 << ' ' << it.value() << '\n';
-
- f.close();
+ int m = static_cast<int>(u.size());
+ int n = static_cast<int>(v.size());
+ sdpmVectorXcd w = sdpmVectorXcd::Zero(n + m - 1);
+ for (int i = 0; i < w.size(); ++i)
+ {
+ int jmin = (i - n + 1 >= 0) ? i - n + 1 : 0;
+ int jmax = (i - m + 1 < 0) ? i : m - 1;
+ for (int j = jmin; j <= jmax; ++j)
+ w(i) += u(j) * v(i - j);
+ }
+ return w;
}
} // namespace sdpm
diff --git a/sdpm/tests/lann.py b/sdpm/tests/lann.py
new file mode 100644
index 0000000..031b183
--- /dev/null
+++ b/sdpm/tests/lann.py
@@ -0,0 +1,86 @@
+# -*- coding: utf-8 -*-
+
+# Copyright 2023 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+# LANN wing
+# Adrien Crovato
+
+import sdpm
+from sdpm.gmsh import MeshLoader
+from sdpm.utils import build_fpath
+from sdpm.testing import *
+import math
+import cmath
+
+def main():
+ # geometry parameters
+ wnc = 20 # number of chordwise panels
+ pars = {'nC': wnc, 'bC': 0.25, 'rS': 1}
+ # flow and time parameters
+ minf = 0.621 # freestream Mach number
+ aoa = 0.59 * math.pi / 180
+ omega = 2 * math.pi * 24 # circular frequency
+ amp = 0.5 * 0.25 * math.pi / 180 # semi pitching amplitude
+ uinf = 0.5 * 0.361 * omega / 0.133; # freestream velocity
+
+ # start timer
+ tms = sdpm.Timers()
+ tms['total'].start()
+ # mesh
+ tms['msh'].start()
+ msh = MeshLoader(build_fpath('models/lann.geo', __file__, 1)).run(pars)
+ tms['msh'].stop()
+ # problem
+ tms['pbl'].start()
+ pbl = sdpm.Problem(msh, 1, [omega])
+ pbl.setGeometry(0.253, 0.361, 0.224, 0., 0., True)
+ pbl.updateFreestream(aoa, 0., minf, uinf)
+ wing = sdpm.Body(msh, 'wing', 'te', 0.361, wnc, 1, 1, uinf*math.cos(aoa))
+ wing.addMotion(uinf*math.cos(aoa), amp, 0.224/math.sqrt(1-minf*minf))
+ pbl.addBody(wing)
+ tms['pbl'].stop()
+ ## solver
+ tms['sol'].start()
+ sol = sdpm.Solver(pbl)
+ sol.update()
+ sol.run()
+ sol.runUnsteady()
+ sol.save(sdpm.GmshExport(msh))
+ tms['sol'].stop()
+ # end timer
+ tms['total'].stop()
+ print(tms)
+
+ # test
+ tests = CTests()
+ tests.add(CTest('CL0', sol.getLiftCoeff(), 0.258, 5e-2)) # MATLAB: 0.265
+ tests.add(CTest('CD0', sol.getDragCoeff(), -0.0029, 5e-2)) # MATLAB: -0.0033
+ tests.add(CTest('CS0', sol.getSideforceCoeff(), 0.008, 5e-2))
+ tests.add(CTest('CM0', sol.getMomentCoeff(), -0.079, 5e-2))
+ tests.add(CTest('Abs(CL1)', abs(sol.getUnsteadyLiftCoeff(0, 0)) / amp / math.pi, 1.54, 5e-2)) # MATLAB: 1.62
+ tests.add(CTest('Ang(CL1)', cmath.phase(sol.getUnsteadyLiftCoeff(0, 0)) * 180 / math.pi, 0.66, 5e-2)) # MATLAB: 2.6
+ tests.add(CTest('Abs(CD1)', abs(sol.getUnsteadyDragCoeff(0, 0)) / amp / math.pi, 0.038, 5e-2)) # MATLAB: 0.041
+ tests.add(CTest('Ang(CD1)', cmath.phase(sol.getUnsteadyDragCoeff(0, 0)) * 180 / math.pi, 14.6, 5e-2)) # MATLAB: 16.8
+ tests.add(CTest('Abs(CS1)', abs(sol.getUnsteadySideforceCoeff(0, 0)) / amp / math.pi, 0.037, 5e-2))
+ tests.add(CTest('Ang(CS1)', cmath.phase(sol.getUnsteadySideforceCoeff(0, 0)) * 180 / math.pi, -10.0, 5e-2))
+ tests.add(CTest('Abs(CM1)', abs(sol.getUnsteadyMomentCoeff(0, 0)) / amp / math.pi, 0.33, 5e-2))
+ tests.add(CTest('Ang(CM1)', cmath.phase(sol.getUnsteadyMomentCoeff(0, 0)) * 180 / math.pi, -170.5, 5e-2))
+ tests.run()
+
+ # eof
+ print('')
+
+if __name__ == '__main__':
+ main()
diff --git a/sdpm/tests/m6.py b/sdpm/tests/m6.py
index eff0489..528d4b8 100644
--- a/sdpm/tests/m6.py
+++ b/sdpm/tests/m6.py
@@ -25,9 +25,10 @@ import math
def main():
# geometry parameters
- pars = {'wNc': 100, 'wNs': 10, 'bC': 0.25, 'pS': 1.0}
+ pars = {'nC': 100, 'nS': 10, 'bC': 0.25, 'pS': 1.0}
# flow parameters
aoa = 3.06 * math.pi / 180
+ minf = 0.
# start timer
tms = sdpm.Timers()
@@ -38,8 +39,10 @@ def main():
tms['msh'].stop()
# problem
tms['pbl'].start()
- pbl = sdpm.Problem(msh, aoa, 0, 0, 0.7528, 0.64, 0, 0, 0, True)
- wing = sdpm.Body(msh, 'wing', 'wTe', 10)
+ pbl = sdpm.Problem(msh)
+ pbl.setGeometry(0.7528, 0.64, 0., 0., 0., True)
+ pbl.updateFreestream(aoa, 0., minf, 1.)
+ wing = sdpm.Body(msh, 'wing', 'te', 0.805)
pbl.addBody(wing)
tms['pbl'].stop()
## solver
diff --git a/sdpm/tests/n12.py b/sdpm/tests/n12.py
index d6fc9c4..25bb02d 100644
--- a/sdpm/tests/n12.py
+++ b/sdpm/tests/n12.py
@@ -26,7 +26,7 @@ import math
def main():
# geometry parameters
sym = 0
- span = 5 # true span if sym=0, half span if sym=1
+ span = 5 if sym==0 else 2.5 # span
pars = {'symY': sym, 'wSpn': span, 'wNc': 100, 'wNs': 10}
# flow parameters
aoa = 3 * math.pi / 180
@@ -40,8 +40,10 @@ def main():
tms['msh'].stop()
# problem
tms['pbl'].start()
- pbl = sdpm.Problem(msh, aoa, 0, 0, span, 1, 0, 0, 0, bool(sym))
- wing = sdpm.Body(msh, 'wing', 'wTe', 10)
+ pbl = sdpm.Problem(msh)
+ pbl.setGeometry(span, 1., 0., 0., 0., bool(sym))
+ pbl.updateFreestream(aoa, 0., 0., 1.)
+ wing = sdpm.Body(msh, 'wing', 'wTe', 1.)
pbl.addBody(wing)
tms['pbl'].stop()
## solver
diff --git a/sdpm/tests/n12tail.py b/sdpm/tests/n12tail.py
index aee9152..d76c7ed 100644
--- a/sdpm/tests/n12tail.py
+++ b/sdpm/tests/n12tail.py
@@ -26,9 +26,9 @@ import math
def main():
# geometry parameters
sym = 0
- span = 5 # wing (true span if sym=0, half span if sym=1)
- tspan = 3 # tail (true span if sym=0, half span if sym=1)
- pars = {'symY': sym, 'wSpn': span, 'wNc': 100, 'wNs': 10, 'tSpn': tspan, 'tNc': 100, 'tNs': 6}
+ wspan = 5 if sym==0 else 2.5 # wing span
+ tspan = 3 if sym==0 else 1.5 # tail span
+ pars = {'symY': sym, 'wSpn': wspan, 'wNc': 100, 'wNs': 10, 'tSpn': tspan, 'tNc': 100, 'tNs': 6}
# flow parameters
aoa = 3 * math.pi / 180
@@ -41,9 +41,11 @@ def main():
tms['msh'].stop()
# problem
tms['pbl'].start()
- pbl = sdpm.Problem(msh, aoa, 0, 0, span, 1, 0, 0, 0, bool(sym))
- wing = sdpm.Body(msh, 'wing', 'wTe', 10)
- tail = sdpm.Body(msh, 'tail', 'tTe', 15)
+ pbl = sdpm.Problem(msh)
+ pbl.setGeometry(wspan, 1., 0., 0., 0., bool(sym))
+ pbl.updateFreestream(aoa, 0., 0., 1.)
+ wing = sdpm.Body(msh, 'wing', 'wTe', 1.0)
+ tail = sdpm.Body(msh, 'tail', 'tTe', 0.5)
pbl.addBody(wing)
pbl.addBody(tail)
tms['pbl'].stop()
@@ -60,14 +62,14 @@ def main():
# test
tests = CTests()
- tests.add(CTest('Wing CL', wing.cl, 0.230, 5e-2))
- tests.add(CTest('Wing CD', wing.cd, 0.003, 5e-2))
- tests.add(CTest('Wing CS', wing.cs, 0.000, 5e-2))
- tests.add(CTest('Wing CM', wing.cm, -0.056, 5e-2))
- tests.add(CTest('Tail CL', tail.cl, 0.045, 5e-2))
- tests.add(CTest('Tail CD', tail.cd, 0.001, 5e-2))
- tests.add(CTest('Tail CS', tail.cs, 0.000, 5e-2))
- tests.add(CTest('Tail CM', tail.cm, -0.229, 5e-2))
+ tests.add(CTest('Wing CL', wing.getLiftCoeff(), 0.230, 5e-2))
+ tests.add(CTest('Wing CD', wing.getDragCoeff(), 0.003, 5e-2))
+ tests.add(CTest('Wing CS', wing.getSideforceCoeff(), 0.000, 5e-2))
+ tests.add(CTest('Wing CM', wing.getMomentCoeff(), -0.056, 5e-2))
+ tests.add(CTest('Tail CL', tail.getLiftCoeff(), 0.045, 5e-2))
+ tests.add(CTest('Tail CD', tail.getDragCoeff(), 0.001, 5e-2))
+ tests.add(CTest('Tail CS', tail.getSideforceCoeff(), 0.000, 5e-2))
+ tests.add(CTest('Tail CM', tail.getMomentCoeff(), -0.229, 5e-2))
tests.run()
# eof
--
GitLab