simpler geometry + comments

SlidingSurface3D/rfpm.geo

 Include "rfpm_common.pro";
 
+SetFactory("OpenCASCADE");
 
-// Tolerance for the geometric identification. 
-tol = 1.1*SlidingSurfaceShift;
-
-// Auxiliary Gmsh functions for OpenCascade solid modelling
-
-Function CheckNumbering
-// Check whether BooleanFragments has preserved numbering as expected
-For num In {0:#Volumes()-1}
-   If( Volumes[num] != vv(num) )
-     Error("BooleanFragment changed volume tag %g -> %g", 
-	    Volumes[num], vv(num)); 
-     Abort;
-   EndIf
-EndFor
-/* EXPLAIN
-For num In {0:#Surfaces()-1}
-  If( Surfaces(num) != vv(num+#Volumes()) )
-    Error("BooleanFragment changed surface tag %g -> %g",
-		  Surfaces(num), vv(num+#Volumes()));
-    Abort;
-  EndIf
-EndFor
-*/
-Return
+Mesh.CharacteristicLengthMin = lc;
+Mesh.CharacteristicLengthMax = lc;
+tol = 1e-6;
+
+AnglePole = (ModelAngleMax - ModelAngleMin) * deg;
+AngleMagnet = 0.8 * AnglePole;
+
+// radial dimensions
+lRotor = 44 * mm;
+lMagnet = 10 * mm;
+lGap = 20 * mm;
+lAirRing = 1 * mm;
+lStator = 30 * mm;
+lInf = 30 * mm;
+
+// heights
+hRotor = 80 * mm;
+hMagnet = 40 * mm;
+hTot = 130 * mm;
+
+// initial rotor bloc
+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};
+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;
+
+// rotor
+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};
+BooleanDifference(7) = { Volume{6}; Delete; }{ Volume{4, 5}; };
+Physical Volume("AirRotor", 4) = 7;
+
+// make rotor parts conformal
+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};
+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};
+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};
+BooleanDifference(16) = { Volume{15}; Delete; }{ Volume{14}; Delete; };
+BooleanDifference(17) = { Volume{16}; Delete; }{ Volume{13}; };
+Physical Volume("AirStator", 5) = 17;
+
+// make stator parts and air ring conformal - but don't glue them to the rotor:
+// we need a "double" surface to act as "SlidingMaster" and "SlidingSlave",
+// which will be matched programmatically with a Link constraint in getdp after
+// rotating the rotor alone
+VolumesStator() = {10,13,17};
+BooleanFragments{ Volume{VolumesStator()}; Delete; }{}
+
+// We're done with the volumes. Due to internal relabelling during
+// BooleanFragments that cannot be controlled by the user, the consistent way to
+// identify the surface tags is to proceed geometrically. The surface tags that
+// are needed are in general pieces of the outer boundary of the solid where
+// specific boundary conditions are to be imposed.  Very often, the volume
+// occupied by a FE model is a parallelepipedic block in some coordinate system
+// (Cartesian, cylindrical, ...).  The boundary conditions are thus associated
+// with the faces of that block and can therefore be identified by means of thin
+// bounding boxes in the coordinate space.
 
 Function ChangeCoordinates
-  If( CoordinateSystem == 1) // cylindrical coordinates
-    coord() = {Hypot(point(0),point(1)),Atan2(point(1),point(0)),point(2)};
+  // return coordinates of a point depending on CoordinateSystem
+  //
+  // input: point() and CoordinateSystem
+  // output: coord()
+  If(CoordinateSystem == 1) // cylindrical coordinates
+    coord() = { Hypot(point(0), point(1)), Atan2(point(1), point(0)), point(2) };
   Else // Cartesian coordinates
     coord() = point();
   EndIf
 Return
 
-/*
-  Due to internal relabelling that cannot be controlled by the user,
-  the consistent way to identify the surface labels (tags) 
-  of a solid modelled with the OpenCascade factory is to proceed
-  geometrically.
-  The surface tags that are needed are in general 
-  pieces of the outer boundary of the solid 
-  where specific boundary conditions are to be imposed.
-  Very often, the volume occupied by a FE model 
-  is a parallelepipedic block in some coordinate system
-  (Cartesian, cylindrical, ...)
-  The boundary conditions are thus associated 
-  with the faces of that block and can therefore be identified
-  by means of thin bounding boxes in the coordinate space.
-  The function 'Volume2Cube' returns 
-  the parallelepipedic coordinate block of a list of volumes
-  'Volumes'.
-  The function 'Cube2Face' returns the list of faces of the model
-  that have all their points located inside a thin box encosing
-  the 'FaceId'th face of that coordinate block. 
-*/
 Function Volume2Cube
-  // 
+  // return the parallelepipedic coordinate block of a list of volumes
+  //
+  // input: Volumes()
+  // output: cube()
   points() = PointsOf{ Volume{ Volumes() }; };
-  point() = Point { points[0] };
+  point() = Point { points(0) };
   Call ChangeCoordinates;
-  cube() = {coord(0),coord(1),coord(2),coord(0),coord(1),coord(2)};
-  For p In { 1:#points()-1 }
-    point() = Point { points(p) };
+  cube() = {coord(0), coord(1), coord(2), coord(0), coord(1), coord(2)};
+  For p In {1 : #points() - 1}
+    point() = Point{ points(p) };
     Call ChangeCoordinates;
-    cube(0) = (coord(0) < cube(0) ? coord(0) : cube(0) );
-    cube(1) = (coord(1) < cube(1) ? coord(1) : cube(1) );
-    cube(2) = (coord(2) < cube(2) ? coord(2) : cube(2) );
-    cube(3) = (coord(0) > cube(3) ? coord(0) : cube(3) );
-    cube(4) = (coord(1) > cube(4) ? coord(1) : cube(4) );
-    cube(5) = (coord(2) > cube(5) ? coord(2) : cube(5) );
+    cube(0) = (coord(0) < cube(0) ? coord(0) : cube(0));
+    cube(1) = (coord(1) < cube(1) ? coord(1) : cube(1));
+    cube(2) = (coord(2) < cube(2) ? coord(2) : cube(2));
+    cube(3) = (coord(0) > cube(3) ? coord(0) : cube(3));
+    cube(4) = (coord(1) > cube(4) ? coord(1) : cube(4));
+    cube(5) = (coord(2) > cube(5) ? coord(2) : cube(5));
   EndFor
 Return
 
 Function Cube2Face
+  // return the list of surfaces of the model that have all their points located
+  // inside a thin box encosing the FaceId'th face of that coordinate block
+  //
+  // input: cube() and FaceId
+  // output: face()
   bbox() = cube();
-  bbox(Modulo(FaceId+3,6)) = bbox(FaceId);
-  //Printf("bbox=", bbox());
-  face()={};
+  bbox(Modulo(FaceId + 3, 6)) = bbox(FaceId);
+  face() = {};
   surfaces() = Boundary{ Volume{ Volumes() }; };
   For s In { 0:#surfaces()-1 }
     NbIn = 0;
     points() = PointsOf { Surface { surfaces(s) } ; };
-    For p In { 0:#points()-1 }
+    For p In { 0 : #points() - 1 }
       point() = Point { points(p) };
       Call ChangeCoordinates;
-      If( (CoordinateSystem == 1) && (coord(0)==0) )
+      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) )
+      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() )
+	If(NbIn == #points())
 	  face() += surfaces(s);
     EndIf
   EndFor
 Return
 
-
-SetFactory("OpenCASCADE");
-
-Solver.AutoMesh = 2;
-Mesh.CharacteristicLengthMin = lc;
-Mesh.CharacteristicLengthMax = lc;
-// Hide mesh in post-processing views
-Mesh.VolumeEdges = 0;
-Mesh.SurfaceEdges = 0;
-
-L_r = 44*mm;
-L_m = 40*mm;
-L_s = 44*mm;
-pp = 3;
-anglePole = Pi/pp; 
-e_in = 0;
-e_r = 3.5*mm;
-b_depth = 0.5*mm;
-e_m = 5.3*mm;
-e_sl = 0;
-e_ar = 1*mm;
-e_ss = 0.1*mm; // air ring width
-e_a = 1*mm; // width stator air gap
-e_s = 5*mm; // width stator core
-sRM = 0.75;
-A_infinity = 66.75*mm;
-R_infinity = 21.75*mm;
-
-// Radii
-R_in = e_in;
-R_r = R_in + e_r;
-R_m_in = R_r - b_depth;
-R_m = R_m_in + e_m;
-R_sl = R_m + e_sl;
-R_ar = R_sl + e_ar;
-R_ss = R_ar + e_ss;
-R_s = R_ar + e_a;
-R_o = R_s + e_s;
-
-// R O T O R
-
-p1r = newp; Point(p1r) = {R_in,0,0};
-p2r = newp; Point(p2r) = {R_in,0,L_r/2};
-p3r = newp; Point(p3r) = {R_r,0,L_r/2};
-p4r = newp; Point(p4r) = {R_r,0,0};
-
-l1r = newl; Line(l1r) = {p1r,p2r};
-l2r = newl; Line(l2r) = {p2r,p3r};
-l3r = newl; Line(l3r) = {p3r,p4r};
-l4r = newl; Line(l4r) = {p4r,p1r};
-
-ll_r = newll; Line Loop(ll_r) = {l1r, l2r, l3r, l4r};
-sRotor = news; Plane Surface(sRotor) = {ll_r};
-vRotorExt() = Extrude {{0, 0, 1}, {0, 0, 0}, anglePole} {
-  Surface{sRotor}; 
-};
-
-
-// Magnets
-
-p1m = newp; Point(p1m) = {R_m_in,0,0};
-p2m = newp; Point(p2m) = {R_m_in,0,L_m/2};
-p3m = newp; Point(p3m) = {R_m,0,L_m/2};
-p4m = newp; Point(p4m) = {R_m,0,0};
-
-l1m = newl; Line(l1m) = {p1m,p2m};
-l2m = newl; Line(l2m) = {p2m,p3m};
-l3m = newl; Line(l3m) = {p3m,p4m};
-l4m = newl; Line(l4m) = {p4m,p1m};
-
-ll_m = newll; Line Loop(ll_m) = {l1m, l2m, l3m, l4m};
-sMagnet = news; Plane Surface(sMagnet) = {ll_m};
-
-Rotate{ {0,0,1}, {0,0,0}, anglePole*(1-sRM)/2 } { 
-  Surface{sMagnet}; 
-}
-vMagnetExt() = Extrude {{0,0,1}, {0,0,0}, anglePole*sRM} {
-  Surface{sMagnet}; 
-};
-
-// Air rotor
-
-p1ar = newp; Point(p1ar) = {R_in,0,0};
-p2ar = newp; Point(p2ar) = {R_in,0,A_infinity/2};
-p3ar = newp; Point(p3ar) = {R_ar,0,A_infinity/2};
-p4ar = newp; Point(p4ar) = {R_ar,0,0};
-  
-l1ar = newl; Line(l1ar) = {p1ar,p2ar};
-l2ar = newl; Line(l2ar) = {p2ar,p3ar};
-l3ar = newl; Line(l3ar) = {p3ar,p4ar};
-l4ar = newl; Line(l4ar) = {p4ar,p1ar};
-
-ll_ar = newll; Line Loop(ll_ar) = {l1ar, l2ar, l3ar, l4ar};
-sAirRotor = news; Plane Surface(sAirRotor) = {ll_ar};
-vAirRotorExt() = Extrude {{0, 0, 1}, {0, 0, 0}, anglePole} {
-  Surface{sAirRotor};
-};
-
-// Create volumes with boolean operations
-
-vMagnet = vMagnetExt[1];
-vRotor = BooleanDifference{
-  Volume{vRotorExt[1]}; Delete; 
-}{ 
-  Volume{vMagnet}; 
-};
-vAirRotor = BooleanDifference{
-  Volume{vAirRotorExt[1]}; Delete; 
-}{ 
-  Volume{vRotor}; Volume{vMagnet};
-};
-
-// Glue up rotor pieces with 'BooleanFragment'. 
-// This operation is necessary, as it connects the different 
-// parts of the model together,
-// but it entails an internal renumbering of surfaces, lines, etc...
-// Hence the need for a geometrical identification afterwards (see below). 
-
-Volumes() = { vMagnet, vRotor, vAirRotor };
-VolumesRotor() = BooleanFragments{
-  Volume { Volumes() }; Delete; }{};
-
-vv() = VolumesRotor(); Call CheckNumbering;
-
-// S T A T O R
-
-// Air ring 
-// This region is specific to the sliding surface mechanism
-// to represent rotor motion.
-// It is a thin volumic region of width 'e_ss' in contact 
-// with the 'SlidingMaster' surface. 
-
-p1ss = newp; Point(p1ss) = {R_ar,0,SlidingSurfaceShift};
-p2ss = newp; Point(p2ss) = {R_ar,0,SlidingSurfaceShift+A_infinity/2};
-p3ss = newp; Point(p3ss) = {R_ss,0,A_infinity/2};
-p4ss = newp; Point(p4ss) = {R_ss,0,0};
-  
-l1ss = newl; Line(l1ss) = {p1ss,p2ss};
-l2ss = newl; Line(l2ss) = {p2ss,p3ss};
-l3ss = newl; Line(l3ss) = {p3ss,p4ss};
-l4ss = newl; Line(l4ss) = {p4ss,p1ss};
-
-ll_ss = newll; Line Loop(ll_ss) = {l1ss, l2ss, l3ss, l4ss};
-sAirRing = news; Plane Surface(sAirRing) = {ll_ss};
-vAirRingExt() = Extrude {{0, 0, 1}, {0, 0, 0}, anglePole} {
-  Surface{sAirRing};
-};
-
-vAirRing = vAirRingExt[1] ;
-
-// Air stator
-
-p1as = newp; Point(p1as) = {R_ss,0,0};
-p2as = newp; Point(p2as) = {R_ss,0,A_infinity/2};
-p3as = newp; Point(p3as) = {R_infinity,0,A_infinity/2};
-p4as = newp; Point(p4as) = {R_infinity,0,0};
-  
-l1as = newl; Line(l1as) = {p1as,p2as};
-l2as = newl; Line(l2as) = {p2as,p3as};
-l3as = newl; Line(l3as) = {p3as,p4as};
-l4as = newl; Line(l4as) = {p4as,p1as};
-
-ll_as = newll; Line Loop(ll_as) = {l1as, l2as, l3as, l4as};
-sAirStator = news; Plane Surface(sAirStator) = {ll_as};
-vAirStatorExt() = Extrude {{0, 0, 1}, {0, 0, 0}, anglePole} {
-  Surface{sAirStator};
-};
-
-// Stator core
-
-p1s = newp; Point(p1s) = {R_s,0,0};
-p2s = newp; Point(p2s) = {R_s,0,L_s/2};
-p3s = newp; Point(p3s) = {R_o,0,L_s/2};
-p4s = newp; Point(p4s) = {R_o,0,0};
-
-l1s = newl; Line(l1s) = {p1s,p2s};
-l2s = newl; Line(l2s) = {p2s,p3s};
-l3s = newl; Line(l3s) = {p3s,p4s};
-l4s = newl; Line(l4s) = {p4s,p1s};
-
-ll_s = newll; Line Loop(ll_s) = {l1s, l2s, l3s, l4s};
-sStator = news; Plane Surface(sStator) = {ll_s};
-vStatorExt() = Extrude {{0, 0, 1}, {0, 0, 0}, anglePole} {
-  Surface{sStator}; 
-};
-
-// Create volumes by boolean operations
-
-vStator = vStatorExt[1];
-vAirStator = BooleanDifference{
-  Volume{vAirStatorExt[1]}; Delete; 
-}{ 
-  Volume{vStator}; 
-};
-
-// Glue up stator pieces
-Volumes() = { vAirRing, vStator, vAirStator };
-VolumesStator() = BooleanFragments{
-  Volume { Volumes() }; Delete; }{};
-
-Printf("VolumesStator=", VolumesStator());
-Printf("VolumesRotor=", VolumesRotor());
-
-// Recover Physical Surfaces geometrically
-
 CoordinateSystem = 1; // cylindrical coordinates
 
+// get useful rotor surfaces
 Volumes() = VolumesRotor(); Call Volume2Cube;
-Printf("cube rotor", cube());
-FaceId = 1; Call Cube2Face; // no statement here FIXME 
+Printf("cube rotor: ", cube());
+FaceId = 1; Call Cube2Face; // FIXME: cannot have statement here
 RotorPerMaster() = face();
-FaceId = 2; Call Cube2Face; 
+FaceId = 2; Call Cube2Face;
 RotorBottom() = face();
-FaceId = 3; Call Cube2Face; 
+FaceId = 3; Call Cube2Face;
 SlidingSlave() = face();
-FaceId = 4; Call Cube2Face; 
+FaceId = 4; Call Cube2Face;
 RotorPerSlave() = face();
-FaceId = 5; Call Cube2Face; 
+FaceId = 5; Call Cube2Face;
 RotorTop() = face();
 
-/* control if the geometrical identification succeeded 
-Physical Surface("Test1", 101) = { RotorPerMaster() };
-Physical Surface("Test2", 102) = { RotorTop() };
-Physical Surface("Test3", 103) = { SlidingSlave() };
-Physical Surface("Test4", 104) = { RotorBottom() };
-Physical Surface("Test5", 105) = { RotorPerSlave() };
-*/
-
+// get useful stator surfaces
 Volumes() = VolumesStator(); Call Volume2Cube;
-Printf("cube stator", cube());
-FaceId = 0; Call Cube2Face; 
+Printf("cube stator: ", cube());
+FaceId = 0; Call Cube2Face;
 SlidingMaster() = face();
-FaceId = 1; Call Cube2Face; 
+FaceId = 1; Call Cube2Face;
 StatorPerMaster() = face();
-FaceId = 2; Call Cube2Face; 
+FaceId = 2; Call Cube2Face;
 StatorBottom() = face();
-FaceId = 3; Call Cube2Face; 
+FaceId = 3; Call Cube2Face;
 StatorOuter() = face();
-FaceId = 4; Call Cube2Face; 
+FaceId = 4; Call Cube2Face;
 StatorPerSlave() = face();
-FaceId = 5; Call Cube2Face; 
+FaceId = 5; Call Cube2Face;
 StatorTop() = face();
 
-/* control if the geometrical identification succeeded 
-Physical Surface("Test11", 111) = { SlidingMaster() };
-Physical Surface("Test12", 112) = { StatorPerMaster() };
-Physical Surface("Test13", 113) = { StatorBottom() };
-Physical Surface("Test14", 114) = { StatorOuter() };
-Physical Surface("Test15", 115) = { StatorPerSlave() };
-Physical Surface("Test16", 116) = { StatorTop() };
-*/
+// Make AirRing a transfinite region.  The width of the air ring is arbitrary,
+// but it is in general thin, especially in machine models where it has to fit
+// into the air gap.  To improve the quality of the mesh, it is thus appropriate
+// to have it transfinite in practice.
 
-
-// Make AirRing a transfinite region
-// The width of the air ring is arbitrary, but it is in general thin,
-// especially in machine models where it has to fit into the air gap.
-// To improve the quality of the mesh, 
-// it is thus appropriate to have it transfinite in practice. 
-
-surfaces() = Boundary { Volume{ vAirRing }; };
-//Printf("%g %g", vAirRing, #surfaces());
+surfaces() = Boundary { Volume{10}; };
 lines() = Boundary { Surface{ SlidingSlave() }; };
 For num In { 0:#surfaces()-1 }
-  lines() += Boundary { Surface{ surfaces[num] }; };
-  Printf("%g %g", surfaces[num], #lines());
+  lines() += Boundary { Surface{ surfaces(num) }; };
 EndFor
 For num In { 0:#lines()-1 }
   points() = PointsOf { Line { lines(num) }; };
-  If( #points() != 2 )
+  If(#points() != 2)
     Error("This is unexpected"); Abort;
   EndIf
   ptA() = Point { points(0) };
   ptB() = Point { points(1) };
-  // Printf("%f %f", Atan2[ ptA[1], ptA[0]]/deg, Atan2[ ptB[1], ptB[0] ]/deg );
-  If( Fabs( Atan2( ptA(1), ptA(0) ] - Atan2[ ptB(1), ptB(0) ) ) > tol )
+  If(Fabs(Atan2(ptA(1), ptA(0)) - Atan2(ptB(1), ptB(0))) > tol)
     // line along theta-axis
-    Transfinite Line{ lines[num] } = NbrDiv+1;
-  ElseIf( Fabs[ ptA(2) - ptB(2) ] > tol ) 
+    Transfinite Line{ lines(num) } = NbrDiv+1;
+  ElseIf(Fabs(ptA(2) - ptB(2)) > tol)
     // line along z-axis
-    Transfinite Line{ lines(num) } = Fabs( ptA(2) - ptB(2) )/lc-1;
-  Else 
+    Transfinite Line{ lines(num) } = Fabs(ptA(2) - ptB(2))/lc-1;
+  Else
     // line along r-axis
     Transfinite Line{ lines(num) } = 2;
   EndIf
 EndFor
 
-Transfinite Surface{ surfaces() } ; 
-Transfinite Surface{ SlidingSlave() } Right;
-Transfinite Volume{ vAirRing };
+Transfinite Surface{ surfaces() } ;
+Transfinite Surface{ SlidingSlave() };
+Transfinite Volume{ 10 };
 
 // Identify 'SlidingSubmaster' Curve
 
-AngleSubSurfaces = (ModelAngleMax-ModelAngleMin)*deg;
-
 lines() = Boundary { Surface{ SlidingMaster() }; };
 For l In { 0:#lines()-1 }
   points() = PointsOf { Line { lines(l) }; };
-  If( #points() != 2 )
+  If(#points() != 2)
     Error("This is unexpected"); Abort;
   EndIf
   ptA() = Point { points(0) };
   ptB() = Point { points(1) };
-  // Printf("%f %f", Atan2[ ptA[1], ptA[0]]/deg, Atan2[ ptB[1], ptB[0] ]/deg );
-  If( Fabs( Atan2(ptA(1), ptA(0)) - AngleSubSurfaces) < tol &&
- 	  Fabs( Atan2(ptB(1), ptB(0)) - AngleSubSurfaces) < tol )
+  If(Fabs(Atan2(ptA(1), ptA(0)) - AnglePole) < tol &&
+     Fabs(Atan2(ptB(1), ptB(0)) - AnglePole) < tol)
     SlidingSubmaster =  lines(l);
   EndIf
 EndFor
@@ -417,30 +234,20 @@ EndFor
 lines() = Boundary { Surface{ SlidingSlave() }; };
 For l In { 0:#lines()-1 }
   points() = PointsOf { Line { lines(l) }; };
-  If( #points() != 2 )
+  If(#points() != 2)
   Error("This is unexpected"); Abort;
   EndIf
   ptA() = Point { points(0) };
   ptB() = Point { points(1) };
-  If( Fabs( Atan2(ptA(1), ptA(0)) - AngleSubSurfaces ) < tol &&
-	  Fabs( Atan2(ptB(1), ptB(0)) - AngleSubSurfaces ) < tol )
+  If(Fabs(Atan2(ptA(1), ptA(0)) - AnglePole) < tol &&
+     Fabs(Atan2(ptB(1), ptB(0)) - AnglePole) < tol)
     SlidingSubslave =  lines(l);
   EndIf
 EndFor
 
-
-// P H Y S I C A L   R E G I O N S
-
-Physical Volume("Rotor"    , 1) = {vRotor};
-Physical Volume("Magnet"   , 2) = {vMagnet};
-Physical Volume("Stator"   , 3) = {vStator};
-Physical Volume("AirRotor" , 4) = {vAirRotor};
-Physical Volume("AirStator", 5) = {vAirStator};
-Physical Volume("AirRing"  , 6) = {vAirRing};
-
 Physical Surface("Outer"          , 10) = {StatorOuter()};
-Physical Surface("Top"            , 11) = {RotorTop(),StatorTop()};
-Physical Surface("Bottom"         , 12) = {RotorBottom(),StatorBottom()};
+Physical Surface("Top"            , 11) = {RotorTop(), StatorTop()};
+Physical Surface("Bottom"         , 12) = {RotorBottom(), StatorBottom()};
 Physical Surface("RotorPerMaster" , 13) = {RotorPerMaster()};
 Physical Surface("RotorPerSlave"  , 14) = {RotorPerSlave()};
 Physical Surface("StatorPerMaster", 15) = {StatorPerMaster()};
@@ -448,24 +255,19 @@ Physical Surface("StatorPerSlave" , 16) = {StatorPerSlave()};
 Physical Surface("SlidingMaster"  , 17) = {SlidingMaster()};
 Physical Surface("SlidingSlave"   , 18) = {SlidingSlave()};
 
-Physical Line("Sliding_Submaster", 20) = {SlidingSubmaster};
-Physical Line("Sliding_Subslave",  21) = {SlidingSubslave};
-
+Physical Line("Sliding_Submaster" , 20) = {SlidingSubmaster};
+Physical Line("Sliding_Subslave"  , 21) = {SlidingSubslave};
 
 // Ensure identical meshes on the periodic faces
-// The 'Air ring' being meshed transfinite with prisms,
-// the periodic surface routine issues a harmless warning  
-// Info    : No periodic vertices in surface xxx 
 For num In {0:#RotorPerMaster()-1}
-  Periodic Surface { RotorPerSlave(num) } 
-                 = { RotorPerMaster(num) } Rotate { {0,0,1}, {0,0,0}, anglePole };
+  Periodic Surface { RotorPerSlave(num) } =
+    { RotorPerMaster(num) } Rotate { {0,0,1}, {0,0,0}, AnglePole };
 EndFor
 For num In {0:#StatorPerMaster()-1}
-  Periodic Surface { StatorPerSlave(num) } 
-                 = { StatorPerMaster(num) } Rotate { {0,0,1}, {0,0,0}, anglePole };
+  Periodic Surface { StatorPerSlave(num) } =
+    { StatorPerMaster(num) } Rotate { {0,0,1}, {0,0,0}, AnglePole };
 EndFor
 
-
 treeLines()  = CombinedBoundary { Physical Surface { 10 }; } ;
 treeLines() += CombinedBoundary { Physical Surface { 11 }; } ;
 treeLines() += CombinedBoundary { Physical Surface { 12 }; } ;
@@ -477,4 +279,3 @@ treeLines() += CombinedBoundary { Physical Surface { 17 }; } ;
 treeLines() += CombinedBoundary { Physical Surface { 18 }; } ;
 
 Physical Line("treeLines", 22) = { treeLines() };
-

SlidingSurface3D/rfpm.pro

@@ -18,7 +18,7 @@ DefineConstant[
 
 Flag_QuasiStatic = !StrCmp[ R_,"QuasiStatic"] ;
 RotorPosition = deg*
-  DefineNumber[0, Name "Options/Rotor position [deg]",
+  DefineNumber[20, Name "Options/Rotor position [deg]",
                Visible !Flag_QuasiStatic];
 AngularStep = deg*
   DefineNumber[1, Name "Options/1Angular step [deg]",

SlidingSurface3D/rfpm_common.pro

 mm = 1.e-3 ;
 deg = Pi/180.;
 
-lc = mm * DefineNumber[2, Name "Options/Global mesh size (mm)"];
-If(lc == 0)
-  Error("Glabal mesh size cannot be zero. Reset to default."); lc=2*mm;
+lc = mm * DefineNumber[10, Name "Options/Global mesh size (mm)"];
+If(lc <= 0)
+  Error("Glabal mesh size cannot be zero. Reset to default.");
+  lc = 10 * mm;
 EndIf
 
 // Number of elements on the sliding surface in rotation direction
-NbrDiv = 30*mm/lc; 
+NbrDiv = 50*mm/lc;
 
-Gauge = DefineNumber[1, Name "Options/Gauging",
-					 Choices{ 0="MUMPS", 1="Tree"}];
+Gauge = DefineNumber[1, Name "Options/Gauging", Choices{ 0="MUMPS", 1="Tree"}];
 
-ModelAngleMin = 0. ; 
+ModelAngleMin = 0. ;
 ModelAngleMax = 60. ;
-SlidingSurfaceShift = 1.e-5;
 
+//SlidingSurfaceShift = 1.e-5;
+SlidingSurfaceShift = 0;