Merge branch 'master' into 'master'

Master Thesis - Dechamps - Implementation of the linear part

See merge request !2

fpm/models/agard445.geo

+/* AGARD 445 wing */
+
+// Mesh parameters
+// domain and mesh
+DefineConstant[ wNc = { 50, Name "Number of panel along chord" }
+                wNs= { 10, Name "Number of panel along span" }
+                bC = { 0.05, Name "Progression along chord" }
+                pS = { 1.0, Name "Progression along span" } ];
+
+//// GEOMETRY
+/// Points
+// Airfoil 1: agard445, 51 points 
+Point(1) = {0.559000,0.000000,0.000000};
+Point(2) = {0.531050,0.000000,0.001398};
+Point(3) = {0.503100,0.000000,0.002739};
+Point(4) = {0.475150,0.000000,0.004075};
+Point(5) = {0.447200,0.000000,0.005406};
+Point(6) = {0.419250,0.000000,0.006680};
+Point(7) = {0.391300,0.000000,0.007837};
+Point(8) = {0.363350,0.000000,0.008866};
+Point(9) = {0.335400,0.000000,0.009743};
+Point(10) = {0.307450,0.000000,0.010442};
+Point(11) = {0.279500,0.000000,0.010923};
+Point(12) = {0.251550,0.000000,0.011158};
+Point(13) = {0.223600,0.000000,0.011163};
+Point(14) = {0.195650,0.000000,0.010968};
+Point(15) = {0.167700,0.000000,0.010576};
+Point(16) = {0.139750,0.000000,0.009995};
+Point(17) = {0.111800,0.000000,0.009196};
+Point(18) = {0.083850,0.000000,0.008156};
+Point(19) = {0.055900,0.000000,0.006781};
+Point(20) = {0.041925,0.000000,0.005920};
+Point(21) = {0.027950,0.000000,0.004891};
+Point(22) = {0.013975,0.000000,0.003617};
+Point(23) = {0.006988,0.000000,0.002622};
+Point(24) = {0.004193,0.000000,0.002057};
+Point(25) = {0.002795,0.000000,0.001699};
+Point(26) = {0.000000,0.000000,0.000000};
+Point(27) = {0.002795,0.000000,-0.001699};
+Point(28) = {0.004193,0.000000,-0.002057};
+Point(29) = {0.006988,0.000000,-0.002622};
+Point(30) = {0.013975,0.000000,-0.003617};
+Point(31) = {0.027950,0.000000,-0.004891};
+Point(32) = {0.041925,0.000000,-0.005920};
+Point(33) = {0.055900,0.000000,-0.006781};
+Point(34) = {0.083850,0.000000,-0.008156};
+Point(35) = {0.111800,0.000000,-0.009196};
+Point(36) = {0.139750,0.000000,-0.009995};
+Point(37) = {0.167700,0.000000,-0.010576};
+Point(38) = {0.195650,0.000000,-0.010968};
+Point(39) = {0.223600,0.000000,-0.011163};
+Point(40) = {0.251550,0.000000,-0.011158};
+Point(41) = {0.279500,0.000000,-0.010923};
+Point(42) = {0.307450,0.000000,-0.010442};
+Point(43) = {0.335400,0.000000,-0.009743};
+Point(44) = {0.363350,0.000000,-0.008866};
+Point(45) = {0.391300,0.000000,-0.007837};
+Point(46) = {0.419250,0.000000,-0.006680};
+Point(47) = {0.447200,0.000000,-0.005406};
+Point(48) = {0.475150,0.000000,-0.004075};
+Point(49) = {0.503100,0.000000,-0.002739};
+Point(50) = {0.531050,0.000000,-0.001398};
+// Airfoil 2: agard445, 51 points 
+Point(52) = {1.178128,0.762000,0.000000};
+Point(53) = {1.159709,0.762000,0.000921};
+Point(54) = {1.141290,0.762000,0.001805};
+Point(55) = {1.122870,0.762000,0.002685};
+Point(56) = {1.104451,0.762000,0.003562};
+Point(57) = {1.086032,0.762000,0.004402};
+Point(58) = {1.067613,0.762000,0.005165};
+Point(59) = {1.049194,0.762000,0.005843};
+Point(60) = {1.030775,0.762000,0.006421};
+Point(61) = {1.012356,0.762000,0.006881};
+Point(62) = {0.993937,0.762000,0.007198};
+Point(63) = {0.975518,0.762000,0.007353};
+Point(64) = {0.957099,0.762000,0.007357};
+Point(65) = {0.938680,0.762000,0.007228};
+Point(66) = {0.920261,0.762000,0.006970};
+Point(67) = {0.901842,0.762000,0.006587};
+Point(68) = {0.883423,0.762000,0.006060};
+Point(69) = {0.865004,0.762000,0.005375};
+Point(70) = {0.846585,0.762000,0.004468};
+Point(71) = {0.837375,0.762000,0.003901};
+Point(72) = {0.828166,0.762000,0.003223};
+Point(73) = {0.818956,0.762000,0.002383};
+Point(74) = {0.814351,0.762000,0.001728};
+Point(75) = {0.812509,0.762000,0.001356};
+Point(76) = {0.811589,0.762000,0.001120};
+Point(77) = {0.809747,0.762000,0.000000};
+Point(78) = {0.811589,0.762000,-0.001120};
+Point(79) = {0.812509,0.762000,-0.001356};
+Point(80) = {0.814351,0.762000,-0.001728};
+Point(81) = {0.818956,0.762000,-0.002383};
+Point(82) = {0.828166,0.762000,-0.003223};
+Point(83) = {0.837375,0.762000,-0.003901};
+Point(84) = {0.846585,0.762000,-0.004468};
+Point(85) = {0.865004,0.762000,-0.005375};
+Point(86) = {0.883423,0.762000,-0.006060};
+Point(87) = {0.901842,0.762000,-0.006587};
+Point(88) = {0.920261,0.762000,-0.006970};
+Point(89) = {0.938680,0.762000,-0.007228};
+Point(90) = {0.957099,0.762000,-0.007357};
+Point(91) = {0.975518,0.762000,-0.007353};
+Point(92) = {0.993937,0.762000,-0.007198};
+Point(93) = {1.012356,0.762000,-0.006881};
+Point(94) = {1.030775,0.762000,-0.006421};
+Point(95) = {1.049194,0.762000,-0.005843};
+Point(96) = {1.067613,0.762000,-0.005165};
+Point(97) = {1.086032,0.762000,-0.004402};
+Point(98) = {1.104451,0.762000,-0.003562};
+Point(99) = {1.122870,0.762000,-0.002685};
+Point(100) = {1.141290,0.762000,-0.001805};
+Point(101) = {1.159709,0.762000,-0.000921};
+
+/// Lines
+// Airfoil 1:
+Spline(1) = {1:26};
+Spline(2) = {26:50,1};
+// Airfoil 2:
+Spline(3) = {52:77};
+Spline(4) = {77:101,52};
+// Airfoil 1 to airfoil 2:
+Line(61) = {1,52};
+Line(62) = {26,77};
+
+
+/// Line loops & Surfaces
+// Wing 1:
+Line Loop(11) = {1,62,-3,-61};
+Line Loop(12) = {2,61,-4,-62};
+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 Surface{11,12};
+Recombine Surface{11,12};
+
+//// PHYSICAL GROUPS
+// Trailing edge
+Physical Line("wTe") = {61};
+// Wing:
+Physical Surface("wing") = {11,12};

fpm/models/n0012.py

@@ -31,7 +31,7 @@ def main(p, tail=False, field=False):
     print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
     msh = gmsh.MeshLoader('n0012.geo', __file__).execute(**p['pars'])
     # problem
-    pbl = fpm.Problem(msh, p['aoa'], p['aos'], p['mach'], p['spn'], 1, 0, 0, 0)
+    pbl = fpm.Problem(msh, p['aoa'], p['aos'], p['mach'], p['spn'], 1, 0, 0, 0, bool(p['sym']))
     wing = fpm.Body(msh, 'wing', ['wTe'], 10)
     pbl.add(wing)
     if tail:
@@ -52,6 +52,8 @@ def main(p, tail=False, field=False):
     tms['total'].stop()
     print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
     print(tms)
+    # return
+    return slv
 
 if __name__ == '__main__':
     main()

fpm/models/oneraM6.geo

+/* Onera M6 wing */
+
+// Mesh parameters
+// domain and mesh
+DefineConstant[ wNc = { 50, Name "Number of panel along chord" }
+                wNs = { 10, Name "Number of panel along span" }
+                bC = { 0.05, Name "Progression along chord" }
+                pS = { 1.0, Name "Progression along span" } ];
+
+//// GEOMETRY
+/// Points
+// Airfoil 1: oneraM6, 143 points 
+Point(1) = {0.805900,0.000000,0.000000};
+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};
+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};
+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};
+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};
+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};
+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};
+// Airfoil 2: oneraM6, 143 points 
+Point(144) = {1.143427,1.196000,0.000000};
+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};
+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};
+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};
+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};
+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};
+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};
+
+/// Lines
+// Airfoil 1:
+Spline(1) = {1:72};
+Spline(2) = {72:142,1};
+// Airfoil 2:
+Spline(3) = {144:215};
+Spline(4) = {215:285,144};
+// Airfoil 1 to airfoil 2:
+Line(61) = {1,144};
+Line(62) = {72,215};
+
+/// Line loops & Surfaces
+// Wing 1:
+Line Loop(11) = {1,62,-3,-61};
+Line Loop(12) = {2,61,-4,-62};
+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 Surface{11,12};
+Recombine Surface{11,12};
+
+//// PHYSICAL GROUPS
+// Trailing edge
+Physical Line("wTe") = {61};
+// Wing:
+Physical Surface("wing") = {11,12};

fpm/src/fBody.cpp

@@ -30,17 +30,6 @@ using namespace fpm;
 Body::Body(std::shared_ptr<MshData> _msh, std::string const &id,
            std::vector<std::string> const &teIds, double xF) : Group(_msh, id), Cl(0), Cd(0), Cs(0), Cm(0)
 {
-    // Get nodes
-    for (auto e : tag->elems)
-        for (auto n : e->nodes)
-            nodes.push_back(n);
-    std::sort(nodes.begin(), nodes.end());
-    nodes.erase(std::unique(nodes.begin(), nodes.end()), nodes.end());
-    // Associate each node to its adjacent elements
-    for (auto e : tag->elems)
-        for (auto n : e->nodes)
-            neMap[n].push_back(e);
-
     // If wakes are requested, check if they are already available, otherwise create them
     bool hasChanged = false;
     size_t j = 0;
@@ -85,14 +74,36 @@ Body::Body(std::shared_ptr<MshData> _msh, std::string const &id,
                 wkNodes.push_back(nN);
                 msh->nodes.push_back(nN);
             }
-            // Create wake elements
+            // Create wake elements and wake
             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->elems.push_back(new Quad4(msh->elems.back()->no + 1, tagp, tage, tag->elems[0]->parts, qnodes));
             }
-            // Create wake
             wakes.push_back(new Wake(_msh, tag->name + "Wake" + std::to_string(j), tag));
+            // Duplicate TE nodes and swap thos of lower wing elements
+            std::map<Node *, Node *> teMap;
+            for (auto n : teNodes)
+            {
+                Node *nN = new Node(msh->nodes.back()->no + 1, n->pos);
+                teMap[n] = nN;
+                msh->nodes.push_back(nN);
+            }
+            for (auto p : wakes.back()->wkMap)
+            {
+                Element *lwe = p.second.second;
+                for (size_t i = 0; i < lwe->nodes.size(); ++i)
+                {
+                    try
+                    {
+                        lwe->nodes[i] = teMap.at(lwe->nodes[i]);
+                    }
+                    catch (const std::out_of_range &)
+                    {
+                        //std::cout << lwe->nodes[i]->no << "not found in map!\n";
+                    }
+                }
+            }
             std::cout << *wakes.back() << " created. " << std::flush;
         }
         ++j;
@@ -104,6 +115,17 @@ Body::Body(std::shared_ptr<MshData> _msh, std::string const &id,
         writer.save(_msh->name);
     }
 
+    // Get nodes (now that they have been duplicated at lower TE)
+    for (auto e : tag->elems)
+        for (auto n : e->nodes)
+            nodes.push_back(n);
+    std::sort(nodes.begin(), nodes.end());
+    nodes.erase(std::unique(nodes.begin(), nodes.end()), nodes.end());
+    // Associate each node to its adjacent elements
+    for (auto e : tag->elems)
+        for (auto n : e->nodes)
+            neMap[n].push_back(e);
+
     // Size load vectors
     cLoadX.resize(nodes.size());
     cLoadY.resize(nodes.size());

fpm/src/fBuilder.h

@@ -24,18 +24,21 @@
 #include <vector>
 #include <memory>
 #include <Eigen/Dense>
+#include <map>
 
 namespace fpm
 {
 
 /**
  * @brief Aerodynamic Influence Coefficients builder
- * @authors Adrien Crovato
- * @todo add symmetry support
- * @todo implement AIC computation
+ * @authors Adrien Crovato and Axel Dechamps
+ * @todo implement field AIC computation
  */
 class FPM_API Builder : public fwk::wSharedObject
 {
+
+std::map<tbox::Element *, int> rows; ///< map linking an element to a matrix cell
+
 public:
     std::shared_ptr<Problem> pbl; ///< problem definition
 
@@ -59,8 +62,11 @@ public:
 #endif
 
 private:
-    double mu(tbox::Element *ei, tbox::Element *ej);
-    double tau(tbox::Element *ei, tbox::Element *ej);
+    Eigen::Matrix3d Glob2loc(bool symy, tbox::Element *ej);
+    Eigen::MatrixXd TrsfCoordSrc(bool symy, tbox::Element *ej, Eigen::Matrix3d const &glob2loc);
+    Eigen::Vector3d TrsfCoordTgt(bool symy, tbox::Element *ei, tbox::Element *ej, Eigen::Matrix3d const &glob2loc);
+    Eigen::Vector3d Distance(bool symy, tbox::Element *ei, tbox::Element *ej, Eigen::Matrix3d const &glob2loc);
+    Eigen::Vector2d MuTau(bool symy, tbox::Element *ei, tbox::Element *ej, Eigen::MatrixXd const &trsfCoordSrc, Eigen::Vector3d const &trsfCoordTgt, Eigen::Vector3d const &distance);
     double sigma(tbox::Element *ei, tbox::Element *ej);
     double sigmaX(tbox::Element *ei, tbox::Element *ej);
     double sigmaY(tbox::Element *ei, tbox::Element *ej);

fpm/src/fProblem.cpp

@@ -15,17 +15,20 @@
  */
 
 #include "fProblem.h"
+#include "fBuilder.h"
 #include "fField.h"
 #include "fBody.h"
+#include "fWake.h"
 #include "wTag.h"
 #include "wElement.h"
+#include "wNode.h"
 #include "wMem.h"
 #include "wCache.h"
 using namespace tbox;
 using namespace fpm;
 
 Problem::Problem(std::shared_ptr<MshData> _msh, double aoa, double aos, double mch,
-                 double sref, double cref, double xref, double yref, double zref) : msh(_msh), alpha(aoa), beta(aos), mach(mch), sR(sref), cR(cref)
+                 double sref, double cref, double xref, double yref, double zref, bool symy) : msh(_msh), alpha(aoa), beta(aos), mach(mch), sR(sref), cR(cref), symY(symy)
 {
     xR(0) = xref;
     xR(1) = yref;
@@ -39,6 +42,7 @@ void Problem::set(std::shared_ptr<Field> f)
 {
     field = f;
 }
+
 /**
  * @brief Add a non-lifting body
  */
@@ -50,21 +54,129 @@ void Problem::add(std::shared_ptr<Body> b)
 /**
  * @brief Compute velocity at element
  */
-Eigen::Vector3d Problem::U(tbox::Element &e, std::vector<double> const &phi)
+Eigen::Vector3d Problem::U(tbox::Element &e, std::vector<double> const &mu, double const &tau)
 {
-    Eigen::Vector3d v(0, 0, 0);
-    // @todo to be implemented if usefull
-    return v;
-    std::cout << std::endl;
+    Eigen::Vector3d U(0, 0, 0);
+    size_t nVertices = e.nodes.size();                          // Number of nodes per panel
+    size_t nGP = e.getVCache().getVGauss().getN();              // Number of Gauss points on the element
+    Eigen::MatrixXd coords(3,nVertices);
+    for (int i = 0; i < 3; ++i)
+        for (int j = 0; j < nVertices; ++j)
+            coords(i,j) = e.nodes[(nVertices - 1) - j]->pos[i]; // Coordinates of the panel
+        
+    // Import of the shape functions at the Gauss points
+    Eigen::MatrixXd ff(nVertices,nGP);
+    for(int i = 0; i < nVertices; ++i)
+        for(int j = 0; j < nGP; ++j)
+            ff(i,j) = e.getVCache().getSf(j)[i];
+
+    // Import of the derivatives of the shape functions at the Gauss points
+    Eigen::MatrixXd dff(2 * nGP,nVertices);
+    for(int i = 0; i < nGP; ++i)
+        for(int j = 0; j < 2; ++j)
+            dff.row(j + (2 * i)) = e.getVCache().getDsf(i).row(j);
+    
+    // Computation of the Jacobian matrix - xi-derivative (first row of the Jacobian)
+    Eigen::MatrixXd dXi(3,nGP);
+    for(int i = 0; i < nGP; ++i)
+    {
+        dXi.col(i) << 0, 0, 0;
+
+        for(int j = 0; j < nVertices; ++j)
+        {
+            dXi(0,i) += dff((2 * i),j) * coords(0,j);
+            dXi(1,i) += dff((2 * i),j) * coords(1,j);
+            dXi(2,i) += dff((2 * i),j) * coords(2,j);
+        }
+    }
+
+    // Computation of the Jacobian matrix - eta-derivative (second row of the Jacobian)
+    Eigen::MatrixXd dEta(3,nGP);
+    for(int i = 0; i < nGP; ++i)
+    {
+        dEta.col(i) << 0, 0, 0;
+
+        for(int j = 0; j < nVertices; ++j)
+        {
+            dEta(0,i) += dff(1 + (2 * i),j) * coords(0,j);
+            dEta(1,i) += dff(1 + (2 * i),j) * coords(1,j);
+            dEta(2,i) += dff(1 + (2 * i),j) * coords(2,j);
+        }
+    }
+
+    // Computation of the Jacobian matrix - cross-product (third row of the Jacobian)
+    Eigen::MatrixXd dZeta(3,nGP);
+
+    for(int i = 0; i < nGP; ++i)
+    {
+        Eigen::Vector3d dxi = dXi.col(i);
+        Eigen::Vector3d deta = dEta.col(i);
+        Eigen::Vector3d dzeta;
+        dzeta = (dxi.cross(deta)).normalized();
+        dZeta.col(i) = dzeta; 
+    }
+
+    // Assembly of the Jacobian and partial computation of the velocity
+    // Done for each Gauss point -> For n Gauss points, there are n values of the velocity -> Average
+    Eigen::MatrixXd v(3,nGP);
+    for(int i = 0; i < nGP; ++i)
+    {
+        // Creation of the Jacobian
+        Eigen::Matrix3d J0;
+        J0.row(0) = dXi.col(i).transpose();
+        J0.row(1) = dEta.col(i).transpose();
+        J0.row(2) = dZeta.col(i).transpose();
+
+        // Inversion of the Jacobian
+        Eigen::Matrix3d Jinv = J0.inverse();
+
+        // Removal of the third column 
+        Eigen::MatrixXd J(3,2);
+        for(int j = 0; j < 2; ++j)
+        {
+            J.col(j) = Jinv.col(j);
+        }
+
+        // Import of the doublets at the nodes of the panel
+        Eigen::VectorXd mu_element(nVertices);
+        for (int j = 0; j < nVertices; ++j)
+        {
+            mu_element(j) = mu[e.nodes[j]->no - 1];
+        }
+
+        // Computation of the velocity
+        Eigen::MatrixXd dff0(2,4);
+        dff0.row(0) = dff.row(2 * i);
+        dff0.row(1) = dff.row(1 + (2 * i));
+        v.col(i) = J * dff0 * mu_element;
+    }
+
+    // Average of the velocity
+    for(int i = 0; i < nGP; ++i)
+    {
+        U += v.col(i);
+    }
+    U = U/nGP;
+
+    // Import of the unit normal vector - Target panel
+    Eigen::Vector3d n = e.normal(); 
+
+    // Contribution of the normal velocity
+    U += (-tau)*n;
+
+    // Contribution of the freestream velocity
+    U += Ui();
+    
+    return U;
 }
 
 /**
  * @brief Compute density at element
  */
-double Problem::Rho(tbox::Element &e, std::vector<double> const &phi)
+double Problem::Rho(Eigen::Vector3d U)
 {
     // particularized: pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), 1 / (gamma - 1));
-    Eigen::Vector3d u = U(e, phi);
+    Eigen::Vector3d u = U;
     double a = 1 + 0.2 * mach * mach * (1 - u.dot(u));
     return sqrt(a * a * a * a * a);
 }
@@ -72,14 +184,14 @@ double Problem::Rho(tbox::Element &e, std::vector<double> const &phi)
 /**
  * @brief Compute Mach number at element
  */
-double Problem::Mach(tbox::Element &e, std::vector<double> const &phi)
+double Problem::Mach(Eigen::Vector3d U)
 {
     if (mach == 0)
         return 0;
     else
     {
         // particularized: pow(1 + 0.5 * (gamma - 1) * mInf * mInf * (1 - gradU2), 1 / (gamma - 1));
-        Eigen::Vector3d u = U(e, phi);
+        Eigen::Vector3d u = U;
         double a2 = 1 / (mach * mach) + 0.2 * (1 - u.dot(u)); // speed of sound squared
         return u.norm() / sqrt(a2);
     }
@@ -88,9 +200,9 @@ double Problem::Mach(tbox::Element &e, std::vector<double> const &phi)
 /**
  * @brief Compute pressure coefficient at element
  */
-double Problem::Cp(tbox::Element &e, std::vector<double> const &phi)
+double Problem::Cp(Eigen::Vector3d U)
 {
-    Eigen::Vector3d u = U(e, phi);
+    Eigen::Vector3d u = U;
     if (mach == 0)
         return 1 - u.dot(u);
     else
@@ -112,8 +224,15 @@ void Problem::updateMem()
             e->getVMem().update(true);
     // Update surface Jacobian
     for (auto s : bodies)
+    {
         for (auto e : s->tag->elems)
+        {
             e->getVMem().update(false);
+        }
+        for(auto w : s-> wakes)
+            for (auto e : w->tag->elems)
+                e->getVMem().update(false);
+    }
 }
 
 /**

fpm/src/fProblem.h

@@ -27,8 +27,7 @@ namespace fpm
 
 /**
  * @brief Manage the problem
- * @authors Adrien Crovato
- * @todo check and implement velocity computation
+ * @authors Adrien Crovato & Axel Dechamps
  */
 class FPM_API Problem : public fwk::wSharedObject
 {
@@ -45,8 +44,9 @@ public:
     double sR;          ///< Reference surface
     double cR;          ///< Reference chord
     Eigen::Vector3d xR; ///< Reference center point (for moment computation)
+    bool symY;          ///< True if mesh is symmetric with respect to y-plane
 
-    Problem(std::shared_ptr<tbox::MshData> _msh, double aoa, double aos, double mch, double sref, double cref, double xref, double yref, double zref);
+    Problem(std::shared_ptr<tbox::MshData> _msh, double aoa, double aos, double mch, double sref, double cref, double xref, double yref, double zref, bool symy);
     ~Problem() { std::cout << "~Problem()\n"; }
 
     void set(std::shared_ptr<Field> f);
@@ -55,10 +55,10 @@ public:
     // Functions
     inline double PhiI(Eigen::Vector3d const &pos);
     inline Eigen::Vector3d Ui();
-    Eigen::Vector3d U(tbox::Element &e, std::vector<double> const &phi);
-    double Rho(tbox::Element &e, std::vector<double> const &phi);
-    double Mach(tbox::Element &e, std::vector<double> const &phi);
-    double Cp(tbox::Element &e, std::vector<double> const &phi);
+    Eigen::Vector3d U(tbox::Element &e, std::vector<double> const &mu, double const &tau);
+    double Rho(Eigen::Vector3d U);
+    double Mach(Eigen::Vector3d U);
+    double Cp(Eigen::Vector3d U);
 
 #ifndef SWIG
     void check() const;

fpm/src/fSolver.cpp

@@ -36,7 +36,7 @@ using namespace fpm;
 
 /**
  * @brief Initialize the solver and perform sanity checks
- * @authors Adrien Crovato
+ * @authors Adrien Crovato & Axel Dechamps
  */
 Solver::Solver(std::shared_ptr<Builder> _aic) : aic(_aic), Cl(0), Cd(0), Cs(0), Cm(0)
 {
@@ -56,8 +56,10 @@ Solver::Solver(std::shared_ptr<Builder> _aic) : aic(_aic), Cl(0), Cd(0), Cs(0),
     // Setup variables
     phi.resize(aic->pbl->msh->nodes.size(), 0.);
     vPhi.resize(aic->pbl->msh->nodes.size(), 0.);
+    muNodes.resize(aic->pbl->msh->nodes.size(), 0.);
     U.resize(aic->pbl->msh->elems.size(), Eigen::Vector3d::Zero());
     rho.resize(aic->pbl->msh->elems.size(), 0.);
+    taue.resize(aic->pbl->msh->elems.size(), 0.);
     M.resize(aic->pbl->msh->elems.size(), 0.);
     Cp.resize(aic->pbl->msh->elems.size(), 0.);
 }
@@ -84,12 +86,18 @@ void Solver::run()
     Eigen::VectorXd mu(aic->nP), tau(aic->nP);
     int ig = 0;
     for (auto body : aic->pbl->bodies)
+    {
         for (int i = 0; i < body->tag->elems.size(); ++i, ++ig)
+        {
             tau(ig) = (aic->pbl->Ui() + Eigen::Vector3d(uX(ig), uY(ig), uZ(ig))).dot(body->tag->elems[i]->normal());
+            taue[body->tag->elems[i]->no-1] = aic->pbl->Ui().dot(body->tag->elems[i]->normal());
+        }
+    }
+
     std::cout << "done!" << std::endl;
 
     std::cout << "Solving for body doublets... " << std::flush;
-    Eigen::PartialPivLU<Eigen::MatrixXd>(aic->A).solve(aic->B * tau);
+    mu = Eigen::PartialPivLU<Eigen::MatrixXd>(aic->A).solve(aic->B * tau);
     std::cout << "done!" << std::endl;
 
     std::cout << "Computing flow and loads on bodies... " << std::flush;
@@ -110,6 +118,8 @@ void Solver::save(std::shared_ptr<MshExport> mshWriter, int n)
     Results results;
     results.scalars_at_nodes["phi"] = &phi;
     results.scalars_at_nodes["varPhi"] = &vPhi;
+    results.scalars_at_elems["mu"] = &mue;
+    results.scalars_at_elems["tau"] = &taue;
     results.vectors_at_elems["U"] = &U;
     results.scalars_at_elems["rho"] = &rho;
     results.scalars_at_elems["Mach"] = &M;
@@ -134,43 +144,57 @@ void Solver::save(std::shared_ptr<MshExport> mshWriter, int n)
  */
 void Solver::computeFlow(const Eigen::VectorXd &mu, const Eigen::VectorXd &tau, const Eigen::VectorXd &sigma)
 {
-    // @todo should we
-    // - compute velocity from phi at nodes or mu at nodes? if mu, then store phi at elements only
+    // Compute phi and mu at elements and transfer at nodes
+    Eigen::VectorXd phie = -aic->A * mu - aic->B * tau - aic->C * sigma;
+
+    mue.resize(mu.size());
+    for(size_t i = 0 ; i < mue.size() ; ++i)
+        mue[i] = mu(i);
 
-    // Compute phi at elements and transfer at nodes
-    Eigen::VectorXd phie = aic->A * mu + aic->B * tau + aic->C * sigma;
     auto pbl = aic->pbl;
     // reset all values
     for (auto n : pbl->msh->nodes)
+    {
         vPhi[n->no - 1] = 0;
+        muNodes[n->no - 1] = 0;
+    }
     int ig = 0;
     for (auto body : pbl->bodies)
     {
-        // associate the potential value to its element
+        // associate the potential and doublet value to its element
         std::map<Element *, double> mapPhi;
+        std::map<Element *, double> mapMu;
         for (auto e : body->tag->elems)
         {
             mapPhi[e] = phie[ig];
+            mapMu[e] = mu[ig];
             ig++;
         }
         // average at nodes
         for (auto n : body->nodes)
-            for (auto e : body->neMap.at(n))
-                vPhi[n->no - 1] += mapPhi.at(e) / body->neMap.at(n).size(); // maybe use area average?
+        {
+            for (auto e : body->neMap.at(n)) 
+            {
+                vPhi[n->no - 1] += mapPhi.at(e) / body->neMap.at(n).size(); 
+                muNodes[n->no - 1] += mapMu.at(e) / body->neMap.at(n).size(); 
+            }
+        }
     }
 
     // Compute surface flow
     for (auto body : pbl->bodies)
     {
         for (auto n : body->nodes)
-            phi[n->no - 1] = pbl->PhiI(n->pos) + vPhi[n->row]; // total potential
+            phi[n->no - 1] = pbl->PhiI(n->pos) + vPhi[n->no - 1]; // total potential
 
+        int ig = 0;
         for (auto e : body->tag->elems)
         {
-            U[e->no - 1] = pbl->U(*e, phi);     // velocity
-            rho[e->no - 1] = pbl->Rho(*e, phi); // density
-            M[e->no - 1] = pbl->Mach(*e, phi);  // Mach number
-            Cp[e->no - 1] = pbl->Cp(*e, phi);   // pressure coefficient
+            U[e->no - 1] = pbl->U(*e, muNodes, tau(ig));    // velocity
+            rho[e->no - 1] = pbl->Rho(U[e->no - 1]);        // density
+            M[e->no - 1] = pbl->Mach(U[e->no - 1]);         // Mach number
+            Cp[e->no - 1] = pbl->Cp(U[e->no - 1]);          // pressure coefficient
+            ig += 1;
         }
     }
 }
@@ -205,7 +229,7 @@ void Solver::computeLoad()
         // Compute aerodynamic load coefficients (i.e. Load / pressure / surface)
         // compute coefficients along x and vertical (y in 2D, z in 3D) directions and pitching moment (along -z in 2D, y in 3D)
         body->Cm = 0;
-        Eigen::Vector3d Cf(0, 0, 0);
+        Eigen::Vector3d Cf(0, 0, 0);        
         for (auto e : body->tag->elems)
         {
             Eigen::Vector3d l = e->getVMem().getCg() - pbl->xR;                       // lever arm

fpm/src/fSolver.h

@@ -30,8 +30,7 @@ namespace fpm
 
 /**
  * @brief Field panel solver
- * @authors Adrien Crovato
- * @todo check and implement velocity computation
+ * @authors Adrien Crovato & Axel Dechamps
  */
 class FPM_API Solver : public fwk::wSharedObject
 {
@@ -40,6 +39,9 @@ public:
 
     std::vector<double> phi;        ///< full potential
     std::vector<double> vPhi;       ///< perturbation potential
+    std::vector<double> taue;       ///< body source at the element
+    std::vector<double> muNodes;    ///< body doublet averaged at the nodes
+    std::vector<double> mue;        ///< body doublet averaged at the elements
     std::vector<Eigen::Vector3d> U; ///< velocity
     std::vector<double> rho;        ///< density
     std::vector<double> M;          ///< mach number

fpm/src/fWake.cpp

@@ -65,7 +65,7 @@ Wake::Wake(std::shared_ptr<MshData> _msh, std::string const &name, Tag *const &w
         if (ed->el1 == NULL || ed->el2 == NULL)
         {
             std::stringstream err;
-            err << "fpm::Lifting: no wing element from tag " << *wTag << "could be associated to wake element " << *(ed->el0) << "!\n";
+            err << "fpm::Wake: no wing element from tag " << *wTag << "could be associated to wake element " << *(ed->el0) << "!\n";
             throw std::runtime_error(err.str());
         }
         wkMap[ed->el0] = std::pair<Element *, Element *>(ed->el1, ed->el2);

fpm/tests/M6.py

+#!/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.
+
+# Onera M6 wing
+# Axel Dechamps
+
+from fwk.testing import *
+from fwk.coloring import ccolors
+import math
+import fpm
+import tbox
+import tbox.gmsh as gmsh
+
+def main():
+    # geometry parameters
+    S_ref = 1.5057      # Reference surface (full if symY = 0, half if symY = 1)
+    c_ref = 0.805       # Reference chord
+    wNc = 100           # Number of panels along the chord
+    wNs = 10            # Number of panels along the span
+    bC = 1.0            # Progression along chord
+    pS = 1.0            # Progression along span
+
+    # flow parameters
+    AoA = 3.06 * math.pi / 180  # Angle of attack (numeric value given in degrees)
+    Mach = 0                    # Mach number
+
+    # timer
+    tms = fpm.Timers()
+    tms['total'].start()
+
+    # mesh
+    print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+    pars = {'wNc' : wNc, 'wNs' : wNs, 'bC' : bC, 'pS' : pS}
+    msh = gmsh.MeshLoader('../models/oneraM6.geo', __file__).execute(**pars)
+
+    # problem
+    pbl = fpm.Problem(msh, AoA, 0., Mach, S_ref, c_ref, 0, 0, 0, True)
+    wing = fpm.Body(msh, 'wing', ['wTe'], 10)
+    pbl.add(wing)
+
+    # AIC builder
+    aic = fpm.Builder(pbl)
+    aic.run()
+
+    # solver
+    slv = fpm.Solver(aic)
+    slv.run()
+    slv.save(tbox.GmshExport(msh))
+
+    # end timer
+    tms['total'].stop()
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print(tms)
+
+    # test
+    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+    tests = CTests()
+    if AoA == 3.06 * math.pi / 180 and Mach == 0:
+        tests.add(CTest('CL', slv.Cl, 0.0998, 5e-2))
+        tests.add(CTest('CD', slv.Cd, 0.0020, 5e-2))
+        tests.add(CTest('CS', slv.Cs, 0.0037, 5e-2))
+        tests.add(CTest('CM', slv.Cm, -0.0580, 5e-2))
+    tests.run()
+
+    # eof
+    print('')
+
+if __name__ == '__main__':
+    main()

fpm/tests/a445.py

+#!/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.
+
+# Agard 445 wing
+# Axel Dechamps
+
+from fwk.testing import *
+from fwk.coloring import ccolors
+import math
+import fpm
+import tbox
+import tbox.gmsh as gmsh
+
+def main():    
+    # geometry parameters
+    S_ref = 0.35        # Reference surface (full if symY = 0, half if symY = 1)
+    c_ref = 0.559       # Reference chord
+    wNc = 100           # Number of panels along the chord
+    wNs = 10            # Number of panels along the span
+    bC = 1.0            # Progression along chord
+    pS = 1.0            # Progression along span
+    
+    # flow parameters
+    AoA = 1 * math.pi / 180  # Angle of attack (numeric value given in degrees)
+    Mach = 0                 # Mach number
+
+    # timer
+    tms = fpm.Timers()
+    tms['total'].start()
+
+    # mesh
+    print(ccolors.ANSI_BLUE + 'PyMeshing...' + ccolors.ANSI_RESET)
+    pars = {'wNc' : wNc, 'wNs' : wNs, 'bC' : bC, 'pS' : pS}
+    msh = gmsh.MeshLoader('../models/agard445.geo', __file__).execute(**pars)
+
+    # problem
+    pbl = fpm.Problem(msh, AoA, 0., Mach, S_ref, c_ref, 0, 0, 0, True)
+    wing = fpm.Body(msh, 'wing', ['wTe'], 10)
+    pbl.add(wing)
+
+    # AIC builder
+    aic = fpm.Builder(pbl)
+    aic.run()
+
+    # solver
+    slv = fpm.Solver(aic)
+    slv.run()
+    slv.save(tbox.GmshExport(msh))
+
+    # end timer
+    tms['total'].stop()
+    print(ccolors.ANSI_BLUE + 'PyTiming...' + ccolors.ANSI_RESET)
+    print(tms)
+
+    # test
+    print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
+    tests = CTests()
+    if AoA == 1 * math.pi / 180 and Mach == 0:
+        tests.add(CTest('CL', slv.Cl, 0.054790, 5e-2)) 
+        tests.add(CTest('CD', slv.Cd, 0.000221, 5e-2)) 
+        tests.add(CTest('CS', slv.Cs, 0.001184, 5e-2))
+        tests.add(CTest('CM', slv.Cm, -0.048214, 5e-2))
+    tests.run()
+    
+    # eof
+    print('')
+
+if __name__ == '__main__':
+    main()
+

fpm/tests/n12.py

@@ -19,22 +19,34 @@
 # Adrien Crovato
 
 from fwk.testing import *
+import math
 
 def main():
     p = {}
     # geometry parameters
-    p['pars'] = {'wSpn' : 5, 'wNc' : 80, 'wNs' : 10}
-    p['spn'] = 5
+    p['sym'] = 0
+    p['spn'] = 5 # true span if sym=0, half span if sym=1
+    p['pars'] = {'symY' : p['sym'],'wSpn' : p['spn'], 'wNc' : 100, 'wNs' : 10} 
     # flow parameters
-    p['aoa'] = 0
+    p['aoa'] = 3 * math.pi / 180 
     p['aos'] = 0
     p['mach'] = 0
     # run
     import fpm.models.n0012 as n12
-    n12.main(p)
+    solver = n12.main(p)
     # test
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
     tests = CTests()
+    if(p['aoa'] == 0 * math.pi / 180 and p['mach'] == 0):
+        tests.add(CTest('CL', solver.Cl, 0.00000, 5e-2)) # aero (old code) = 0.00000
+        tests.add(CTest('CD', solver.Cd, 0.00021, 5e-2)) # aero (old code) = 0.00013
+        tests.add(CTest('CS', solver.Cs, -0.00000, 5e-2))
+        tests.add(CTest('CM', solver.Cm, -0.00000, 5e-2))
+    if(p['aoa'] == 3 * math.pi / 180 and p['mach'] == 0):
+        tests.add(CTest('CL', solver.Cl, 0.229, 5e-2)) # aero (old code) = 0.226
+        tests.add(CTest('CD', solver.Cd, 0.003, 5e-2)) # aero (old code) = 0.003
+        tests.add(CTest('CS', solver.Cs, -0.000, 5e-2))
+        tests.add(CTest('CM', solver.Cm, -0.056, 5e-2))
     tests.run()
     
     # eof

fpm/tests/n12f.py

@@ -19,14 +19,16 @@
 # Adrien Crovato
 
 from fwk.testing import *
+import math
 
 def main():
     p = {}
     # geometry parameters
-    p['pars'] = {'wSpn' : 5, 'wNc' : 80, 'wNs' : 10, 'fWdt' : 10, 'nX' : 10, 'nY' : 10, 'nZ' : 4}
-    p['spn'] = 5
+    p['sym'] = 0
+    p['spn'] = 5 # true span if sym=0, half span if sym=1
+    p['pars'] = {'symY' : p['sym'], 'wSpn' : p['spn'], 'wNc' : 100, 'wNs' : 10, 'fWdt' : 10, 'nX' : 10, 'nY' : 10, 'nZ' : 4}
     # flow parameters
-    p['aoa'] = 0
+    p['aoa'] = 3 * math.pi / 180 
     p['aos'] = 0
     p['mach'] = 0.8
     # run

fpm/tests/n12t.py

@@ -19,22 +19,30 @@
 # Adrien Crovato
 
 from fwk.testing import *
+import math
 
 def main():
     p = {}
     # geometry parameters
-    p['pars'] = {'wSpn' : 5, 'wNc' : 80, 'wNs' : 10, 'tSpn' : 3, 'tNc' : 40, 'tNs' : 5}
-    p['spn'] = 5
+    p['sym'] = 0
+    p['spn'] = 5 # true span if sym=0, half span if sym=1
+    p['tspn'] = 3 # true span if sym=0, half span if sym=1
+    p['pars'] = {'symY' : p['sym'], 'wSpn' : p['spn'], 'wNc' : 100, 'wNs' : 10, 'tSpn' : p['tspn'], 'tNc' : 100, 'tNs' : 10}
     # flow parameters
-    p['aoa'] = 0
+    p['aoa'] = 3 * math.pi / 180 
     p['aos'] = 0
     p['mach'] = 0
     # run
     import fpm.models.n0012 as n12
-    n12.main(p, tail=True)
+    solver = n12.main(p, tail=True)
     # test
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)
     tests = CTests()
+    if(p['aoa'] == 3 * math.pi / 180 and p['mach'] == 0):
+        tests.add(CTest('CL', solver.Cl, 0.273, 5e-2))
+        tests.add(CTest('CD', solver.Cd, 0.004, 5e-2))
+        tests.add(CTest('CS', solver.Cs, -0.000, 5e-2))
+        tests.add(CTest('CM', solver.Cm, -0.277, 5e-2))
     tests.run()
     
     # eof