tutorial PotentialFlow

PotentialFlow/generateOptFile.py

+#!/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()
+
+
+
+
+
+
+
+
+

PotentialFlow/magnus.geo

+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 };
+

PotentialFlow/magnus.pro

+/* -------------------------------------------------------------------
+   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}
+];
+

PotentialFlow/magnus_common.pro

+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]"];
+

PotentialFlow/nacaAirfoil.geo

+
+
+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}; */