diff --git a/DiffractionGratings/grating3D.geo b/DiffractionGratings/grating3D.geo
index b596ca1f40e00b8d22ca9b6619118303a45262eb..5047c4e81edec4702f6320eba8504061bff05217 100644
--- a/DiffractionGratings/grating3D.geo
+++ b/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};
-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 };
+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)
+ 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
+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}; };
@@ -223,4 +258,4 @@ If (tag_geom==3) // Split torus weird otherwise
EndIf
// Mesh.SurfaceEdges = 0;
Mesh.VolumeEdges = 0;
-Mesh.ElementOrder = og;
+Mesh.ElementOrder = og;
\ No newline at end of file