From 2765e1d4e67dc9f16585989b89d2355151ea1bb6 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Mon, 25 Nov 2019 23:36:52 +0100
Subject: [PATCH] better geo parameters + cleanup

---
 SlidingSurface3D/rfpm.geo | 48 ++++++++++++------------
 SlidingSurface3D/rfpm.pro | 78 ++++++++++++++++++---------------------
 2 files changed, 61 insertions(+), 65 deletions(-)

diff --git a/SlidingSurface3D/rfpm.geo b/SlidingSurface3D/rfpm.geo
index a168ace..abad140 100644
--- a/SlidingSurface3D/rfpm.geo
+++ b/SlidingSurface3D/rfpm.geo
@@ -11,8 +11,10 @@ AngleMagnet = 0.8 * AnglePole;
 
 // radial dimensions
 lRotor = 44 * mm;
-lMagnet = 10 * mm;
-lGap = 20 * mm;
+lMagnetIn = 10 * mm;
+lMagnetOut = 5 * mm;
+lGapRotor = 7 * mm;
+lGapStator = 3 * mm;
 lAirRing = 1 * mm;
 lStator = 30 * mm;
 lInf = 30 * mm;
@@ -26,8 +28,8 @@ hTot = 130 * mm;
 Cylinder(1) = {0,0,0, 0,0,hRotor, lRotor, AnglePole};
 
 // magnet
-Cylinder(2) = {0,0,0, 0,0,hMagnet, lRotor + lMagnet / 2, AngleMagnet};
-Cylinder(3) = {0,0,0, 0,0,hMagnet, lRotor - lMagnet / 2, AngleMagnet};
+Cylinder(2) = {0,0,0, 0,0,hMagnet, lRotor + lMagnetOut, AngleMagnet};
+Cylinder(3) = {0,0,0, 0,0,hMagnet, lRotor - lMagnetIn, AngleMagnet};
 BooleanDifference(4) = { Volume{2}; Delete; }{ Volume{3}; Delete; };
 Rotate{ {0,0,1}, {0,0,0}, (AnglePole - AngleMagnet) / 2 }{ Volume{4}; }
 Physical Volume("Magnet", 2) = 4;
@@ -37,7 +39,7 @@ BooleanDifference(5) = { Volume{1}; Delete; }{ Volume{4}; };
 Physical Volume("Rotor", 1) = 5;
 
 // rotor air
-Cylinder(6) = {0,0,0, 0,0,hTot, lRotor + lGap / 2, AnglePole};
+Cylinder(6) = {0,0,0, 0,0,hTot, lRotor + lGapRotor, AnglePole};
 BooleanDifference(7) = { Volume{6}; Delete; }{ Volume{4, 5}; };
 Physical Volume("AirRotor", 4) = 7;
 
@@ -46,20 +48,20 @@ VolumesRotor() = {4,5,7};
 BooleanFragments{ Volume{VolumesRotor()}; Delete; }{}
 
 // air ring
-Cylinder(8) = {0,0,0, 0,0,hTot, lRotor + lGap / 2 , AnglePole};
-Cylinder(9) = {0,0,0, 0,0,hTot, lRotor + lGap / 2 + lAirRing, AnglePole};
+Cylinder(8) = {0,0,0, 0,0,hTot, lRotor + lGapRotor , AnglePole};
+Cylinder(9) = {0,0,0, 0,0,hTot, lRotor + lGapRotor + lAirRing, AnglePole};
 BooleanDifference(10) = { Volume{9}; Delete; }{ Volume{8}; Delete; };
 Physical Volume("AirRing", 6) = 10;
 
 // stator
-Cylinder(11) = {0,0,0, 0,0,hRotor, lRotor + lGap, AnglePole};
-Cylinder(12) = {0,0,0, 0,0,hRotor, lRotor + lGap + lStator, AnglePole};
+Cylinder(11) = {0,0,0, 0,0,hRotor, lRotor + lGapRotor + lGapStator, AnglePole};
+Cylinder(12) = {0,0,0, 0,0,hRotor, lRotor + lGapRotor + lGapStator + lStator, AnglePole};
 BooleanDifference(13) = { Volume{12}; Delete; }{ Volume{11}; Delete; };
 Physical Volume("Stator", 3) = 13;
 
 // stator air
-Cylinder(14) = {0,0,0, 0,0,hTot, lRotor + lGap / 2 + lAirRing, AnglePole};
-Cylinder(15) = {0,0,0, 0,0,hTot, lRotor + lGap + lStator + lInf, AnglePole};
+Cylinder(14) = {0,0,0, 0,0,hTot, lRotor + lGapRotor + lAirRing, AnglePole};
+Cylinder(15) = {0,0,0, 0,0,hTot, lRotor + lGapRotor + lGapStator + lStator + lInf, AnglePole};
 BooleanDifference(16) = { Volume{15}; Delete; }{ Volume{14}; Delete; };
 BooleanDifference(17) = { Volume{16}; Delete; }{ Volume{13}; };
 Physical Volume("AirStator", 5) = 17;
@@ -130,20 +132,20 @@ Function Cube2Face
     For p In { 0 : #points() - 1 }
       point() = Point { points(p) };
       Call ChangeCoordinates;
-      If((CoordinateSystem == 1) && (coord(0)==0))
-  	    coord(1) = bbox(1);
+      If((CoordinateSystem == 1) && (coord(0) == 0))
+        coord(1) = bbox(1);
       EndIf
-      If((coord(0)>bbox(0)-tol) &&
-		  (coord(1)>bbox(1)-tol) &&
-		  (coord(2)>bbox(2)-tol) &&
-          (coord(0)<bbox(3)+tol) &&
-		  (coord(1)<bbox(4)+tol) &&
-		  (coord(2)<bbox(5)+tol))
-         NbIn++;
+      If((coord(0) > bbox(0) - tol) &&
+         (coord(1) > bbox(1) - tol) &&
+         (coord(2) > bbox(2) - tol) &&
+         (coord(0) < bbox(3) + tol) &&
+         (coord(1) < bbox(4) + tol) &&
+         (coord(2) < bbox(5) + tol))
+        NbIn++;
       EndIf
     EndFor
-	If(NbIn == #points())
-	  face() += surfaces(s);
+    If(NbIn == #points())
+      face() += surfaces(s);
     EndIf
   EndFor
 Return
@@ -278,4 +280,4 @@ treeLines() += CombinedBoundary { Physical Surface { 16 }; } ;
 treeLines() += CombinedBoundary { Physical Surface { 17 }; } ;
 treeLines() += CombinedBoundary { Physical Surface { 18 }; } ;
 
-Physical Line("treeLines", 22) = { treeLines() };
+Physical Line("TreeLines", 22) = { treeLines() };
diff --git a/SlidingSurface3D/rfpm.pro b/SlidingSurface3D/rfpm.pro
index cbc46fb..eef5946 100644
--- a/SlidingSurface3D/rfpm.pro
+++ b/SlidingSurface3D/rfpm.pro
@@ -10,31 +10,29 @@ Include "rfpm_common.pro";
 
 DefineConstant[
   R_ = {"Static", Choices {"Static", "QuasiStatic"}, Visible 1,
-        Name "GetDP/1ResolutionChoices"},
-  C_ = {"-solve -v 4 -pos", Name "GetDP/9ComputeCommand", Visible 0}
+    Name "GetDP/1ResolutionChoices"},
+  C_ = {"-solve -v 4 -pos -v2", Visible 0,
+    Name "GetDP/9ComputeCommand"},
   P_ = { "Check_Periodicity", Choices {"Fields", "Check_Periodicity"}, Visible 0,
-         Name "GetDP/2PostOperationChoices"}
+    Name "GetDP/2PostOperationChoices"}
 ];
 
-Flag_QuasiStatic = !StrCmp[ R_,"QuasiStatic"] ;
-RotorPosition = deg*
-  DefineNumber[20, Name "Options/Rotor position [deg]",
-               Visible !Flag_QuasiStatic];
-AngularStep = deg*
-  DefineNumber[1, Name "Options/1Angular step [deg]",
-               Visible Flag_QuasiStatic];
-NumStep =
-  DefineNumber[20, Name "Options/1Number of steps",
-               Visible Flag_QuasiStatic];
+Flag_QuasiStatic = !StrCmp[ R_, "QuasiStatic"];
 
-rpm = 2.*Pi/60. ;
-w1 = 1000*rpm ;
-Mag_DTime = AngularStep/Fabs[w1] ;
-Mag_TimeMax = Mag_DTime*NumStep;
+RotorPosition = deg * DefineNumber[0, Name "Options/Rotor position [deg]",
+  Visible !Flag_QuasiStatic];
+AngularStep = deg * DefineNumber[2, Name "Options/1Angular step [deg]",
+  Visible Flag_QuasiStatic];
+NumStep = DefineNumber[20, Name "Options/1Number of steps",
+  Visible Flag_QuasiStatic];
 
+rpm = 2. * Pi / 60. ;
+w1 = 1000 * rpm ;
+Mag_DTime = AngularStep / Fabs[w1] ;
+Mag_TimeMax = Mag_DTime * NumStep;
 
 Group {
- // Geometrical regions
+  // Geometrical regions
   Rotor     = Region[1];
   Magnet    = Region[2];
   Stator    = Region[3];
@@ -42,9 +40,9 @@ Group {
   AirStator = Region[5];
   AirRing   = Region[6];
 
-  Outer           = Region[10];
-  Top             = Region[11];
-  Bottom          = Region[12];
+  Outer               = Region[10];
+  Top                 = Region[11];
+  Bottom              = Region[12];
   Sur_RotorPerMaster  = Region[13];
   Sur_RotorPerSlave   = Region[14];
   Sur_StatorPerMaster = Region[15];
@@ -72,7 +70,6 @@ Group {
                       Sur_RotorPerMaster, Sur_RotorPerSlave} ] ;
   Dom_Hcurl_a = Region[ { Vol_Mag, Sur_Link } ]; // No Neumann surface
 
-
   Sur_Tree_Sliding = Region[ { Sur_SlidingMaster, Sur_SlidingSlave } ];
   Lin_Tree = Region[ Lin_TreeLines ];
   Sur_Tree = Region[ { Sur_Link, Sur_Dirichlet } ];
@@ -92,7 +89,6 @@ Function{
   js[] = 0;
 }
 
-
 Function {
   // Sliding surface:
   // a1 is the angular span of Region and RegionRef.
@@ -122,11 +118,11 @@ Function {
   Coef[] = ( fFloor[ ((Atan2[ Y[], X[] ] - a2[] )/a1) ] % 2 ) ? Periodicity : 1. ;
 }
 
-
 Group {
   DefineGroup[ DomainInf ];
   DefineVariable[ Val_Rint, Val_Rext ];
 }
+
 Jacobian {
   { Name Vol ;
     Case { { Region DomainInf ; Jacobian VolSphShell {Val_Rint, Val_Rext} ; }
@@ -195,12 +191,11 @@ FunctionSpace {
 	Support Dom_Hcurl_a ; Entity EdgesOf[ All ]; }
     }
     Constraint {
-      { NameOfCoef ae;  EntityType EdgesOf ;
-      	NameOfConstraint a; }
-	  If( Gauge==1 ) // spanning tree gauging
-		{ NameOfCoef ae;  EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
-		  NameOfConstraint GaugeCondition_a ; }
-	  EndIf
+      { NameOfCoef ae;  EntityType EdgesOf ; NameOfConstraint a; }
+      If(Gauge == 1) // spanning tree gauging
+        { NameOfCoef ae;  EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
+          NameOfConstraint GaugeCondition_a ; }
+      EndIf
     }
   }
 }
@@ -221,7 +216,6 @@ Formulation {
   }
 }
 
-
 Resolution {
   { Name Static ;
     System {
@@ -272,9 +266,9 @@ PostProcessing {
   { Name MagSta_a ; NameOfFormulation MagSta_a ;
     PostQuantity {
       { Name b ; Value { Local { [ {d a} ];
-	    In Vol_Mag ; Jacobian Vol; } } }
+            In Vol_Mag ; Jacobian Vol; } } }
       { Name bsurf ; Value { Local { [ {d a} ];
-	    In Sur_Link ; Jacobian Sur; } } }
+            In Sur_Link ; Jacobian Sur; } } }
       { Name asurf ; Value { Local { [ {a} ];
 	    In Sur_Link ; Jacobian Sur; } } }
       { Name br ; Value { Local { [ br[] ];
@@ -302,8 +296,6 @@ PostOperation Fields UsingPost MagSta_a {
   	File "tmp.geo", LastTimeStepOnly] ;
 }
 
-
-
 PostOperation {
   { Name pos; NameOfPostProcessing MagSta_a;
     Operation {
@@ -318,13 +310,15 @@ PostOperation {
 }
 
 PostOperation Check_Periodicity UsingPost MagSta_a {
-  // Print the tree for debugging purposes (the group name
-  PrintGroup[ _CO_Entity_39 , In Vol_Tree, File "Tree.pos"];
-
-  Echo[ Str["l=PostProcessing.NbViews-1;",
-			"View[l].ColorTable = { DarkRed };",
-			"View[l].LineWidth = 5;"] ,
-		File "tmp.geo", LastTimeStepOnly] ;
+  If(Gauge == 1)
+    // Print the tree for debugging purposes (the group name
+    PrintGroup[ _CO_Entity_39 , In Vol_Tree, File "Tree.pos"];
+
+    Echo[ Str["l=PostProcessing.NbViews-1;",
+              "View[l].ColorTable = { DarkRed };",
+              "View[l].LineWidth = 5;"] ,
+            File "tmp.geo", LastTimeStepOnly] ;
+  EndIf
 
   Print[ bsurf, OnElementsOf Region[ {Sur_SlidingMaster} ],
 	 LastTimeStepOnly, File "bsm.pos"] ;
-- 
GitLab