From 0fe2c76cab1287716e646e8c25d41f5f53416159 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Thu, 6 May 2004 15:38:53 +0000
Subject: [PATCH] added a couple of nice examples

---
 benchmarks/2d/machine3/Machine.geo |  37 ++++
 benchmarks/2d/machine3/Rotor.geo   | 319 +++++++++++++++++++++++++++++
 benchmarks/2d/machine3/Stator.geo  | 303 +++++++++++++++++++++++++++
 benchmarks/3d/periodic.geo         |  61 ++++++
 4 files changed, 720 insertions(+)
 create mode 100644 benchmarks/2d/machine3/Machine.geo
 create mode 100644 benchmarks/2d/machine3/Rotor.geo
 create mode 100644 benchmarks/2d/machine3/Stator.geo
 create mode 100644 benchmarks/3d/periodic.geo

diff --git a/benchmarks/2d/machine3/Machine.geo b/benchmarks/2d/machine3/Machine.geo
new file mode 100644
index 0000000000..a2bf788a5e
--- /dev/null
+++ b/benchmarks/2d/machine3/Machine.geo
@@ -0,0 +1,37 @@
+
+Include "Rotor.geo" ;
+Include "Stator.geo" ;
+
+MB = 9999;
+If (MB)
+  dH=newreg; 
+  Line Loop(dH) = {OuterMB_[]};
+  Line Loop(dH+1) = {InnerMB_[]};
+  Plane Surface(dH) = {dH,dH+1};
+  Physical Surface(MB) = {dH};
+  Color Black {Surface{dH};}
+EndIf
+
+
+
+
+Mesh.CharacteristicLengthFactor = 1.2;
+
+Mesh.Smoothing = 5;
+//Mesh.Algorithm = 2;
+Mesh.ElementOrder = 2;
+
+Coherence;
+
+dP=newp-1;
+Point(dP+1)  = {0.0258, 0.0045, 0, p*2};  //stator
+//Point(dP+2)  = {0.0248, 0.0043, 0, p*2};  
+Point(dP+2)  = {0.0242, 0.007, 0, p*2};  
+
+
+
+
+
+
+
+
diff --git a/benchmarks/2d/machine3/Rotor.geo b/benchmarks/2d/machine3/Rotor.geo
new file mode 100644
index 0000000000..dad85f9d38
--- /dev/null
+++ b/benchmarks/2d/machine3/Rotor.geo
@@ -0,0 +1,319 @@
+
+Geometry.AutoCoherence = 0;      
+Mesh.CharacteristicLengthFactor = 1.;
+
+NbrPoles = 1; // number of rotor poles in FE model
+NbrPolesT = 4; // number of poles in complete cross-section 
+NbrSectT = 40; // number of "rotor teeth"
+NbrSect = NbrSectT*NbrPoles/NbrPolesT; // number of "rotor teeth" in FE model
+
+RotorAngle_R = Pi/NbrSectT -Pi/2; // initial rotor angle (radians)
+RotorAngle_S = RotorAngle_R;
+
+
+/* physical rotor numbers (for GetDP and Gmsh) */
+
+RotorBars= 17900; gg=RotorBars;
+  RotorBar01 = gg+01; RotorBar02 = gg+02; RotorBar03 = gg+03; RotorBar04 = gg+04; RotorBar05 = gg+05; RotorBar06 = gg+06;
+  RotorBar07 = gg+07; RotorBar08 = gg+08; RotorBar09 = gg+09; RotorBar10 = gg+10; RotorBar11 = gg+11; RotorBar12 = gg+12;  
+  RotorBar13 = gg+13; RotorBar14 = gg+14; RotorBar15 = gg+15; RotorBar16 = gg+16; RotorBar17 = gg+17; RotorBar18 = gg+18;
+  RotorBar19 = gg+19; RotorBar20 = gg+20; RotorBar21 = gg+21; RotorBar22 = gg+22; RotorBar23 = gg+23; RotorBar24 = gg+24;  
+  RotorBar25 = gg+25; RotorBar26 = gg+26; RotorBar27 = gg+27; RotorBar28 = gg+28; RotorBar29 = gg+29; RotorBar30 = gg+30; 
+  RotorBar31 = gg+31; RotorBar32 = gg+32; RotorBar33 = gg+33; RotorBar34 = gg+34; RotorBar35 = gg+35; RotorBar36 = gg+36;
+  RotorBar37 = gg+37; RotorBar38 = gg+38; RotorBar39 = gg+39; RotorBar40 = gg+40; RotorBar41 = gg+41; RotorBar42 = gg+42;
+
+
+RotorIron = 5000;
+RotorShaft = 5555;
+RotorAirgapLayer = 1999;
+RotorSlotOpening = 2999;
+
+RotorPeriod_Reference = 200000; 
+RotorPeriod_Dependent = 200001;
+RotorBoundary = 900;
+OuterShaft = 6666;
+InnerMB = 16000; MB_R1 = InnerMB+1; MB_R2 = InnerMB+2; MB_R3 = InnerMB+3; MB_R4 = InnerMB+4;
+
+
+u=1e-3; // scale factor
+
+AG = u* 0.47;       // airgap width
+R2 = u* 92/2 - AG;      // outer rotor radius
+//R2 = u* 92/2;      // outer rotor radius
+R3 = u* 31.75/2;      // shaft radius
+//R1 = u* 94/2;        // inner radius of moving band
+R1 = R2+AG/3;        // inner radius of moving band
+
+/*
+R1 = u* 88/2;       // inner radius of moving band
+R2 = u* 86/2;     // outer rotor radius
+R3 = u* 31.75/2;      // shaft radius
+*/
+
+
+// parameters for conductor and slot opening
+h1  = u* 1;
+h2  = u* 14.25;
+d1  = u* 2;
+Rsl = u* 4.26/2;
+      
+// characteristic lengths
+uc = u* 1.3;
+
+pslo = uc* 0.3; // slot opening
+psl  = uc* 0.6; // upper part slot 
+pslu = uc* 1; // lower part slot 
+psha = uc* 2; // shaft radius
+pMB  = AG; // MB
+//pMB  = uc* 0.5; // MB
+p  = uc* 2; // 
+
+Y1 = Sqrt(R2*R2-d1*d1/4) ;
+Y2 = Sqrt(Rsl*Rsl-d1*d1/4) ;
+Y3 = Sqrt(R1*R1-d1*d1/4) ;
+RX = Rsl*Cos(Pi/NbrSectT) ;
+RY = Rsl*Sin(Pi/NbrSectT) ;
+RR = (h2-Rsl*(1+1/Sin(Pi/NbrSectT)))/(1-1/Sin(Pi/NbrSectT));
+RX2 = RR*Cos(Pi/NbrSectT) ;
+RY2 = RR*Sin(Pi/NbrSectT) ;
+
+
+lc_1=0.0005;
+lc_2=0.0015;
+lc_3=0.003;
+lc_4=0.005;
+
+a_x=0.000784;
+a_y=0.04755;
+b_x=0.002319;
+b_y=0.067049;
+c_x=0.000700;
+c_y=0.070521;
+d_x=0.0007;
+d_y=0.075021;
+e_x=0.00145;
+e_y=0.075650;
+f_x=0.001450;
+f_y=0.077333;
+g_x=0;
+g_y=0.078550;
+h_x=(0.053/2)*Sin(4.5*(Pi/180.));
+h_y=(0.053/2)*Cos(4.5*(Pi/180.));
+hh_x=(0.053/2)*Sin(0*(Pi/180.));
+hh_y=(0.053/2)*Cos(0*(Pi/180.));
+R_i=0.080-0.00065;
+i_x=R_i*Sin(4.5*Pi/180);
+i_y=R_i*Cos(4.5*Pi/180);
+
+
+k_y=0.080-0.000650;
+l_y=0.080-(0.000650*(2/3));
+m_x=l_y*Sin(4.5*(Pi/180.));
+m_y=l_y*Cos(4.5*(Pi/180.));
+o_x=R_i*Sin(1.5*(Pi/180.));
+o_y=R_i*Cos(1.5*(Pi/180.));
+o_x=1.45e-3;
+o_y=Sqrt(R_i^2-o_x^2);
+
+oo_x=1.45e-3;
+oo_y=Sqrt(l_y^2-o_x^2);
+
+
+ii_x=l_y*Sin(4.5*Pi/180);
+ii_y=l_y*Cos(4.5*Pi/180);
+
+Point(0)  = {0,0,0,p};
+
+For i In {0:NbrSect-1}
+For j In {0:1}
+
+dP=newp-1;
+
+Point(dP+1) = {h_x,h_y,0,lc_4};
+Point(dP+2) = {hh_x,hh_y,0,lc_4};
+Point(dP+3) = {i_x,i_y,0, 0.7e-3};
+Point(dP+4) = {ii_x,ii_y,0, 0.7e-3};
+Point(dP+5) = {a_x,a_y,0, 2e-3};
+Point(dP+6) = {0,a_y,0,1e-3};
+Point(dP+7) = {b_x,b_y,0,lc_3*0.8};
+Point(dP+8) = {c_x,c_y,0,lc_2*0.8};
+Point(dP+9) = {0,c_y,0,lc_2};
+Point(dP+10) = {d_x,d_y,0,lc_2};
+Point(dP+11) = {0,d_y,0,lc_2};
+Point(dP+12) = {e_x,e_y,0,lc_2};
+Point(dP+13) = {f_x,f_y,0,lc_1};
+Point(dP+14) = {g_x,g_y,0,lc_1};
+Point(dP+15) = {0,R_i,0,0.3e-3};
+Point(dP+16) = {0,l_y,0,0.3e-3};
+Point(dP+17) = {o_x,o_y,0, 0.6e-3};
+Point(dP+18) = {oo_x,oo_y,0, 0.6e-3};
+
+For t In {dP+1:dP+18}
+  Rotate {{0,0,1},{0,0,0}, RotorAngle_R+2*Pi*i/NbrSectT} {Point{t};}
+EndFor
+
+If (j==1)
+  For t In {dP+1:dP+18}
+    Symmetry {Cos(RotorAngle_S+2*Pi*i/NbrSectT),Sin(RotorAngle_S+2*Pi*i/NbrSectT),0,0} { Point{t}; }  
+  EndFor
+EndIf
+
+dR=newreg-1;
+Line(dR+1) = {dP+16,dP+15};
+Line(dR+2) = {dP+15,dP+14};
+Line(dR+3) = {dP+14,dP+11};
+Line(dR+4) = {dP+14,dP+13};
+Line(dR+5) = {dP+13,dP+12};
+Line(dR+6) = {dP+12,dP+10};
+Line(dR+7) = {dP+4,dP+3};
+Line(dR+8) = {dP+11,dP+9};
+Line(dR+9) = {dP+10,dP+8};
+Line(dR+10) = {dP+8,dP+7};
+Line(dR+11) = {dP+7,dP+5};
+Line(dR+12) = {dP+5,dP+6};
+Line(dR+13) = {dP+6,dP+9};
+Line(dR+14) = {dP+6,dP+2};
+Line(dR+15) = {dP+3,dP+1};
+Circle(dR+16) = {dP+3,0,dP+17};
+Circle(dR+17) = {dP+17,0,dP+15};
+Circle(dR+18) = {dP+1,0,dP+2};
+Circle(dR+19) = {dP+4,0,dP+18};
+Circle(dR+20) = {dP+18,0,dP+16};
+Line(dR+21) = {dP+13,dP+17};
+
+
+// physical lines
+
+OuterShaft_[{2*i+j}] = dR+18;
+//RotorBoundary_[{14*i+7*j:14*i+7*j+6}] = {dR+12,dR+15,dR+13,dR+1,dR+16,dR+17,dR+14};
+RotorBoundary_[{20*i+10*j:20*i+10*j+9}] = {dR+4,dR+5,dR+6,dR+9,dR+10,dR+11,dR+12,dR+18,dR+16,dR+17};
+If (j==0)
+  InnerMB_[{4*i+2*j:4*i+2*j+1}] = {dR+19,dR+20};
+EndIf
+If (j==1)
+  InnerMB_[{4*i+2*j:4*i+2*j+1}] = {-dR-19,-dR-20};
+EndIf
+
+If (NbrSectT != NbrSect)
+  If (i==0 && j==0)
+    Physical Line(RotorPeriod_Reference) =  {dR+7,dR+15};
+    RotorPeriod_Reference_[] =  {dR+7,dR+15};
+  EndIf
+  If (i == NbrSect-1  && j==1)
+    Physical Line(RotorPeriod_Dependent) = {dR+7,dR+15};
+    RotorPeriod_Dependent_[] =  {dR+7,dR+15};
+  EndIf
+EndIf
+
+
+dH=newreg;  // rotor conductor
+Line Loop(dH) = {-dR-12,-dR-11,-dR-10,-dR-9,-dR-6,-dR-5,-dR-4,dR+3,dR+8,-dR-13};                    
+Plane Surface(dH) ={dH};
+RotorConductor_[2*i+j] = dH;
+
+
+dH=newreg;  // rotor iron
+Line Loop(dH) = {dR+16,-dR-21,dR+5,dR+6,dR+9,dR+10,dR+11,dR+12,dR+14,-dR-18,-dR-15};
+Plane Surface(dH) = {dH};
+RotorIron_[2*i+j] = dH;
+
+/*
+dH=newreg;   // rotor shaft
+Line Loop(dH) = {-dR-5,dR+2,dR+12};
+Plane Surface(dH) = {dH};
+RotorShaft_[2*i+j] = dH;
+*/
+
+dH=newreg;   // rotor slot opening  
+Line Loop(dH) = {dR+21,dR+17,dR+2,dR+4};
+Plane Surface(dH) ={dH};
+RotorSlotOpening_[2*i+j] = dH;
+
+dH=newreg;  // rotor airgap layer
+Line Loop(dH) = {-dR-1,-dR-20,-dR-19,dR+7,dR+16,dR+17};
+Plane Surface(dH) = {dH};
+RotorAirgapLayer_[2*i+j] = dH;
+
+
+EndFor
+EndFor
+
+
+
+For i In {0:NbrSect-1}
+   Physical Surface(RotorBars+1+i) = RotorConductor_[{2*i:2*i+1}];
+   Color Orange {Surface{RotorConductor_[{2*i}]};}
+   Color Orange {Surface{RotorConductor_[{2*i+1}]};}
+EndFor
+
+
+Physical Surface(RotorIron) = {RotorIron_[{0:NbrSect*2-1}]};
+Physical Surface(RotorSlotOpening) = {RotorSlotOpening_[{0:NbrSect*2-1}]};
+Physical Surface(RotorAirgapLayer) = {RotorAirgapLayer_[{0:NbrSect*2-1}]};
+
+Color Red {Surface{RotorIron_[{0:NbrSect*2-1}]};}
+Color Black {Surface{RotorSlotOpening_[{0:NbrSect*2-1}]};}
+Color Black {Surface{RotorAirgapLayer_[{0:NbrSect*2-1}]};}
+//Color Black {Surface{RotorShaft_[{0:NbrSect*2-1}]};}
+
+Physical Line(OuterShaft) = {OuterShaft_[]};
+
+If (NbrSectT != NbrSect)
+  Physical Line(RotorBoundary) = {RotorBoundary_[],RotorPeriod_Reference_[],RotorPeriod_Dependent_[]};
+EndIf
+If (NbrSectT == NbrSect)
+  Physical Line(RotorBoundary) = {RotorBoundary_[]};
+EndIf
+
+
+
+// moving band
+
+For i In {NbrSect:NbrSectT-1}
+If (i < NbrSectT)
+  For j In {0:1}
+    dP=newp-1;
+
+
+    Point(dP+4) = {ii_x,ii_y,0,0.7e-3};
+    Point(dP+16) = {0,l_y,0,0.3e-3};
+    Point(dP+18) = {oo_x,oo_y,0,0.6e-3};
+    Rotate {{0,0,1},{0,0,0}, RotorAngle_R+2*Pi*i/NbrSectT} { Point{dP+4}; Point{dP+16}; Point{dP+18}; }
+    If (j==1)
+      Symmetry {Cos(RotorAngle_S+2*Pi*i/NbrSectT),Sin(RotorAngle_S+2*Pi*i/NbrSectT),0,0} 
+                { Point{dP+4}; Point{dP+16}; Point{dP+18}; }  
+    EndIf
+
+    dR=newreg-1;
+    Circle(dR+19) = {dP+4,0,dP+18};
+    Circle(dR+20) = {dP+18,0,dP+16};
+
+    If (j==0)
+      InnerMB_[{4*i+2*j:4*i+2*j+1}] = {dR+19,dR+20};
+    EndIf
+    If (j==1)
+      InnerMB_[{4*i+2*j:4*i+2*j+1}] = {-dR-19,-dR-20};
+    EndIf
+  EndFor
+EndIf
+EndFor
+
+For i In {0:NbrPolesT-1}
+  Physical Line(InnerMB+i+1) = {InnerMB_[{i*4*NbrSect/NbrPoles:(i*4+4)*NbrSect/NbrPoles-1}]};
+EndFor
+  
+
+
+Coherence;
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/benchmarks/2d/machine3/Stator.geo b/benchmarks/2d/machine3/Stator.geo
new file mode 100644
index 0000000000..62bb04b6fd
--- /dev/null
+++ b/benchmarks/2d/machine3/Stator.geo
@@ -0,0 +1,303 @@
+
+Mesh.CharacteristicLengthFactor = 1.;
+
+Geometry.AutoCoherence = 0;      
+
+NbrPoles = 1; // number of poles in FE model
+NbrPolesT = 4; // number of poles in complete cross-section 
+NbrSectT = 48; // number of stator teeth
+NbrSect = NbrSectT*NbrPoles/NbrPolesT; // number of stator teeth in FE model
+
+StatorAngle_ = Pi/NbrSectT-Pi/2; // initial stator angle (radians)
+//StatorAngle_ = -2*Pi/NbrSectT; // initial stator angle (radians)
+StatorAngle_S = StatorAngle_;
+
+// physical stator numbers (for GetDP and Gmsh)
+
+StatorConductor= 1000; //12000
+ Stat_Up = StatorConductor+1; Stat_Wm = StatorConductor+2; Stat_Vp = StatorConductor+3; 
+ Stat_Um = StatorConductor+4; Stat_Wp = StatorConductor+5; Stat_Vm = StatorConductor+6; 
+StatorIron        = 10000;
+StatorSlotOpening = 14000;
+StatorAirgapLayer = 11000;
+
+OuterStator = 40000;
+StatorPeriod_Reference = 100000; 
+StatorPeriod_Dependent = 100001;
+StatorBoundary  = 800;
+OuterMB = 17000; MB_S1 = OuterMB+1; MB_S2 = OuterMB+2; MB_S3 = OuterMB+3; MB_S4 = OuterMB+4;
+      
+
+u=1e-3;
+AG = u* 0.47;       // airgap width
+R2 = u* 92/2;     // inner stator radius
+R3 = u* 150/2;      // outer stator radius
+//R1 = u* 90/2;       // outer radius of moving band
+R1 = R2-AG/3;       // outer radius of moving band
+
+// parameters for conductor and slot opening
+h1  = u* 1;
+h2  = u* 15.3;
+d1  = u* 2.5;
+Rsl = u* 6.36/2;
+ss = 0.05;
+
+RR = (h2-Rsl*(1+1/Sin(Pi/NbrSectT)))/(1-1/Sin(Pi/NbrSectT));
+Y1 = Sqrt(R2*R2-d1*d1/4) ;
+Y2 = Sqrt(RR*RR-d1*d1/4) ;
+Y3 = Sqrt(R1*R1-d1*d1/4) ;
+RX = Rsl*Cos(Pi/NbrSectT) ;
+RY = Rsl*Sin(Pi/NbrSectT) ;
+RX2 = RR*Cos(Pi/NbrSectT) ;
+RY2 = RR*Sin(Pi/NbrSectT) ;
+
+
+
+// characteristic lengths
+uc = u* 1.4 ;
+
+pslo = uc* 0.3; // slot opening
+psl  = uc* 0.6; // upper part slot 
+pslu = uc* 1; // lower part slot 
+pout = uc* 2; // outer radius
+pMB  = AG*1; // MB 
+//pMB  = uc* 0.5; // MB
+p  = uc* 2; // 
+
+
+lc_1=0.0004;
+lc_2=0.001;
+lc_3=0.003;
+a_x=0.001500;
+a_y=0.079986;
+b_x=0.001500;
+b_y=0.080900;
+c_x=0.002658;
+c_y=0.081380;
+d_x=0.003664;
+d_y=0.097100;
+e_x=0.000006;
+e_y=0.101000;
+f_x=0.080*Sin(3.75*(Pi/180.));
+f_y=0.080*Cos(3.75*(Pi/180.));
+g_x=0.120*Sin((3.75*(Pi/180.)));
+g_y=0.120*Cos((3.75*(Pi/180.)));
+h_y=0.097334;
+R_sup = Sqrt(d_x^2+(h_y-d_y)^2);
+i_x=0.000000;
+i_y=0.101000;
+j_x=0.002762;
+j_y=0.083000;
+k_x=((0.080-(0.00065/3))*Sin(3.75*(Pi/180.)));
+k_y=((0.080-(0.00065/3))*Cos(3.75*(Pi/180.)));
+l_x=0.001500;
+l_y=0.079986-0.00065/3;
+m_x=((0.080-(0.00065/3))*Sin(0*(Pi/180.)));
+m_y=((0.080-(0.00065/3))*Cos(0*(Pi/180.)));
+n_x=0.080*Sin(0*(Pi/180.));
+n_y=0.080*Cos(0*(Pi/180.));
+o_x=0.120*Sin((0*(Pi/180.)));
+o_y=0.120*Cos((0*(Pi/180.)));
+
+
+
+//Point(0)  = {0,0,0,p};
+
+
+For i In {0:NbrSect-1}
+For j In {0:1}
+
+dP=newp-1;
+Point(dP+1) = {a_x,a_y,0,lc_1};
+Point(dP+2) = {b_x,b_y,0,lc_1};
+Point(dP+3) = {c_x,c_y,0,lc_2*0.6};
+Point(dP+4) = {d_x,d_y,0,lc_3};
+Point(dP+5) = {0,h_y+R_sup,0,lc_3*0.5};
+Point(dP+6) = {g_x,g_y,0,lc_3*1.3};
+Point(dP+7) = {f_x,f_y,0,lc_1*1.4};
+Point(dP+8) = {0,h_y,0,lc_3};
+Point(dP+9) = {k_x,k_y,0,lc_1*1.4};
+Point(dP+10) = {l_x,l_y,0,lc_1};
+Point(dP+11) = {m_x,m_y,0,lc_1};
+Point(dP+12) = {n_x,n_y,0,lc_1};
+Point(dP+13) = {o_x,o_y,0,lc_3*1.3};
+Point(dP+14) = {j_x,j_y,0,lc_2};
+Point(dP+15) = {0,j_y,0,lc_2};
+
+
+For t In {dP+1:dP+15}
+  Rotate {{0,0,1},{0,0,0}, StatorAngle_+2*Pi*i/NbrSectT} {Point{t};}
+EndFor
+
+If (j==1)
+  For t In {dP+1:dP+15}
+    Symmetry {Cos(StatorAngle_S+2*Pi*i/NbrSectT),Sin(StatorAngle_S+2*Pi*i/NbrSectT),0,0} { Point{t}; }  
+  EndFor
+EndIf
+
+
+dR=newreg-1;
+Line(dR+1) = {dP+7,dP+6};
+Line(dR+2) = {dP+10,dP+1};
+Line(dR+3) = {dP+1,dP+2};
+Line(dR+4) = {dP+2,dP+3};
+Line(dR+5) = {dP+3,dP+14};
+Line(dR+6) = {dP+14,dP+4};
+Line(dR+7) = {dP+11,dP+12};
+Line(dR+8) = {dP+12,dP+15};
+Line(dR+9) = {dP+15,dP+8};
+Line(dR+10) = {dP+8,dP+5};
+Line(dR+11) = {dP+5,dP+13};
+Line(dR+12) = {dP+15,dP+14};
+Circle(dR+13) = {dP+6,0,dP+13};
+Circle(dR+14) = {dP+4,dP+8,dP+5};
+Circle(dR+15) = {dP+9,0,dP+10};
+Circle(dR+16) = {dP+10,0,dP+11};
+Circle(dR+17) = {dP+7,0,dP+1};
+Circle(dR+18) = {dP+1,0,dP+12};
+Line(dR+19) = {dP+9,dP+7};
+
+
+// physical lines
+
+OuterStator_[{2*i+j}] = dR+13;
+StatorBoundary_[{16*i+8*j:16*i+8*j+7}] = {dR+3,dR+4,dR+5,dR+6,dR+13,dR+14,dR+17,dR+12};
+
+If (j==0)
+  OuterMB_[{4*i+2*j:4*i+2*j+1}] = {dR+15,dR+16};
+EndIf
+If (j==1)
+  OuterMB_[{4*i+2*j:4*i+2*j+1}] = {-dR-15,-dR-16};
+EndIf
+
+
+If (NbrSectT != NbrSect)
+  If (i==0 && j==0)
+    Physical Line(StatorPeriod_Reference) = {dR+1,dR+19};
+    StatorPeriod_Reference_[] = {dR+1,dR+19};
+  EndIf
+  If (i == NbrSect-1  && j==1)
+    Physical Line(StatorPeriod_Dependent) = {dR+1,dR+19};
+    StatorPeriod_Dependent_[] = {dR+1,dR+19};
+  EndIf
+EndIf
+
+dH=newreg; 
+Line Loop(dH) = {dR+6,dR+14,-dR-10,-dR-9,dR+12};                    
+Plane Surface(dH) ={dH};
+StatorConductor_[2*i+j] = dH;
+
+dH=newreg;
+Line Loop(dH) = {dR+3,dR+4,dR+5,dR+6,dR+14,dR+11,-dR-13,-dR-1,dR+17};
+Plane Surface(dH) = {dH};
+StatorIron_[2*i+j] = dH;
+
+dH=newreg;     
+Line Loop(dH) = {-dR-12,-dR-8,-dR-18,dR+3,dR+4,dR+5};
+Plane Surface(dH) ={dH};
+StatorSlotOpening_[2*i+j] = dH;
+
+dH=newreg;
+Line Loop(dH) = {-dR-7,dR+17,dR+18,-dR-16,-dR-15,dR+19};
+Plane Surface(dH) = {dH};
+StatorAirgapLayer_[2*i+j] = dH;
+
+EndFor
+EndFor
+
+
+
+qq=4;
+
+For f In {0:5}
+  nCon=0;
+  For i In {0:NbrSect/qq-1}
+    If (Fmod(i,6) == f)
+      For j In {0:qq-1}
+        nCon+=2; Con[{nCon-2,nCon-1}] = {StatorConductor_[2*i*qq+2*j],StatorConductor_[2*i*qq+2*j+1]};
+      EndFor
+    EndIf
+  EndFor
+  If (nCon > 0)
+    Physical Surface(StatorConductor+1+f) = {Con[{0:nCon-1}]};
+// Color {R[f],G[f],B[f],1} {Surface{Con[{0:nCon-1}]};}
+    If (f == 0)
+      Color Pink {Surface{Con[{0:nCon-1}]};}
+    EndIf
+    If (f == 1)
+      Color Cyan {Surface{Con[{0:nCon-1}]};}
+    EndIf
+    If (f == 2)
+      Color Orange {Surface{Con[{0:nCon-1}]};}
+    EndIf
+  EndIf
+EndFor
+
+
+
+Physical Surface(StatorIron) = {StatorIron_[{0:NbrSect*2-1}]};
+Physical Surface(StatorSlotOpening) = {StatorSlotOpening_[{0:NbrSect*2-1}]};
+Physical Surface(StatorAirgapLayer) = {StatorAirgapLayer_[{0:NbrSect*2-1}]};
+
+Color Red {Surface{StatorIron_[{0:NbrSect*2-1}]};}
+Color Black {Surface{StatorSlotOpening_[{0:NbrSect*2-1}]};}
+Color Black {Surface{StatorAirgapLayer_[{0:NbrSect*2-1}]};}
+
+Physical Line(OuterStator) = {OuterStator_[]};
+
+If (NbrSectT != NbrSect)
+  Physical Line(StatorBoundary) = {StatorBoundary_[],StatorPeriod_Reference_[],StatorPeriod_Dependent_[]};
+EndIf
+If (NbrSectT == NbrSect)
+  Physical Line(StatorBoundary) = {StatorBoundary_[]};
+EndIf
+
+
+
+
+// moving band
+
+For i In {NbrSect:NbrSectT-1}
+If (i < NbrSectT)
+  For j In {0:1}
+    dP=newp-1;
+
+    Point(dP+9) = {k_x,k_y,0,lc_1*1.4};
+    Point(dP+10) = {l_x,l_y,0,lc_1};
+    Point(dP+11) = {m_x,m_y,0,lc_1};
+   // Point(dP+12) = {n_x,n_y,0,lc_1};
+
+    Rotate {{0,0,1},{0,0,0}, StatorAngle_+2*Pi*i/NbrSectT} { Point{dP+9}; Point{dP+10}; Point{dP+11}; }
+    If (j==1)
+      Symmetry {Cos(StatorAngle_S+2*Pi*i/NbrSectT),Sin(StatorAngle_S+2*Pi*i/NbrSectT),0,0} 
+                { Point{dP+9}; Point{dP+10}; Point{dP+11}; }  
+    EndIf
+    dR=newreg-1;
+    Circle(dR+15) = {dP+9,0,dP+10};
+    Circle(dR+16) = {dP+10,0,dP+11};
+
+    If (j==0)
+      OuterMB_[{4*i+2*j:4*i+2*j+1}] = {dR+15,dR+16};
+    EndIf
+    If (j==1)
+      OuterMB_[{4*i+2*j:4*i+2*j+1}] = {-dR-15,-dR-16};
+    EndIf
+  EndFor
+EndIf
+EndFor
+
+For i In {0:NbrPolesT-1}
+  Physical Line(OuterMB+i+1) = {OuterMB_[{i*4*NbrSect/NbrPoles:(i*4+4)*NbrSect/NbrPoles-1}]};
+EndFor
+  
+
+Coherence;
+
+
+
+
+
+
+
+
+
diff --git a/benchmarks/3d/periodic.geo b/benchmarks/3d/periodic.geo
new file mode 100644
index 0000000000..6196155fce
--- /dev/null
+++ b/benchmarks/3d/periodic.geo
@@ -0,0 +1,61 @@
+// original file from benarafa@alpes.cea.fr
+
+Lc = 0.003 ;
+nb = 12;
+use_prisms = 1;
+
+Point(1) = {0.,0.,0.,Lc};
+L = 0.045;
+R = 0.01085;
+Point(2) = {R,0.,0.,Lc};
+Point(3) = {-R,0.,0.,Lc};
+Point(4) = {0.,R,0.,Lc};
+Point(5) = {0.,-R,0.,Lc};
+Point(6) = {0.5*L,-0.5*L,0.,Lc};
+Point(7) = {-0.5*L,-0.5*L,0.,Lc};
+Point(8) = {-0.5*L,0.5*L,0.,Lc};
+Point(9) = {0.5*L,0.5*L,0.,Lc};
+Point(10) = {0.5*L-R,0.5*L,0.,Lc};
+Point(11) = {0.5*L,0.5*L-R,0.,Lc};
+Point(12) = {-0.5*L+R,0.5*L,0.,Lc};
+Point(13) = {-0.5*L+R,-0.5*L,0.,Lc};
+Point(14) = {-0.5*L,-0.5*L+R,0.,Lc};
+Point(15) = {0.5*L,-0.5*L+R,0.,Lc};
+Point(16) = {-0.5*L,0.5*L-R,0.,Lc};
+
+Point(17) = {0.5*L-R,-0.5*L,0.,Lc};
+Line(1) = {12,10};
+Line(2) = {16,14};
+Line(3) = {11,15};
+Line(4) = {13,17};
+Circle(5) = {10,9,11};
+Circle(6) = {15,6,17};
+Circle(7) = {13,7,14};
+Circle(8) = {16,8,12};
+Circle(9) = {4,1,2};
+Circle(10) = {2,1,5};
+Circle(11) = {5,1,3};
+
+Circle(12) = {3,1,4};
+Line Loop(13) = {3,6,-4,7,-2,8,1,5};
+Line Loop(14) = {10,11,12,9};
+Plane Surface(15) = {13,14};
+Physical Line(16) = {5,6,7,8,12,9,10,11};
+Physical Line(17) = {1,4};
+Physical Line(18) = {2,3};
+Physical Surface(19) = {15};
+
+If(use_prisms)
+  Extrude Surface {15, {0.,0.,2.*R}}{Layers{nb,83,1}; Recombine; };
+EndIf
+
+If(!use_prisms)
+  Extrude Surface {15, {0.,0.,2.*R}};
+  Transfinite Line {27,55,1,59,23,43,4,39,21,34,3,35,25,51,2,47} = nb Using Progression 1.;
+  Transfinite Surface {52} = {16,14,33,37};
+  Transfinite Surface {36} = {11,15,19,18};
+  Transfinite Surface {44} = {13,17,24,28};
+  Transfinite Surface {60} = {12,10,46,42};
+  Surface Loop(82) = {81,36,15,40,44,48,52,56,60,64,68,72,76,80};
+  Volume(83) = {82};
+EndIf
-- 
GitLab