some improvements, got rid of Python

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.pro

@@ -2,14 +2,13 @@
    Tutorial 6 : Potential flow, lift and Magnus effect
 
    Features:
-   - Discontinuous scalar field 
+   - Potential flow, irrotational flow, lift and Magnus effect
    - Multivalued scalar field
-   - Potential flow
    - Non-linear iteration on an Associated global quantity
-   - Using a run-time variable in a Post-Operation (via Python) 
+   - Use a run-time variable in a Post-Operation 
 
    To compute the solution in a terminal: 
-     getdp magnus.pro -solve LPE -pos LPE
+     getdp magnus.pro -solve PotentialFlow -pos PotentialFlow
      gmsh magnus.geo velocity.pos
 
    To compute the solution interactively from the Gmsh GUI:
@@ -77,19 +76,21 @@
    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 
+   is located at the trailing edge of the airfoil,.
+   This is true 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.
+   close to the trailing edge of the airfoil. 
    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";
@@ -196,6 +197,9 @@ FunctionSpace{
       	  NameOfConstraint DeltaPhi;}
       Else
         If( Flag_Object == 0 ) // if Cylinder only
+      // In case of the Airfoil, the contraint on Dmdt is omitted
+      // and replaced by an equation "Dmdt=$newDmdt" 
+      // See the "Resolution" section.
           {NameOfCoef Dmdt; EntityType GroupsOfNodesOf; 
       	  NameOfConstraint MassFlowRate;}
         EndIf
@@ -205,7 +209,7 @@ FunctionSpace{
 }
 
 Formulation{
-  {Name LPE; Type FemEquation;
+  {Name PotentialFlow; Type FemEquation;
     Quantity{
       {Name phi; Type Local; NameOfSpace Vh;}
       { Name phiCont ; Type Local ; NameOfSpace Vh[phiCont] ; }
@@ -228,9 +232,9 @@ Formulation{
 }
 
 Resolution{
-  {Name LPE;
+  {Name PotentialFlow;
     System{
-      {Name A; NameOfFormulation LPE;}
+      {Name A; NameOfFormulation PotentialFlow;}
     }
     Operation{
       InitSolution[A];
@@ -241,9 +245,11 @@ Resolution{
 
       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]
+	// The run-time variable $newDmdt is used in Generate[A]
+        // whereas $circ, $dmdt, $argV and $phiTrailing are evaluated 
+	// by the PostOperation[Trailing].
+	// A file containing "PhiTrailing = $phiTrailing;" is after all
+	// written on disk to be used for formatting the 
 
 	DeleteFile["KJiter.txt"];
 
@@ -256,7 +262,7 @@ Resolution{
         PostOperation[Trailing];
 
         Print[{$syscount, $circ, $dmdt, $argV}, 
-	    Format "iter = %3g Circ = %7.5e Dmdt = %7.5e argV = %7.5e"];
+	    Format "iter = %3g Circ = %5.2f Dmdt = %5.2f argV = %5.2e"];
 
         While[ Norm[ $argV ] > 1e-3 && $syscount < 50] {
 
@@ -273,20 +279,17 @@ Resolution{
 	  Evaluate[$argVp = $argV];
 	  PostOperation[Trailing];
 
-	  Print[{$syscount, $circ, $dmdt, $argV}, 
-		Format "iter = %3g Circ = %7.5e Dmdt = %7.5e argV = %7.5e"];
+	  Print[{$syscount, $circ, $dmdt, $jac, $argV}, 
+		Format "iter = %3g Circ = %5.2f Dmdt = %5.2f jac=%5.2f argV = %5.3e"];
 	}
-	// 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"} ];
+	Print[{$phiTrailing}, Format "PhiTrailing=%g;", File "PhiTrailing.txt"];
 	EndIf
     }
   }
 }
 
 PostProcessing{
-  {Name LPE; NameOfFormulation LPE;
+  {Name PotentialFlow; NameOfFormulation PotentialFlow;
     Quantity{
       {Name phi; Value { 
 	  Local{ [ {phi} ] ; In Dom_Vh; Jacobian Vol; } } }
@@ -296,6 +299,8 @@ PostProcessing{
  	  Local { [ { phiDisc } ] ; In Dom_Vh ; Jacobian Vol ; } } }
       {Name velocity; Value { 
 	  Local { [ {d phi} ]; In Dom_Vh; Jacobian Vol; } } }
+      {Name normVelocity; Value { 
+	  Local { [ Norm[{d phi}] ]; In Dom_Vh; Jacobian Vol; } } }
       {Name pressure; Value { 
 	  Local { [-0.5*rho[]*SquNorm[ {d phi} ]]; 
 	    In Dom_Vh; Jacobian Vol; } } }
@@ -331,7 +336,7 @@ PostProcessing{
 }
 
 PostOperation{
-  {Name LPE; NameOfPostProcessing LPE;
+  {Name PotentialFlow; NameOfPostProcessing PotentialFlow;
     Operation{
 
       Print[Circ, OnRegion Sur_Cut, File > "output.txt", Color "Ivory",
@@ -354,35 +359,41 @@ PostOperation{
       //Print[pressure, OnElementsOf Vol_rho, File "p.pos"];
 
       Print[ velocity, OnElementsOf Vol_rho, File "velocity.pos"];
-      Echo[ StrCat["l=PostProcessing.NbViews-1;",
+      Echo[ Str["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;"],
+      // Stagnation points are points where Norm[{d phi}] is zero
+      // This is better seen in log scale. 
+      Print[normVelocity, OnElementsOf Vol_rho, File "p.pos"];
+      Echo[ Str["l=PostProcessing.NbViews-1;",
+		"View[l].Name = 'stagnation points';",
+		"View[l].ScaleType = 2;"],
 	    File "tmp.geo", LastTimeStepOnly] ;
 
-      Else
+      // Show the isovalue phi=$PhiTrailing
+      // which should be perpendicular to the airfoil
+      // if the Kutta condition is fulfilled
       Print[phi, OnElementsOf Vol_rho, File "KJ.pos"];
-        Echo[ StrCat["l=PostProcessing.NbViews-1;",
+      Echo[Str["Merge 'PhiTrailing.txt';",
+	       "l=PostProcessing.NbViews-1;",
+	       "View[l].Name = 'isovalue phiTrailing';",
 	       "View[l].IntervalsType = 3;",
-		     "View[l].NbIso = 40;"],
+	       "Mesh.SurfaceEdges = 0;",
+	       "View[l].RangeType = 2;",
+	       "View[l].CustomMin = PhiTrailing;",
+	       "View[l].CustomMax = 1.001*PhiTrailing;",
+	       "View[l].NbIso = 10;"], 
 	   File "tmp.geo", LastTimeStepOnly] ;
-      EndIf
     }
   }
 }
 
 PostOperation{ // for the Airfoil model
-  {Name Trailing; NameOfPostProcessing LPE;
+  {Name Trailing; NameOfPostProcessing PotentialFlow;
     Operation{
        Print[Circ, OnRegion Sur_Cut, File > "KJiter.txt", Format Table,
      	     StoreInVariable $circ, 
@@ -403,8 +414,8 @@ PostOperation{ // for the Airfoil model
 
 
 DefineConstant[
-  R_ = {"LPE", Name "GetDP/1ResolutionChoices", Visible 0},
+  R_ = {"PotentialFlow", Name "GetDP/1ResolutionChoices", Visible 0},
   C_ = {"-solve -pos", Name "GetDP/9ComputeCommand", Visible 0},
-  P_ = {"LPE", Name "GetDP/2PostOperationChoices", Visible 0}
+  P_ = {"PotentialFlow", Name "GetDP/2PostOperationChoices", Visible 0}
 ];