diff --git a/PotentialFlow/generateOptFile.py b/PotentialFlow/generateOptFile.py
new file mode 100644
index 0000000000000000000000000000000000000000..7f121d051cbe7be23549e32571fc38bc8a57d24a
--- /dev/null
+++ b/PotentialFlow/generateOptFile.py
@@ -0,0 +1,26 @@
+#!/usr/bin/env python
+# This Python script generates a formay file for the "phi.pos"
+# map generated by the model.
+
+#print input[0]
+phiTrail = input[0]
+
+fout = open('KJ.pos.opt', 'w')
+fout.write('l=PostProcessing.NbViews-1;\n')
+fout.write('View[l].IntervalsType = 3;\n')
+fout.write('Mesh.SurfaceEdges = 0;\n')
+fout.write('View[l].RangeType = 2;\n')
+fout.write('View[l].CustomMin = ' + str(phiTrail) + ';\n')
+fout.write('View[l].CustomMax = ' + str(1.001*phiTrail) + ';\n')
+fout.write('View[l].NbIso = 10;\n\n')
+
+fout.close()
+
+
+
+
+
+
+
+
+
diff --git a/PotentialFlow/magnus.geo b/PotentialFlow/magnus.geo
new file mode 100644
index 0000000000000000000000000000000000000000..c78c6ac81c8d10598121a1cbcaccbb01cffb466e
--- /dev/null
+++ b/PotentialFlow/magnus.geo
@@ -0,0 +1,71 @@
+Include "magnus_common.pro";
+
+A = (BoxSize-1)/2.;
+B = A;
+R = 0.5;
+lc = A/10;
+
+
+If( Flag_Object == 0 ) // Cylinder
+
+lcr = 0.1;
+Point(1) = { 2*R, 0.0, 0.0, lcr};
+Point(2) = { R, -R, 0.0, lcr};
+Point(3) = { 0.0, 0.0, 0.0, lcr};
+Point(4) = { R, R, 0.0 , lcr};
+Point(5) = { R, 0.0, 0.0 , lcr};
+Circle(1) = {1,5,2};
+Circle(2) = {2,5,3};
+Circle(3) = {3,5,4};
+Circle(4) = {4,5,1};
+
+// Points to be connected with the outer boundary
+PtA = 4;
+PtB = 2;
+
+Else // naca airfoil
+
+lca = 0.03;
+Include "nacaAirFoil.geo";
+PtA = 121;
+PtB = 81;
+
+EndIf
+
+Point(306) = { 0, B, 0, lc};
+Point(307) = { 0,-B, 0, lc};
+
+Line(5) = { PtA, 306 };
+Line(6) = { PtB, 307 };
+
+Point(308) = {-A,-B, 0, lc};
+Point(309) = {-A, B, 0, lc};
+Point(310) = { A+1, B, 0, lc};
+Point(311) = { A+1,-B, 0, lc};
+
+Line( 7) = { 306, 309 };
+Line( 8) = { 309, 308 };
+Line( 9) = { 308, 307 };
+Line(10) = { 307, 311 };
+Line(11) = { 311, 310 };
+Line(12) = { 310, 306 };
+
+Line Loop(21) = { 7, 8, 9, -6, 2, 3, 5 }; // aire a gauche
+Line Loop(22) = { 10, 11, 12, -5, 4, 1, 6 };
+
+Plane Surface(24) = { 21 };
+Plane Surface(25) = { 22 };
+
+Physical Surface("Fluid", 2) = { 24, 25 };
+Physical Line("UpStream", 10) = { 8 };
+Physical Line("DownStream", 11) = { 11 };
+Physical Line("Airfoil", 12) = { 1 ... 4 };
+Physical Line("Wake", 13) = { 5 };
+
+
+
+//Line Loop(20) = {1,2,3,4};
+//Plane Surface(23) = {20};
+//Physical Surface("Cylinder", 1) = {23};
+//Physical Line("Outer", 10) = { 7 ... 12 };
+
diff --git a/PotentialFlow/magnus.pro b/PotentialFlow/magnus.pro
new file mode 100644
index 0000000000000000000000000000000000000000..ab45cee046c9fedff1f4b85122b3638876a14209
--- /dev/null
+++ b/PotentialFlow/magnus.pro
@@ -0,0 +1,410 @@
+/* -------------------------------------------------------------------
+ Tutorial 6 : Potential flow, lift and Magnus effect
+
+ Features:
+ - Discontinuous scalar field
+ - Multivalued scalar field
+ - Potential flow
+ - Non-linear iteration on an Associated global quantity
+ - Using a run-time variable in a Post-Operation (via Python)
+
+ To compute the solution in a terminal:
+ getdp magnus.pro -solve LPE -pos LPE
+ gmsh magnus.geo velocity.pos
+
+ To compute the solution interactively from the Gmsh GUI:
+ File > Open > magnus.pro
+ Run (button at the bottom of the left panel)
+ ------------------------------------------------------------------- */
+/*
+ This model solves a 2D potential flow around a cylinder or a naca airfoil
+ placed in a uniform "V_infinity" flow.
+ Potential flows are defined by
+
+ V = grad phi => curl V = 0
+
+ where the multivalued scalar potential "phi"
+ presents a discontinuity of magnitude "deltaPhi"
+ across a cut "Sur_Cut".
+ In consequence, the velocity field "V" is curl-free
+ over the wole domain of analysis "Vol_rho"
+ but its circulation over any closed curve circling around the object,
+ a quantity often noted "Gamma" in the literature,
+ gives "deltaPhi" as a result.
+ The circulation of "V" over closed curves
+ not circling around the object is zero.
+
+ In the getDP model, the discontinuity "deltaPhi" [m^2/s]
+ is defined as a global quantity in the Functional space of "phi".
+ The associated (dual) global quantity is the mass flow density,
+ noted "Dmdt", which has [kg/s] as a unit.
+
+ Governing equations are
+
+ div ( rho[] grad phi ) = 0 in Vol_rho
+ rho[] grad phi . n = rho[] Vn = 0 on Sur_Neu
+
+ whereas the uniform velocity field "V_infinity" is imposed
+ by means of Dirichlet boundary conditions
+ on the upstream and downstream surfaces "GammaUp" and "GammaDown".
+
+ Momentum equation is decoupled in case of stationary potential flows.
+ The velocity filed can be solved first, and the pressure
+ is obtained afterwards by invoking Bernoulli's theorem:
+
+ p = -0.5 * rho[] * SquNorm [ {d phi} ]
+
+ The "Lift" force in "Y" direction is then evaluated
+ either by integrating "p * CompY[Normal[]]"
+ over the contour of the object,
+ or by the Kutta-Jukowski theorem
+
+ Lift = - L_z rho[] V_infinity Gamma
+
+ which is a first order approximation of the former.
+
+ "Cylinder" case:
+ Either the circulation "deltaPhi" or the mass flow rate "Dmdt"
+ can be imposed, according to the "Flag_Circ" flag.
+ The Lift is evaluated by both the integration of pressure
+ and the Kutta-Jukowski approximation.
+
+ "Airfoil" case:
+ If "Flag_Circ == 1", the model works similar to the "Cylinder" case.
+ If "Flag_Circ == 0", the mass flow rate is not directly imposed
+ in this case.
+ It is determined so that the Kutta condition is verified.
+ This condition states that the actual value of "Dmdt"
+ is the one for which the stagnation point
+ (singularity of the velocity field)
+ is located at the trailing edge of the airfoil,
+ i.e., when the function
+
+ argVTrail[] = ( Atan2[CompY[$1], CompX[$1] ] - Incidence ) / deg ;
+
+ returns zero when evaluated for the velocity at a point "P_edge"
+ close to the trailing edge.
+ A non linear iteration is thus done by means of a pseudo-newton scheme,
+ updating the value of the imposed "Dmdt" until the norm of the residual
+
+ f(V) = argVTrail[ {d phi } ] OnPoint {1.0001,0,0}
+
+ comes below a fixed tolerance.
+ */
+
+Include "magnus_common.pro";
+
+Group{
+ // Physical regions
+ Fluid = Region[ 2 ];
+ GammaUp = Region[ 10 ];
+ GammaDown = Region[ 11 ];
+ GammaAirf = Region[ 12 ];
+ GammaCut = Region[ 13 ];
+
+ // Abstract regions
+ Vol_rho = Region[ { Fluid } ];
+ Sur_Dir = Region[ { GammaUp, GammaDown } ];
+ Sur_Neu = Region[ { GammaAirf } ];
+ Sur_Cut = Region[ { GammaCut } ];
+}
+
+Function{
+ rho[] = 1.225; // kg/m^3
+ MassFlowRate[] = MassFlowRate; // kg/s
+ DeltaPhi[] = DeltaPhi; // m^2/s
+ argVTrail[] = ( Atan2[CompY[$1], CompX[$1] ] - Incidence ) / deg ;
+}
+
+Jacobian {
+ { Name Vol ;
+ Case {
+ { Region All ; Jacobian Vol ; }
+ }
+ }
+ { Name Sur ;
+ Case {
+ { Region All ; Jacobian Sur ; }
+ }
+ }
+}
+
+Integration {
+ { Name Int ;
+ Case {
+ { Type Gauss ;
+ Case {
+ { GeoElement Point ; NumberOfPoints 1 ; }
+ { GeoElement Line ; NumberOfPoints 4 ; }
+ { GeoElement Triangle ; NumberOfPoints 6 ; }
+ { GeoElement Quadrangle ; NumberOfPoints 7 ; }
+ { GeoElement Tetrahedron ; NumberOfPoints 15 ; }
+ { GeoElement Hexahedron ; NumberOfPoints 34 ; }
+ }
+ }
+ }
+ }
+}
+
+Constraint{
+ { Name Dirichlet; Type Assign;
+ Case{
+ // Boundary conditions for the V_infinity uniform flow
+ { Region GammaDown; Value Velocity*BoxSize; }
+ { Region GammaUp ; Value 0; }
+ }
+ }
+ { Name DeltaPhi; Type Assign;
+ Case{
+ { Region Sur_Cut; Value DeltaPhi[]; }
+ }
+ }
+ { Name MassFlowRate; Type Assign;
+ Case{
+ { Region Sur_Cut; Value MassFlowRate[]; }
+ }
+ }
+}
+
+// Domains of definition used in the description of the function space
+// These groups contain both volume and surface regions/elements.
+Group{
+ Dom_Vh = Region[ { Vol_rho, Sur_Dir, Sur_Neu, Sur_Cut } ];
+ Dom_Cut = ElementsOf[ Dom_Vh, OnPositiveSideOf Sur_Cut ];
+}
+FunctionSpace{
+ { Name Vh; Type Form0;
+ BasisFunction{
+ { Name vn; NameOfCoef phin; Function BF_Node;
+ Support Dom_Vh; Entity NodesOf[All]; }
+ { Name vc; NameOfCoef dphi; Function BF_GroupOfNodes;
+ Support Dom_Cut; Entity GroupsOfNodesOf[ Sur_Cut ];}
+ }
+ SubSpace {
+ { Name phiCont ; NameOfBasisFunction { vn } ; }
+ { Name phiDisc ; NameOfBasisFunction { vc } ; }
+ }
+ GlobalQuantity {
+ { Name Circ ; Type AliasOf ; NameOfCoef dphi ; }
+ { Name Dmdt ; Type AssociatedWith ; NameOfCoef dphi ; }
+ }
+ Constraint{
+ {NameOfCoef phin; EntityType NodesOf;
+ NameOfConstraint Dirichlet;}
+ If( Flag_Circ )
+ {NameOfCoef Circ; EntityType GroupsOfNodesOf;
+ NameOfConstraint DeltaPhi;}
+ Else
+ If( Flag_Object == 0 ) // if Cylinder only
+ {NameOfCoef Dmdt; EntityType GroupsOfNodesOf;
+ NameOfConstraint MassFlowRate;}
+ EndIf
+ EndIf
+ }
+ }
+}
+
+Formulation{
+ {Name LPE; Type FemEquation;
+ Quantity{
+ {Name phi; Type Local; NameOfSpace Vh;}
+ { Name phiCont ; Type Local ; NameOfSpace Vh[phiCont] ; }
+ { Name phiDisc ; Type Local ; NameOfSpace Vh[phiDisc] ; }
+ { Name Circ ; Type Global ; NameOfSpace Vh[Circ] ; }
+ { Name Dmdt ; Type Global ; NameOfSpace Vh[Dmdt] ; }
+ }
+ Equation{
+ Galerkin{[ rho[] * Dof{d phi}, {d phi}];
+ In Vol_rho; Jacobian Vol; Integration Int;}
+ If( !Flag_Circ )
+ GlobalTerm { [ -Dof{Dmdt} , {Circ} ] ; In Sur_Cut ; }
+ If( Flag_Object == 1 )
+ GlobalTerm { [ -Dof{Dmdt} , {Dmdt} ] ; In Sur_Cut ; }
+ GlobalTerm { [ $newDmdt , {Dmdt} ] ; In Sur_Cut ; }
+ EndIf
+ EndIf
+ }
+ }
+}
+
+Resolution{
+ {Name LPE;
+ System{
+ {Name A; NameOfFormulation LPE;}
+ }
+ Operation{
+ InitSolution[A];
+
+ If( Flag_Circ || ( Flag_Object == 0 ) )
+
+ Generate[A]; Solve[A]; SaveSolution[A];
+
+ Else
+ // Pseudo-Newton iteration to determine Dmdt in the Airfoil case.
+ // The run-time variable $newDmdt used in Generate[A]
+ // whereas $circ, $dmdt and argV are evaluated
+ // by PostOperation[Trailing]
+
+ DeleteFile["KJiter.txt"];
+
+ Evaluate[$newDmdt = MassFlowRate];
+ Evaluate[ $syscount = 0 ];
+ Generate[A]; Solve[A]; SaveSolution[A];
+
+ Evaluate[$dmdtp = MassFlowRate];
+ Evaluate[$argVp = DeltaPhi];
+ PostOperation[Trailing];
+
+ Print[{$syscount, $circ, $dmdt, $argV},
+ Format "iter = %3g Circ = %7.5e Dmdt = %7.5e argV = %7.5e"];
+
+ While[ Norm[ $argV ] > 1e-3 && $syscount < 50] {
+
+ Test[ $syscount && Norm[$argV-$argVp] > 1e-3 ]
+ {Evaluate[$jac = Min[ ($dmdt-$dmdtp)/($argV-$argVp), 0.2 ] ] ;}
+ {Evaluate[$jac = 0.2];}
+ Print[{$jac}, Format "Jac = %7.5e"];
+
+ Evaluate[$newDmdt = $dmdt - $jac * $argV];
+ Evaluate[ $syscount = $syscount + 1 ];
+ Generate[A]; Solve[A]; SaveSolution[A];
+
+ Evaluate[$dmdtp = $dmdt];
+ Evaluate[$argVp = $argV];
+ PostOperation[Trailing];
+
+ Print[{$syscount, $circ, $dmdt, $argV},
+ Format "iter = %3g Circ = %7.5e Dmdt = %7.5e argV = %7.5e"];
+ }
+ // A python script is used to generate a format file "KJ.pos.opt"
+ // for the color map "KJ.pos" that depends
+ // on the run-time variable "phiTrailing".
+ Evaluate[ Python[$phiTrailing]{"generateOptFile.py"} ];
+ EndIf
+ }
+ }
+}
+
+PostProcessing{
+ {Name LPE; NameOfFormulation LPE;
+ Quantity{
+ {Name phi; Value {
+ Local{ [ {phi} ] ; In Dom_Vh; Jacobian Vol; } } }
+ { Name phiCont ; Value {
+ Local { [ { phiCont } ] ; In Dom_Vh ; Jacobian Vol ; } } }
+ { Name phiDisc ; Value {
+ Local { [ { phiDisc } ] ; In Dom_Vh ; Jacobian Vol ; } } }
+ {Name velocity; Value {
+ Local { [ {d phi} ]; In Dom_Vh; Jacobian Vol; } } }
+ {Name pressure; Value {
+ Local { [ -0.5*rho[]*SquNorm[ {d phi} ]];
+ In Dom_Vh; Jacobian Vol; } } }
+ {Name Angle; Value {
+ Local{ [ argVTrail[{d phi}] ];
+ In Dom_Vh; Jacobian Vol; } } }
+
+ { Name Circ; Value { Local { [ {Circ} ]; In Sur_Cut; } } }
+
+ { Name Dmdt; Value { Local { [ {Dmdt} ]; In Sur_Cut; } } }
+
+ // Kutta-Jukowski approximation for Lift
+ { Name LiftKJ; Value { Local { [ -rho[]*{Circ}*Velocity ];
+ In Sur_Cut; } } }
+
+ // Lift computed with the real pressure field
+ { Name Lift;
+ Value {
+ Integral { [ -0.5*rho[]*SquNorm[{d phi}]*CompY[Normal[]] ];
+ In Dom_Vh ; Jacobian Sur ; Integration Int; }
+ }
+ }
+
+ { Name circulation;
+ Value {
+ Integral { [ {d phi} * Tangent[] ] ;
+ In Dom_Vh ; Jacobian Sur ; Integration Int; }
+ }
+ }
+
+ }
+ }
+}
+
+PostOperation{
+ {Name LPE; NameOfPostProcessing LPE;
+ Operation{
+
+ Print[Circ, OnRegion Sur_Cut, File > "output.txt", Color "Ivory",
+ Format Table, SendToServer "Output/Circ" ];
+
+ Print[Dmdt, OnRegion Sur_Cut, File > "output.txt", Color "Ivory",
+ Format Table, SendToServer "Output/Dmdt" ];
+
+ Print[LiftKJ, OnRegion Sur_Cut, File > "output.txt", Color "Ivory",
+ Format Table, SendToServer "Output/LiftKJ" ];
+
+ Print[Lift[GammaAirf], OnGlobal, Format Table,
+ File > "output.txt", Color "Ivory",
+ SendToServer "Output/Lift"];
+
+ // Uncomment these lines to have more color maps
+ //Print[phi, OnElementsOf Vol_rho, File "phi.pos"];
+ //Print[phiDisc, OnElementsOf Vol_rho, File "phiDisc.pos"];
+ //Print[phiCont, OnElementsOf Vol_rho, File "phiCont.pos"];
+ //Print[pressure, OnElementsOf Vol_rho, File "p.pos"];
+
+ Print[ velocity, OnElementsOf Vol_rho, File "velocity.pos"];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].VectorType = 1;",
+ "View[l].LineWidth = 2;",
+ "View[l].ArrowSizeMax = 100;",
+ "View[l].CenterGlyphs = 1;"],
+ File "tmp.geo", LastTimeStepOnly] ;
+
+ If( Flag_Circ || ( Flag_Object == 0 ) )
+ //same condition as in the Resolution section
+
+ Print[phi, OnElementsOf Vol_rho, File "phi.pos"];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].IntervalsType = 3;",
+ "View[l].NbIso = 40;"],
+ File "tmp.geo", LastTimeStepOnly] ;
+
+ Else
+ Print[phi, OnElementsOf Vol_rho, File "KJ.pos"];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].IntervalsType = 3;",
+ "View[l].NbIso = 40;"],
+ File "tmp.geo", LastTimeStepOnly] ;
+ EndIf
+ }
+ }
+}
+
+PostOperation{ // for the Airfoil model
+ {Name Trailing; NameOfPostProcessing LPE;
+ Operation{
+ Print[Circ, OnRegion Sur_Cut, File > "KJiter.txt", Format Table,
+ StoreInVariable $circ,
+ SendToServer "Output/Circ" ];
+ Print[Dmdt, OnRegion Sur_Cut, File > "KJiter.txt", Format Table,
+ StoreInVariable $dmdt,
+ SendToServer "Output/Dmdt" ];
+ Print[phi, OnPoint{1.0001,0,0}, File > "KJiter.txt", Color "Ivory",
+ StoreInVariable $phiTrailing,
+ Format Table, SendToServer "Output/PhiTrailing"];
+
+ Print[Angle, OnPoint{1.0001,0,0}, File > "KJiter.txt", Color "Ivory",
+ StoreInVariable $argV,
+ Format Table, SendToServer "Output/argVTrailing"];
+ }
+ }
+}
+
+
+DefineConstant[
+ R_ = {"LPE", Name "GetDP/1ResolutionChoices", Visible 0},
+ C_ = {"-solve -pos", Name "GetDP/9ComputeCommand", Visible 0},
+ P_ = {"LPE", Name "GetDP/2PostOperationChoices", Visible 0}
+];
+
diff --git a/PotentialFlow/magnus_common.pro b/PotentialFlow/magnus_common.pro
new file mode 100644
index 0000000000000000000000000000000000000000..b9c659b20ef1c87b8df1d90fdef029ab3a4034b2
--- /dev/null
+++ b/PotentialFlow/magnus_common.pro
@@ -0,0 +1,27 @@
+cm = 1e-2;
+deg = Pi/180;
+kmh = 10./36.;
+
+Flag_Circ =
+ DefineNumber[0, Name "Model/Impose circulation", Choices {0, 1} ];
+Flag_Object =
+ DefineNumber[1, Name "Model/Object",
+ Choices {0="Cylinder", 1="Airfoil"} ];
+
+BoxSize =
+ DefineNumber[7, Name "Model/Air box size [m]", Help "In flow direction."];
+Velocity = kmh*
+ DefineNumber[100, Name "Model/V", Label "Model/Flow velocity [km/h]"];
+Incidence = -deg*
+ DefineNumber[10, Name "Model/Angle of attack [deg]", Visible Flag_Object];
+
+// Negative circulation for a positive lift.
+DeltaPhi = (-1)*
+ DefineNumber[10, Name "Model/Circ", Visible Flag_Circ,
+ Min 0, Max 20, Step 2,
+ Label "Circulation around foil [m^2/s]"];
+
+MassFlowRate =
+ DefineNumber[-100, Name "Model/Dmdt", Visible !Flag_Circ,
+ Label "Mass flow rate [kg/s]"];
+
diff --git a/PotentialFlow/nacaAirfoil.geo b/PotentialFlow/nacaAirfoil.geo
new file mode 100644
index 0000000000000000000000000000000000000000..9ef61996e0759624d206570cf17e658117423d36
--- /dev/null
+++ b/PotentialFlow/nacaAirfoil.geo
@@ -0,0 +1,222 @@
+
+
+Point(1) = {1.000000,0.000000,0,lca};
+Point(2) = {0.999753,-0.000035,0,lca};
+Point(3) = {0.999013,-0.000141,0,lca};
+Point(4) = {0.997781,-0.000317,0,lca};
+Point(5) = {0.996057,-0.000562,0,lca};
+Point(6) = {0.993844,-0.000876,0,lca};
+Point(7) = {0.991144,-0.001258,0,lca};
+Point(8) = {0.987958,-0.001707,0,lca};
+Point(9) = {0.984292,-0.002222,0,lca};
+Point(10) = {0.980147,-0.002801,0,lca};
+Point(11) = {0.975528,-0.003443,0,lca};
+Point(12) = {0.970440,-0.004147,0,lca};
+Point(13) = {0.964888,-0.004909,0,lca};
+Point(14) = {0.958877,-0.005729,0,lca};
+Point(15) = {0.952414,-0.006603,0,lca};
+Point(16) = {0.945503,-0.007531,0,lca};
+Point(17) = {0.938153,-0.008510,0,lca};
+Point(18) = {0.930371,-0.009537,0,lca};
+Point(19) = {0.922164,-0.010610,0,lca};
+Point(20) = {0.913540,-0.011726,0,lca};
+Point(21) = {0.904508,-0.012883,0,lca};
+Point(22) = {0.895078,-0.014079,0,lca};
+Point(23) = {0.885257,-0.015310,0,lca};
+Point(24) = {0.875056,-0.016574,0,lca};
+Point(25) = {0.864484,-0.017868,0,lca};
+Point(26) = {0.853553,-0.019189,0,lca};
+Point(27) = {0.842274,-0.020535,0,lca};
+Point(28) = {0.830656,-0.021904,0,lca};
+Point(29) = {0.818712,-0.023291,0,lca};
+Point(30) = {0.806454,-0.024694,0,lca};
+Point(31) = {0.793893,-0.026111,0,lca};
+Point(32) = {0.781042,-0.027539,0,lca};
+Point(33) = {0.767913,-0.028974,0,lca};
+Point(34) = {0.754521,-0.030414,0,lca};
+Point(35) = {0.740877,-0.031856,0,lca};
+Point(36) = {0.726995,-0.033296,0,lca};
+Point(37) = {0.712890,-0.034733,0,lca};
+Point(38) = {0.698574,-0.036163,0,lca};
+Point(39) = {0.684062,-0.037582,0,lca};
+Point(40) = {0.669369,-0.038988,0,lca};
+Point(41) = {0.654508,-0.040378,0,lca};
+Point(42) = {0.639496,-0.041747,0,lca};
+Point(43) = {0.624345,-0.043094,0,lca};
+Point(44) = {0.609072,-0.044414,0,lca};
+Point(45) = {0.593691,-0.045705,0,lca};
+Point(46) = {0.578217,-0.046962,0,lca};
+Point(47) = {0.562667,-0.048182,0,lca};
+Point(48) = {0.547054,-0.049362,0,lca};
+Point(49) = {0.531395,-0.050499,0,lca};
+Point(50) = {0.515705,-0.051587,0,lca};
+Point(51) = {0.500000,-0.052625,0,lca};
+Point(52) = {0.484295,-0.053608,0,lca};
+Point(53) = {0.468605,-0.054534,0,lca};
+Point(54) = {0.452946,-0.055397,0,lca};
+Point(55) = {0.437333,-0.056195,0,lca};
+Point(56) = {0.421783,-0.056924,0,lca};
+Point(57) = {0.406309,-0.057581,0,lca};
+Point(58) = {0.390928,-0.058163,0,lca};
+Point(59) = {0.375655,-0.058666,0,lca};
+Point(60) = {0.360504,-0.059087,0,lca};
+Point(61) = {0.345492,-0.059424,0,lca};
+Point(62) = {0.330631,-0.059674,0,lca};
+Point(63) = {0.315938,-0.059834,0,lca};
+Point(64) = {0.301426,-0.059902,0,lca};
+Point(65) = {0.287110,-0.059876,0,lca};
+Point(66) = {0.273005,-0.059754,0,lca};
+Point(67) = {0.259123,-0.059535,0,lca};
+Point(68) = {0.245479,-0.059217,0,lca};
+Point(69) = {0.232087,-0.058799,0,lca};
+Point(70) = {0.218958,-0.058280,0,lca};
+Point(71) = {0.206107,-0.057661,0,lca};
+Point(72) = {0.193546,-0.056940,0,lca};
+Point(73) = {0.181288,-0.056119,0,lca};
+Point(74) = {0.169344,-0.055197,0,lca};
+Point(75) = {0.157726,-0.054176,0,lca};
+Point(76) = {0.146447,-0.053056,0,lca};
+Point(77) = {0.135516,-0.051839,0,lca};
+Point(78) = {0.124944,-0.050527,0,lca};
+Point(79) = {0.114743,-0.049121,0,lca};
+Point(80) = {0.104922,-0.047624,0,lca};
+Point(81) = {0.095492,-0.046037,0,lca};
+Point(82) = {0.086460,-0.044364,0,lca};
+Point(83) = {0.077836,-0.042608,0,lca};
+Point(84) = {0.069629,-0.040770,0,lca};
+Point(85) = {0.061847,-0.038854,0,lca};
+Point(86) = {0.054497,-0.036863,0,lca};
+Point(87) = {0.047586,-0.034800,0,lca};
+Point(88) = {0.041123,-0.032668,0,lca};
+Point(89) = {0.035112,-0.030471,0,lca};
+Point(90) = {0.029560,-0.028212,0,lca};
+Point(91) = {0.024472,-0.025893,0,lca};
+Point(92) = {0.019853,-0.023517,0,lca};
+Point(93) = {0.015708,-0.021088,0,lca};
+Point(94) = {0.012042,-0.018607,0,lca};
+Point(95) = {0.008856,-0.016078,0,lca};
+Point(96) = {0.006156,-0.013503,0,lca};
+Point(97) = {0.003943,-0.010884,0,lca};
+Point(98) = {0.002219,-0.008223,0,lca};
+Point(99) = {0.000987,-0.005521,0,lca};
+Point(100) = {0.000247,-0.002779,0,lca};
+Point(101) = {0.000000,0.000000,0,lca};
+Point(102) = {0.000247,0.002779,0,lca};
+Point(103) = {0.000987,0.005521,0,lca};
+Point(104) = {0.002219,0.008223,0,lca};
+Point(105) = {0.003943,0.010884,0,lca};
+Point(106) = {0.006156,0.013503,0,lca};
+Point(107) = {0.008856,0.016078,0,lca};
+Point(108) = {0.012042,0.018607,0,lca};
+Point(109) = {0.015708,0.021088,0,lca};
+Point(110) = {0.019853,0.023517,0,lca};
+Point(111) = {0.024472,0.025893,0,lca};
+Point(112) = {0.029560,0.028212,0,lca};
+Point(113) = {0.035112,0.030471,0,lca};
+Point(114) = {0.041123,0.032668,0,lca};
+Point(115) = {0.047586,0.034800,0,lca};
+Point(116) = {0.054497,0.036863,0,lca};
+Point(117) = {0.061847,0.038854,0,lca};
+Point(118) = {0.069629,0.040770,0,lca};
+Point(119) = {0.077836,0.042608,0,lca};
+Point(120) = {0.086460,0.044364,0,lca};
+Point(121) = {0.095492,0.046037,0,lca};
+Point(122) = {0.104922,0.047624,0,lca};
+Point(123) = {0.114743,0.049121,0,lca};
+Point(124) = {0.124944,0.050527,0,lca};
+Point(125) = {0.135516,0.051839,0,lca};
+Point(126) = {0.146447,0.053056,0,lca};
+Point(127) = {0.157726,0.054176,0,lca};
+Point(128) = {0.169344,0.055197,0,lca};
+Point(129) = {0.181288,0.056119,0,lca};
+Point(130) = {0.193546,0.056940,0,lca};
+Point(131) = {0.206107,0.057661,0,lca};
+Point(132) = {0.218958,0.058280,0,lca};
+Point(133) = {0.232087,0.058799,0,lca};
+Point(134) = {0.245479,0.059217,0,lca};
+Point(135) = {0.259123,0.059535,0,lca};
+Point(136) = {0.273005,0.059754,0,lca};
+Point(137) = {0.287110,0.059876,0,lca};
+Point(138) = {0.301426,0.059902,0,lca};
+Point(139) = {0.315938,0.059834,0,lca};
+Point(140) = {0.330631,0.059674,0,lca};
+Point(141) = {0.345492,0.059424,0,lca};
+Point(142) = {0.360504,0.059087,0,lca};
+Point(143) = {0.375655,0.058666,0,lca};
+Point(144) = {0.390928,0.058163,0,lca};
+Point(145) = {0.406309,0.057581,0,lca};
+Point(146) = {0.421783,0.056924,0,lca};
+Point(147) = {0.437333,0.056195,0,lca};
+Point(148) = {0.452946,0.055397,0,lca};
+Point(149) = {0.468605,0.054534,0,lca};
+Point(150) = {0.484295,0.053608,0,lca};
+Point(151) = {0.500000,0.052625,0,lca};
+Point(152) = {0.515705,0.051587,0,lca};
+Point(153) = {0.531395,0.050499,0,lca};
+Point(154) = {0.547054,0.049362,0,lca};
+Point(155) = {0.562667,0.048182,0,lca};
+Point(156) = {0.578217,0.046962,0,lca};
+Point(157) = {0.593691,0.045705,0,lca};
+Point(158) = {0.609072,0.044414,0,lca};
+Point(159) = {0.624345,0.043094,0,lca};
+Point(160) = {0.639496,0.041747,0,lca};
+Point(161) = {0.654508,0.040378,0,lca};
+Point(162) = {0.669369,0.038988,0,lca};
+Point(163) = {0.684062,0.037582,0,lca};
+Point(164) = {0.698574,0.036163,0,lca};
+Point(165) = {0.712890,0.034733,0,lca};
+Point(166) = {0.726995,0.033296,0,lca};
+Point(167) = {0.740877,0.031856,0,lca};
+Point(168) = {0.754521,0.030414,0,lca};
+Point(169) = {0.767913,0.028974,0,lca};
+Point(170) = {0.781042,0.027539,0,lca};
+Point(171) = {0.793893,0.026111,0,lca};
+Point(172) = {0.806454,0.024694,0,lca};
+Point(173) = {0.818712,0.023291,0,lca};
+Point(174) = {0.830656,0.021904,0,lca};
+Point(175) = {0.842274,0.020535,0,lca};
+Point(176) = {0.853553,0.019189,0,lca};
+Point(177) = {0.864484,0.017868,0,lca};
+Point(178) = {0.875056,0.016574,0,lca};
+Point(179) = {0.885257,0.015310,0,lca};
+Point(180) = {0.895078,0.014079,0,lca};
+Point(181) = {0.904508,0.012883,0,lca};
+Point(182) = {0.913540,0.011726,0,lca};
+Point(183) = {0.922164,0.010610,0,lca};
+Point(184) = {0.930371,0.009537,0,lca};
+Point(185) = {0.938153,0.008510,0,lca};
+Point(186) = {0.945503,0.007531,0,lca};
+Point(187) = {0.952414,0.006603,0,lca};
+Point(188) = {0.958877,0.005729,0,lca};
+Point(189) = {0.964888,0.004909,0,lca};
+Point(190) = {0.970440,0.004147,0,lca};
+Point(191) = {0.975528,0.003443,0,lca};
+Point(192) = {0.980147,0.002801,0,lca};
+Point(193) = {0.984292,0.002222,0,lca};
+Point(194) = {0.987958,0.001707,0,lca};
+Point(195) = {0.991144,0.001258,0,lca};
+Point(196) = {0.993844,0.000876,0,lca};
+Point(197) = {0.996057,0.000562,0,lca};
+Point(198) = {0.997781,0.000317,0,lca};
+Point(199) = {0.999013,0.000141,0,lca};
+Point(200) = {0.999753,0.000035,0,lca};
+
+//Points numerotation:
+numberLeadingEgde = 101 ;
+numberLowerSurface = 81 ;
+numberUpperSurface = 121 ;
+
+//Distances:
+distanceTrailingLowerPoint = 0.908262 ;
+distanceLowerUpperPoint = 0.222905 ;
+distanceUpperTrailingPoint = 0.908262 ;
+
+Line(1) = { 1 ... 81 };
+Line(2) = { 81 ... 101 };
+Line(3) = { 101 ... 121 };
+Line(4) = { 121 ... 200, 1 };
+
+Rotate { {0,0,1}, {1,0,0}, Incidence } { Line { 1 ... 4 } ; }
+
+/* Line Loop(1) = {1,2,3,4}; */
+/* Plane Surface(1) = {1}; */
diff --git a/PotentialFlow/screenshot1.png b/PotentialFlow/screenshot1.png
new file mode 100644
index 0000000000000000000000000000000000000000..1d57354aae0328c97969aa1172957ee3c0cb71ac
Binary files /dev/null and b/PotentialFlow/screenshot1.png differ
diff --git a/PotentialFlow/screenshot1_512.png b/PotentialFlow/screenshot1_512.png
new file mode 100644
index 0000000000000000000000000000000000000000..dd9835b9a054e89193c321e3d6aa447779f6f509
Binary files /dev/null and b/PotentialFlow/screenshot1_512.png differ
diff --git a/PotentialFlow/screenshot2.png b/PotentialFlow/screenshot2.png
new file mode 100644
index 0000000000000000000000000000000000000000..f6d2095515c7249100035fda62937111963170d1
Binary files /dev/null and b/PotentialFlow/screenshot2.png differ
diff --git a/PotentialFlow/screenshot2_512.png b/PotentialFlow/screenshot2_512.png
new file mode 100644
index 0000000000000000000000000000000000000000..18f399d7bcccf30811dafbecef2353c594692693
Binary files /dev/null and b/PotentialFlow/screenshot2_512.png differ