skew geometry

DiffractionGratings/grating3D.geo

@@ -7,13 +7,26 @@
 Include StrCat["grating3D_data_",test_case,".geo"];
 
 SetFactory("OpenCASCADE");
+e = 5*nm;
+E = 10*(period_x+period_x/2);
 
 Macro SetPBCs
-  e = 5*nm;
-  masterX() = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, PML_bot_hh-e,-period_x/2+e, period_y/2+e, PML_top_hh+PML_top+e};
-  slaveX()  = Surface In BoundingBox{ period_x/2-e,-period_y/2-e, PML_bot_hh-e, period_x/2+e, period_y/2+e, PML_top_hh+PML_top+e};
-  masterY() = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, PML_bot_hh-e, period_x/2+e,-period_y/2+e, PML_top_hh+PML_top+e};
-  slaveY()  = Surface In BoundingBox{-period_x/2-e, period_y/2-e, PML_bot_hh-e, period_x/2+e, period_y/2+e, PML_top_hh+PML_top+e};
+  masterX() = Surface In BoundingBox{.5*(-period_x-dys)-e,-dyc/2-e, PML_bot_hh-e,(-period_x+dys)/2+e, dyc/2+e, PML_top_hh+PML_top+e};
+  slaveX()  = Surface In BoundingBox{.5*( period_x-dys)-e,-dyc/2-e, PML_bot_hh-e,( period_x+dys)/2+e, dyc/2+e, PML_top_hh+PML_top+e};
+  masterY() = Surface In BoundingBox{.5*(-period_x-dys)-e,-dyc/2-e, PML_bot_hh-e,( period_x+dys)/2+e,-dyc/2+e, PML_top_hh+PML_top+e};
+  slaveY()  = Surface In BoundingBox{.5*(-period_x-dys)-e, dyc/2-e, PML_bot_hh-e,( period_x+dys)/2+e, dyc/2+e, PML_top_hh+PML_top+e};
+  // For k In {1:#masterX()-1}
+  //   Printf("masterX surf %g",masterX(k));
+  // EndFor
+  // For k In {1:#masterY()-1}
+  //   Printf("masterY surf %g",masterY(k));
+  // EndFor
+  // For k In {1:#slaveX()-1}
+  //   Printf("slaveX surf %g",slaveX(k));
+  // EndFor
+  // For k In {1:#slaveY()-1}
+  //   Printf("slaveY surf %g",slaveY(k));
+  // EndFor
   If (tag_geom==6) // bi-sin : BoundingBox does not catch BSpline Surfaces?
     masterX()+={48,52};
     slaveX()+={50,54};
@@ -21,7 +34,7 @@ Macro SetPBCs
     slaveY()+={49,53};
   EndIf
   Periodic Surface{slaveX()} = {masterX()} Translate{period_x,        0,        0};
-  Periodic Surface{slaveY()} = {masterY()} Translate{       0, period_y,        0};
+  Periodic Surface{slaveY()} = {masterY()} Translate{    dys, dyc,        0};
 Return
 
 If (tag_geom==6)
@@ -99,16 +112,38 @@ If (tag_geom==6)
   Volume(10) = {2};
 EndIf
 
-Box(1) = {-period_x/2,-period_y/2, PML_bot_hh,period_x,period_y, PML_bot  };
-Box(2) = {-period_x/2,-period_y/2,     hh_L_6,period_x,period_y, thick_L_6};
-Box(3) = {-period_x/2,-period_y/2,     hh_L_5,period_x,period_y, thick_L_5};
-Box(4) = {-period_x/2,-period_y/2,     hh_L_4,period_x,period_y, thick_L_4};
+dyc = period_y*Cos(xsi);
+dys = period_y*Sin(xsi);
+hh_L_7 = PML_bot_hh;
+hh_L_0 = PML_top_hh;
+thick_L_7 = PML_bot;
+thick_L_0 = PML_top;
+
+For k In {7:0:-1}
   If (tag_geom!=6)
-  Box(5) = {-period_x/2,-period_y/2,     hh_L_3,period_x,period_y, thick_L_3};
+    p=newp; l=newl; ll=newll; s=news;
+    Point(p)   = {0.5*(-period_x-dys), -dyc/2, hh_L~{k}};
+    Point(p+1) = {0.5*( period_x-dys), -dyc/2, hh_L~{k}};
+    Point(p+2) = {0.5*( period_x+dys),  dyc/2, hh_L~{k}};
+    Point(p+3) = {0.5*(-period_x+dys),  dyc/2, hh_L~{k}};
+    Line(l)={p,p+1};Line(l+1)={p+1,p+2};Line(l+2)={p+2,p+3};Line(l+3)={p+3,p};
+    Curve Loop(ll) = {l,l+1,l+2,l+3};
+    Surface(s) = {ll};
+    Extrude {0, 0, thick_L~{k}} {Surface{s};}
+  Else
+    v=newv;
   EndIf
-Box(6) = {-period_x/2,-period_y/2,     hh_L_2,period_x,period_y, thick_L_2};
-Box(7) = {-period_x/2,-period_y/2,     hh_L_1,period_x,period_y, thick_L_1};
-Box(8) = {-period_x/2,-period_y/2, PML_top_hh,period_x,period_y, PML_top  };
+EndFor
+// Box(1) = {-period_x/2,-period_y/2, PML_bot_hh,period_x,period_y, PML_bot  };
+// Box(2) = {-period_x/2,-period_y/2,     hh_L_6,period_x,period_y, thick_L_6};
+// Box(3) = {-period_x/2,-period_y/2,     hh_L_5,period_x,period_y, thick_L_5};
+// Box(4) = {-period_x/2,-period_y/2,     hh_L_4,period_x,period_y, thick_L_4};
+// If (tag_geom!=6)
+//   Box(5) = {-period_x/2,-period_y/2,     hh_L_3,period_x,period_y, thick_L_3};
+// EndIf
+// Box(6) = {-period_x/2,-period_y/2,     hh_L_2,period_x,period_y, thick_L_2};
+// Box(7) = {-period_x/2,-period_y/2,     hh_L_1,period_x,period_y, thick_L_1};
+// Box(8) = {-period_x/2,-period_y/2, PML_top_hh,period_x,period_y, PML_top  };
 
 If (tag_geom==1)
   p=newp;
@@ -158,7 +193,7 @@ EndIf
 
 
 Coherence;
-// // // // Physical Surface("ENDPML",770) = {69,37};
+// // // Physical Surface("ENDPML",770) = {69,37};
 
 Call SetPBCs;
 
@@ -176,10 +211,10 @@ Physical Surface("BXM" ,101 ) = masterX();
 Physical Surface("BXP" ,102 ) = slaveX();
 Physical Surface("BYM" ,201 ) = masterY();
 Physical Surface("BYP" ,202 ) = slaveY();
-Physical Surface("STOP",301 ) = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, PML_top_hh-e,period_x/2+e, period_y/2+e, PML_top_hh+e};
-Physical Surface("SBOT",302 ) = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, hh_L_6-e,period_x/2+e, period_y/2+e, hh_L_6+e};
-Physical Surface("SPMLTOP",401 ) = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, PML_top_hh+PML_top-e,period_x/2+e, period_y/2+e, PML_top_hh+PML_top+e};
-Physical Surface("SPMLBOT",402 ) = Surface In BoundingBox{-period_x/2-e,-period_y/2-e, PML_bot_hh-e,period_x/2+e, period_y/2+e, PML_bot_hh+e};
+Physical Surface("STOP",301 ) = Surface In BoundingBox{-period_x/2-E,-period_y/2-E, PML_top_hh-e,period_x/2+E, period_y/2+E, PML_top_hh+e};
+Physical Surface("SBOT",302 ) = Surface In BoundingBox{-period_x/2-E,-period_y/2-E, hh_L_6-e,period_x/2+E    , period_y/2+E, hh_L_6+e};
+Physical Surface("SPMLTOP",401 ) = Surface In BoundingBox{-period_x/2-E,-period_y/2-E, PML_top_hh+PML_top-e,period_x/2+E, period_y/2+E, PML_top_hh+PML_top+e};
+Physical Surface("SPMLBOT",402 ) = Surface In BoundingBox{-period_x/2-E,-period_y/2-E, PML_bot_hh-e,period_x/2+E, period_y/2+E, PML_bot_hh+e};
 
 pts_PMLBOT()   = PointsOf{ Physical Volume{1}; };
 pts_LAYER_L6() = PointsOf{ Physical Volume{2}; };