From bb3f0e0c2f18f6b0df2d6a59af40d106fda9c7dd Mon Sep 17 00:00:00 2001
From: Dechamps Axel <axel.dechamps@student.uliege.be>
Date: Tue, 1 Feb 2022 08:13:44 +0000
Subject: [PATCH] Master Thesis - Dechamps - Implementation of the linear part
---
fpm/models/agard445.geo | 144 ++++++++++++
fpm/models/n0012.geo | 445 +++++++++++++++++++------------------
fpm/models/n0012.py | 4 +-
fpm/models/oneraM6.geo | 327 ++++++++++++++++++++++++++++
fpm/src/fBody.cpp | 48 ++--
fpm/src/fBuilder.cpp | 470 +++++++++++++++++++++++++++++++++++++---
fpm/src/fBuilder.h | 16 +-
fpm/src/fProblem.cpp | 143 +++++++++++-
fpm/src/fProblem.h | 14 +-
fpm/src/fSolver.cpp | 54 +++--
fpm/src/fSolver.h | 6 +-
fpm/src/fWake.cpp | 2 +-
fpm/tests/M6.py | 83 +++++++
fpm/tests/a445.py | 84 +++++++
fpm/tests/n12.py | 20 +-
fpm/tests/n12f.py | 8 +-
fpm/tests/n12t.py | 16 +-
17 files changed, 1567 insertions(+), 317 deletions(-)
create mode 100644 fpm/models/agard445.geo
create mode 100644 fpm/models/oneraM6.geo
create mode 100644 fpm/tests/M6.py
create mode 100644 fpm/tests/a445.py
diff --git a/fpm/models/agard445.geo b/fpm/models/agard445.geo
new file mode 100644
index 0000000..c419f7a
--- /dev/null
+++ b/fpm/models/agard445.geo
@@ -0,0 +1,144 @@
+/* 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};
diff --git a/fpm/models/n0012.geo b/fpm/models/n0012.geo
index c1c0f80..1f70a5d 100644
--- a/fpm/models/n0012.geo
+++ b/fpm/models/n0012.geo
@@ -2,122 +2,139 @@
// Parameters
// domain and mesh
-DefineConstant[ wSpn = { 2.0, Name "Wing span" } ];
-DefineConstant[ wNc = { 80, Name "Number of cells along wing chord" } ];
-DefineConstant[ wNs = { 10, Name "Number of cells along wing span" } ];
-DefineConstant[ tSpn = { 0.0, Name "Tail span" } ];
-DefineConstant[ tNc = { 40, Name "Number of cells along tail chord" } ];
-DefineConstant[ tNs = { 6, Name "Number of cells along tail span" } ];
-DefineConstant[ fWdt = { 0.0, Name "Field width" } ];
-DefineConstant[ nX = { 10, Name "Number of field cells along X" } ];
-DefineConstant[ nY = { 10, Name "Number of field cells along Y" } ];
-DefineConstant[ nZ = { 10, Name "Number of field cells along Z" } ];
+DefineConstant[
+ symY = { 1, Choices{0, 1}, Name "Y-plane symmetry" }
+ wSpn = { 2.0, Name "Wing span" }
+ wNc = { 80, Name "Number of cells along wing chord" }
+ wNs = { 10, Name "Number of cells along wing span" }
+ tSpn = { 0.0, Name "Tail span" }
+ tNc = { 40, Name "Number of cells along tail chord" }
+ tNs = { 6, Name "Number of cells along tail span" }
+ fWdt = { 0.0, Name "Field width" }
+ nX = { 10, Name "Number of field cells along X" }
+ nY = { 10, Name "Number of field cells along Y" }
+ nZ = { 10, Name "Number of field cells along Z" }
+];
//// GEOMETRY
+If(!symY)
+ _wSpn = -wSpn/2;
+ _tSpn = -tSpn/2;
+ _fWdt = -fWdt/2;
+EndIf
+If(symY)
+ _wSpn = 0;
+ _tSpn = 0;
+ _fWdt = 0;
+ wSpn = 2*wSpn;
+ tSpn = 2*tSpn;
+ fWdt = 2*fWdt;
+EndIf
+
/// Points
// Wing
-Point(1) = {1.000000,-wSpn/2,0.000000};
-Point(2) = {0.999013,-wSpn/2,0.000143};
-Point(3) = {0.996057,-wSpn/2,0.000572};
-Point(4) = {0.991144,-wSpn/2,0.001280};
-Point(5) = {0.984292,-wSpn/2,0.002260};
-Point(6) = {0.975528,-wSpn/2,0.003501};
-Point(7) = {0.964888,-wSpn/2,0.004990};
-Point(8) = {0.952414,-wSpn/2,0.006710};
-Point(9) = {0.938153,-wSpn/2,0.008643};
-Point(10) = {0.922164,-wSpn/2,0.010770};
-Point(11) = {0.904508,-wSpn/2,0.013071};
-Point(12) = {0.885257,-wSpn/2,0.015523};
-Point(13) = {0.864484,-wSpn/2,0.018106};
-Point(14) = {0.842274,-wSpn/2,0.020795};
-Point(15) = {0.818712,-wSpn/2,0.023569};
-Point(16) = {0.793893,-wSpn/2,0.026405};
-Point(17) = {0.767913,-wSpn/2,0.029279};
-Point(18) = {0.740877,-wSpn/2,0.032168};
-Point(19) = {0.712890,-wSpn/2,0.035048};
-Point(20) = {0.684062,-wSpn/2,0.037896};
-Point(21) = {0.654508,-wSpn/2,0.040686};
-Point(22) = {0.624345,-wSpn/2,0.043394};
-Point(23) = {0.593691,-wSpn/2,0.045992};
-Point(24) = {0.562667,-wSpn/2,0.048455};
-Point(25) = {0.531395,-wSpn/2,0.050754};
-Point(26) = {0.500000,-wSpn/2,0.052862};
-Point(27) = {0.468605,-wSpn/2,0.054749};
-Point(28) = {0.437333,-wSpn/2,0.056390};
-Point(29) = {0.406309,-wSpn/2,0.057755};
-Point(30) = {0.375655,-wSpn/2,0.058819};
-Point(31) = {0.345492,-wSpn/2,0.059557};
-Point(32) = {0.315938,-wSpn/2,0.059947};
-Point(33) = {0.287110,-wSpn/2,0.059971};
-Point(34) = {0.259123,-wSpn/2,0.059614};
-Point(35) = {0.232087,-wSpn/2,0.058863};
-Point(36) = {0.206107,-wSpn/2,0.057712};
-Point(37) = {0.181288,-wSpn/2,0.056159};
-Point(38) = {0.157726,-wSpn/2,0.054206};
-Point(39) = {0.135516,-wSpn/2,0.051862};
-Point(40) = {0.114743,-wSpn/2,0.049138};
-Point(41) = {0.095492,-wSpn/2,0.046049};
-Point(42) = {0.077836,-wSpn/2,0.042615};
-Point(43) = {0.061847,-wSpn/2,0.038859};
-Point(44) = {0.047586,-wSpn/2,0.034803};
-Point(45) = {0.035112,-wSpn/2,0.030473};
-Point(46) = {0.024472,-wSpn/2,0.025893};
-Point(47) = {0.015708,-wSpn/2,0.021088};
-Point(48) = {0.008856,-wSpn/2,0.016078};
-Point(49) = {0.003943,-wSpn/2,0.010884};
-Point(50) = {0.000987,-wSpn/2,0.005521};
-Point(51) = {0.000000,-wSpn/2,0.000000};
-Point(52) = {0.000987,-wSpn/2,-0.005521};
-Point(53) = {0.003943,-wSpn/2,-0.010884};
-Point(54) = {0.008856,-wSpn/2,-0.016078};
-Point(55) = {0.015708,-wSpn/2,-0.021088};
-Point(56) = {0.024472,-wSpn/2,-0.025893};
-Point(57) = {0.035112,-wSpn/2,-0.030473};
-Point(58) = {0.047586,-wSpn/2,-0.034803};
-Point(59) = {0.061847,-wSpn/2,-0.038859};
-Point(60) = {0.077836,-wSpn/2,-0.042615};
-Point(61) = {0.095492,-wSpn/2,-0.046049};
-Point(62) = {0.114743,-wSpn/2,-0.049138};
-Point(63) = {0.135516,-wSpn/2,-0.051862};
-Point(64) = {0.157726,-wSpn/2,-0.054206};
-Point(65) = {0.181288,-wSpn/2,-0.056159};
-Point(66) = {0.206107,-wSpn/2,-0.057712};
-Point(67) = {0.232087,-wSpn/2,-0.058863};
-Point(68) = {0.259123,-wSpn/2,-0.059614};
-Point(69) = {0.287110,-wSpn/2,-0.059971};
-Point(70) = {0.315938,-wSpn/2,-0.059947};
-Point(71) = {0.345492,-wSpn/2,-0.059557};
-Point(72) = {0.375655,-wSpn/2,-0.058819};
-Point(73) = {0.406309,-wSpn/2,-0.057755};
-Point(74) = {0.437333,-wSpn/2,-0.056390};
-Point(75) = {0.468605,-wSpn/2,-0.054749};
-Point(76) = {0.500000,-wSpn/2,-0.052862};
-Point(77) = {0.531395,-wSpn/2,-0.050754};
-Point(78) = {0.562667,-wSpn/2,-0.048455};
-Point(79) = {0.593691,-wSpn/2,-0.045992};
-Point(80) = {0.624345,-wSpn/2,-0.043394};
-Point(81) = {0.654508,-wSpn/2,-0.040686};
-Point(82) = {0.684062,-wSpn/2,-0.037896};
-Point(83) = {0.712890,-wSpn/2,-0.035048};
-Point(84) = {0.740877,-wSpn/2,-0.032168};
-Point(85) = {0.767913,-wSpn/2,-0.029279};
-Point(86) = {0.793893,-wSpn/2,-0.026405};
-Point(87) = {0.818712,-wSpn/2,-0.023569};
-Point(88) = {0.842274,-wSpn/2,-0.020795};
-Point(89) = {0.864484,-wSpn/2,-0.018106};
-Point(90) = {0.885257,-wSpn/2,-0.015523};
-Point(91) = {0.904508,-wSpn/2,-0.013071};
-Point(92) = {0.922164,-wSpn/2,-0.010770};
-Point(93) = {0.938153,-wSpn/2,-0.008643};
-Point(94) = {0.952414,-wSpn/2,-0.006710};
-Point(95) = {0.964888,-wSpn/2,-0.004990};
-Point(96) = {0.975528,-wSpn/2,-0.003501};
-Point(97) = {0.984292,-wSpn/2,-0.002260};
-Point(98) = {0.991144,-wSpn/2,-0.001280};
-Point(99) = {0.996057,-wSpn/2,-0.000572};
-Point(100) = {0.999013,-wSpn/2,-0.000143};
+Point(1) = {1.000000,_wSpn,0.000000};
+Point(2) = {0.999013,_wSpn,0.000143};
+Point(3) = {0.996057,_wSpn,0.000572};
+Point(4) = {0.991144,_wSpn,0.001280};
+Point(5) = {0.984292,_wSpn,0.002260};
+Point(6) = {0.975528,_wSpn,0.003501};
+Point(7) = {0.964888,_wSpn,0.004990};
+Point(8) = {0.952414,_wSpn,0.006710};
+Point(9) = {0.938153,_wSpn,0.008643};
+Point(10) = {0.922164,_wSpn,0.010770};
+Point(11) = {0.904508,_wSpn,0.013071};
+Point(12) = {0.885257,_wSpn,0.015523};
+Point(13) = {0.864484,_wSpn,0.018106};
+Point(14) = {0.842274,_wSpn,0.020795};
+Point(15) = {0.818712,_wSpn,0.023569};
+Point(16) = {0.793893,_wSpn,0.026405};
+Point(17) = {0.767913,_wSpn,0.029279};
+Point(18) = {0.740877,_wSpn,0.032168};
+Point(19) = {0.712890,_wSpn,0.035048};
+Point(20) = {0.684062,_wSpn,0.037896};
+Point(21) = {0.654508,_wSpn,0.040686};
+Point(22) = {0.624345,_wSpn,0.043394};
+Point(23) = {0.593691,_wSpn,0.045992};
+Point(24) = {0.562667,_wSpn,0.048455};
+Point(25) = {0.531395,_wSpn,0.050754};
+Point(26) = {0.500000,_wSpn,0.052862};
+Point(27) = {0.468605,_wSpn,0.054749};
+Point(28) = {0.437333,_wSpn,0.056390};
+Point(29) = {0.406309,_wSpn,0.057755};
+Point(30) = {0.375655,_wSpn,0.058819};
+Point(31) = {0.345492,_wSpn,0.059557};
+Point(32) = {0.315938,_wSpn,0.059947};
+Point(33) = {0.287110,_wSpn,0.059971};
+Point(34) = {0.259123,_wSpn,0.059614};
+Point(35) = {0.232087,_wSpn,0.058863};
+Point(36) = {0.206107,_wSpn,0.057712};
+Point(37) = {0.181288,_wSpn,0.056159};
+Point(38) = {0.157726,_wSpn,0.054206};
+Point(39) = {0.135516,_wSpn,0.051862};
+Point(40) = {0.114743,_wSpn,0.049138};
+Point(41) = {0.095492,_wSpn,0.046049};
+Point(42) = {0.077836,_wSpn,0.042615};
+Point(43) = {0.061847,_wSpn,0.038859};
+Point(44) = {0.047586,_wSpn,0.034803};
+Point(45) = {0.035112,_wSpn,0.030473};
+Point(46) = {0.024472,_wSpn,0.025893};
+Point(47) = {0.015708,_wSpn,0.021088};
+Point(48) = {0.008856,_wSpn,0.016078};
+Point(49) = {0.003943,_wSpn,0.010884};
+Point(50) = {0.000987,_wSpn,0.005521};
+Point(51) = {0.000000,_wSpn,0.000000};
+Point(52) = {0.000987,_wSpn,-0.005521};
+Point(53) = {0.003943,_wSpn,-0.010884};
+Point(54) = {0.008856,_wSpn,-0.016078};
+Point(55) = {0.015708,_wSpn,-0.021088};
+Point(56) = {0.024472,_wSpn,-0.025893};
+Point(57) = {0.035112,_wSpn,-0.030473};
+Point(58) = {0.047586,_wSpn,-0.034803};
+Point(59) = {0.061847,_wSpn,-0.038859};
+Point(60) = {0.077836,_wSpn,-0.042615};
+Point(61) = {0.095492,_wSpn,-0.046049};
+Point(62) = {0.114743,_wSpn,-0.049138};
+Point(63) = {0.135516,_wSpn,-0.051862};
+Point(64) = {0.157726,_wSpn,-0.054206};
+Point(65) = {0.181288,_wSpn,-0.056159};
+Point(66) = {0.206107,_wSpn,-0.057712};
+Point(67) = {0.232087,_wSpn,-0.058863};
+Point(68) = {0.259123,_wSpn,-0.059614};
+Point(69) = {0.287110,_wSpn,-0.059971};
+Point(70) = {0.315938,_wSpn,-0.059947};
+Point(71) = {0.345492,_wSpn,-0.059557};
+Point(72) = {0.375655,_wSpn,-0.058819};
+Point(73) = {0.406309,_wSpn,-0.057755};
+Point(74) = {0.437333,_wSpn,-0.056390};
+Point(75) = {0.468605,_wSpn,-0.054749};
+Point(76) = {0.500000,_wSpn,-0.052862};
+Point(77) = {0.531395,_wSpn,-0.050754};
+Point(78) = {0.562667,_wSpn,-0.048455};
+Point(79) = {0.593691,_wSpn,-0.045992};
+Point(80) = {0.624345,_wSpn,-0.043394};
+Point(81) = {0.654508,_wSpn,-0.040686};
+Point(82) = {0.684062,_wSpn,-0.037896};
+Point(83) = {0.712890,_wSpn,-0.035048};
+Point(84) = {0.740877,_wSpn,-0.032168};
+Point(85) = {0.767913,_wSpn,-0.029279};
+Point(86) = {0.793893,_wSpn,-0.026405};
+Point(87) = {0.818712,_wSpn,-0.023569};
+Point(88) = {0.842274,_wSpn,-0.020795};
+Point(89) = {0.864484,_wSpn,-0.018106};
+Point(90) = {0.885257,_wSpn,-0.015523};
+Point(91) = {0.904508,_wSpn,-0.013071};
+Point(92) = {0.922164,_wSpn,-0.010770};
+Point(93) = {0.938153,_wSpn,-0.008643};
+Point(94) = {0.952414,_wSpn,-0.006710};
+Point(95) = {0.964888,_wSpn,-0.004990};
+Point(96) = {0.975528,_wSpn,-0.003501};
+Point(97) = {0.984292,_wSpn,-0.002260};
+Point(98) = {0.991144,_wSpn,-0.001280};
+Point(99) = {0.996057,_wSpn,-0.000572};
+Point(100) = {0.999013,_wSpn,-0.000143};
Point(102) = {1.000000,wSpn/2,0.000000};
Point(103) = {0.999013,wSpn/2,0.000143};
@@ -222,106 +239,106 @@ Point(201) = {0.999013,wSpn/2,-0.000143};
// Tail
If(tSpn != 0)
- Point(301) = {5+1.000000/2,-tSpn/2,0.5+0.000000/2};
- Point(302) = {5+0.999013/2,-tSpn/2,0.5+0.000143/2};
- Point(303) = {5+0.996057/2,-tSpn/2,0.5+0.000572/2};
- Point(304) = {5+0.991144/2,-tSpn/2,0.5+0.001280/2};
- Point(305) = {5+0.984292/2,-tSpn/2,0.5+0.002260/2};
- Point(306) = {5+0.975528/2,-tSpn/2,0.5+0.003501/2};
- Point(307) = {5+0.964888/2,-tSpn/2,0.5+0.004990/2};
- Point(308) = {5+0.952414/2,-tSpn/2,0.5+0.006710/2};
- Point(309) = {5+0.938153/2,-tSpn/2,0.5+0.008643/2};
- Point(310) = {5+0.922164/2,-tSpn/2,0.5+0.010770/2};
- Point(311) = {5+0.904508/2,-tSpn/2,0.5+0.013071/2};
- Point(312) = {5+0.885257/2,-tSpn/2,0.5+0.015523/2};
- Point(313) = {5+0.864484/2,-tSpn/2,0.5+0.018106/2};
- Point(314) = {5+0.842274/2,-tSpn/2,0.5+0.020795/2};
- Point(315) = {5+0.818712/2,-tSpn/2,0.5+0.023569/2};
- Point(316) = {5+0.793893/2,-tSpn/2,0.5+0.026405/2};
- Point(317) = {5+0.767913/2,-tSpn/2,0.5+0.029279/2};
- Point(318) = {5+0.740877/2,-tSpn/2,0.5+0.032168/2};
- Point(319) = {5+0.712890/2,-tSpn/2,0.5+0.035048/2};
- Point(320) = {5+0.684062/2,-tSpn/2,0.5+0.037896/2};
- Point(321) = {5+0.654508/2,-tSpn/2,0.5+0.040686/2};
- Point(322) = {5+0.624345/2,-tSpn/2,0.5+0.043394/2};
- Point(323) = {5+0.593691/2,-tSpn/2,0.5+0.045992/2};
- Point(324) = {5+0.562667/2,-tSpn/2,0.5+0.048455/2};
- Point(325) = {5+0.531395/2,-tSpn/2,0.5+0.050754/2};
- Point(326) = {5+0.500000/2,-tSpn/2,0.5+0.052862/2};
- Point(327) = {5+0.468605/2,-tSpn/2,0.5+0.054749/2};
- Point(328) = {5+0.437333/2,-tSpn/2,0.5+0.056390/2};
- Point(329) = {5+0.406309/2,-tSpn/2,0.5+0.057755/2};
- Point(330) = {5+0.375655/2,-tSpn/2,0.5+0.058819/2};
- Point(331) = {5+0.345492/2,-tSpn/2,0.5+0.059557/2};
- Point(332) = {5+0.315938/2,-tSpn/2,0.5+0.059947/2};
- Point(333) = {5+0.287110/2,-tSpn/2,0.5+0.059971/2};
- Point(334) = {5+0.259123/2,-tSpn/2,0.5+0.059614/2};
- Point(335) = {5+0.232087/2,-tSpn/2,0.5+0.058863/2};
- Point(336) = {5+0.206107/2,-tSpn/2,0.5+0.057712/2};
- Point(337) = {5+0.181288/2,-tSpn/2,0.5+0.056159/2};
- Point(338) = {5+0.157726/2,-tSpn/2,0.5+0.054206/2};
- Point(339) = {5+0.135516/2,-tSpn/2,0.5+0.051862/2};
- Point(340) = {5+0.114743/2,-tSpn/2,0.5+0.049138/2};
- Point(341) = {5+0.095492/2,-tSpn/2,0.5+0.046049/2};
- Point(342) = {5+0.077836/2,-tSpn/2,0.5+0.042615/2};
- Point(343) = {5+0.061847/2,-tSpn/2,0.5+0.038859/2};
- Point(344) = {5+0.047586/2,-tSpn/2,0.5+0.034803/2};
- Point(345) = {5+0.035112/2,-tSpn/2,0.5+0.030473/2};
- Point(346) = {5+0.024472/2,-tSpn/2,0.5+0.025893/2};
- Point(347) = {5+0.015708/2,-tSpn/2,0.5+0.021088/2};
- Point(348) = {5+0.008856/2,-tSpn/2,0.5+0.016078/2};
- Point(349) = {5+0.003943/2,-tSpn/2,0.5+0.010884/2};
- Point(350) = {5+0.000987/2,-tSpn/2,0.5+0.005521/2};
- Point(351) = {5+0.000000/2,-tSpn/2,0.5+0.000000/2};
- Point(352) = {5+0.000987/2,-tSpn/2,0.5-0.005521/2};
- Point(353) = {5+0.003943/2,-tSpn/2,0.5-0.010884/2};
- Point(354) = {5+0.008856/2,-tSpn/2,0.5-0.016078/2};
- Point(355) = {5+0.015708/2,-tSpn/2,0.5-0.021088/2};
- Point(356) = {5+0.024472/2,-tSpn/2,0.5-0.025893/2};
- Point(357) = {5+0.035112/2,-tSpn/2,0.5-0.030473/2};
- Point(358) = {5+0.047586/2,-tSpn/2,0.5-0.034803/2};
- Point(359) = {5+0.061847/2,-tSpn/2,0.5-0.038859/2};
- Point(360) = {5+0.077836/2,-tSpn/2,0.5-0.042615/2};
- Point(361) = {5+0.095492/2,-tSpn/2,0.5-0.046049/2};
- Point(362) = {5+0.114743/2,-tSpn/2,0.5-0.049138/2};
- Point(363) = {5+0.135516/2,-tSpn/2,0.5-0.051862/2};
- Point(364) = {5+0.157726/2,-tSpn/2,0.5-0.054206/2};
- Point(365) = {5+0.181288/2,-tSpn/2,0.5-0.056159/2};
- Point(366) = {5+0.206107/2,-tSpn/2,0.5-0.057712/2};
- Point(367) = {5+0.232087/2,-tSpn/2,0.5-0.058863/2};
- Point(368) = {5+0.259123/2,-tSpn/2,0.5-0.059614/2};
- Point(369) = {5+0.287110/2,-tSpn/2,0.5-0.059971/2};
- Point(370) = {5+0.315938/2,-tSpn/2,0.5-0.059947/2};
- Point(371) = {5+0.345492/2,-tSpn/2,0.5-0.059557/2};
- Point(372) = {5+0.375655/2,-tSpn/2,0.5-0.058819/2};
- Point(373) = {5+0.406309/2,-tSpn/2,0.5-0.057755/2};
- Point(374) = {5+0.437333/2,-tSpn/2,0.5-0.056390/2};
- Point(375) = {5+0.468605/2,-tSpn/2,0.5-0.054749/2};
- Point(376) = {5+0.500000/2,-tSpn/2,0.5-0.052862/2};
- Point(377) = {5+0.531395/2,-tSpn/2,0.5-0.050754/2};
- Point(378) = {5+0.562667/2,-tSpn/2,0.5-0.048455/2};
- Point(379) = {5+0.593691/2,-tSpn/2,0.5-0.045992/2};
- Point(380) = {5+0.624345/2,-tSpn/2,0.5-0.043394/2};
- Point(381) = {5+0.654508/2,-tSpn/2,0.5-0.040686/2};
- Point(382) = {5+0.684062/2,-tSpn/2,0.5-0.037896/2};
- Point(383) = {5+0.712890/2,-tSpn/2,0.5-0.035048/2};
- Point(384) = {5+0.740877/2,-tSpn/2,0.5-0.032168/2};
- Point(385) = {5+0.767913/2,-tSpn/2,0.5-0.029279/2};
- Point(386) = {5+0.793893/2,-tSpn/2,0.5-0.026405/2};
- Point(387) = {5+0.818712/2,-tSpn/2,0.5-0.023569/2};
- Point(388) = {5+0.842274/2,-tSpn/2,0.5-0.020795/2};
- Point(389) = {5+0.864484/2,-tSpn/2,0.5-0.018106/2};
- Point(390) = {5+0.885257/2,-tSpn/2,0.5-0.015523/2};
- Point(391) = {5+0.904508/2,-tSpn/2,0.5-0.013071/2};
- Point(392) = {5+0.922164/2,-tSpn/2,0.5-0.010770/2};
- Point(393) = {5+0.938153/2,-tSpn/2,0.5-0.008643/2};
- Point(394) = {5+0.952414/2,-tSpn/2,0.5-0.006710/2};
- Point(395) = {5+0.964888/2,-tSpn/2,0.5-0.004990/2};
- Point(396) = {5+0.975528/2,-tSpn/2,0.5-0.003501/2};
- Point(397) = {5+0.984292/2,-tSpn/2,0.5-0.002260/2};
- Point(398) = {5+0.991144/2,-tSpn/2,0.5-0.001280/2};
- Point(399) = {5+0.996057/2,-tSpn/2,0.5-0.000572/2};
- Point(400) = {5+0.999013/2,-tSpn/2,0.5-0.000143/2};
+ Point(301) = {5+1.000000/2,_tSpn,0.5+0.000000/2};
+ Point(302) = {5+0.999013/2,_tSpn,0.5+0.000143/2};
+ Point(303) = {5+0.996057/2,_tSpn,0.5+0.000572/2};
+ Point(304) = {5+0.991144/2,_tSpn,0.5+0.001280/2};
+ Point(305) = {5+0.984292/2,_tSpn,0.5+0.002260/2};
+ Point(306) = {5+0.975528/2,_tSpn,0.5+0.003501/2};
+ Point(307) = {5+0.964888/2,_tSpn,0.5+0.004990/2};
+ Point(308) = {5+0.952414/2,_tSpn,0.5+0.006710/2};
+ Point(309) = {5+0.938153/2,_tSpn,0.5+0.008643/2};
+ Point(310) = {5+0.922164/2,_tSpn,0.5+0.010770/2};
+ Point(311) = {5+0.904508/2,_tSpn,0.5+0.013071/2};
+ Point(312) = {5+0.885257/2,_tSpn,0.5+0.015523/2};
+ Point(313) = {5+0.864484/2,_tSpn,0.5+0.018106/2};
+ Point(314) = {5+0.842274/2,_tSpn,0.5+0.020795/2};
+ Point(315) = {5+0.818712/2,_tSpn,0.5+0.023569/2};
+ Point(316) = {5+0.793893/2,_tSpn,0.5+0.026405/2};
+ Point(317) = {5+0.767913/2,_tSpn,0.5+0.029279/2};
+ Point(318) = {5+0.740877/2,_tSpn,0.5+0.032168/2};
+ Point(319) = {5+0.712890/2,_tSpn,0.5+0.035048/2};
+ Point(320) = {5+0.684062/2,_tSpn,0.5+0.037896/2};
+ Point(321) = {5+0.654508/2,_tSpn,0.5+0.040686/2};
+ Point(322) = {5+0.624345/2,_tSpn,0.5+0.043394/2};
+ Point(323) = {5+0.593691/2,_tSpn,0.5+0.045992/2};
+ Point(324) = {5+0.562667/2,_tSpn,0.5+0.048455/2};
+ Point(325) = {5+0.531395/2,_tSpn,0.5+0.050754/2};
+ Point(326) = {5+0.500000/2,_tSpn,0.5+0.052862/2};
+ Point(327) = {5+0.468605/2,_tSpn,0.5+0.054749/2};
+ Point(328) = {5+0.437333/2,_tSpn,0.5+0.056390/2};
+ Point(329) = {5+0.406309/2,_tSpn,0.5+0.057755/2};
+ Point(330) = {5+0.375655/2,_tSpn,0.5+0.058819/2};
+ Point(331) = {5+0.345492/2,_tSpn,0.5+0.059557/2};
+ Point(332) = {5+0.315938/2,_tSpn,0.5+0.059947/2};
+ Point(333) = {5+0.287110/2,_tSpn,0.5+0.059971/2};
+ Point(334) = {5+0.259123/2,_tSpn,0.5+0.059614/2};
+ Point(335) = {5+0.232087/2,_tSpn,0.5+0.058863/2};
+ Point(336) = {5+0.206107/2,_tSpn,0.5+0.057712/2};
+ Point(337) = {5+0.181288/2,_tSpn,0.5+0.056159/2};
+ Point(338) = {5+0.157726/2,_tSpn,0.5+0.054206/2};
+ Point(339) = {5+0.135516/2,_tSpn,0.5+0.051862/2};
+ Point(340) = {5+0.114743/2,_tSpn,0.5+0.049138/2};
+ Point(341) = {5+0.095492/2,_tSpn,0.5+0.046049/2};
+ Point(342) = {5+0.077836/2,_tSpn,0.5+0.042615/2};
+ Point(343) = {5+0.061847/2,_tSpn,0.5+0.038859/2};
+ Point(344) = {5+0.047586/2,_tSpn,0.5+0.034803/2};
+ Point(345) = {5+0.035112/2,_tSpn,0.5+0.030473/2};
+ Point(346) = {5+0.024472/2,_tSpn,0.5+0.025893/2};
+ Point(347) = {5+0.015708/2,_tSpn,0.5+0.021088/2};
+ Point(348) = {5+0.008856/2,_tSpn,0.5+0.016078/2};
+ Point(349) = {5+0.003943/2,_tSpn,0.5+0.010884/2};
+ Point(350) = {5+0.000987/2,_tSpn,0.5+0.005521/2};
+ Point(351) = {5+0.000000/2,_tSpn,0.5+0.000000/2};
+ Point(352) = {5+0.000987/2,_tSpn,0.5-0.005521/2};
+ Point(353) = {5+0.003943/2,_tSpn,0.5-0.010884/2};
+ Point(354) = {5+0.008856/2,_tSpn,0.5-0.016078/2};
+ Point(355) = {5+0.015708/2,_tSpn,0.5-0.021088/2};
+ Point(356) = {5+0.024472/2,_tSpn,0.5-0.025893/2};
+ Point(357) = {5+0.035112/2,_tSpn,0.5-0.030473/2};
+ Point(358) = {5+0.047586/2,_tSpn,0.5-0.034803/2};
+ Point(359) = {5+0.061847/2,_tSpn,0.5-0.038859/2};
+ Point(360) = {5+0.077836/2,_tSpn,0.5-0.042615/2};
+ Point(361) = {5+0.095492/2,_tSpn,0.5-0.046049/2};
+ Point(362) = {5+0.114743/2,_tSpn,0.5-0.049138/2};
+ Point(363) = {5+0.135516/2,_tSpn,0.5-0.051862/2};
+ Point(364) = {5+0.157726/2,_tSpn,0.5-0.054206/2};
+ Point(365) = {5+0.181288/2,_tSpn,0.5-0.056159/2};
+ Point(366) = {5+0.206107/2,_tSpn,0.5-0.057712/2};
+ Point(367) = {5+0.232087/2,_tSpn,0.5-0.058863/2};
+ Point(368) = {5+0.259123/2,_tSpn,0.5-0.059614/2};
+ Point(369) = {5+0.287110/2,_tSpn,0.5-0.059971/2};
+ Point(370) = {5+0.315938/2,_tSpn,0.5-0.059947/2};
+ Point(371) = {5+0.345492/2,_tSpn,0.5-0.059557/2};
+ Point(372) = {5+0.375655/2,_tSpn,0.5-0.058819/2};
+ Point(373) = {5+0.406309/2,_tSpn,0.5-0.057755/2};
+ Point(374) = {5+0.437333/2,_tSpn,0.5-0.056390/2};
+ Point(375) = {5+0.468605/2,_tSpn,0.5-0.054749/2};
+ Point(376) = {5+0.500000/2,_tSpn,0.5-0.052862/2};
+ Point(377) = {5+0.531395/2,_tSpn,0.5-0.050754/2};
+ Point(378) = {5+0.562667/2,_tSpn,0.5-0.048455/2};
+ Point(379) = {5+0.593691/2,_tSpn,0.5-0.045992/2};
+ Point(380) = {5+0.624345/2,_tSpn,0.5-0.043394/2};
+ Point(381) = {5+0.654508/2,_tSpn,0.5-0.040686/2};
+ Point(382) = {5+0.684062/2,_tSpn,0.5-0.037896/2};
+ Point(383) = {5+0.712890/2,_tSpn,0.5-0.035048/2};
+ Point(384) = {5+0.740877/2,_tSpn,0.5-0.032168/2};
+ Point(385) = {5+0.767913/2,_tSpn,0.5-0.029279/2};
+ Point(386) = {5+0.793893/2,_tSpn,0.5-0.026405/2};
+ Point(387) = {5+0.818712/2,_tSpn,0.5-0.023569/2};
+ Point(388) = {5+0.842274/2,_tSpn,0.5-0.020795/2};
+ Point(389) = {5+0.864484/2,_tSpn,0.5-0.018106/2};
+ Point(390) = {5+0.885257/2,_tSpn,0.5-0.015523/2};
+ Point(391) = {5+0.904508/2,_tSpn,0.5-0.013071/2};
+ Point(392) = {5+0.922164/2,_tSpn,0.5-0.010770/2};
+ Point(393) = {5+0.938153/2,_tSpn,0.5-0.008643/2};
+ Point(394) = {5+0.952414/2,_tSpn,0.5-0.006710/2};
+ Point(395) = {5+0.964888/2,_tSpn,0.5-0.004990/2};
+ Point(396) = {5+0.975528/2,_tSpn,0.5-0.003501/2};
+ Point(397) = {5+0.984292/2,_tSpn,0.5-0.002260/2};
+ Point(398) = {5+0.991144/2,_tSpn,0.5-0.001280/2};
+ Point(399) = {5+0.996057/2,_tSpn,0.5-0.000572/2};
+ Point(400) = {5+0.999013/2,_tSpn,0.5-0.000143/2};
Point(402) = {5+1.000000/2,tSpn/2,0.5+0.000000/2};
Point(403) = {5+0.999013/2,tSpn/2,0.5+0.000143/2};
@@ -427,10 +444,10 @@ EndIf
// Field
If(fWdt != 0)
- Point(1000001) = {-2, -fWdt/2, -2};
- Point(1000002) = {8, -fWdt/2, -2};
- Point(1000003) = {8, -fWdt/2, 2};
- Point(1000004) = {-2, -fWdt/2, 2};
+ Point(1000001) = {-2, _fWdt, -2};
+ Point(1000002) = {8, _fWdt, -2};
+ Point(1000003) = {8, _fWdt, 2};
+ Point(1000004) = {-2, _fWdt, 2};
Point(1000005) = {-2, fWdt/2, -2};
Point(1000006) = {8, fWdt/2, -2};
Point(1000007) = {8, fWdt/2, 2};
diff --git a/fpm/models/n0012.py b/fpm/models/n0012.py
index aecc964..bb51b25 100644
--- a/fpm/models/n0012.py
+++ b/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()
diff --git a/fpm/models/oneraM6.geo b/fpm/models/oneraM6.geo
new file mode 100644
index 0000000..418d940
--- /dev/null
+++ b/fpm/models/oneraM6.geo
@@ -0,0 +1,327 @@
+/* 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};
diff --git a/fpm/src/fBody.cpp b/fpm/src/fBody.cpp
index 571fb58..2b9bf3e 100644
--- a/fpm/src/fBody.cpp
+++ b/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());
diff --git a/fpm/src/fBuilder.cpp b/fpm/src/fBuilder.cpp
index 720b3b8..fe5a573 100644
--- a/fpm/src/fBuilder.cpp
+++ b/fpm/src/fBuilder.cpp
@@ -22,6 +22,7 @@
#include "wTag.h"
#include "wElement.h"
+#include "wNode.h"
#include "wMem.h"
#define PI 3.14159
@@ -31,13 +32,12 @@ using namespace fpm;
/**
* @brief Initialize the solver and perform sanity checks
- * @authors Adrien Crovato
+ * @authors Adrien Crovato & Axel Dechamps
*/
Builder::Builder(std::shared_ptr<Problem> _pbl) : pbl(_pbl), valid(false)
{
- // Check the problem and update element memory
+ // Check the problem
pbl->check();
- pbl->updateMem();
// Setup AIC matrices
pbl->field ? nF = pbl->field->tag->elems.size() : nF = 0;
@@ -50,15 +50,8 @@ Builder::Builder(std::shared_ptr<Problem> _pbl) : pbl(_pbl), valid(false)
Cx = Eigen::MatrixXd::Zero(nP, nF);
Cy = Eigen::MatrixXd::Zero(nP, nF);
Cz = Eigen::MatrixXd::Zero(nP, nF);
-}
-/**
- * @brief Fill AIC matrices
- */
-void Builder::run()
-{
// global->local id
- std::map<Element *, int> rows;
int i = 0;
for (auto body : pbl->bodies)
for (auto e : body->tag->elems)
@@ -66,28 +59,175 @@ void Builder::run()
rows[e] = i;
i++;
}
+}
- // body to body
- for (auto ibody : pbl->bodies)
- for (auto ei : ibody->tag->elems)
+/**
+ * @brief Fill AIC matrices
+ */
+void Builder::run()
+{
+ // Update element memory
+ pbl->updateMem();
+
+ // body to body - body panels
+ for (auto jbody : pbl->bodies)
+ {
+ for (auto ej : jbody->tag->elems)
{
+ // Initialisation of storage variables for the iteration
+ Eigen::Matrix3d glob2loc;
+ size_t nVertices = ej->nodes.size();
+ Eigen::MatrixXd trsfCoordSrc(3, nVertices);
+
+ // Computations related to source panels
+ glob2loc = Glob2loc(false, ej); // Transformation matrix
+ trsfCoordSrc = TrsfCoordSrc(false, ej, glob2loc); // Local coordinates of the source panel
+
// body panels
- for (auto jbody : pbl->bodies)
- for (auto ej : jbody->tag->elems)
+ for (auto ibody : pbl->bodies)
+ {
+ for (auto ei : ibody->tag->elems)
+ {
+ // Initialisation of storage variables for the iteration
+ Eigen::Vector3d trsfCoordTgt;
+ Eigen::Vector3d distance;
+ Eigen::Vector2d mutau;
+
+ // Computations related to target panels
+ trsfCoordTgt = TrsfCoordTgt(false, ei, ej, glob2loc); // Local coordinates of the target panel
+ distance = Distance(false, ei, ej, glob2loc); // Distances between source and target panels
+ mutau = MuTau(false, ei, ej, trsfCoordSrc, trsfCoordTgt, distance); // Surface AIC
+
+ // Assembly of the matrices
+ A(rows.at(ei), rows.at(ej)) = mutau(0);
+ B(rows.at(ei), rows.at(ej)) = mutau(1);
+ }
+ }
+ }
+ }
+
+ // body to body - body panels - symmetry
+ if (pbl->symY)
+ {
+ for (auto jbody : pbl->bodies)
+ {
+ for (auto ej : jbody->tag->elems)
+ {
+ // Initialisation of storage variables for the iteration
+ Eigen::Matrix3d glob2loc;
+ size_t nVertices = ej->nodes.size();
+ Eigen::MatrixXd trsfCoordSrc(3, nVertices);
+
+ // Computations related to source panels
+ glob2loc = Glob2loc(true, ej); // Transformation matrix
+ trsfCoordSrc = TrsfCoordSrc(true, ej, glob2loc); // Local coordinates of the source panel
+
+ // body panels
+ for (auto ibody : pbl->bodies)
+ {
+ for (auto ei : ibody->tag->elems)
+ {
+ // Initialisation of storage variables for the iteration
+ Eigen::Vector3d trsfCoordTgt;
+ Eigen::Vector3d distance;
+ Eigen::Vector2d mutau;
+
+ // Computations related to target panels
+ trsfCoordTgt = TrsfCoordTgt(true, ei, ej, glob2loc); // Local coordinates of the target panel
+ distance = Distance(true, ei, ej, glob2loc); // Distances between source and target panels
+ mutau = MuTau(true, ei, ej, trsfCoordSrc, trsfCoordTgt, distance); // Surface AIC
+
+ // Assembly of the matrices
+ A(rows.at(ei), rows.at(ej)) -= mutau(0);
+ B(rows.at(ei), rows.at(ej)) -= mutau(1);
+ }
+ }
+ }
+ }
+ }
+
+ // body to body - wake panels
+ for (auto jbody : pbl->bodies)
+ {
+ for (auto wake : jbody->wakes)
+ {
+ for (auto ew : wake->tag->elems)
+ {
+ // Initialisation of storage variables for the iteration
+ Eigen::Matrix3d glob2loc;
+ size_t nVertices = ew->nodes.size();
+ Eigen::MatrixXd trsfCoordSrc(3, nVertices);
+
+ // Computations related to source panels
+ glob2loc = Glob2loc(false, ew); // Transformation matrix
+ trsfCoordSrc = TrsfCoordSrc(false, ew, glob2loc); // Local coordinates of the source panel
+
+ // wake panels
+ for (auto ibody : pbl->bodies)
{
- A(rows.at(ei), rows.at(ej)) = mu(ei, ej);
- B(rows.at(ei), rows.at(ej)) = tau(ei, ej);
+ for (auto ei : ibody->tag->elems)
+ {
+ // Initialisation of storage variables for the iteration
+ Eigen::Vector3d trsfCoordTgt;
+ Eigen::Vector3d distance;
+ Eigen::Vector2d wmutau;
+
+ // Computations related to target panels
+ trsfCoordTgt = TrsfCoordTgt(false, ei, ew, glob2loc); // Local coordinates of the target panel
+ distance = Distance(false, ei, ew, glob2loc); // Distances between source and target panels
+ wmutau = MuTau(false, ei, ew, trsfCoordSrc, trsfCoordTgt, distance); // Surface AIC
+
+ // Assembly of the matrices
+ A(rows.at(ei), rows.at(wake->wkMap.at(ew).first)) += wmutau(0);
+ A(rows.at(ei), rows.at(wake->wkMap.at(ew).second)) -= wmutau(0);
+ }
}
- // wake panels
- for (auto jbody : pbl->bodies)
- for (auto wake : jbody->wakes)
- for (auto ew : wake->tag->elems)
+ }
+ }
+ }
+
+ // body to body - wake panels - symmetry
+ if (pbl->symY)
+ {
+ for (auto jbody : pbl->bodies)
+ {
+ for (auto wake : jbody->wakes)
+ {
+ for (auto ew : wake->tag->elems)
+ {
+ // Initialisation of storage variables for the iteration
+ Eigen::Matrix3d glob2loc;
+ size_t nVertices = ew->nodes.size();
+ Eigen::MatrixXd trsfCoordSrc(3, nVertices);
+
+ // Computations related to source panels
+ glob2loc = Glob2loc(true, ew); // Transformation matrix
+ trsfCoordSrc = TrsfCoordSrc(true, ew, glob2loc); // Local coordinates of the source panel
+
+ // wake panels
+ for (auto ibody : pbl->bodies)
{
- double wmu = mu(ei, ew);
- A(rows.at(ei), rows.at(wake->wkMap.at(ew).first)) += wmu;
- A(rows.at(ei), rows.at(wake->wkMap.at(ew).second)) -= wmu;
+ for (auto ei : ibody->tag->elems)
+ {
+ // Initialisation of storage variables for the iteration
+ Eigen::Vector3d trsfCoordTgt;
+ Eigen::Vector3d distance;
+ Eigen::Vector2d wmutau;
+
+ // Computations related to target panels
+ trsfCoordTgt = TrsfCoordTgt(true, ei, ew, glob2loc); // Local coordinates of the target panel
+ distance = Distance(true, ei, ew, glob2loc); // Distances between source and target panels
+ wmutau = MuTau(true, ei, ew, trsfCoordSrc, trsfCoordTgt, distance); // Surface AIC
+
+ // Assembly of the matrices
+ A(rows.at(ei), rows.at(wake->wkMap.at(ew).first)) -= wmutau(0);
+ A(rows.at(ei), rows.at(wake->wkMap.at(ew).second)) += wmutau(0);
+ }
}
+ }
+ }
}
+ }
// field to body
if (pbl->field)
@@ -98,58 +238,316 @@ void Builder::run()
C(rows.at(e), j) = sigma(e, pbl->field->tag->elems[j]);
Cx(rows.at(e), j) = sigmaX(e, pbl->field->tag->elems[j]);
Cy(rows.at(e), j) = sigmaY(e, pbl->field->tag->elems[j]);
- Cz(rows.at(e), j) = sigmaY(e, pbl->field->tag->elems[j]);
+ Cz(rows.at(e), j) = sigmaZ(e, pbl->field->tag->elems[j]);
}
valid = true;
}
/**
- * @brief Compute doublet AIC from body panel ej to body panel ei
+ * @brief Compute unit axes of the source panel
*/
-double Builder::mu(Element *ei, Element *ej)
+Eigen::Matrix3d Builder::Glob2loc(bool symy, Element *ej)
{
- if (ei == ej)
- return 0.5;
- else
- return 0;
+ Eigen::Matrix3d glob2loc; // Transformation matrix
+ Eigen::Vector3d n = ej->normal(); // Unit normal vector of the source panel
+
+ // Computation of the unit longitudinal vector - Source panel
+ Eigen::Vector3d cg = ej->getVMem().getCg();
+ Eigen::Vector3d x0 = ej->nodes[0]->pos; // Coordinates of node 1
+ Eigen::Vector3d l = (x0 - cg).normalized();
+
+ // Computation of the unit perpendicular vector - Source panel
+ Eigen::Vector3d p = (n.cross(l)).normalized();
+
+ // Handling of the outputs
+ glob2loc.row(0) = l;
+ glob2loc.row(1) = p;
+ glob2loc.row(2) = n;
+
+ if(symy)
+ {
+ glob2loc.row(1) = -glob2loc.row(1);
+ glob2loc.col(1) = -glob2loc.col(1);
+ }
+
+ return glob2loc;
+}
+
+/**
+ * @brief Compute coordinates of the source panel from global to local axes
+ */
+Eigen::MatrixXd Builder::TrsfCoordSrc(bool symy, Element *ej, Eigen::Matrix3d const &glob2loc)
+{
+ size_t nVertices = ej->nodes.size();
+ Eigen::MatrixXd trsfCoordSrc(3, nVertices);
+ Eigen::MatrixXd coords(3, nVertices);
+
+ for (int i = 0; i < 3; ++i)
+ for (int j = 0; j < nVertices; ++j)
+ coords(i, j) = ej->nodes[(nVertices - 1) - j]->pos[i]; // Global coordinates of the nodes of the source panel
+
+ if(symy)
+ coords.row(1) = -coords.row(1);
+
+ for (int i = 0; i < nVertices; ++i)
+ trsfCoordSrc.col(i) = glob2loc * coords.col(i); // Local coordinates of the nodes of the source panel
+
+ return trsfCoordSrc;
}
+
/**
- * @brief Compute source AIC from body panel ej to body panel ei
+ * @brief Compute coordinates of the target panel from global to local axes
*/
-double Builder::tau(Element *ei, Element *ej)
+Eigen::Vector3d Builder::TrsfCoordTgt(bool symy, Element *ei, Element *ej, Eigen::Matrix3d const &glob2loc)
{
- if (ei == ej)
- return 0;
+ Eigen::Vector3d trsfCoordTgt;
+ Eigen::Vector3d cg = ei->getVMem().getCg(); // Global coordinates of the CG of the target panel
+ trsfCoordTgt = glob2loc * cg; // Local coordinates of the target panel
+
+ return trsfCoordTgt;
+}
+
+/**
+ * @brief Compute distances between two panels
+ */
+Eigen::Vector3d Builder::Distance(bool symy, Element *ei, Element *ej, Eigen::Matrix3d const &glob2loc)
+{
+ Eigen::Vector3d distance; // Distance between the target and source panels
+ Eigen::Vector3d cgSrc = ej->getVMem().getCg();
+ Eigen::Vector3d cgTgt = ei->getVMem().getCg();
+
+ for (int i = 0; i < 3; ++i)
+ distance(i) = cgTgt(i) - cgSrc(i);
+
+ if(symy)
+ distance(1) = cgTgt(1) + 2 * cgSrc(1);
+
+ distance = glob2loc * distance;
+
+ return distance;
+}
+
+/**
+ * @brief Compute doublet AIC and source AIC from body panel ej to body panel ei
+ */
+Eigen::Vector2d Builder::MuTau(bool symy, Element *ei, Element *ej, Eigen::MatrixXd const &trsfCoordSrc, Eigen::Vector3d const &trsfCoordTgt, Eigen::Vector3d const &distance)
+{
+ double mu = 0;
+ double tau = 0;
+ Eigen::Vector2d mu_tau;
+ Eigen::Vector4d d, e, h, r, F, G;
+ size_t nVertices = ei->nodes.size();
+
+ // Computation of the doublet AIC - mu
+ if (ei->no == ej->no && !symy)
+ mu = 0.5; // Doublet term - Influence of a panel on itself
else
- return 0;
+ {
+ // Coefficients
+ for (int i = 1; i <= nVertices; ++i)
+ {
+ e(i - 1) = (trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) + distance(2) * distance(2);
+ h(i - 1) = (trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1));
+ r(i - 1) = sqrt((trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) + (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) * (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) + distance(2) * distance(2));
+ }
+ for (int i = 1; i <= nVertices; ++i)
+ {
+ int j = i % nVertices + 1;
+
+ F(i - 1) = (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)) * e(i - 1) - (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * h(i - 1);
+ G(i - 1) = (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)) * e(j - 1) - (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * h(j - 1);
+ }
+ // Doublet term - Influence of a panel on another panel
+ for (int i = 1; i <= nVertices; ++i)
+ {
+ int j = i % nVertices + 1;
+
+ mu = mu + atan2(distance(2) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * (F(i - 1) * r(j - 1) - G(i - 1) * r(i - 1)),
+ distance(2) * distance(2) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * r(i - 1) * r(j - 1) + F(i - 1) * G(i - 1));
+ }
+ mu = -1 / (4 * PI) * mu;
+ }
+
+ // Computation of the source AIC - tau
+ // Coefficients
+ for (int i = 1; i <= nVertices; ++i)
+ {
+ int j = i % nVertices + 1;
+
+ d(i - 1) = sqrt((trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1)) + (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)) * (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)));
+ r(i - 1) = sqrt((trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) + (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) * (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) + distance(2) * distance(2));
+ }
+
+ if (ei->no == ej->no && !symy)
+ {
+ // Source term - Influence of a panel on itself
+ for (int i = 1; i <= nVertices; ++i)
+ {
+ int j = i % nVertices + 1;
+
+ tau = tau + ((trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)) - (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1))) / d(i - 1) * log((r(i - 1) + r(j - 1) + d(i - 1)) / (r(i - 1) + r(j - 1) - d(i - 1)));
+ }
+ tau = -1 / (4 * PI) * tau;
+ }
+ else
+ {
+ // Source term - Influence of a panel on another panel
+ for (int i = 1; i <= nVertices; ++i)
+ {
+ int j = i % nVertices + 1;
+
+ tau = tau + ((trsfCoordTgt(0) - trsfCoordSrc(0,i - 1)) * (trsfCoordSrc(1,j - 1) - trsfCoordSrc(1,i - 1)) - (trsfCoordTgt(1) - trsfCoordSrc(1,i - 1)) * (trsfCoordSrc(0,j - 1) - trsfCoordSrc(0,i - 1))) / d(i - 1) * log((r(i - 1) + r(j - 1) + d(i - 1)) / (r(i - 1) + r(j - 1) - d(i - 1)));
+ }
+
+ tau = -1 / (4 * PI) * tau - distance(2) * mu;
+ }
+
+ mu_tau(0) = mu;
+ mu_tau(1) = tau;
+
+ return mu_tau;
}
+
/**
* @brief Compute source AIC from field panel ej to body panel ei
*/
double Builder::sigma(Element *ei, Element *ej)
{
+ /* Initialisation of the parameters */
+ // double sigma = 0;
+ // double A, B, C, I, II, III, R;
+
+ /* Import of the coordinates of the influencing field cell vertices */
+ // x = 0; y = 0; z = 0;
+ // xC = 0; yC = 0; zC = 0;
+
+ /* Field source potential influence between two field panels */
+ // for (int i = 1; i <= 2; ++i) {
+ // for (int j = 1; j <= 2; ++j) {
+ // for (int k = 1; k <= 2; ++k) {
+ // A = x - xC[i-1];
+ // B = y - yC[j-1];
+ // C = z - zC[k-1];
+ // R = sqrt(A * A + B * B + C * C);
+ // I = B * C * log(A + R) - A * A / 2 * atan((B * C) / (A * R));
+ // II = C * A * log(B + R) - B * B / 2 * atan((C * A) / (B * R));
+ // III = A * B * log(C + R) - C * C / 2 * atan((A * B) / (C * R));
+
+ // sigma -= pow((-1),(double)(i + j + k)) * (I + II + III);
+ // }
+ // }
+ // }
+
+ // sigma = - 1 / (4 * PI) * sigma;
+
+ // return sigma;
+
return 0;
}
+
/**
* @brief Compute source x-velocity AIC from field panel ej to body panel ei
*/
double Builder::sigmaX(Element *ei, Element *ej)
{
+ /* Initialisation of the parameters */
+ // double sigmaX = 0;
+ // double A, B, C, R, t0;
+
+ /* Import of the coordinates of the influencing field cell vertices */
+ // x = 0; y = 0; z = 0;
+ // xC = 0; yC = 0; zC = 0;
+
+ /* Field source velocity influence between a field panel and a body panel */
+ // for (int i = 1; i <= 2; ++i) {
+ // for (int j = 1; j <= 2; ++j) {
+ // for (int k = 1; k <= 2; ++k) {
+ // A = (x - xC[i-1]);
+ // B = (y - yC[j-1]);
+ // C = (z - zC[k-1]);
+ // R = sqrt(A * A + B * B + C * C);
+ // t0 = B / 2 * log((R + C) / (R - C)) + C / 2 * log((R + B) / (R - B)) - A * atan((B * C) / (A * R));
+
+ // sigmaX += pow((-1),(double)(i + j + k)) * t0;
+ // }
+ // }
+ // }
+
+ // sigmaX = 1 / (4 * PI) * sigmaX;
+
+ // return sigmaX;
+
return 0;
}
+
/**
* @brief Compute source y-velocity AIC from field panel ej to body panel ei
*/
double Builder::sigmaY(Element *ei, Element *ej)
{
+ /* Initialisation of the parameters */
+ // double sigmaY = 0;
+ // double A, B, C, R, t1;
+
+ /* Import of the coordinates of the influencing field cell vertices */
+ // x = 0; y = 0; z = 0;
+ // xC = 0; yC = 0; zC = 0;
+
+ /* Field source velocity influence between a field panel and a body panel */
+ // for (int i = 1; i <= 2; ++i) {
+ // for (int j = 1; j <= 2; ++j) {
+ // for (int k = 1; k <= 2; ++k) {
+ // A = (x - xC[i-1]);
+ // B = (y - yC[j-1]);
+ // C = (z - zC[k-1]);
+ // R = sqrt(A * A + B * B + C * C);
+ // t1 = C / 2 * log((R + A) / (R - A)) + A / 2 * log((R + C) / (R - C)) - B * atan((C * A) / (B * R));
+
+ // sigmaY += pow((-1),(double)(i + j + k)) * t1;
+ // }
+ // }
+ // }
+
+ // sigmaY = 1 / (4 * PI) * sigmaY;
+
+ // return sigmaY;
+
return 0;
}
+
/**
* @brief Compute source z-velocity AIC from field panel ej to body panel ei
*/
double Builder::sigmaZ(Element *ei, Element *ej)
{
+ /* Initialisation of the parameters */
+ // double sigmaZ = 0;
+ // double A, B, C, R, t2;
+
+ /* Import of the coordinates of the influencing field cell vertices */
+ // x = 0; y = 0; z = 0;
+ // xC = 0; yC = 0; zC = 0;
+
+ /* Field source velocity influence between a field panel and a body panel */
+ // for (int i = 1; i <= 2; ++i) {
+ // for (int j = 1; j <= 2; ++j) {
+ // for (int k = 1; k <= 2; ++k) {
+ // A = (x - xC[i-1]);
+ // B = (y - yC[j-1]);
+ // C = (z - zC[k-1]);
+ // R = sqrt(A * A + B * B + C * C);
+ // t2 = A / 2 * log((R + B) / (R - B)) + B / 2 * log((R + A) / (R - A)) - C * atan((A * B) / (C * R));
+
+ // sigmaZ += pow((-1),(double)(i + j + k)) * t2;
+ // }
+ // }
+ // }
+
+ // sigmaZ = 1 / (4 * PI) * sigmaZ;
+
+ // return sigmaZ;
+
return 0;
}
diff --git a/fpm/src/fBuilder.h b/fpm/src/fBuilder.h
index bb290cc..0e721a8 100644
--- a/fpm/src/fBuilder.h
+++ b/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);
diff --git a/fpm/src/fProblem.cpp b/fpm/src/fProblem.cpp
index ef7843f..5d1b840 100644
--- a/fpm/src/fProblem.cpp
+++ b/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);
+ }
}
/**
diff --git a/fpm/src/fProblem.h b/fpm/src/fProblem.h
index 4bf8c33..3cdd15c 100644
--- a/fpm/src/fProblem.h
+++ b/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;
diff --git a/fpm/src/fSolver.cpp b/fpm/src/fSolver.cpp
index 51066d4..6e49acc 100644
--- a/fpm/src/fSolver.cpp
+++ b/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"] = φ
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"] = ρ
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
diff --git a/fpm/src/fSolver.h b/fpm/src/fSolver.h
index a7c9084..c2cdb1a 100644
--- a/fpm/src/fSolver.h
+++ b/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
diff --git a/fpm/src/fWake.cpp b/fpm/src/fWake.cpp
index 2b95ce7..0b40761 100644
--- a/fpm/src/fWake.cpp
+++ b/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);
diff --git a/fpm/tests/M6.py b/fpm/tests/M6.py
new file mode 100644
index 0000000..702110d
--- /dev/null
+++ b/fpm/tests/M6.py
@@ -0,0 +1,83 @@
+#!/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()
diff --git a/fpm/tests/a445.py b/fpm/tests/a445.py
new file mode 100644
index 0000000..8cd41f6
--- /dev/null
+++ b/fpm/tests/a445.py
@@ -0,0 +1,84 @@
+#!/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()
+
diff --git a/fpm/tests/n12.py b/fpm/tests/n12.py
index 663229f..59010c7 100644
--- a/fpm/tests/n12.py
+++ b/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
diff --git a/fpm/tests/n12f.py b/fpm/tests/n12f.py
index edd633a..086b301 100644
--- a/fpm/tests/n12f.py
+++ b/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
diff --git a/fpm/tests/n12t.py b/fpm/tests/n12t.py
index 7834735..fa1d843 100644
--- a/fpm/tests/n12t.py
+++ b/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
--
GitLab