diff --git a/dart/models/oneraM6_VII.geo b/dart/models/oneraM6_VII.geo
new file mode 100644
index 0000000000000000000000000000000000000000..6c14627beea6c5588ab593509a2b95cf405377c6
--- /dev/null
+++ b/dart/models/oneraM6_VII.geo
@@ -0,0 +1,590 @@
+/* Onera M6 wing */
+// Initially generated by unsRgridWingGen.m
+// For gmsh 4, use line 425 instead of line 426 (Bezier <- BSpline)
+
+// Parameters
+// domain and mesh
+DefineConstant[ lgt = { 3.0, Name "Channel half length" } ];
+DefineConstant[ wdt = { 2.0, Name "Channel width" } ];
+DefineConstant[ hgt = { 3.0, Name "Channel half height" } ];
+DefineConstant[ msLeRt = { 0.01, Name "Root leading edge mesh size" } ];
+DefineConstant[ msTeRt = { 0.01, Name "Root trailing edge mesh size" } ];
+DefineConstant[ msLeTp = { 0.01, Name "Tip leading edge mesh size" } ];
+DefineConstant[ msTeTp = { 0.01, Name "Tip trailing edge mesh size" } ];
+DefineConstant[ msF = { 1.0, Name "Farfield mesh size" } ];
+
+
+//// GEOMETRY
+
+
+/// Points
+
+// Airfoil 1: oneraM6, 143 points
+Point(1) = {0.805900,0.000000,0.000000,msTeRt};
+Point(2) = {0.804129,0.000000,0.000232};
+Point(3) = {0.802038,0.000000,0.000505};
+Point(4) = {0.799569,0.000000,0.000828};
+Point(5) = {0.796652,0.000000,0.001210};
+Point(6) = {0.793208,0.000000,0.001661};
+Point(7) = {0.789139,0.000000,0.002193};
+Point(8) = {0.784331,0.000000,0.002822};
+Point(9) = {0.778649,0.000000,0.003565};
+Point(10) = {0.771932,0.000000,0.004444};
+Point(11) = {0.763989,0.000000,0.005483};
+Point(12) = {0.754592,0.000000,0.006710};
+Point(13) = {0.743470,0.000000,0.008151};
+Point(14) = {0.730299,0.000000,0.009828,msTeRt};
+Point(15) = {0.714691,0.000000,0.011690};
+Point(16) = {0.696182,0.000000,0.013774};
+Point(17) = {0.677546,0.000000,0.015783};
+Point(18) = {0.658809,0.000000,0.017717};
+Point(19) = {0.639970,0.000000,0.019586};
+Point(20) = {0.621026,0.000000,0.021397};
+Point(21) = {0.601980,0.000000,0.023159};
+Point(22) = {0.582826,0.000000,0.024875};
+Point(23) = {0.563568,0.000000,0.026547};
+Point(24) = {0.544202,0.000000,0.028170};
+Point(25) = {0.524729,0.000000,0.029737};
+Point(26) = {0.505147,0.000000,0.031238};
+Point(27) = {0.485455,0.000000,0.032658};
+Point(28) = {0.465652,0.000000,0.033984};
+Point(29) = {0.445738,0.000000,0.035197};
+Point(30) = {0.425711,0.000000,0.036282};
+Point(31) = {0.405570,0.000000,0.037225};
+Point(32) = {0.385317,0.000000,0.038011};
+Point(33) = {0.364947,0.000000,0.038631};
+Point(34) = {0.344460,0.000000,0.039073};
+Point(35) = {0.323856,0.000000,0.039344};
+Point(36) = {0.303135,0.000000,0.039432};
+Point(37) = {0.282293,0.000000,0.039343};
+Point(38) = {0.261331,0.000000,0.039078};
+Point(39) = {0.240248,0.000000,0.038642,1.5*msLeRt};
+Point(40) = {0.219042,0.000000,0.038037};
+Point(41) = {0.197712,0.000000,0.037261};
+Point(42) = {0.176258,0.000000,0.036306};
+Point(43) = {0.154679,0.000000,0.035154};
+Point(44) = {0.132972,0.000000,0.033774};
+Point(45) = {0.111131,0.000000,0.032117};
+Point(46) = {0.092662,0.000000,0.030457};
+Point(47) = {0.077073,0.000000,0.028830};
+Point(48) = {0.063937,0.000000,0.027269};
+Point(49) = {0.052891,0.000000,0.025800};
+Point(50) = {0.043619,0.000000,0.024441};
+Point(51) = {0.035851,0.000000,0.023203};
+Point(52) = {0.029356,0.000000,0.022027};
+Point(53) = {0.023936,0.000000,0.020812};
+Point(54) = {0.019428,0.000000,0.019503};
+Point(55) = {0.015684,0.000000,0.018096};
+Point(56) = {0.012586,0.000000,0.016619};
+Point(57) = {0.010032,0.000000,0.015114};
+Point(58) = {0.007931,0.000000,0.013618};
+Point(59) = {0.006214,0.000000,0.012165};
+Point(60) = {0.004815,0.000000,0.010776};
+Point(61) = {0.003683,0.000000,0.009463};
+Point(62) = {0.002775,0.000000,0.008233};
+Point(63) = {0.002050,0.000000,0.007089};
+Point(64) = {0.001480,0.000000,0.006027};
+Point(65) = {0.001037,0.000000,0.005045};
+Point(66) = {0.000698,0.000000,0.004138};
+Point(67) = {0.000444,0.000000,0.003301};
+Point(68) = {0.000260,0.000000,0.002529};
+Point(69) = {0.000135,0.000000,0.001818};
+Point(70) = {0.000056,0.000000,0.001162};
+Point(71) = {0.000013,0.000000,0.000557};
+Point(72) = {0.000000,0.000000,0.000000,msLeRt};
+Point(73) = {0.000013,0.000000,-0.000557};
+Point(74) = {0.000056,0.000000,-0.001162};
+Point(75) = {0.000135,0.000000,-0.001818};
+Point(76) = {0.000260,0.000000,-0.002529};
+Point(77) = {0.000444,0.000000,-0.003301};
+Point(78) = {0.000698,0.000000,-0.004138};
+Point(79) = {0.001037,0.000000,-0.005045};
+Point(80) = {0.001480,0.000000,-0.006027};
+Point(81) = {0.002050,0.000000,-0.007089};
+Point(82) = {0.002775,0.000000,-0.008233};
+Point(83) = {0.003683,0.000000,-0.009463};
+Point(84) = {0.004815,0.000000,-0.010776};
+Point(85) = {0.006214,0.000000,-0.012165};
+Point(86) = {0.007931,0.000000,-0.013618};
+Point(87) = {0.010032,0.000000,-0.015114};
+Point(88) = {0.012586,0.000000,-0.016619};
+Point(89) = {0.015684,0.000000,-0.018096};
+Point(90) = {0.019428,0.000000,-0.019503};
+Point(91) = {0.023936,0.000000,-0.020812};
+Point(92) = {0.029356,0.000000,-0.022027};
+Point(93) = {0.035851,0.000000,-0.023203};
+Point(94) = {0.043619,0.000000,-0.024441};
+Point(95) = {0.052891,0.000000,-0.025800};
+Point(96) = {0.063937,0.000000,-0.027269};
+Point(97) = {0.077073,0.000000,-0.028830};
+Point(98) = {0.092662,0.000000,-0.030457};
+Point(99) = {0.111131,0.000000,-0.032117};
+Point(100) = {0.132972,0.000000,-0.033774};
+Point(101) = {0.154679,0.000000,-0.035154};
+Point(102) = {0.176258,0.000000,-0.036306};
+Point(103) = {0.197712,0.000000,-0.037261};
+Point(104) = {0.219042,0.000000,-0.038037};
+Point(105) = {0.240248,0.000000,-0.038642,1.5*msLeRt};
+Point(106) = {0.261331,0.000000,-0.039078};
+Point(107) = {0.282293,0.000000,-0.039343};
+Point(108) = {0.303135,0.000000,-0.039432};
+Point(109) = {0.323856,0.000000,-0.039344};
+Point(110) = {0.344460,0.000000,-0.039073};
+Point(111) = {0.364947,0.000000,-0.038631};
+Point(112) = {0.385317,0.000000,-0.038011};
+Point(113) = {0.405570,0.000000,-0.037225};
+Point(114) = {0.425711,0.000000,-0.036282};
+Point(115) = {0.445738,0.000000,-0.035197};
+Point(116) = {0.465652,0.000000,-0.033984};
+Point(117) = {0.485455,0.000000,-0.032658};
+Point(118) = {0.505147,0.000000,-0.031238};
+Point(119) = {0.524729,0.000000,-0.029737};
+Point(120) = {0.544202,0.000000,-0.028170};
+Point(121) = {0.563568,0.000000,-0.026547};
+Point(122) = {0.582826,0.000000,-0.024875};
+Point(123) = {0.601980,0.000000,-0.023159};
+Point(124) = {0.621026,0.000000,-0.021397};
+Point(125) = {0.639970,0.000000,-0.019586};
+Point(126) = {0.658809,0.000000,-0.017717};
+Point(127) = {0.677546,0.000000,-0.015783};
+Point(128) = {0.696182,0.000000,-0.013774};
+Point(129) = {0.714691,0.000000,-0.011690};
+Point(130) = {0.730299,0.000000,-0.009828,msTeRt};
+Point(131) = {0.743470,0.000000,-0.008151};
+Point(132) = {0.754592,0.000000,-0.006710};
+Point(133) = {0.763989,0.000000,-0.005483};
+Point(134) = {0.771932,0.000000,-0.004444};
+Point(135) = {0.778649,0.000000,-0.003565};
+Point(136) = {0.784331,0.000000,-0.002822};
+Point(137) = {0.789139,0.000000,-0.002193};
+Point(138) = {0.793208,0.000000,-0.001661};
+Point(139) = {0.796652,0.000000,-0.001210};
+Point(140) = {0.799569,0.000000,-0.000828};
+Point(141) = {0.802038,0.000000,-0.000505};
+Point(142) = {0.804129,0.000000,-0.000232,msTeRt};
+
+// Airfoil 2: oneraM6, 143 points
+Point(144) = {1.143427,1.196000,0.000000,msTeTp};
+Point(145) = {1.142432,1.196000,0.000130};
+Point(146) = {1.141256,1.196000,0.000284};
+Point(147) = {1.139869,1.196000,0.000466};
+Point(148) = {1.138230,1.196000,0.000680};
+Point(149) = {1.136294,1.196000,0.000933};
+Point(150) = {1.134007,1.196000,0.001232};
+Point(151) = {1.131305,1.196000,0.001586};
+Point(152) = {1.128112,1.196000,0.002004};
+Point(153) = {1.124337,1.196000,0.002498};
+Point(154) = {1.119873,1.196000,0.003082};
+Point(155) = {1.114592,1.196000,0.003771};
+Point(156) = {1.108341,1.196000,0.004581};
+Point(157) = {1.100939,1.196000,0.005523,1.5*msTeTp};
+Point(158) = {1.092167,1.196000,0.006570};
+Point(159) = {1.081765,1.196000,0.007741};
+Point(160) = {1.071292,1.196000,0.008870};
+Point(161) = {1.060762,1.196000,0.009957};
+Point(162) = {1.050174,1.196000,0.011007};
+Point(163) = {1.039528,1.196000,0.012025};
+Point(164) = {1.028824,1.196000,0.013015};
+Point(165) = {1.018059,1.196000,0.013980};
+Point(166) = {1.007236,1.196000,0.014919};
+Point(167) = {0.996353,1.196000,0.015831};
+Point(168) = {0.985409,1.196000,0.016712};
+Point(169) = {0.974403,1.196000,0.017556};
+Point(170) = {0.963336,1.196000,0.018354};
+Point(171) = {0.952208,1.196000,0.019099};
+Point(172) = {0.941016,1.196000,0.019781};
+Point(173) = {0.929760,1.196000,0.020391};
+Point(174) = {0.918441,1.196000,0.020920};
+Point(175) = {0.907059,1.196000,0.021362};
+Point(176) = {0.895611,1.196000,0.021711};
+Point(177) = {0.884097,1.196000,0.021959};
+Point(178) = {0.872518,1.196000,0.022111};
+Point(179) = {0.860873,1.196000,0.022161};
+Point(180) = {0.849160,1.196000,0.022111};
+Point(181) = {0.837379,1.196000,0.021962};
+Point(182) = {0.825530,1.196000,0.021717,1.5*msLeTp};
+Point(183) = {0.813612,1.196000,0.021377};
+Point(184) = {0.801625,1.196000,0.020941};
+Point(185) = {0.789568,1.196000,0.020404};
+Point(186) = {0.777440,1.196000,0.019757};
+Point(187) = {0.765241,1.196000,0.018981};
+Point(188) = {0.752966,1.196000,0.018050};
+Point(189) = {0.742587,1.196000,0.017117};
+Point(190) = {0.733826,1.196000,0.016203};
+Point(191) = {0.726444,1.196000,0.015325};
+Point(192) = {0.720236,1.196000,0.014500};
+Point(193) = {0.715025,1.196000,0.013736};
+Point(194) = {0.710659,1.196000,0.013040};
+Point(195) = {0.707009,1.196000,0.012379};
+Point(196) = {0.703963,1.196000,0.011696};
+Point(197) = {0.701429,1.196000,0.010961};
+Point(198) = {0.699325,1.196000,0.010170};
+Point(199) = {0.697584,1.196000,0.009340};
+Point(200) = {0.696149,1.196000,0.008494};
+Point(201) = {0.694968,1.196000,0.007654};
+Point(202) = {0.694003,1.196000,0.006837};
+Point(203) = {0.693217,1.196000,0.006056};
+Point(204) = {0.692581,1.196000,0.005318};
+Point(205) = {0.692070,1.196000,0.004627};
+Point(206) = {0.691663,1.196000,0.003984};
+Point(207) = {0.691343,1.196000,0.003387};
+Point(208) = {0.691094,1.196000,0.002835};
+Point(209) = {0.690903,1.196000,0.002325};
+Point(210) = {0.690760,1.196000,0.001855};
+Point(211) = {0.690657,1.196000,0.001421};
+Point(212) = {0.690587,1.196000,0.001022};
+Point(213) = {0.690542,1.196000,0.000653};
+Point(214) = {0.690518,1.196000,0.000313};
+Point(215) = {0.690511,1.196000,0.000000,msLeTp};
+Point(216) = {0.690518,1.196000,-0.000313};
+Point(217) = {0.690542,1.196000,-0.000653};
+Point(218) = {0.690587,1.196000,-0.001022};
+Point(219) = {0.690657,1.196000,-0.001421};
+Point(220) = {0.690760,1.196000,-0.001855};
+Point(221) = {0.690903,1.196000,-0.002325};
+Point(222) = {0.691094,1.196000,-0.002835};
+Point(223) = {0.691343,1.196000,-0.003387};
+Point(224) = {0.691663,1.196000,-0.003984};
+Point(225) = {0.692070,1.196000,-0.004627};
+Point(226) = {0.692581,1.196000,-0.005318};
+Point(227) = {0.693217,1.196000,-0.006056};
+Point(228) = {0.694003,1.196000,-0.006837};
+Point(229) = {0.694968,1.196000,-0.007654};
+Point(230) = {0.696149,1.196000,-0.008494};
+Point(231) = {0.697584,1.196000,-0.009340};
+Point(232) = {0.699325,1.196000,-0.010170};
+Point(233) = {0.701429,1.196000,-0.010961};
+Point(234) = {0.703963,1.196000,-0.011696};
+Point(235) = {0.707009,1.196000,-0.012379};
+Point(236) = {0.710659,1.196000,-0.013040};
+Point(237) = {0.715025,1.196000,-0.013736};
+Point(238) = {0.720236,1.196000,-0.014500};
+Point(239) = {0.726444,1.196000,-0.015325};
+Point(240) = {0.733826,1.196000,-0.016203};
+Point(241) = {0.742587,1.196000,-0.017117};
+Point(242) = {0.752966,1.196000,-0.018050};
+Point(243) = {0.765241,1.196000,-0.018981};
+Point(244) = {0.777440,1.196000,-0.019757};
+Point(245) = {0.789568,1.196000,-0.020404};
+Point(246) = {0.801625,1.196000,-0.020941};
+Point(247) = {0.813612,1.196000,-0.021377};
+Point(248) = {0.825530,1.196000,-0.021717,1.5*msLeTp};
+Point(249) = {0.837379,1.196000,-0.021962};
+Point(250) = {0.849160,1.196000,-0.022111};
+Point(251) = {0.860873,1.196000,-0.022161};
+Point(252) = {0.872518,1.196000,-0.022111};
+Point(253) = {0.884097,1.196000,-0.021959};
+Point(254) = {0.895611,1.196000,-0.021711};
+Point(255) = {0.907059,1.196000,-0.021362};
+Point(256) = {0.918441,1.196000,-0.020920};
+Point(257) = {0.929760,1.196000,-0.020391};
+Point(258) = {0.941016,1.196000,-0.019781};
+Point(259) = {0.952208,1.196000,-0.019099};
+Point(260) = {0.963336,1.196000,-0.018354};
+Point(261) = {0.974403,1.196000,-0.017556};
+Point(262) = {0.985409,1.196000,-0.016712};
+Point(263) = {0.996353,1.196000,-0.015831};
+Point(264) = {1.007236,1.196000,-0.014919};
+Point(265) = {1.018059,1.196000,-0.013980};
+Point(266) = {1.028824,1.196000,-0.013015};
+Point(267) = {1.039528,1.196000,-0.012025};
+Point(268) = {1.050174,1.196000,-0.011007};
+Point(269) = {1.060762,1.196000,-0.009957};
+Point(270) = {1.071292,1.196000,-0.008870};
+Point(271) = {1.081765,1.196000,-0.007741};
+Point(272) = {1.092167,1.196000,-0.006570};
+Point(273) = {1.100939,1.196000,-0.005523,1.5*msTeTp};
+Point(274) = {1.108341,1.196000,-0.004581};
+Point(275) = {1.114592,1.196000,-0.003771};
+Point(276) = {1.119873,1.196000,-0.003082};
+Point(277) = {1.124337,1.196000,-0.002498};
+Point(278) = {1.128112,1.196000,-0.002004};
+Point(279) = {1.131305,1.196000,-0.001586};
+Point(280) = {1.134007,1.196000,-0.001232};
+Point(281) = {1.136294,1.196000,-0.000933};
+Point(282) = {1.138230,1.196000,-0.000680};
+Point(283) = {1.139869,1.196000,-0.000466};
+Point(284) = {1.141256,1.196000,-0.000284};
+Point(285) = {1.142432,1.196000,-0.000130,msTeTp};
+
+// Box:
+Point(5001) = {-lgt,0.000000,-hgt,msF};
+Point(5002) = {1+lgt,0.000000,-hgt,msF};
+Point(5003) = {-lgt,0.000000,hgt,msF};
+Point(5004) = {1+lgt,0.000000,hgt,msF};
+Point(5005) = {-lgt,wdt,-hgt,msF};
+Point(5006) = {1+lgt,wdt,-hgt,msF};
+Point(5007) = {-lgt,wdt,hgt,msF};
+Point(5008) = {1+lgt,wdt,hgt,msF};
+
+// Tip:
+Point(5101) = {1.142432,1.196000,0.000000};
+Point(5102) = {1.141256,1.196040,0.000000};
+Point(5103) = {1.139869,1.196088,0.000000};
+Point(5104) = {1.138230,1.196144,0.000000};
+Point(5105) = {1.136294,1.196210,0.000000};
+Point(5106) = {1.134007,1.196289,0.000000};
+Point(5107) = {1.131305,1.196381,0.000000};
+Point(5108) = {1.128112,1.196491,0.000000};
+Point(5109) = {1.124337,1.196620,0.000000};
+Point(5110) = {1.119873,1.196773,0.000000};
+Point(5111) = {1.114592,1.196954,0.000000};
+Point(5112) = {1.108341,1.197168,0.000000};
+Point(5113) = {1.100939,1.197422,0.000000,1.5*msTeTp};
+Point(5114) = {1.092167,1.197722,0.000000};
+Point(5115) = {1.081765,1.198079,0.000000};
+Point(5116) = {1.071292,1.198438,0.000000};
+Point(5117) = {1.060762,1.198798,0.000000};
+Point(5118) = {1.050174,1.199161,0.000000};
+Point(5119) = {1.039528,1.199526,0.000000};
+Point(5120) = {1.028824,1.199893,0.000000};
+Point(5121) = {1.018059,1.200262,0.000000};
+Point(5122) = {1.007236,1.200632,0.000000};
+Point(5123) = {0.996353,1.201005,0.000000};
+Point(5124) = {0.985409,1.201380,0.000000};
+Point(5125) = {0.974403,1.201757,0.000000};
+Point(5126) = {0.963336,1.202137,0.000000};
+Point(5127) = {0.952208,1.202518,0.000000};
+Point(5128) = {0.941016,1.202901,0.000000};
+Point(5129) = {0.929760,1.203287,0.000000};
+Point(5130) = {0.918441,1.203675,0.000000};
+Point(5131) = {0.907059,1.204065,0.000000};
+Point(5132) = {0.895611,1.204457,0.000000};
+Point(5133) = {0.884097,1.204852,0.000000};
+Point(5134) = {0.872518,1.205248,0.000000};
+Point(5135) = {0.860873,1.205648,0.000000};
+Point(5136) = {0.849160,1.206049,0.000000};
+Point(5137) = {0.837379,1.206453,0.000000};
+Point(5138) = {0.825530,1.206859,0.000000,1.5*msLeTp};
+Point(5139) = {0.699294,1.211213,0.000000};
+
+// Dummy tip center:
+Point(5349) = {1.100939,1.196000,0.000000};
+Point(5350) = {0.825530,1.196000,0.000000};
+
+
+// Midplane:
+Point(5351) = {1+lgt,0.000000,0.000000,msF};
+Point(5352) = {1+lgt,1.196000,0.000000,msF};
+Point(5353) = {1+lgt,wdt,0.000000,msF};
+Point(5354) = {1.143427,wdt,0.000000,msF};
+Point(5355) = {1.100939,wdt,0.000000,msF};
+Point(5356) = {0.825530,wdt,0.000000,msF};
+Point(5357) = {0.690511,wdt,0.000000,msF};
+Point(5358) = {-lgt,wdt,0.000000,msF};
+Point(5359) = {-lgt,1.196000,0.000000,msF};
+Point(5360) = {-lgt,0.000000,0.000000,msF};
+
+/// Lines
+
+// Airfoil 1:
+Spline(1) = {1,2,3,4,5,6,7,8,9,10,11,12,13,14};
+Spline(2) = {14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39};
+Spline(3) = {39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72};
+Spline(4) = {72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105};
+Spline(5) = {105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130};
+Spline(6) = {130,131,132,133,134,135,136,137,138,139,140,141,142,1};
+
+// Airfoil 2:
+Spline(7) = {144,145,146,147,148,149,150,151,152,153,154,155,156,157};
+Spline(8) = {157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182};
+Spline(9) = {182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215};
+Spline(10) = {215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248};
+Spline(11) = {248,249,250,251,252,253,254,255,256,257,258,259,260,261,262,263,264,265,266,267,268,269,270,271,272,273};
+Spline(12) = {273,274,275,276,277,278,279,280,281,282,283,284,285,144};
+
+
+// Box:
+Line(41) = {5001,5002};
+Line(42) = {5003,5004};
+Line(43) = {5005,5006};
+Line(44) = {5007,5008};
+Line(45) = {5001,5005};
+Line(46) = {5002,5006};
+Line(47) = {5003,5007};
+Line(48) = {5004,5008};
+Line(49) = {5001,5360};
+Line(50) = {5002,5351};
+Line(51) = {5003,5360};
+Line(52) = {5004,5351};
+Line(53) = {5005,5358};
+Line(54) = {5006,5353};
+Line(55) = {5007,5358};
+Line(56) = {5008,5353};
+
+// Airfoil 1 to airfoil 2:
+Line(61) = {1,144};
+Line(62) = {14,157};
+Line(63) = {39,182};
+Line(64) = {72,215};
+Line(65) = {105,248};
+Line(66) = {130,273};
+
+
+// Tip:
+Spline(121) = {144,5101,5102,5103,5104,5105,5106,5107,5108,5109,5110,5111,5112,5113};
+Spline(122) = {5113,5114,5115,5116,5117,5118,5119,5120,5121,5122,5123,5124,5125,5126,5127,5128,5129,5130,5131,5132,5133,5134,5135,5136,5137,5138};
+If(GMSH_MAJOR_VERSION >= 4)
+ Bezier(123) = {5138,5139,215};
+Else
+ BSpline(123) = {5138,5139,215};
+EndIf
+Ellipse(124) = {157,5349,157,5113};
+Ellipse(125) = {182,5350,182,5138};
+Ellipse(126) = {248,5350,248,5138};
+Ellipse(127) = {273,5349,273,5113};
+
+// Midplane:
+Line(131) = {5351,5352};
+Line(132) = {5352,5353};
+Line(133) = {5353,5354};
+Line(134) = {5354,5355};
+Line(135) = {5355,5356};
+Line(136) = {5356,5357};
+Line(137) = {5357,5358};
+Line(138) = {5358,5359};
+Line(139) = {5359,5360};
+
+// Wing to midplane:
+Line(161) = {1,5351};
+Line(162) = {144,5352};
+Line(163) = {144,5354};
+Line(164) = {5113,5355};
+Line(165) = {5138,5356};
+Line(166) = {215,5357};
+Line(167) = {215,5359};
+Line(168) = {72,5360};
+
+/// Line loops & Surfaces
+
+// Box:
+Line Loop(1) = {1,2,3,168,-51,42,52,-161};
+Line Loop(2) = {48,56,-132,-131,-52};
+Line Loop(3) = {44,56,133,134,135,136,137,-55};
+Line Loop(4) = {47,55,138,139,-51};
+Line Loop(5) = {42,48,-44,-47};
+Line Loop(6) = {-6,-5,-4,168,-49,41,50,-161};
+Line Loop(7) = {46,54,-132,-131,-50};
+Line Loop(8) = {43,54,133,134,135,136,137,-53};
+Line Loop(9) = {45,53,138,139,-49};
+Line Loop(10) = {41,46,-43,-45};
+Plane Surface(1) = {1};
+Plane Surface(2) = {2};
+Plane Surface(3) = {-3};
+Plane Surface(4) = {-4};
+Plane Surface(5) = {-5};
+Plane Surface(6) = {-6};
+Plane Surface(7) = {-7};
+Plane Surface(8) = {8};
+Plane Surface(9) = {9};
+Plane Surface(10) = {10};
+
+// Wing 1:
+Line Loop(11) = {1,62,-7,-61};
+Line Loop(12) = {2,63,-8,-62};
+Line Loop(13) = {3,64,-9,-63};
+Line Loop(14) = {4,65,-10,-64};
+Line Loop(15) = {5,66,-11,-65};
+Line Loop(16) = {6,61,-12,-66};
+Surface(11) = {-11};
+Surface(12) = {-12};
+Surface(13) = {-13};
+Surface(14) = {-14};
+Surface(15) = {-15};
+Surface(16) = {-16};
+
+// Wingtip:
+Line Loop(71) = {7,124,-121};
+Line Loop(72) = {8,125,-122,-124};
+Line Loop(73) = {9,-123,-125};
+Line Loop(74) = {10,126,123};
+Line Loop(75) = {11,127,122,-126};
+Line Loop(76) = {12,121,-127};
+Surface(71) = {-71};
+Surface(72) = {-72};
+Surface(73) = {-73};
+Surface(74) = {-74};
+Surface(75) = {-75};
+Surface(76) = {-76};
+
+// Midplane:
+Line Loop(81) = {161,131,-162,-61};
+Line Loop(82) = {162,132,133,-163};
+Line Loop(83) = {163,134,-164,-121};
+Line Loop(84) = {164,135,-165,-122};
+Line Loop(85) = {165,136,-166,-123};
+Line Loop(86) = {167,-138,-137,-166};
+Line Loop(87) = {167,139,-168,64};
+Surface(81) = {81};
+Surface(82) = {82};
+Surface(83) = {83};
+Surface(84) = {84};
+Surface(85) = {85};
+Surface(86) = {-86};
+Surface(87) = {87};
+
+/// Surface loops & Volumes
+
+// Upper:
+Surface Loop(1) = {11,12,13,71,72,73,1,2,3,4,5,81,82,83,84,85,86,87};
+Volume(1) = {1};
+// Lower:
+Surface Loop(2) = {14,15,16,74,75,76,6,7,8,9,10,81,82,83,84,85,86,87};
+Volume(2) = {2};
+
+
+
+//// MESHING ALGORITHM
+
+
+/// 2D:
+
+///Wing, farfield and symmetry plane:
+MeshAlgorithm Surface{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,81,82,83,84,85,86,87} = 5;
+
+///Tip:
+MeshAlgorithm Surface{73,74} = 1;
+MeshAlgorithm Surface{71,72,75,76} = 5;
+
+/// 3D:
+
+Mesh.Algorithm3D = 2;
+Mesh.OptimizeNetgen = 1;
+Mesh.Smoothing = 10;
+Mesh.SmoothNormals = 1;
+
+
+
+//// PHYSICAL GROUPS
+
+// Trailing edge and wake tip
+Physical Line("wakeTip") = {162};
+Physical Line("teTip") = {61, 162};
+
+// Internal Field:
+Physical Volume("field") = {1};
+Physical Volume("field") += {2};
+
+// Wing:
+Physical Surface("wing") = {11,12,13,71,72,73};
+Physical Surface("wing_") = {14,15,16,74,75,76};
+
+// Tip
+//Physical Surface("tip") = {71,72,73,74,75,76};
+
+// Symmetry:
+Physical Surface("symmetry") = {1};
+Physical Surface("symmetry") += {6};
+
+// Farfield:
+Physical Surface("upstream") = {10};
+Physical Surface("farfield") = {3,4,5,8,9};
+
+// Downstream:
+Physical Surface("downstream") = {2};
+Physical Surface("downstream") += {7};
+
+// Wake:
+Physical Surface("wake") = {81};
+
+Coherence;
diff --git a/dart/pyVII/vAPI.py b/dart/pyVII/vAPI.py
deleted file mode 100644
index 59a9599bfc30cc1cd990a0c31ae80d759304a76c..0000000000000000000000000000000000000000
--- a/dart/pyVII/vAPI.py
+++ /dev/null
@@ -1,102 +0,0 @@
-import numpy as np
-import dart.pyVII.vStruct as vStruct
-from matplotlib import pyplot as plt
-from scipy import interpolate
-
-class dartAPI:
-
- def __init__(self, _dartSolver, bnd_airfoil, bnd_wake, _nSections):
- self.nSections = _nSections
- self.iSolver = _dartSolver
- self.invStruct = [[vStruct.Airfoil(bnd_airfoil), vStruct.Wake(bnd_wake)] for _ in range(self.nSections)] # Parameters passing protocols
- self.viscStruct = [[vStruct.viscousIPS("upper"), vStruct.viscousIPS("lower"), vStruct.viscousIPS("wake")] for _ in range(self.nSections)]
-
- def GetInvBC(self, vSolver, couplIter):
-
- # Extract parameters
- for iSec in range(self.nSections):
- dataIn = [None, None]
- for n in range(0, len(self.invStruct[iSec])):
- for i in range(0, len(self.invStruct[iSec][n].boundary.nodes)):
-
- self.invStruct[iSec][n].v[i,0] = self.iSolver.U[self.invStruct[iSec][n].boundary.nodes[i].row][0]
- self.invStruct[iSec][n].v[i,1] = self.iSolver.U[self.invStruct[iSec][n].boundary.nodes[i].row][1]
- self.invStruct[iSec][n].v[i,2] = self.iSolver.U[self.invStruct[iSec][n].boundary.nodes[i].row][2]
- self.invStruct[iSec][n].Me[i] = self.iSolver.M[self.invStruct[iSec][n].boundary.nodes[i].row]
- self.invStruct[iSec][n].rhoe[i] = self.iSolver.rho[self.invStruct[iSec][n].boundary.nodes[i].row]
-
- self.invStruct[iSec][n].connectListNodes, self.invStruct[iSec][n].connectListElems, dataIn[n] = self.invStruct[iSec][n].connectList()
-
- if len(dataIn[n]) == 2:
- self.viscStruct[iSec][0].SetValues(dataIn[n][0])
- self.viscStruct[iSec][1].SetValues(dataIn[n][1])
- else:
- self.viscStruct[iSec][2].SetValues(dataIn[n])
-
- for iRegion in range(len(self.viscStruct[iSec])):
- if couplIter==0:
- self.viscStruct[iSec][iRegion].xPrev = self.viscStruct[iSec][iRegion].x
- vSolver.SetGlobMesh(iSec, iRegion, self.viscStruct[iSec][iRegion].x,
- self.viscStruct[iSec][iRegion].y,
- self.viscStruct[iSec][iRegion].z)
-
- vSolver.SetVelocities(iSec, iRegion, self.viscStruct[iSec][iRegion].vx,
- self.viscStruct[iSec][iRegion].vy,
- self.viscStruct[iSec][iRegion].vz)
-
- vSolver.SetMe(iSec, iRegion, self.viscStruct[iSec][iRegion].Me)
- vSolver.SetRhoe(iSec, iRegion, self.viscStruct[iSec][iRegion].Rhoe)
- vSolver.SetxxExt(iSec, iRegion, self.viscStruct[iSec][iRegion].xxExt)
- if len(self.viscStruct[iSec][iRegion].x) != len(self.viscStruct[iSec][iRegion].xPrev):
- f = interpolate.interp1d(self.viscStruct[iSec][iRegion].xPrev, self.viscStruct[iSec][iRegion].DeltaStarExt, bounds_error=False, fill_value="extrapolate")
- newDS = f(self.viscStruct[iSec][iRegion].x)
- plt.plot(self.viscStruct[iSec][iRegion].xPrev, self.viscStruct[iSec][iRegion].DeltaStarExt, '--')
- plt.plot(self.viscStruct[iSec][iRegion].x, newDS)
- plt.show()
- self.viscStruct[iSec][iRegion].DeltaStarExt = newDS
- vSolver.SetDeltaStarExt(iSec, iRegion, self.viscStruct[iSec][iRegion].DeltaStarExt)
- self.viscStruct[iSec][iRegion].xPrev = self.viscStruct[iSec][iRegion].x
-
- def SetBlowingVel(self, vSolver):
- """ Collects blowing velocities from inviscid solver, maps them to inviscid mesh and """
-
- for iSec in range(self.nSections):
- for iRegion in range(len(self.viscStruct[iSec])):
- self.viscStruct[iSec][iRegion].BlowingVel = np.asarray(vSolver.ExtractBlowingVel(iSec, iRegion))
- self.viscStruct[iSec][iRegion].DeltaStarExt = vSolver.ExtractDeltaStar(iSec, iRegion)
- self.viscStruct[iSec][iRegion].xxExt = vSolver.Extractxx(iSec, iRegion)
-
- blwVelAF = np.concatenate((self.viscStruct[iSec][0].BlowingVel, self.viscStruct[iSec][1].BlowingVel))
- blwVelWK = self.viscStruct[iSec][2].BlowingVel
-
- self.invStruct[iSec][0].u = blwVelAF[self.invStruct[iSec][0].connectListElems.argsort()]
- self.invStruct[iSec][1].u = blwVelWK[self.invStruct[iSec][1].connectListElems.argsort()]
-
- for n in range(0, len(self.invStruct[iSec])):
- for i in range(0, self.invStruct[iSec][n].nE):
- self.invStruct[iSec][n].boundary.setU(i, self.invStruct[iSec][n].u[i])
-
- def RunSolver(self):
- self.iSolver.run()
-
- def ResetSolver(self):
- self.iSolver.reset()
-
- def GetCp(self,bnd):
- import dart.utils as floU
-
- Cp = floU.extract(bnd, self.iSolver.Cp)
- return Cp
-
- def GetAlpha(self):
- return self.iSolver.pbl.alpha
-
- def GetMinf(self):
- return self.iSolver.pbl.M_inf
-
- def GetCl(self):
- return self.iSolver.Cl
-
- def GetCm(self):
- return self.iSolver.Cm
-
\ No newline at end of file
diff --git a/dart/pyVII/vAPI3.py b/dart/pyVII/vAPI3.py
deleted file mode 100644
index b5718a536de65cc63a8cd36364b851f344969108..0000000000000000000000000000000000000000
--- a/dart/pyVII/vAPI3.py
+++ /dev/null
@@ -1,107 +0,0 @@
-import numpy as np
-import dart.pyVII.vStruct as vStruct
-from matplotlib import pyplot as plt
-from scipy import interpolate
-
-class dartAPI:
-
- def __init__(self, _dartSolver, bnd_wing, bnd_wake, _nSections, _interp):
- self.nSections = _nSections
- self.iSolver = _dartSolver
- self.bndWing = bnd_wing
- self.bndWk = bnd_wake
- #self.invStruct = [[vStruct.Airfoil(bnd_airfoil), vStruct.Wake(bnd_wake)] for _ in range(self.nSections)] # Parameters passing protocols
- self.viscStruct = [[vStruct.viscousIPS("upper"), vStruct.viscousIPS("lower"), vStruct.viscousIPS("wake")] for _ in range(self.nSections)]
- self.interp = _interp
-
- def GetInvBC(self, vSolver, couplIter):
- Ux = np.zeros(len(self.bndWing.nodes))
- Uy = np.zeros(len(self.bndWing.nodes))
- Uz = np.zeros(len(self.bndWing.nodes))
- Me = np.zeros(len(self.bndWing.nodes))
- Rho = np.zeros(len(self.bndWing.nodes))
- UxWk = np.zeros(len(self.bndWk.nodMap))
- UyWk = np.zeros(len(self.bndWk.nodMap))
- UzWk = np.zeros(len(self.bndWk.nodMap))
- MeWk = np.zeros(len(self.bndWk.nodMap))
- RhoWk = np.zeros(len(self.bndWk.nodMap))
-
- # From wing
- for node in range(len(self.bndWing.nodes)):
- Ux[node] = self.iSolver.U[self.bndWing.nodes[node].row][0]
- Uy[node] = self.iSolver.U[self.bndWing.nodes[node].row][1]
- Uz[node] = self.iSolver.U[self.bndWing.nodes[node].row][2]
-
- Me[node] = self.iSolver.M[self.bndWing.nodes[node].row]
- Rho[node] = self.iSolver.rho[self.bndWing.nodes[node].row]
-
- # From wake
- for iNodeWk, n in enumerate(self.bndWk.nodMap):
- UxWk[iNodeWk] = self.iSolver.U[n.row][0]
- UyWk[iNodeWk] = self.iSolver.U[n.row][1]
- UzWk[iNodeWk] = self.iSolver.U[n.row][2]
-
- MeWk[iNodeWk] = self.iSolver.M[n.row]
- RhoWk[iNodeWk] = self.iSolver.rho[n.row]
-
- mapedUx = self.interp.InterpToSections(Ux, UxWk)
- mapedUy = self.interp.InterpToSections(Uy, UyWk)
- mapedUz = self.interp.InterpToSections(Uz, UzWk)
-
- """plt.plot(self.interp.Sections[3][0][:,0], np.sqrt(mapedUx[3][0]**2 + mapedUy[3][0]**2 +mapedUz[3][0]**2))
- plt.plot(self.interp.Sections[3][1][:,0], np.sqrt(mapedUx[3][1]**2 + mapedUy[3][1]**2 +mapedUz[3][1]**2))
- plt.plot(self.interp.Sections[3][2][:,0], np.sqrt(mapedUx[3][2]**2 + mapedUy[3][2]**2 +mapedUz[3][2]**2))
- plt.show()"""
-
- mapedMe = self.interp.InterpToSections(Me, MeWk)
- mapedRho = self.interp.InterpToSections(Rho, RhoWk)
-
- for iSec in range(len(self.interp.Sections)):
- for iReg in range(3):
- vSolver.SetGlobMesh(iSec, iReg, self.interp.Sections[iSec][iReg][:,0],
- self.interp.Sections[iSec][iReg][:,1],
- self.interp.Sections[iSec][iReg][:,2])
- vSolver.SetVelocities(iSec, iReg, mapedUx[iSec][iReg], mapedUy[iSec][iReg], mapedUz[iSec][iReg])
- vSolver.SetMe(iSec, iReg, mapedMe[iSec][iReg])
- vSolver.SetRhoe(iSec, iReg, mapedRho[iSec][iReg])
- if couplIter == 0:
- self.viscStruct[iSec][iReg].xxExt = np.zeros(len(self.interp.Sections[iSec][iReg][:,0]))
- self.viscStruct[iSec][iReg].DeltaStarExt = np.zeros(len(self.interp.Sections[iSec][iReg][:,0]))
- vSolver.SetxxExt(iSec, iReg, self.viscStruct[iSec][iReg].xxExt)
- vSolver.SetDeltaStarExt(iSec, iReg, self.viscStruct[iSec][iReg].DeltaStarExt)
-
- def SetBlowingVel(self, vSolver):
- """ Collects blowing velocities from inviscid solver, maps them to inviscid mesh and """
-
- for iSec in range(self.nSections):
- for iRegion in range(len(self.viscStruct[iSec])):
- self.viscStruct[iSec][iRegion].BlowingVel = np.asarray(vSolver.ExtractBlowingVel(iSec, iRegion))
- self.viscStruct[iSec][iRegion].DeltaStarExt = vSolver.ExtractDeltaStar(iSec, iRegion)
- self.viscStruct[iSec][iRegion].xxExt = vSolver.Extractxx(iSec, iRegion)
-
- self.interp.InterpToElements(self.viscStruct)
-
- def RunSolver(self):
- self.iSolver.run()
-
- def ResetSolver(self):
- self.iSolver.reset()
-
- def GetCp(self,bnd):
- import dart.utils as floU
-
- Cp = floU.extract(bnd, self.iSolver.Cp)
- return Cp
-
- def GetAlpha(self):
- return self.iSolver.pbl.alpha
-
- def GetMinf(self):
- return self.iSolver.pbl.M_inf
-
- def GetCl(self):
- return self.iSolver.Cl
-
- def GetCm(self):
- return self.iSolver.Cm
-
\ No newline at end of file
diff --git a/dart/pyVII/vCoupler.py b/dart/pyVII/vCoupler.py
index 8c623cab11c7408867949ee15df5b479450b1bb4..7f0b8c97109498a8534218903db521294ef4ff17 100644
--- a/dart/pyVII/vCoupler.py
+++ b/dart/pyVII/vCoupler.py
@@ -38,10 +38,11 @@ class Coupler:
self.couplTol = 1e-4
self.tms=fwk.Timers()
- self.resetInv = False
+ self.resetInv = True
pass
def Run(self):
+ self.tms['tot'].start()
# Initilize meshes in c++ objects
@@ -82,14 +83,15 @@ class Coupler:
print(ccolors.ANSI_GREEN,'{:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>8.2f}\n'.format(couplIter,
self.iSolverAPI.GetCl(), self.vSolver.Cdt, self.vSolver.Getxtr(0,0), self.vSolver.Getxtr(0,1), np.log10(error)),ccolors.ANSI_RESET)
+ self.tms['tot'].stop()
+
return 0
cdPrev = self.vSolver.Cdt
- if couplIter%5==0:
- print(' {:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>8.2f}'.format(couplIter,
- self.iSolverAPI.GetCl(), self.vSolver.Cdt, self.vSolver.Getxtr(0,0), self.vSolver.Getxtr(0,1), np.log10(error)))
+ print('- {:>4.0f}| {:>10.5f} {:>10.5f} | {:>8.2f} {:>6.2f} {:>8.2f}'.format(couplIter,
+ self.iSolverAPI.GetCl(), self.vSolver.Cdt, self.vSolver.Getxtr(0,0), self.vSolver.Getxtr(0,1), np.log10(error)))
couplIter += 1
"""if couplIter == 5:
@@ -97,6 +99,8 @@ class Coupler:
self.tms['processing'].stop()
else:
print(ccolors.ANSI_RED, 'Maximum number of {:<.0f} coupling iterations reached'.format(self.maxCouplIter), ccolors.ANSI_RESET)
+ self.tms['tot'].stop()
+
return couplIter
diff --git a/dart/pyVII/vInterpolator.py b/dart/pyVII/vInterpolator.py
index 01518f4e15fe148f518c932ec342b0ecee15f6b1..9dc992089c4ef81e4151820d9a420f20dd23d637 100644
--- a/dart/pyVII/vInterpolator.py
+++ b/dart/pyVII/vInterpolator.py
@@ -1,55 +1,375 @@
+import fwk
+
import numpy as np
import math as m
from matplotlib import pyplot as plt
-from scipy import interpolate
+from scipy.interpolate import bisplrep
+from scipy.interpolate import bisplev
+from scipy.interpolate import interp1d
+from scipy.interpolate import RBFInterpolator
class Interpolator:
- def __init__(self, bnd) -> None:
- arf, up, lw = self.GetArf(bnd)
+ def __init__(self, _dartSolver, bnd_wing, bnd_wake, nSections, nDim, nPoints=[200, 25], eps=0.01, k=[5, 5]) -> None:
+ self.iSolver = _dartSolver
+
+ self.bndWing = bnd_wing
+ self.bndWake = bnd_wake
+
+ self.nDim = nDim
+
+ self.Sections, self.Wing, self.Wake = self.__GetWing(nSections, nDim, nPoints, eps, k)
+ self.xxPrev = [[np.zeros(len(self.Sections[iSec][iReg][:,0])) for iReg in range(len(self.Sections[iSec]))] for iSec in range(len(self.Sections))]
+ self.DeltaStarPrev = [[np.zeros(len(self.Sections[iSec][iReg][:,0])) for iReg in range(len(self.Sections[iSec]))] for iSec in range(len(self.Sections))]
+ self.tms = fwk.Timers()
+
+ def GetInvBC(self, vSolver, couplIter):
+ self.tms['InvToVisc'].start()
+
+ # Extract parameters from the inviscid solver
+ # Particular care must be taken with space directions; In the inviscid, the airfoil is in the x,z plane.
+ # While setting up the viscous solver, y and z must be inverted.
+
+ # Wing
+
+ nElm = len(self.bndWing.nodes)
+ U_wg, M_wg, Rho_wg = np.zeros((nElm, 3)), np.zeros(nElm), np.zeros(nElm)
+ for iNode, n in enumerate(self.bndWing.nodes):
+ U_wg[iNode,0] = self.iSolver.U[n.row][0]
+ U_wg[iNode,1] = self.iSolver.U[n.row][1]
+ U_wg[iNode,2] = self.iSolver.U[n.row][2]
+ M_wg[iNode] = self.iSolver.M[n.row]
+ Rho_wg[iNode] = self.iSolver.rho[n.row]
+
+ # Wake
+
+ nElm = len(self.bndWake.nodes)
+ U_wk, M_wk, Rho_wk = np.zeros((nElm, 3)), np.zeros(nElm), np.zeros(nElm)
+ for iNode, n in enumerate(self.bndWake.nodes):
+ U_wk[iNode,0] = self.iSolver.U[n.row][0]
+ U_wk[iNode,1] = self.iSolver.U[n.row][1]
+ U_wk[iNode,2] = self.iSolver.U[n.row][2]
+ M_wk[iNode] = self.iSolver.M[n.row]
+ Rho_wk[iNode] = self.iSolver.rho[n.row]
+
+ self.tms['Rbf'].start()
+ Ux = self.__InterpToSections(U_wg[:,0], U_wk[:,0])
+ Uy = self.__InterpToSections(U_wg[:,1], U_wk[:,1])
+ Uz = self.__InterpToSections(U_wg[:,2], U_wk[:,2])
+
+ Me = self.__InterpToSections(M_wg, M_wk)
+ Rho = self.__InterpToSections(Rho_wg, Rho_wk)
+ self.tms['Rbf'].stop()
+
+ # Find stagnation point
+ for iSec in range(len(self.Sections)):
+
+ idxStgUp = np.argmin(np.sqrt(Ux[iSec][0]**2+Uy[iSec][0]**2))
+ idxStgLw = np.argmin(np.sqrt(Ux[iSec][1]**2+Uy[iSec][1]**2))
+
+ if np.sqrt(Ux[iSec][0][idxStgUp]**2+Uy[iSec][0][idxStgUp]**2) >= np.sqrt(Ux[iSec][1][idxStgLw]**2+Uy[iSec][1][idxStgLw]**2):
+ reg = 1
+ idx = idxStgLw
+ else:
+ reg = 0
+ idx = idxStgUp
+
+ xUp, xLw = self.__Remesh(self.Sections[iSec][0][:,0], self.Sections[iSec][1][:,0], idx, reg)
+ yUp, yLw = self.__Remesh(self.Sections[iSec][0][:,1], self.Sections[iSec][1][:,1], idx, reg)
+ zUp, zLw = self.__Remesh(self.Sections[iSec][0][:,2], self.Sections[iSec][1][:,2], idx, reg)
+
+ vSolver.SetGlobMesh(iSec, 0, xUp, zUp, yUp)
+ vSolver.SetGlobMesh(iSec, 1, xLw, zLw, yLw)
+ vSolver.SetGlobMesh(iSec, 2, self.Sections[iSec][2][:,0],
+ self.Sections[iSec][2][:,1],
+ self.Sections[iSec][2][:,2])
+
+ UxUp, UxLw = self.__Remesh(Ux[iSec][0], Ux[iSec][1], idx, reg)
+ UyUp, UyLw = self.__Remesh(Uy[iSec][0], Uy[iSec][1], idx, reg)
+ UzUp, UzLw = self.__Remesh(Uz[iSec][0], Uz[iSec][1], idx, reg)
+
+ vSolver.SetVelocities(iSec, 0, UxUp, UzUp, UyUp)
+ vSolver.SetVelocities(iSec, 1, UxLw, UzLw, UyLw)
+ vSolver.SetVelocities(iSec, 2, Ux[iSec][2], Uz[iSec][2], Uy[iSec][2])
+
+ MeUp, MeLw = self.__Remesh(Me[iSec][0], Me[iSec][1], idx, reg)
+ vSolver.SetMe(iSec, 0, MeUp)
+ vSolver.SetMe(iSec, 1, MeLw)
+ vSolver.SetMe(iSec, 2, Me[iSec][2])
+
+ RhoUp, RhoLw = self.__Remesh(Rho[iSec][0], Rho[iSec][1], idx, reg)
+ vSolver.SetRhoe(iSec, 0, RhoUp)
+ vSolver.SetRhoe(iSec, 1, RhoLw)
+ vSolver.SetRhoe(iSec, 2, Rho[iSec][2])
+
+ for iReg in range(3):
+ vSolver.SetxxExt(iSec, iReg, self.xxPrev[iSec][iReg])
+ vSolver.SetDeltaStarExt(iSec, iReg, self.DeltaStarPrev[iSec][iReg])
+ self.tms['InvToVisc'].stop()
+
+ def SetBlowingVel(self, vSolver):
+ self.tms['ViscToInv'].start()
+ BlowingVel = [[np.zeros(0) for iReg in range(len(self.Sections[iSec]))] for iSec in range(len(self.Sections))]
+ xCoord = [[np.zeros(0) for iReg in range(len(self.Sections[iSec]))] for iSec in range(len(self.Sections))]
+ for iSec in range(len(self.Sections)):
+ for iRegion in range(len(self.Sections[iSec])):
+ xCoord[iSec][iRegion] = np.asarray(vSolver.ExtractxCoord(iSec, iRegion))
+ BlowingVel[iSec][iRegion] = np.asarray(vSolver.ExtractBlowingVel(iSec, iRegion))
+ self.DeltaStarPrev[iSec][iRegion] = vSolver.ExtractDeltaStar(iSec, iRegion)
+ self.xxPrev[iSec][iRegion] = vSolver.Extractxx(iSec, iRegion)
+
+ if self.nDim == 2:
+
+ # Compute cell centroid x positions (in the viscous mesh)
+
+ xCellsUp = np.zeros(len(xCoord[0][0])-1)
+ xCellsUp[0] = (xCoord[0][0][1] - xCoord[0][0][0])/2
+ cumDx = xCoord[0][0][1] - xCoord[0][0][0]
+ for iCell in range(1, len(xCoord[0][0])-1):
+ xCellsUp[iCell] = cumDx + (xCoord[0][0][iCell+1] - xCoord[0][0][iCell])/2
+ cumDx += (xCoord[0][0][iCell] - xCoord[0][0][iCell-1]) # Cumulative space between cells
+
+ xCellsLw = np.zeros(len(xCoord[0][1])-1)
+ xCellsLw[0] = (xCoord[0][1][1] - xCoord[0][1][0])/2
+ cumDx = xCoord[0][1][1] - xCoord[0][0][0]
+ for iCell in range(1, len(xCoord[0][1])-1):
+ xCellsLw[iCell] = cumDx + (xCoord[0][1][iCell+1] - xCoord[0][1][iCell])/2
+ cumDx += (xCoord[0][1][iCell] - xCoord[0][1][iCell-1]) # Cumulative space between cells
+
+ xCellsWk = np.zeros(len(xCoord[0][2])-1)
+ xCellsWk[0] = (xCoord[0][2][1] - xCoord[0][2][0])/2
+ cumDx = xCoord[0][2][1] - xCoord[0][2][0]
+ for iCell in range(1, len(xCoord[0][2])-1):
+ xCellsWk[iCell] = cumDx + (xCoord[0][2][iCell+1] - xCoord[0][2][iCell])/2
+ cumDx += (xCoord[0][2][iCell] - xCoord[0][2][iCell-1]) # Cumulative space between cells
+
+ # Create interpolation functions
+
+ blwSplnUp = interp1d(xCellsUp, BlowingVel[0][0], bounds_error=False, fill_value="extrapolate")
+ blwSplnLw = interp1d(xCellsLw, BlowingVel[0][1], bounds_error=False, fill_value="extrapolate")
+ blwSplnWk = interp1d(xCellsWk, BlowingVel[0][2], bounds_error=False, fill_value="extrapolate")
+
+ # Airfoil (Cell centroids are now computed in the inviscid mesh)
+
+ for iNode in range(len(self.bndWing.nodes)-1):
+ if self.bndWing.nodes[iNode].pos[1]>=0:
+ self.bndWing.setU(iNode, float(blwSplnUp(self.bndWing.nodes[iNode].pos[0])))
+ else:
+ self.bndWing.setU(iNode, float(blwSplnLw(self.bndWing.nodes[iNode].pos[0])))
+
+ for iNode in range(len(self.bndWake.nodes)-1):
+ self.bndWake.setU(iNode, float(blwSplnWk(self.bndWake.nodes[iNode].pos[0])))
+
+
+ elif self.nDim == 3:
+ pass
+ self.tms['ViscToInv'].stop()
+
+
+ def RunSolver(self):
+ self.iSolver.run()
+
+ def ResetSolver(self):
+ self.iSolver.reset()
+
+ def GetCp(self,bnd):
+ import dart.utils as floU
+
+ Cp = floU.extract(bnd, self.iSolver.Cp)
+ return Cp
+
+ def GetAlpha(self):
+ return self.iSolver.pbl.alpha
- def GetArf(self, bnd):
+ def GetMinf(self):
+ return self.iSolver.pbl.M_inf
+
+ def GetCl(self):
+ return self.iSolver.Cl
+
+ def GetCm(self):
+ return self.iSolver.Cm
+
+
+ def __InterpToSections(self, var, varWk):
+
+ varWk = np.reshape(varWk, [len(varWk), 1])
+ completeWake = np.hstack((self.Wake, varWk))
+ _, idx = np.unique(completeWake[:,:3], axis=0, return_index=True)
+ completeWake = completeWake[np.sort(idx)]
+
+ var = np.reshape(var, [len(var), 1])
+ completeWing = np.hstack((self.Wing, var))
+
+ wingUp = completeWing[completeWing[:,2]>=0]
+ wingLw = completeWing[completeWing[:,2]<=0]
+
+ mappedVar = [[np.zeros(0) for _ in range(3)] for _ in range(len(self.Sections))]
+
+ for iSec in range(len(mappedVar)):
+ for iReg in range(len(mappedVar[iSec])):
+ if iReg==0:
+ mappedVar[iSec][iReg] = RBFInterpolator(wingUp[:,:3], wingUp[:,3], neighbors=200, smoothing=0.1, kernel='linear', degree=0)(self.Sections[iSec][iReg])
+ if iReg==1:
+ mappedVar[iSec][iReg] = RBFInterpolator(wingLw[:,:3], wingLw[:,3], neighbors=200, smoothing=0.1, kernel='linear', degree=0)(self.Sections[iSec][iReg])
+ if iReg==2:
+ mappedVar[iSec][iReg] = RBFInterpolator(completeWake[:,:2], completeWake[:,3], smoothing=0.1, kernel='linear', degree=0)(self.Sections[iSec][iReg][:,:2])
+ return mappedVar
+
+
+ def __Remesh(self, vectorUp, vectorLw, idx, reg):
+ """Rebuilds input vector to match the stagnation point position.
+
+ args:
+ ----
+ - VectorUp/VectorLw upper/lower side vector of the desired variable from x_glob = 0 to x_glob = 1
+ - idx : Stagnation point position from leading edge point (x_glob=0)
+ - reg : 0/1 indicates on which side the stagnation point is.
+
+ return : vUp/vLw distribution from x_glob=stag to x_glob=1 on upper/lower sides.
+ ------
+ """
+
+ if reg==0:
+ vLw = np.insert(vectorLw, 0, np.flip(vectorUp[1:idx]))
+ vUp = np.delete(vectorUp, np.s_[:idx-1])
+ elif reg==1:
+ vUp = np.insert(vectorUp, 0, np.flip(vectorLw[1:idx]))
+ vLw = np.delete(vectorLw, np.s_[:idx-1])
+ else:
+ raise RuntimeError('Incorrect region specified during remeshing.')
+
+ return vUp, vLw
+
+
+ def __GetWing(self, nSections, nDim, nPoints, eps, k):
"Extract and sort the node coordinates of the airfoil."
- arf = np.zeros((len(self.bnd.nodes),3))
+
+ if nDim == 2 and nSections != 1:
+ raise RuntimeError('Cannot create more than 1 section on 2D geometry. Specified:', nSections)
+
+ wing = np.zeros((len(self.bndWing.nodes),3))
+ wake = np.zeros((len(self.bndWake.nodes),3))
# Extract coordinates.
- for iNode, n in enumerate(bnd.nodes):
- for iDir in range(len(arf[0,:])):
- arf[iNode,iDir] = n.pos[iDir]
+ for iNode, n in enumerate(self.bndWing.nodes):
+ for iDir in range(len(wing[0,:])):
+ wing[iNode,iDir] = n.pos[iDir]
+
+ for iNode, n in enumerate(self.bndWake.nodes):
+ for iDir in range(len(wake[0,:])):
+ wake[iNode, iDir] = n.pos[iDir]
+
+ if nDim == 2: # 2D case
+ wing[:,[1,2]] = wing[:,[2,1]]
+ wake[:,[1,2]] = wake[:,[2,1]]
+
+ # Separate upper and lower sides
+
+ wingUp = wing[wing[:,2]>=0]
+ wingLw = wing[wing[:,2]<=0]
+
+ # Sort upper/lower side and remove duplicates
+ wingUp = np.delete(wingUp, wingUp[:,1]>1, 0)
+ wingLw = np.delete(wingLw, wingLw[:,1]>1, 0)
+
+ _, idxUp = np.unique(wingUp, axis=0, return_index=True)
+ wingUp = wingUp[np.sort(idxUp)]
+ _, idxLw = np.unique(wingLw, axis=0, return_index=True)
+ wingLw = wingLw[np.sort(idxLw)]
+
+ if nDim == 2:
+ arfSplnUp = interp1d(wingUp[:,0], wingUp[:,2], kind='cubic')
+ arfSplnLw = interp1d(wingLw[:,0], wingLw[:,2], kind='cubic')
+
+ Sections = [[np.zeros((nPoints[1] if i==2 else nPoints[0], 3)) for i in range(3)]]
+ for iReg in range(len(Sections[0])):
+ Sections[0][iReg][:,1] = wingUp[0,1]
+ if iReg != 2:
+ Sections[0][iReg][:,0] = np.linspace(min(wing[:,0]), max(wing[:,0]), nPoints[0])
+
+ if iReg == 0:
+ for iPoint, x in enumerate(Sections[0][iReg][:,0]):
+ Sections[0][iReg][iPoint,2] = arfSplnUp(x)
+ elif iReg == 1:
+ for iPoint, x in enumerate(Sections[0][iReg][:,0]):
+ Sections[0][iReg][iPoint,2] = arfSplnLw(x)
+ else:
+ Sections[0][iReg][:,0] = np.linspace(min(wake[:,0]), max(wake[:,0]), nPoints[1])
+
- upper = np.empty((0,3))
- lower = np.empty((0,3))
+ elif nDim == 3:
+ sUp = (len(wingUp[:,0])-np.sqrt(2*len(wingUp[:,0])) + len(wingUp[:,0])+np.sqrt(2*len(wingUp[:,0])))/2
+ sLw = (len(wingLw[:,0])-np.sqrt(2*len(wingLw[:,0])) + len(wingLw[:,0])+np.sqrt(2*len(wingLw[:,0])))/2
- # Split between upper and lower surfaces.
- for i in range(len(arf[:,0])):
- if arf[i,1]>=0:
- upper=np.vstack([upper, arf[i,:]])
- else:
- lower=np.vstack([lower, arf[i,:]])
+ tckUp = bisplrep(wingUp[:,0], wingUp[:,1], wingUp[:,2], kx=k[0], ky=k[1], s=sUp)
+ tckLw = bisplrep(wingLw[:,0], wingLw[:,1], wingLw[:,2], kx=k[0], ky=k[1], s=sLw)
+
+ zDistri = np.linspace(min(wing[:,1]), max(wing[:,1]), nSections)
+
+ zDistri[-1] -= 0.2
+
+ Sections = [[np.zeros((nPoints[1] if i==2 else nPoints[0], 3)) for i in range(3)] for _ in range(nSections)]
+
+ for iSec, y in enumerate(zDistri):
+ for iReg in range(3):
+ print('Sec', iSec, 'y', y, min(wing[abs(wing[:,1]-y)<eps, 0]), max(wing[abs(wing[:,1]-y)<eps, 0]))
+ print(wing[abs(wing[:,1]-y)<0.02, :])
+ if iReg == 0:
+ Sections[iSec][iReg][:,0] = np.linspace(min(wing[abs(wing[:,1]-y)<eps, 0]), max(wing[abs(wing[:,1]-y)<eps, 0]), nPoints[0])
+ elif iReg == 1:
+ Sections[iSec][iReg][:,0] = np.linspace(min(wing[abs(wing[:,1]-y)<eps, 0]), max(wing[abs(wing[:,1]-y)<eps, 0]), nPoints[0])
+ elif iReg == 2:
+ #print('sec', iSec, 'min', min(wake[:, 0]), 'max', max(wake[:, 0]))
+ Sections[iSec][iReg][:,0] = np.linspace(max(wing[abs(wing[:,1]-y)<eps, 0]), max(wake[:, 0]), nPoints[1])
+ Sections[iSec][iReg][:,1] = y
+ for iPt in range(len(Sections[iSec][iReg][:,0])):
+ if iReg==0:
+ Sections[iSec][iReg][iPt,2] = bisplev(Sections[iSec][iReg][iPt, 0], Sections[iSec][iReg][iPt, 1], tckUp)
+ elif iReg==1:
+ Sections[iSec][iReg][iPt,2] = bisplev(Sections[iSec][iReg][iPt, 0], Sections[iSec][iReg][iPt, 1], tckLw)
+ elif iReg == 2:
+ Sections[iSec][iReg][iPt,2] = 0
+ else:
+ raise RuntimeError("dart::vInterpolator3 l.83 0<=iReg<=2")
+
+ # Add leading edge point
+ #Sections.pop()
+
+ for iSec in range(len(Sections)):
+ for iReg in range(2):
+ Sections[iSec][iReg][0,2] = 0
- # Sort according to x.
- upper = upper[upper[:, 0].argsort()]
- lower = lower[lower[:, 0].argsort()]
-
- # Look for doublons.
- for i in range(len(upper[:,0])-1):
- if upper[i+1,0] == 0 and upper[i,0] == 1:
- lower = np.vstack([lower, upper[i,:]])
- upper = np.delete(upper, i, 0)
- break
- return arf, upper, lower
-
-
- def meshGeneration(self):
- print(arf[:,0])
- x_up = np.linspace(0,1,100)
- x_lw = np.linspace(0,1,100)
- func_up = interpolate.interp1d(arf[:self.sep,0], arf[:self.sep,1], bounds_error=False, fill_value="extrapolate")
- func_lw = interpolate.interp1d(arf[self.sep:,0], arf[self.sep:,1], bounds_error=False, fill_value="extrapolate")
- y_up = func_up(x_up)
- y_lw = func_lw(x_lw)
-
- plt.plot(arf[:self.sep-1,0], arf[:self.sep-1,1])
- plt.plot(x_up,y_up, 'o')
- plt.plot(x_lw,y_lw, 'o')
- plt.show()
\ No newline at end of file
+ """fig = plt.figure()
+ ax = fig.add_subplot(projection='3d')
+ ax.scatter(wing[:,0], wing[:,1], wing[:,2], s=0.3)
+ ax.scatter(wake[:,0], wake[:,1], wake[:,2], s=0.3)
+ for iSec in range(len(Sections)):
+ ax.scatter(Sections[iSec][0][:,0], Sections[iSec][0][:,1], Sections[iSec][0][:,2], color='red')
+ ax.scatter(Sections[iSec][1][:,0], Sections[iSec][1][:,1], Sections[iSec][1][:,2], color='red')
+ ax.scatter(Sections[iSec][2][:,0], Sections[iSec][2][:,1], Sections[iSec][2][:,2], color='red')
+ ax.set_zlim([-0.5, 0.5])
+ ax.set_xlabel('x')
+ ax.set_ylabel('y')
+ ax.set_zlabel('z')
+ plt.show()
+
+ plt.plot(Sections[0][0][:,0], Sections[0][0][:,2], 'x-', color='red')
+ plt.plot(Sections[0][1][:,0], Sections[0][1][:,2], 'x-', color='red')
+ plt.plot(wingUp[:,0], wingUp[:,2], 'o', color='blue')
+ plt.plot(wingLw[:,0], wingLw[:,2], 'o', color='blue')
+ plt.ylim([-0.5, 0.5])
+ plt.show()
+ #quit()
+
+ for iSec in range(len(Sections)):
+ plt.plot(Sections[iSec][0][:,0], Sections[iSec][0][:,2], 'x-', color='red')
+ plt.plot(Sections[iSec][1][:,0], Sections[iSec][1][:,2], 'x-', color='red')
+ plt.xlim([-0.5, 1.5])
+ plt.ylim([-0.5,0.5])
+ plt.show()"""
+ return Sections, wing, wake
+
+
+
+
\ No newline at end of file
diff --git a/dart/pyVII/vInterpolator3.py b/dart/pyVII/vInterpolator3.py
deleted file mode 100644
index c02a5384c6ee64bcde555d5d6725c62325021ef2..0000000000000000000000000000000000000000
--- a/dart/pyVII/vInterpolator3.py
+++ /dev/null
@@ -1,195 +0,0 @@
-import numpy as np
-import math as m
-from matplotlib import pyplot as plt
-from scipy.interpolate import bisplrep
-from scipy.interpolate import bisplev
-from scipy import stats
-from scipy.interpolate import RBFInterpolator
-
-import sys
-
-class Interpolator:
-
- def __init__(self, bnd, wk, nSections, blw) -> None:
- self.bndWing = bnd
- self.bndWake = wk
- self.blwWing = blw[0]
- self.blwWake = blw[1]
- self.Sections, self.wing, self.wake = self.GetWing(bnd, wk, nSections)
-
-
- def GetWing(self, bnd, wk, nSections, nPoints=[500, 25], eps=5e-2, k=[5, 5]):
-
- "Extract and sort the node coordinates of the airfoil."
- wing = np.zeros((len(bnd.nodes),3))
- wake = np.zeros((len(wk.nodMap),3))
-
- # Extract coordinates.
- for iNode, n in enumerate(bnd.nodes):
- for iDir in range(len(wing[0,:])):
- wing[iNode,iDir] = n.pos[iDir]
-
- for iNode, n in enumerate(wk.nodMap):
- for iDir in range(len(wake[0,:])):
- wake[iNode, iDir] = n.pos[iDir]
-
- # Separate upper and lower sides
-
- wingUp = wing[wing[:,2]>=0]
- wingLw = wing[wing[:,2]<=0]
-
- wingUp = np.delete(wingUp, wingUp[:,1]>1, 0)
- wingLw = np.delete(wingLw, wingLw[:,1]>1, 0)
-
- xUp = wingUp[:,0].copy()
- xLw = wingLw[:,0].copy()
-
- np.set_printoptions(threshold=sys.maxsize)
-
- sUp = (len(xUp)-np.sqrt(2*len(xUp)) + len(xUp)+np.sqrt(2*len(xUp)))/2
- sLw = (len(xLw)-np.sqrt(2*len(xLw)) + len(xLw)+np.sqrt(2*len(xLw)))/2
-
- tckUp = bisplrep(wingUp[:,0], wingUp[:,1], wingUp[:,2], kx=k[0], ky=k[1], s=sUp)
- tckLw = bisplrep(wingLw[:,0], wingLw[:,1], wingLw[:,2], kx=k[0], ky=k[1], s=sLw)
-
- yDistri = np.linspace(min(wing[:,1]), max(wing[:,1]), nSections)
-
- Sections = [[np.zeros((nPoints[1] if i==2 else nPoints[0], 3)) for i in range(3)] for _ in range(nSections)]
-
- for iSec, y in enumerate(yDistri):
- for iReg in range(3):
- if iReg == 0:
- Sections[iSec][iReg][:,0] = np.linspace(min(wing[abs(wing[:,1]-y)<eps, 0]), max(wing[abs(wing[:,1]-y)<eps, 0]), nPoints[0])
- elif iReg == 1:
- Sections[iSec][iReg][:,0] = np.linspace(min(wing[abs(wing[:,1]-y)<eps, 0]), max(wing[abs(wing[:,1]-y)<eps, 0]), nPoints[0])
- elif iReg == 2:
- print('sec', iSec, 'min', min(wake[:, 0]), 'max', max(wake[:, 0]))
- Sections[iSec][iReg][:,0] = np.linspace(max(wing[abs(wing[:,1]-y)<eps, 0]), max(wake[:, 0]), nPoints[1])
- Sections[iSec][iReg][:,1] = y
- for iPt in range(len(Sections[iSec][iReg][:,0])):
- if iReg==0:
- Sections[iSec][iReg][iPt,2] = bisplev(Sections[iSec][iReg][iPt, 0], Sections[iSec][iReg][iPt, 1], tckUp)
- elif iReg==1:
- Sections[iSec][iReg][iPt,2] = bisplev(Sections[iSec][iReg][iPt, 0], Sections[iSec][iReg][iPt, 1], tckLw)
- elif iReg == 2:
- Sections[iSec][iReg][iPt,2] = 0
- else:
- raise RuntimeError("dart::vInterpolator3 l.83 0<=iReg<=2")
-
- # Add leading edge point
-
- for iSec in range(len(Sections)):
- for iReg in range(2):
- Sections[iSec][iReg][0,2] = 0
-
- """fig = plt.figure()
- ax = fig.add_subplot(projection='3d')
- ax.scatter(wing[:,0], wing[:,1], wing[:,2], s=0.3)
- ax.scatter(wake[:,0], wake[:,1], wake[:,2], s=0.3)
- for iSec in range(len(Sections)):
- ax.scatter(Sections[iSec][0][:,0], Sections[iSec][0][:,1], Sections[iSec][0][:,2], color='red')
- ax.scatter(Sections[iSec][1][:,0], Sections[iSec][1][:,1], Sections[iSec][1][:,2], color='red')
- ax.scatter(Sections[iSec][2][:,0], Sections[iSec][2][:,1], Sections[iSec][2][:,2], color='red')
- ax.set_zlim([-0.5, 0.5])
- ax.set_xlabel('x')
- ax.set_ylabel('y')
- ax.set_zlabel('z')
- plt.show()
-
- plt.plot(Sections[0][0][:,0], Sections[0][0][:,2], 'x-', color='red')
- plt.plot(Sections[0][1][:,0], Sections[0][1][:,2], 'x-', color='red')
- plt.plot(wingUp[:,0], wingUp[:,2], 'o', color='blue')
- plt.plot(wingLw[:,0], wingLw[:,2], 'o', color='blue')
- plt.ylim([-0.5, 0.5])
- plt.show()"""
- return Sections, wing, wake
-
- def InterpToSections(self, var, varWk):
-
- varWk = np.reshape(varWk, [len(varWk), 1])
- completeWake = np.hstack((self.wake, varWk))
- _, idx = np.unique(completeWake[:,:3], axis=0, return_index=True)
- completeWake = completeWake[np.sort(idx)]
-
- var = np.reshape(var, [len(var), 1])
- completeWing = np.hstack((self.wing, var))
-
- wingUp = completeWing[completeWing[:,2]>=0]
- wingLw = completeWing[completeWing[:,2]<=0]
-
- mappedVar = [[np.zeros(0) for _ in range(3)] for _ in range(len(self.Sections))]
-
- for iSec in range(len(mappedVar)):
- for iReg in range(len(mappedVar[iSec])):
- if iReg==0:
- mappedVar[iSec][iReg] = RBFInterpolator(wingUp[:,:3], wingUp[:,3], neighbors=75, smoothing=1.5)(self.Sections[iSec][iReg])
- if iReg==1:
- mappedVar[iSec][iReg] = RBFInterpolator(wingLw[:,:3], wingLw[:,3], neighbors=75, smoothing=1.5)(self.Sections[iSec][iReg])
- if iReg==2:
- mappedVar[iSec][iReg] = RBFInterpolator(completeWake[:,:2], completeWake[:,3], smoothing=1.5)(self.Sections[iSec][iReg][:,:2])
- return mappedVar
-
- def InterpToElements(self, viscStruct):
-
- blwwing = np.zeros((len(self.blwWing.nodes),3))
- blwwake = np.zeros((len(self.blwWake.nodes),3))
-
- # Extract coordinates.
- for iNode, n in enumerate(self.blwWing.nodes):
- for iDir in range(len(blwwing[0,:])):
- blwwing[iNode,iDir] = n.pos[iDir]
-
- for iNode, n in enumerate(self.blwWake.nodes):
- for iDir in range(len(blwwake[0,:])):
- blwwake[iNode, iDir] = n.pos[iDir]
-
- coord_Up=np.zeros((3))
- coord_Lw=np.zeros((3))
- coord_Wk=np.zeros((3))
- BlwTot_Up=np.zeros(1)
- BlwTot_Lw=np.zeros(1)
- BlwTot_Wk=np.zeros(1)
- for iSec in range(len(viscStruct)):
-
- coord_Up = np.vstack((coord_Up, self.Sections[iSec][0]))
- if iSec!=0:
- BlwTot_Up = np.insert(BlwTot_Up, -1, 0)
- BlwTot_Up = np.concatenate([BlwTot_Up, viscStruct[iSec][0].BlowingVel], axis=None)
-
- coord_Lw = np.vstack((coord_Lw, self.Sections[iSec][1]))
- if iSec!=0:
- BlwTot_Lw = np.insert(BlwTot_Lw, -1, 0)
- BlwTot_Lw = np.concatenate([BlwTot_Lw, viscStruct[iSec][1].BlowingVel], axis=None)
-
- coord_Wk = np.vstack((coord_Wk, self.Sections[iSec][2]))
- if iSec!=0:
- BlwTot_Wk = np.insert(BlwTot_Wk, -1, 0)
- BlwTot_Wk = np.concatenate([BlwTot_Wk, viscStruct[iSec][2].BlowingVel], axis=None)
-
- coord_Up = np.delete(coord_Up, 0, axis=0)
- coord_Lw = np.delete(coord_Lw, 0, axis=0)
- coord_Wk = np.delete(coord_Wk, 0, axis=0)
-
- coord = np.vstack((coord_Up, coord_Lw))
- blw = np.concatenate((BlwTot_Up, BlwTot_Lw), axis=0)
-
-
- blw = np.reshape(blw, [len(blw), 1])
- BlwTot_Wk = np.reshape(BlwTot_Wk, [len(BlwTot_Wk), 1])
- completeWing = np.hstack((coord, blw))
- completeWake = np.hstack((coord_Wk, BlwTot_Wk))
-
- _, idx = np.unique(completeWing[:,:3], axis=0, return_index=True)
- completeWing = completeWing[np.sort(idx)]
-
- _, idxWk = np.unique(completeWake[:,:3], axis=0, return_index=True)
- completeWake = completeWake[np.sort(idxWk)]
-
- mappedBlowingWing = RBFInterpolator(completeWing[:,:3], completeWing[:,3], smoothing=1.5)(blwwing)
- mappedBlowingWake = RBFInterpolator(completeWake[:,:2], completeWake[:,3], smoothing=1.5)(blwwake[:,:2])
-
-
- for iNode in range(len(blwwing)):
- self.blwWing.setU(iNode, mappedBlowingWing[iNode])
- for iNode in range(len(blwwake)):
- self.blwWake.setU(iNode, mappedBlowingWake[iNode])
\ No newline at end of file
diff --git a/dart/pyVII/vStruct.py b/dart/pyVII/vStruct.py
deleted file mode 100644
index 65a895cbbb5bbad4ff171e7124cb0d3c7f449818..0000000000000000000000000000000000000000
--- a/dart/pyVII/vStruct.py
+++ /dev/null
@@ -1,193 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2020 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.
-
-
-## Information passing structures (IPS) for communication between solvers.
-
-import numpy as np
-
-class viscousIPS:
- def __init__(self, _name) -> None:
- self.name = _name
- self.DeltaStarExt = np.zeros(1)
- self.xxExt = np.zeros(1)
-
- def SetValues(self, data):
-
- self.x = data[:,0]
- self.y = data[:,1]
- self.z = data[:,2]
- self.vx = data[:,3]
- self.vy = data[:,4]
- self.vz = np.zeros(len(data[:,4]))
- self.Me = data[:,5]
- self.Rhoe = data[:,6]
- #self.DeltaStarExt = data[:,7]
- #self.xxExt = data[:,8]
-
- self.BlowingVel = np.zeros(len(data[:,0]))
-
-class inviscidIPS:
- def __init__(self, _boundary):
- self.boundary = _boundary # boundary of interest
- self.nN = len(self.boundary.nodes) # number of nodes
- self.nE = len(self.boundary.tag.elems) # number of elements
- self.v = np.zeros((self.nN, 3)) # velocity at edge of BL
- self.u = np.zeros(self.nE) # blowing velocity
- self.Me = np.zeros(self.nN) # Mach number at the edge of the boundary layer
- self.rhoe = np.zeros(self.nN) # density at the edge of the boundary layer
- self.connectListnodes = np.zeros(self.nN) # connectivity for nodes
- self.connectListElems = np.zeros(self.nE) # connectivity for elements
- self.deltaStar = np.zeros(self.nN) # displacement thickness
- self.xx = np.zeros(self.nN) # coordinates in the reference of the physical surface. Used to store the values for the interaction law
-
-class Airfoil(inviscidIPS):
- def __init__(self, _boundary):
- inviscidIPS.__init__(self, _boundary)
-
- def connectList(self):
- ''' Sort the value read by the viscous solver/ Create list of connectivity
- '''
- N1 = np.zeros(self.nN, dtype=int) # Node number
- connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
- connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
- data = np.zeros((self.boundary.nodes.size(), 10))
- i = 0
- for n in self.boundary.nodes:
- data[i,0] = n.no
- data[i,1] = n.pos[0]
- data[i,2] = n.pos[1]
- data[i,3] = n.pos[2]
- data[i,4] = self.v[i,0]
- data[i,5] = self.v[i,1]
- data[i,6] = self.Me[i]
- data[i,7] = self.rhoe[i]
- data[i,8] = self.deltaStar[i]
- data[i,9] = self.xx[i]
- i += 1
- # Table containing the element and its nodes
- eData = np.zeros((self.nE,3), dtype=int)
- for i in range(0, self.nE):
- eData[i,0] = self.boundary.tag.elems[i].no
- eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
- eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
- # Find the stagnation point
- idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2))
- globStag = int(data[idxStag,0]) # position of the stagnation point in boundary.nodes
- # Find TE nodes (assuming suction side element is above pressure side element in the y direction)
- idxTE = np.where(data[:,1] == np.max(data[:,1]))[0] # TE nodes
- idxTEe0 = np.where(np.logical_or(eData[:,1] == data[idxTE[0],0], eData[:,2] == data[idxTE[0],0]))[0] # TE element 1
- idxTEe1 = np.where(np.logical_or(eData[:,1] == data[idxTE[1],0], eData[:,2] == data[idxTE[1],0]))[0] # TE element 2
- ye0 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe0[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe0[0],2])[0],2]) # y coordinates of TE element 1 CG
- ye1 = 0.5 * (data[np.where(data[:,0] == eData[idxTEe1[0],1])[0],2] + data[np.where(data[:,0] == eData[idxTEe1[0],2])[0],2]) # y coordinates of TE element 2 CG
- if ye0 - ye1 > 0: # element 1 containing TE node 1 is above element 2
- upperGlobTE = int(data[idxTE[0],0])
- lowerGlobTE = int(data[idxTE[1],0])
- else:
- upperGlobTE = int(data[idxTE[1],0])
- lowerGlobTE = int(data[idxTE[0],0])
- # Connectivity
- connectListElems[0] = np.where(eData[:,1] == globStag)[0]
- N1[0] = eData[connectListElems[0],1] # number of the stag node elems.nodes -> globStag
- connectListNodes[0] = np.where(data[:,0] == N1[0])[0]
- i = 1
- upperTE = 0
- lowerTE = 0
- # Sort the suction part
- while upperTE == 0:
- N1[i] = eData[connectListElems[i-1],2] # Second node of the element
- connectListElems[i] = np.where(eData[:,1] == N1[i])[0] # # Index of the first node of the next element in elems.nodes
- connectListNodes[i] = np.where(data[:,0] == N1[i])[0] # Index of the node in boundary.nodes
- if eData[connectListElems[i],2] == upperGlobTE:
- upperTE = 1
- i += 1
- # Sort the pressure side
- connectListElems[i] = np.where(eData[:,2] == globStag)[0]
- connectListNodes[i] = np.where(data[:,0] == upperGlobTE)[0]
- N1[i] = eData[connectListElems[i],2]
- while lowerTE == 0:
- N1[i+1] = eData[connectListElems[i],1] # First node of the element
- connectListElems[i+1] = np.where(eData[:,2] == N1[i+1])[0] # # Index of the second node of the next element in elems.nodes
- connectListNodes[i+1] = np.where(data[:,0] == N1[i+1])[0] # Index of the node in boundary.nodes
- if eData[connectListElems[i+1],1] == lowerGlobTE:
- lowerTE = 1
- i += 1
- connectListNodes[i+1] = np.where(data[:,0] == lowerGlobTE)[0]
- data[:,0] = data[connectListNodes,0]
- data[:,1] = data[connectListNodes,1]
- data[:,2] = data[connectListNodes,2]
- data[:,3] = data[connectListNodes,3]
- data[:,4] = data[connectListNodes,4]
- data[:,5] = data[connectListNodes,5]
- data[:,6] = data[connectListNodes,6]
- data[:,7] = data[connectListNodes,7]
- data[:,8] = data[connectListNodes,8]
- data[:,9] = data[connectListNodes,9]
- # Separated upper and lower part
- data = np.delete(data,0,1)
- uData = data[0:np.argmax(data[:,0])+1]
- lData = data[np.argmax(data[:,0])+1:None]
- lData = np.insert(lData, 0, uData[0,:], axis = 0) #double the stagnation point
- return connectListNodes, connectListElems,[uData, lData]
-
-class Wake(inviscidIPS):
- def __init__(self, _boundary):
- inviscidIPS.__init__(self, _boundary)
-
- def connectList(self):
- ''' Sort the value read by the viscous solver/ Create list of connectivity
- '''
- N1 = np.zeros(self.nN, dtype=int) # Node number
- connectListNodes = np.zeros(self.nN, dtype=int) # Index in boundary.nodes
- connectListElems = np.zeros(self.nE, dtype=int) # Index in boundary.elems
- data = np.zeros((self.boundary.nodes.size(), 10))
- i = 0
- for n in self.boundary.nodes:
- data[i,0] = n.no
- data[i,1] = n.pos[0]
- data[i,2] = n.pos[1]
- data[i,3] = n.pos[2]
- data[i,4] = self.v[i,0]
- data[i,5] = self.v[i,1]
- data[i,6] = self.Me[i]
- data[i,7] = self.rhoe[i]
- data[i,8] = self.deltaStar[i]
- data[i,9] = self.xx[i]
- i += 1
- # Table containing the element and its nodes
- eData = np.zeros((self.nE,3), dtype=int)
- for i in range(0, self.nE):
- eData[i,0] = self.boundary.tag.elems[i].no
- eData[i,1] = self.boundary.tag.elems[i].nodes[0].no
- eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
- connectListNodes = data[:,1].argsort()
- connectListElems[0] = np.where(eData[:,1] == min(data[:,1]))[0]
- for i in range(1, len(eData[:,0])):
- connectListElems[i] = np.where(eData[:,1] == eData[connectListElems[i-1],2])[0]
- data[:,0] = data[connectListNodes,0]
- data[:,1] = data[connectListNodes,1]
- data[:,2] = data[connectListNodes,2]
- data[:,3] = data[connectListNodes,3]
- data[:,4] = data[connectListNodes,4]
- data[:,5] = data[connectListNodes,5]
- data[:,6] = data[connectListNodes,6]
- data[:,7] = data[connectListNodes,7]
- data[:,8] = data[connectListNodes,8]
- data[:,9] = data[connectListNodes,9]
- # Separated upper and lower part
- data = np.delete(data,0,1)
- return connectListNodes, connectListElems, data
\ No newline at end of file
diff --git a/dart/src/wTimeSolver.cpp b/dart/src/wTimeSolver.cpp
index a3c0c7d65271a53172ef550d5be3b1001f1d390e..49a5d4f0044b695d3bea6aaba21926526db5e2f4 100644
--- a/dart/src/wTimeSolver.cpp
+++ b/dart/src/wTimeSolver.cpp
@@ -57,6 +57,11 @@ void TimeSolver::InitialCondition(size_t iPoint, BLRegion *bl)
{
bl->U[iPoint * nVar + 4] = 0.03;
} // End if
+
+ if (bl->U[iPoint*nVar+3]<0)
+ {
+ std::cout << "dart::TimeSolver negative velocity encountered at point " << iPoint << std::endl;
+ }
} // End InitialCondition
int TimeSolver::Integration(size_t iPoint, BLRegion *bl)
diff --git a/dart/src/wViscSolver.cpp b/dart/src/wViscSolver.cpp
index 3e993035f3b4b33fa7fd1512e04bdf93211ed2bd..5e911afa8923ee602c45fc1e11e379095cc766a4 100644
--- a/dart/src/wViscSolver.cpp
+++ b/dart/src/wViscSolver.cpp
@@ -15,6 +15,7 @@
#include "wTimeSolver.h"
#include "wInterpolator.h"
#include <iostream>
+#include <iomanip>
#include <vector>
#include <cmath>
@@ -95,14 +96,13 @@ ViscSolver::~ViscSolver()
*/
int ViscSolver::Run(unsigned int couplIter)
{
- std::cout << "Starting solve" << std::endl;
bool lockTrans;
int solverExitCode = 0;
int pointExitCode;
for (size_t iSec=0; iSec<Sections.size(); ++iSec)
{
- std::cout << "Section: " << iSec << std::endl;
+ std::cout << "\n - Sec " << iSec << " ";
/* Prepare time integration */
@@ -124,7 +124,7 @@ int ViscSolver::Run(unsigned int couplIter)
std::cout << Sections[iSec][0]->mesh->GetDiffnElm() <<
Sections[iSec][1]->mesh->GetDiffnElm() <<
Sections[iSec][2]->mesh->GetDiffnElm() << std::endl;*/
-
+ resetSolution = true;
if (resetSolution)
{
tSolver->SetCFL0(1);
@@ -156,7 +156,7 @@ int ViscSolver::Run(unsigned int couplIter)
}
if (printOn)
{
- std::cout << "Restart solution" << std::endl;
+ std::cout << "restart. ";
}
}
@@ -168,7 +168,7 @@ int ViscSolver::Run(unsigned int couplIter)
stagPtMvmt[iSec] = false;
if (printOn)
{
- std::cout << "Updating solution" << std::endl;
+ std::cout << "update. ";
}
}
@@ -192,7 +192,7 @@ int ViscSolver::Run(unsigned int couplIter)
{
if (printOn)
{
- std::cout << "- " << Sections[iSec][iRegion]->name << " region";
+ std::cout << Sections[iSec][iRegion]->name << " ";
}
lockTrans = false;
@@ -229,7 +229,9 @@ int ViscSolver::Run(unsigned int couplIter)
tSolver->InitialCondition(iPoint, Sections[iSec][iRegion]);
/* Time integration */
- /*std::cout << " iSec " << iSec
+ /*if (iRegion == 0 && iPoint <4)
+ {
+ std::cout << " iSec " << iSec
<< " iPoint " << iPoint
<< " Loc " << Sections[iSec][iRegion]->mesh->GetLoc(iPoint)
<< " Vt " << Sections[iSec][iRegion]->invBnd->GetVt(iPoint)
@@ -237,7 +239,8 @@ int ViscSolver::Run(unsigned int couplIter)
<< " xxExt " << Sections[iSec][iRegion]->mesh->GetExt(iPoint)
<< " DeltaStarExt " << Sections[iSec][iRegion]->invBnd->GetDeltaStar(iPoint)
<< std::endl;
- Sections[iSec][iRegion]->PrintSolution(iPoint);*/
+ Sections[iSec][iRegion]->PrintSolution(iPoint);
+ }*/
pointExitCode = tSolver->Integration(iPoint, Sections[iSec][iRegion]);
@@ -261,11 +264,11 @@ int ViscSolver::Run(unsigned int couplIter)
{
AverageTransition(iPoint, Sections[iSec][iRegion], 0);
if (printOn)
- {
- std::cout << ", free transition x/c = "
- << Sections[iSec][iRegion]->xtr
- << ". Marker " << Sections[iSec][iRegion]->transMarker
- << std::endl;
+ {
+ std::cout << std::fixed;
+ std::cout << std::setprecision(2);
+ std::cout << Sections[iSec][iRegion]->xtr
+ << " (" << Sections[iSec][iRegion]->transMarker << ") ";
}
lockTrans = true;
} // End if
@@ -389,7 +392,7 @@ void ViscSolver::AverageTransition(size_t iPoint, BLRegion *bl, int forced)
* Cdf = integral(Cf dx)_upper + integral(Cf dx)_lower.
* Cdp = Cdtot - Cdf.
*/
-void ViscSolver::ComputeDrag(std::vector<BLRegion *> bl)
+void ViscSolver::ComputeDrag(std::vector<BLRegion *> const &bl)
{
unsigned int nVar = bl[0]->GetnVar();
@@ -475,11 +478,21 @@ void ViscSolver::SetDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double
Sections[iSec][iRegion]->invBnd->SetDeltaStar(_DeltaStarExt);
}
-std::vector<double> ViscSolver::ExtractBlowingVel(size_t iSec, size_t iRegion)
+std::vector<double> ViscSolver::ExtractxCoord(size_t iSec, size_t iRegion) const
+{
+ std::vector<double> xCoord;
+ for (size_t iPoint=0; iPoint<Sections[iSec][iRegion]->mesh->GetnMarkers(); ++iPoint)
+ {
+ xCoord.push_back(Sections[iSec][iRegion]->mesh->Getx(iPoint));
+ }
+ return xCoord;
+}
+
+std::vector<double> ViscSolver::ExtractBlowingVel(size_t iSec, size_t iRegion) const
{
return Sections[iSec][iRegion]->Blowing;
}
-std::vector<double> ViscSolver::Extractxx(size_t iSec, size_t iRegion)
+std::vector<double> ViscSolver::Extractxx(size_t iSec, size_t iRegion) const
{
std::vector<double> xx(Sections[iSec][iRegion]->mesh->GetnMarkers(),0);
for (size_t iPoint=0; iPoint<Sections[iSec][iRegion]->mesh->GetnMarkers(); ++iPoint)
@@ -488,7 +501,7 @@ std::vector<double> ViscSolver::Extractxx(size_t iSec, size_t iRegion)
}
return xx;
}
-std::vector<double> ViscSolver::ExtractDeltaStar(size_t iSec, size_t iRegion)
+std::vector<double> ViscSolver::ExtractDeltaStar(size_t iSec, size_t iRegion) const
{
return Sections[iSec][iRegion]->DeltaStar;
}
diff --git a/dart/src/wViscSolver.h b/dart/src/wViscSolver.h
index d9486092be26593449bbfb1dfddbdb1d5b9fcbdd..a998acf0ebe9a3bf51fcf31a2daf26139a29ed94 100644
--- a/dart/src/wViscSolver.h
+++ b/dart/src/wViscSolver.h
@@ -38,16 +38,17 @@ public:
void SetxxExt(size_t iSec, size_t iRegion, std::vector<double> _xxExt);
void SetDeltaStarExt(size_t iSec, size_t iRegion, std::vector<double> _DeltaStarExt);
- std::vector<double> ExtractBlowingVel(size_t iSec, size_t iRegion);
- std::vector<double> Extractxx(size_t iSec, size_t iRegion);
- std::vector<double> ExtractDeltaStar(size_t iSec, size_t iRegion);
+ std::vector<double> ExtractBlowingVel(size_t iSec, size_t iRegion) const;
+ std::vector<double> Extractxx(size_t iSec, size_t iRegion) const;
+ std::vector<double> ExtractDeltaStar(size_t iSec, size_t iRegion) const;
+ std::vector<double> ExtractxCoord(size_t iSec, size_t iRegion) const;
double Getxtr(size_t iSec, size_t iRegion) const {return Sections[iSec][iRegion]->xtr;};
private:
void CheckWaves(size_t iPoint, BLRegion *bl);
void AverageTransition(size_t iPoint, BLRegion *bl, int forced);
- void ComputeDrag(std::vector<BLRegion *> bl);
+ void ComputeDrag(std::vector<BLRegion *> const &bl);
void ComputeBlwVel();
void PrintBanner();
};
diff --git a/dart/tests/bli3.py b/dart/tests/bli3.py
index 17e63da7ee8717b753b528fe19310260d0149bd0..5ab3cfe3b37eba033d85e4a003867d9be41d1f81 100644
--- a/dart/tests/bli3.py
+++ b/dart/tests/bli3.py
@@ -30,9 +30,8 @@ import dart.utils as floU
import dart.default as floD
import dart.pyVII.vUtils as viscU
-import dart.pyVII.vAPI3 as viscAPI
import dart.pyVII.vCoupler as viscC
-import dart.pyVII.vInterpolator3 as vInterpol
+import dart.pyVII.vInterpolator as vInterpol
import fwk
from fwk.testing import *
@@ -44,9 +43,9 @@ def main():
tms['total'].start()
# define flow variables
- Re = 1e7
- alpha = 3*math.pi/180
- M_inf = 0.3
+ Re = 11.72e6
+ alpha = 3.06*math.pi/180
+ M_inf = 0.8395
dim = 3
CFL0 = 1
@@ -60,13 +59,13 @@ def main():
S_ref = 1.*spn # reference area
c_ref = 1 # reference chord
fms = 1.0 # farfield mesh size
- nms = 0.02 # nearfield mesh size
+ nms = 0.01 # nearfield mesh size
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
- pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msLe' : nms, 'msTe' : nms, 'msF' : fms}
- msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
+ pars = {'spn' : spn, 'lgt' : lgt, 'wdt' : wdt, 'hgt' : hgt, 'msLe' : 0.001, 'msTe' : nms, 'msF' : fms}
+ msh, gmshWriter = floD.mesh(dim, 'models/oneraM6.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
tms['msh'].stop()
# set the problem and add fluid, initial/boundary conditions
@@ -75,20 +74,26 @@ def main():
tms['pre'].stop()
# solve problem
- vInterp = vInterpol.Interpolator(bnd, wk, nSections, blw)
+ iSolverAPI = vInterpol.Interpolator(floD.newton(pbl), blw[0], blw[1], nSections, dim)
vSolver = viscU.initBL(Re, M_inf, CFL0, nSections)
- iSolverAPI = viscAPI.dartAPI(floD.newton(pbl), bnd, wk, nSections, vInterp)
+ #iSolverAPI = viscAPI.dartAPI(floD.newton(pbl), bnd, wk, nSections, vInterp)
coupler = viscC.Coupler(iSolverAPI, vSolver)
- coupler.Run()
+ try:
+ coupler.Run()
+ except Exception as e:
+ print('VII convergence terminated, ', str(e.args[0]))
# display timers
tms['total'].stop()
print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
print('CPU statistics')
print(tms)
+ print(coupler.tms)
# visualize solution
+ iSolverAPI.GetCp()
+ iSolverAPI.iSolver.save(gmshWriter)
floD.initViewer(pbl)
# check results
diff --git a/dart/tests/bli_lowLift.py b/dart/tests/bli_lowLift.py
index 3c79f984844fbba1040784e1c55a1c13417651b8..7544e5eb81e5226c9b874296c9a72d6a74364d5f 100644
--- a/dart/tests/bli_lowLift.py
+++ b/dart/tests/bli_lowLift.py
@@ -34,11 +34,9 @@
# This test is provided to ensure that the solver works properly.
# Mesh refinement may have to be performed to obtain physical results.
-import dart.utils as floU
import dart.default as floD
import dart.pyVII.vUtils as viscU
-import dart.pyVII.vAPI as viscAPI
import dart.pyVII.vCoupler as viscC
import dart.pyVII.vInterpolator as vInterpol
@@ -60,29 +58,23 @@ def main():
# define flow variables
Re = 1e7
- alpha = 10.*math.pi/180
+ alpha = 5.*math.pi/180
#alphaSeq = np.arange(-5, 5.5, 0.5)
M_inf = 0.
- meshFactor = 1
CFL0 = 1
- plotVar = 'H'
-
# ========================================================================================
c_ref = 1
dim = 2
- if dim == 2:
- nSections = 1
- else:
- nSections = 10
+ nSections = 1
# mesh the geometry
print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
tms['msh'].start()
- pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.01}
+ pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1, 'msTe' : 0.01, 'msLe' : 0.001}
msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
tms['msh'].stop()
@@ -96,9 +88,8 @@ def main():
print(ccolors.ANSI_BLUE + 'PySolving...' + ccolors.ANSI_RESET)
tms['solver'].start()
- vInterp = vInterpol.Interpolator(bnd)
+ iSolverAPI = vInterpol.Interpolator(floD.newton(pbl), blw[0], blw[1], nSections, dim)
vSolver = viscU.initBL(Re, M_inf, CFL0, nSections)
- iSolverAPI = viscAPI.dartAPI(floD.newton(pbl), blw[0], blw[1], nSections)
coupler = viscC.Coupler(iSolverAPI, vSolver)
#coupler.RunPolar(alphaSeq)
@@ -120,6 +111,7 @@ def main():
print(tms)
print('SOLVERS statistics')
print(coupler.tms)
+ print(iSolverAPI.tms)
"""# visualize solution and plot results
#floD.initViewer(pbl)
diff --git a/dart/tests/blipython.py b/dart/tests/blipython.py
index 99e0816e43936b69b96f65126798551e9a59836e..d52d07bc4f7bce543d812c2ad8be37cbe540a159 100755
--- a/dart/tests/blipython.py
+++ b/dart/tests/blipython.py
@@ -60,8 +60,8 @@ def main():
# define flow variables
Re = 1e7
- alpha = 2.*math.pi/180
- M_inf = 0.715
+ alpha = 12.*math.pi/180
+ M_inf = 0.
#user_xtr=[0.41,0.41]
#user_xtr=[0,None]
@@ -72,7 +72,7 @@ def main():
nFactor = 2
# Time solver parameters
- Time_Domain = True # True for unsteady solver, False for steady solver
+ Time_Domain = False # True for unsteady solver, False for steady solver
CFL0 = 1
# ========================================================================================
diff --git a/dart/viscous/Solvers/Steady/StdAirfoil.py b/dart/viscous/Solvers/Steady/StdAirfoil.py
index 514e4681453473a376e1379044e683c7d06d7940..4a63b542de554d75db3f25c41ecfe1cc3257dde0 100755
--- a/dart/viscous/Solvers/Steady/StdAirfoil.py
+++ b/dart/viscous/Solvers/Steady/StdAirfoil.py
@@ -73,6 +73,7 @@ class Airfoil(Boundary):
eData[i,2] = self.boundary.tag.elems[i].nodes[1].no
# Find the stagnation point
idxStag = np.argmin(np.sqrt(data[:,4]**2+data[:,5]**2))
+ print('here',idxStag)
globStag = int(data[idxStag,0]) # position of the stagnation point in boundary.nodes
# Find TE nodes (assuming suction side element is above pressure side element in the y direction)
idxTE = np.where(data[:,1] == np.max(data[:,1]))[0] # TE nodes