diff --git a/SlidingSurface3D/rfpm.geo b/SlidingSurface3D/rfpm.geo
new file mode 100644
index 0000000000000000000000000000000000000000..aa7b9d1cddba8461c0ffbba72c89c69474d3f356
--- /dev/null
+++ b/SlidingSurface3D/rfpm.geo
@@ -0,0 +1,480 @@
+Include "rfpm_common.pro";
+
+
+// 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
+
+Function ChangeCoordinates
+ 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
+ //
+ points() = PointsOf{ Volume{ Volumes() }; };
+ 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) };
+ 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) );
+ EndFor
+Return
+
+Function Cube2Face
+ bbox() = cube();
+ bbox(Modulo(FaceId+3,6)) = bbox(FaceId);
+ //Printf("bbox=", bbox());
+ face()={};
+ surfaces() = Boundary{ Volume{ Volumes() }; };
+ For s In { 0:#surfaces()-1 }
+ NbIn = 0;
+ points() = PointsOf { Surface { surfaces(s) } ; };
+ For p In { 0:#points()-1 }
+ point() = Point { points(p) };
+ Call ChangeCoordinates;
+ 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++;
+ EndIf
+ EndFor
+ 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
+
+Volumes() = VolumesRotor(); Call Volume2Cube;
+Printf("cube rotor", cube());
+FaceId = 1; Call Cube2Face; // no statement here FIXME
+RotorPerMaster() = face();
+FaceId = 2; Call Cube2Face;
+RotorBottom() = face();
+FaceId = 3; Call Cube2Face;
+SlidingSlave() = face();
+FaceId = 4; Call Cube2Face;
+RotorPerSlave() = face();
+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() };
+*/
+
+Volumes() = VolumesStator(); Call Volume2Cube;
+Printf("cube stator", cube());
+FaceId = 0; Call Cube2Face;
+SlidingMaster() = face();
+FaceId = 1; Call Cube2Face;
+StatorPerMaster() = face();
+FaceId = 2; Call Cube2Face;
+StatorBottom() = face();
+FaceId = 3; Call Cube2Face;
+StatorOuter() = face();
+FaceId = 4; Call Cube2Face;
+StatorPerSlave() = face();
+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.
+
+surfaces() = Boundary { Volume{ vAirRing }; };
+//Printf("%g %g", vAirRing, #surfaces());
+lines() = Boundary { Surface{ SlidingSlave() }; };
+For num In { 0:#surfaces()-1 }
+ lines() += Boundary { Surface{ surfaces[num] }; };
+ Printf("%g %g", surfaces[num], #lines());
+EndFor
+For num In { 0:#lines()-1 }
+ points() = PointsOf { Line { lines(num) }; };
+ 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 )
+ // line along theta-axis
+ 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
+ // line along r-axis
+ Transfinite Line{ lines(num) } = 2;
+ EndIf
+EndFor
+
+Transfinite Surface{ surfaces() } ;
+Transfinite Surface{ SlidingSlave() } Right;
+Transfinite Volume{ vAirRing };
+
+// 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 )
+ 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 )
+ SlidingSubmaster = lines(l);
+ EndIf
+EndFor
+
+// Identify 'SlidingSubslave' curve
+
+lines() = Boundary { Surface{ SlidingSlave() }; };
+For l In { 0:#lines()-1 }
+ points() = PointsOf { Line { lines(l) }; };
+ 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 )
+ 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("RotorPerMaster" , 13) = {RotorPerMaster()};
+Physical Surface("RotorPerSlave" , 14) = {RotorPerSlave()};
+Physical Surface("StatorPerMaster", 15) = {StatorPerMaster()};
+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};
+
+
+// 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 };
+EndFor
+For num In {0:#StatorPerMaster()-1}
+ 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 }; } ;
+treeLines() += CombinedBoundary { Physical Surface { 13 }; } ;
+treeLines() += CombinedBoundary { Physical Surface { 14 }; } ;
+treeLines() += CombinedBoundary { Physical Surface { 15 }; } ;
+treeLines() += CombinedBoundary { Physical Surface { 16 }; } ;
+treeLines() += CombinedBoundary { Physical Surface { 17 }; } ;
+treeLines() += CombinedBoundary { Physical Surface { 18 }; } ;
+
+Physical Line("treeLines", 22) = { treeLines() };
+
diff --git a/SlidingSurface3D/rfpm.pro b/SlidingSurface3D/rfpm.pro
new file mode 100644
index 0000000000000000000000000000000000000000..25c8379a03376d586f39158fe49d5893a9c29469
--- /dev/null
+++ b/SlidingSurface3D/rfpm.pro
@@ -0,0 +1,368 @@
+/*
+ 3D Magnetostatics with
+ - OpenCascade solid modelling
+ - geometrical identfication of faces
+ - antiperiodic boundary conditions
+ - sliding surface
+*/
+
+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}
+ P_ = { "Check_Periodicity", Choices {"Fields", "Check_Periodicity"}, Visible 0,
+ Name "GetDP/2PostOperationChoices"}
+];
+
+Flag_QuasiStatic = !StrCmp[ R_,"QuasiStatic"] ;
+RotorPosition = deg*
+ DefineNumber[0, 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];
+
+rpm = 2.*Pi/60. ;
+w1 = 1000*rpm ;
+Mag_DTime = AngularStep/Fabs[w1] ;
+Mag_TimeMax = Mag_DTime*NumStep;
+
+
+Group {
+ // Geometrical regions
+ Rotor = Region[1];
+ Magnet = Region[2];
+ Stator = Region[3];
+ AirRotor = Region[4];
+ AirStator = Region[5];
+ AirRing = Region[6];
+
+ Outer = Region[10];
+ Top = Region[11];
+ Bottom = Region[12];
+ Sur_RotorPerMaster = Region[13];
+ Sur_RotorPerSlave = Region[14];
+ Sur_StatorPerMaster = Region[15];
+ Sur_StatorPerSlave = Region[16];
+ Sur_SlidingMaster = Region[17];
+ Sur_SlidingSlave = Region[18];
+
+ Lin_SlidingSubmaster = Region[20];
+ Lin_SlidingSubslave = Region[21];
+ Lin_TreeLines = Region[22];
+
+ // Abstract regions
+ Iron = Region[ { Rotor, Stator } ];
+ Air = Region[ { AirRotor, AirStator, AirRing } ];
+
+ Vol_Conductors = Region[ {} ] ;
+ Vol_Magnets = Region[ Magnet ];
+ Vol_Mag = Region[ { Air, Magnet, Iron } ];
+ Vol_Moving = Region[ { Rotor, AirRotor, Magnet } ] ;
+
+ Sur_Dirichlet = Region[ { Outer, Top, Bottom } ];
+
+ Sur_Link = Region[ {Sur_SlidingMaster, Sur_SlidingSlave,
+ Sur_StatorPerMaster, Sur_StatorPerSlave,
+ 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 } ];
+ Vol_Tree = Region[ Vol_Mag ];
+}
+
+Function{
+ mu0 = Pi*4e-7;
+ br = 1.2;
+ mur = 3;
+
+ nu[ Iron ] = 1/(mur*mu0);
+ nu[ Air ] = 1/mu0;
+ nu[ Magnet ] = 1/mu0;
+ br[] = br*Unit[ Vector[X[], Y[], 0] ];
+
+ js[] = 0;
+}
+
+
+Function {
+ // Sliding surface:
+ // a1 is the angular span of Region and RegionRef.
+ // a0 is the mesh step of the meshes of
+ // a2[] is the angular position (in radian) of the rotor
+ // relative to its reference position (as in the msh file)
+ // assumed to be aligned with the stator.
+ // a3[] is the shear angle of region AIRBM
+ // (smaller than half the discretization step of the sliding surface
+ // to adapt to a2[] values that are not multiple of this step).
+ // AlignWithMaster[] maps a point of coordinates {X[], Y[], Z[]} onto
+ // its image in the open set RegionRef-SubRegionRef by the symmetry mapping.
+ // Coef[] is evaluated on Master nodes
+ // fFloor[] is a safe version of Floor[] for real represented int variables
+ // fRem[a,b] substracts from a multiple of b so that the result is in [0,1[
+
+ Periodicity = -1. ; // -1 for antiperiodicity, 1 for periodicity
+ RotatePZ[] = Rotate[ XYZ[], 0, 0, $1 ] ;
+ Tol = 1e-8 ; fFloor[] = Floor[ $1 + Tol ] ; fRem[] = $1 - fFloor[ $1 / $2 ] * $2;
+ deg = Pi/180;
+ a1 = (ModelAngleMax-ModelAngleMin)*deg ;
+ a0 = a1/NbrDiv ; // angular span of one moving band element
+ a2[] = $RotorPosition ;
+ a3[] = ( ( fRem[ a2[], a0 ]#2 <= 0.5*a0 ) ? #2 : #2-a0 ) ;
+ AlignWithMaster[] = RotatePZ[ - fFloor[ (Atan2[ Y[], X[] ] - a3[] )/a1 ] * a1 ] ;
+ RestoreRef[] = XYZ[] - Vector[ 0, 0, SlidingSurfaceShift ] ;
+ 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} ; }
+ { Region All ; Jacobian Vol ; }
+ }
+ }
+ { Name Sur ;
+ Case { { Region All ; Jacobian Sur ; }
+ }
+ }
+}
+
+Integration {
+ { Name Int ;
+ Case {
+ { Type Gauss ;
+ Case {
+ { GeoElement Line ; NumberOfPoints 4 ; }
+ { GeoElement Triangle ; NumberOfPoints 4 ; }
+ { GeoElement Quadrangle ; NumberOfPoints 4 ; }
+ { GeoElement Tetrahedron ; NumberOfPoints 4 ; }
+ { GeoElement Hexahedron ; NumberOfPoints 6 ; }
+ { GeoElement Prism ; NumberOfPoints 21 ; }
+ }
+ }
+ }
+ }
+}
+
+Constraint {
+ { Name a ;
+ Case {
+ { Region Sur_Dirichlet ; Value 0. ; }
+
+ // Periodicity condition on lateral faces
+ { Type Link ; Region Sur_StatorPerSlave ; RegionRef Sur_StatorPerMaster ;
+ ToleranceFactor 1e-8;
+ Coefficient Periodicity ; Function RotatePZ[ -a1 ] ; }
+ { Type Link ; Region Sur_RotorPerSlave ; RegionRef Sur_RotorPerMaster ;
+ ToleranceFactor 1e-8;
+ Coefficient Periodicity ; Function RotatePZ[ -a1 ] ; }
+
+ // Sliding surface
+ { Type Link ; Region Sur_SlidingSlave ; SubRegion Lin_SlidingSubslave ;
+ RegionRef Sur_SlidingMaster ; SubRegionRef Lin_SlidingSubmaster ;
+ ToleranceFactor 1e-8;
+ Coefficient Coef[] ; Function AlignWithMaster[] ;
+ FunctionRef RestoreRef[] ;
+ }
+ }
+ }
+ { Name GaugeCondition_a ; Type Assign ;
+ Case {
+ { Region Vol_Tree ; Value 0. ;
+ SubRegion Region[ { Sur_Tree, Lin_Tree} ] ;
+ SubRegion2 Region[ Sur_Tree_Sliding, AlignedWith Z ] ;
+ }
+ }
+ }
+}
+
+FunctionSpace {
+ { Name Hcurl_a; Type Form1;
+ BasisFunction {
+ { Name se; NameOfCoef ae; Function BF_Edge;
+ 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
+ }
+ }
+}
+
+Formulation {
+ { Name MagSta_a; Type FemEquation ;
+ Quantity {
+ { Name a ; Type Local ; NameOfSpace Hcurl_a ; }
+ }
+ Equation {
+ Galerkin { [ nu[] * Dof{d a} , {d a} ] ;
+ In Vol_Mag ; Jacobian Vol ; Integration Int ; }
+ Galerkin { [ -nu[] * br[] , {d a} ] ;
+ In Vol_Magnets ; Jacobian Vol ; Integration Int ; }
+ Galerkin { [ -js[] , {a} ] ; In Vol_Conductors ;
+ Jacobian Vol ; Integration Int ; }
+ }
+ }
+}
+
+
+Resolution {
+ { Name Static ;
+ System {
+ { Name Sys_Mag ; NameOfFormulation MagSta_a ;}
+ }
+ Operation {
+ Evaluate[$RotorPosition = RotorPosition];
+ ChangeOfCoordinates [ NodesOf[ Vol_Moving ], RotatePZ[ a2[] ] ] ;
+ ChangeOfCoordinates [ NodesOf[ Sur_SlidingMaster ], RotatePZ[ a3[] ] ] ;
+ Evaluate[ $a2 = a2[] ];
+ Evaluate[ $a3 = a3[] ];
+ Evaluate[ $aa = Fmod[ a2[], a0 ] ];
+ Evaluate[ $bb = fRem[ a2[], a0 ] ];
+ Print[ {$RotorPosition/deg, $a2/deg, $a3/deg, $aa/deg, $bb/deg},
+ Format "wt=%e a2=%e a3=%e %e %e"] ;
+
+ UpdateConstraint[Sys_Mag] ;
+ Generate Sys_Mag ; Solve Sys_Mag ; SaveSolution Sys_Mag ;
+ PostOperation [Fields];
+ }
+ }
+
+ { Name QuasiStatic ;
+ System {
+ { Name Sys_Mag ; NameOfFormulation MagSta_a ; }
+ }
+ Operation {
+ Evaluate[$RotorPosition = RotorPosition];
+ Generate Sys_Mag ; Solve Sys_Mag ; SaveSolution Sys_Mag ;
+ PostOperation [Fields];
+ TimeLoopTheta{
+ Time0 0 ; TimeMax Mag_TimeMax ; DTime Mag_DTime ; Theta 1. ;
+ Operation {
+ Evaluate[$RotorPosition = RotorPosition + w1*$Time];
+ ChangeOfCoordinates[NodesOf[ Vol_Moving ], RotatePZ[ w1*$DTime ] ] ;
+ ChangeOfCoordinates[NodesOf[ Sur_SlidingMaster ], RotatePZ[ a3[] ] ] ;
+ UpdateConstraint[Sys_Mag] ;
+ Generate Sys_Mag ; Solve Sys_Mag ; SaveSolution Sys_Mag ;
+ PostOperation [Fields];
+ ChangeOfCoordinates [ NodesOf[ Sur_SlidingMaster ], RotatePZ[ -a3[] ] ] ;
+ }
+ }
+ }
+ }
+}
+
+PostProcessing {
+ { Name MagSta_a ; NameOfFormulation MagSta_a ;
+ PostQuantity {
+ { Name b ; Value { Local { [ {d a} ];
+ In Vol_Mag ; Jacobian Vol; } } }
+ { Name bsurf ; Value { Local { [ {d a} ];
+ In Sur_Link ; Jacobian Sur; } } }
+ { Name asurf ; Value { Local { [ {a} ];
+ In Sur_Link ; Jacobian Sur; } } }
+ { Name br ; Value { Local { [ br[] ];
+ In Vol_Magnets ; Jacobian Vol; } } }
+ { Name a ; Value { Local { [ {a} ];
+ In Vol_Mag ; Jacobian Vol; } } }
+ }
+ }
+}
+
+
+PostOperation Fields UsingPost MagSta_a {
+ Print[ b, OnElementsOf Region[ {Vol_Mag} ],
+ LastTimeStepOnly, File "b.pos"] ;
+ Echo[ Str["l=PostProcessing.NbViews-1;",
+ "View[l].ArrowSizeMax = 100;",
+ "View[l].CenterGlyphs = 1;",
+ "View[l].GlyphLocation = 1;",
+ "View[l].RangeType = 1;",
+ "View[l].SaturateValues = 1;",
+ "View[l].CustomMax = 2;",
+ "View[l].CustomMin = 1e-06;",
+ "View[l].LineWidth = 2;",
+ "View[l].VectorType = 4;" ] ,
+ File "tmp.geo", LastTimeStepOnly] ;
+}
+
+
+
+PostOperation {
+ { Name pos; NameOfPostProcessing MagSta_a;
+ Operation {
+ Print[ b, OnElementsOf Dom_Hcurl_a, File "b.pos" ];
+
+ // Print[W_mag_Airgap[AIRGAP], OnGlobal, Format TimeTable,
+ // File StrCat[resPath, "Energy_airgap.pos"]];
+ // Print[Flux[PM], OnGlobal, Format TimeTable,
+ // File >> StrCat[resPath, "Flux.pos"]];
+ }
+ }
+}
+
+PostOperation Check_Periodicity UsingPost MagSta_a {
+ PrintGroup[ EdgesOfTreeIn[ { Vol_Tree }, StartingOn { Sur_Tree } ],
+ In Vol_Tree, File "Tree.pos"];
+ Echo[ Str["l=PostProcessing.NbViews-1;",
+ "View[l].ColorTable = { DarkRed };",
+ "View[l].LineWidth = 5;"] ,
+ File "tmp.geo", LastTimeStepOnly] ;
+
+ Print[ bsurf, OnElementsOf Region[ {Sur_SlidingMaster} ],
+ LastTimeStepOnly, File "bsm.pos"] ;
+ Echo[ Str["l=PostProcessing.NbViews-1;",
+ "View[l].ArrowSizeMax = 100;",
+ "View[l].CenterGlyphs = 0;",
+ "View[l].GlyphLocation = 1;",
+ "View[l].RangeType = 1;",
+ "View[l].SaturateValues = 1;",
+ "View[l].CustomMax = 2;",
+ "View[l].CustomMin = 0;",
+ "View[l].ScaleType = 1;",
+ "View[l].LineWidth = 3;",
+ "View[l].VectorType = 4;" ] ,
+ File "tmp.geo", LastTimeStepOnly] ;
+
+ Print[ bsurf, OnElementsOf Region[ {Sur_SlidingSlave} ],
+ LastTimeStepOnly, File "bss.pos"] ;
+ Echo[ Str["l=PostProcessing.NbViews-1;",
+ "View[l].ArrowSizeMax = 100;",
+ "View[l].CenterGlyphs = 0;",
+ "View[l].GlyphLocation = 1;",
+ "View[l].RangeType = 1;",
+ "View[l].SaturateValues = 1;",
+ "View[l].CustomMax = 2;",
+ "View[l].CustomMin = 0;",
+ "View[l].ScaleType = 1;",
+ "View[l].LineWidth = 3;",
+ "View[l].VectorType = 4;" ] ,
+ File "tmp.geo", LastTimeStepOnly] ;
+
+ Print[ bsurf, OnElementsOf Region[ Sur_Link ],
+ LastTimeStepOnly, File "bsl.pos"] ;
+ Echo[ Str["l=PostProcessing.NbViews-1;",
+ "View[l].ArrowSizeMax = 100;",
+ "View[l].CenterGlyphs = 0;",
+ "View[l].GlyphLocation = 1;",
+ "View[l].ScaleType = 1;",
+ "View[l].LineWidth = 3;",
+ "View[l].VectorType = 4;" ] ,
+ File "tmp.geo", LastTimeStepOnly] ;
+}
diff --git a/SlidingSurface3D/rfpm_common.pro b/SlidingSurface3D/rfpm_common.pro
new file mode 100644
index 0000000000000000000000000000000000000000..c4054fda5c7d0dfdc10ae0ea48b50a248bcf13dc
--- /dev/null
+++ b/SlidingSurface3D/rfpm_common.pro
@@ -0,0 +1,18 @@
+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;
+EndIf
+
+// Number of elements on the sliding surface in rotation direction
+NbrDiv = 30*mm/lc;
+
+Gauge = DefineNumber[1, Name "Options/Gauging",
+ Choices{ 0="MUMPS", 1="Tree"}];
+
+ModelAngleMin = 0. ;
+ModelAngleMax = 60. ;
+SlidingSurfaceShift = 1.e-5;
+