complicated bi-sin solution - periodic surface fails

DiffractionGratings/grating3D.geo

@@ -21,6 +21,78 @@ Macro SetPBCs
     Periodic Surface{slaveY()} = {masterY()} Translate{       0, period_y,        0};
 Return
 
+If (tag_geom==6)
+    nsin = 20;
+    For i In {0:nsin}
+        xx1~{i} = -period_x/2+i*period_x/2/(nsin-1);
+        xx2~{i} =             i*period_x/2/(nsin-1);
+    EndFor
+    pp=newp;//pp=pp-1;
+    For i In {0:nsin-1}
+        yy=-period_y/2;
+        Point(newp) = {xx1~{i},yy,rz/4*(Cos(2*Pi*xx1~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+        Point(newp) = {yy,xx1~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx1~{i}/period_y))+hh_L_3+thick_L_3/2};
+        yy=0;
+        Point(newp) = {xx1~{i},yy,rz/4*(Cos(2*Pi*xx1~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+        Point(newp) = {yy,xx1~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx1~{i}/period_y))+hh_L_3+thick_L_3/2};
+        yy=period_y/2;
+        Point(newp) = {xx1~{i},yy,rz/4*(Cos(2*Pi*xx1~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+        Point(newp) = {yy,xx1~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx1~{i}/period_y))+hh_L_3+thick_L_3/2};
+        yy=-period_y/2;
+        Point(newp) = {xx2~{i},yy,rz/4*(Cos(2*Pi*xx2~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+        Point(newp) = {yy,xx2~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx2~{i}/period_y))+hh_L_3+thick_L_3/2};
+        yy=0;
+        Point(newp) = {xx2~{i},yy,rz/4*(Cos(2*Pi*xx2~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+        Point(newp) = {yy,xx2~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx2~{i}/period_y))+hh_L_3+thick_L_3/2};
+        yy=period_y/2;
+        Point(newp) = {xx2~{i},yy,rz/4*(Cos(2*Pi*xx2~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
+        Point(newp) = {yy,xx2~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx2~{i}/period_y))+hh_L_3+thick_L_3/2};
+    EndFor
+    For k In {1:11}
+        BSpline(newl) = {pp+k:pp+12*nsin:12};
+    EndFor
+    BSpline(newl) = {pp  :pp+12*(nsin-1):12};
+    Coherence;
+    Box(1) = {-period_x/2,-period_y/2,     hh_L_3,period_x,period_y, thick_L_3};
+    Delete{Volume{1}; Surface{1,2,3,4} ; Line{13,15,17,19};}
+    Curve Loop(7) = {21, 20, -23, -16};
+    Curve Loop(8) = {1, 2, -3, -12};
+    Surface(7) = {8};
+    Curve Loop(10) = {6, 5, -8, -3};
+    Surface(8) = {10};
+    Curve Loop(12) = {9, 10, -11, -8};
+    Surface(9) = {12};
+    Curve Loop(14) = {2, 9, -4, -7};
+    Surface(10) = {14};
+    Line(25) = {239, 229};
+    Line(26) = {229, 238};
+    Line(27) = {243, 235};
+    Line(28) = {235, 242};
+    Line(29) = {237, 244};
+    Line(30) = {233, 240};
+    Line(31) = {237, 245};
+    Line(32) = {241, 233};
+    Curve Loop(16) = {25, 12, 6, -27, -21};
+    Surface(11) = {16};
+    Curve Loop(18) = {20, -31, -11, -5, -27};
+    Surface(12) = {18};
+    Curve Loop(20) = {29, -18, -28, 5, 11};
+    Surface(13) = {20};
+    Curve Loop(22) = {10, 29, -24, -30, 4};
+    Surface(14) = {22};
+    Curve Loop(24) = {30, -14, -26, 1, 7};
+    Surface(15) = {24};
+    Curve Loop(26) = {25, 1, 7, -32, -16};
+    Surface(16) = {26};
+    Curve Loop(28) = {26, 22, -28, -6, -12};
+    Surface(17) = {28};
+    Curve Loop(30) = {31, -23, 32, 4, 10};
+    Surface(18) = {30};
+    Surface Loop(3) = {5, 16, 18, 12, 11, 9, 8, 7, 10};
+    Volume(9) = {3};    
+    Surface Loop(2) = {6, 15, 14, 10, 7, 8, 9, 13, 17};
+    Volume(10) = {2};
+EndIf
 
 L1_n     = Sqrt(eps_re_L_1);
 L2_n     = Sqrt(eps_re_L_2);
@@ -95,101 +167,61 @@ If (tag_geom==5)
     Rotate {{0, 0, 1}, {0,0,0}, Pi/4} {Volume{9};}
 EndIf
 
-If (tag_geom==6)
-    nsin = 10;
-    For i In {0:nsin}
-        xx1~{i} = -period_x/2+i*period_x/2/(nsin-1);
-        xx2~{i} =             i*period_x/2/(nsin-1);
-    EndFor
-    pp=newp;//pp=pp-1;
-    For i In {0:nsin-1}
-        yy=-period_y/2;
-        Point(newp) = {xx1~{i},yy,rz/4*(Cos(2*Pi*xx1~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
-        Point(newp) = {yy,xx1~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx1~{i}/period_y))+hh_L_3+thick_L_3/2};
-        yy=0;
-        Point(newp) = {xx1~{i},yy,rz/4*(Cos(2*Pi*xx1~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
-        Point(newp) = {yy,xx1~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx1~{i}/period_y))+hh_L_3+thick_L_3/2};
-        yy=period_y/2;
-        Point(newp) = {xx1~{i},yy,rz/4*(Cos(2*Pi*xx1~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
-        Point(newp) = {yy,xx1~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx1~{i}/period_y))+hh_L_3+thick_L_3/2};
-        yy=-period_y/2;
-        Point(newp) = {xx2~{i},yy,rz/4*(Cos(2*Pi*xx2~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
-        Point(newp) = {yy,xx2~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx2~{i}/period_y))+hh_L_3+thick_L_3/2};
-        yy=0;
-        Point(newp) = {xx2~{i},yy,rz/4*(Cos(2*Pi*xx2~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
-        Point(newp) = {yy,xx2~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx2~{i}/period_y))+hh_L_3+thick_L_3/2};
-        yy=period_y/2;
-        Point(newp) = {xx2~{i},yy,rz/4*(Cos(2*Pi*xx2~{i}/period_x)+Cos(2*Pi*yy/period_y))+hh_L_3+thick_L_3/2};
-        Point(newp) = {yy,xx2~{i},rz/4*(Cos(2*Pi*yy/period_x)+Cos(2*Pi*xx2~{i}/period_y))+hh_L_3+thick_L_3/2};
-    EndFor
-    For k In {1:11}
-        BSpline(newl) = {pp+k:pp+12*nsin:12};
-    EndFor
-    BSpline(newl) = {pp  :pp+12*(nsin-1):12};
-    Line(newl) = {25,pp};
-    Line(newl) = {34,pp};
-    Line(newl) = {29,pp+5};
-    Line(newl) = {38,pp+5};
-    Line(newl) = {27,pp+4};
-    Line(newl) = {36,pp+4};
-    pp=newp;
-    Line(newl) = {40,pp-1};
-    Line(newl) = {31,pp-1};
-    Surface(news)=
-    {85,86,87,96};
-EndIf
 
-// Coherence;
+
+
+Coherence;
 // // // // Physical Surface("ENDPML",770) = {69,37};
 
-// Call SetPBCs;
-
-// Physical Volume("PMLBOT"        ,1) = {1};
-// Physical Volume("LAYER_L6_SUBS" ,2) = {2};
-// Physical Volume("LAYER_L5_LOW"  ,3) = {3};
-// Physical Volume("LAYER_L4_GUI"  ,4) = {4};
-// Physical Volume("LAYER_L3_EMB"  ,5) = {10};
-// Physical Volume("LAYER_L2_TOP"  ,6) = {6};
-// Physical Volume("LAYER_L1_SUPER",7) = {7};
-// Physical Volume("PMLTOP"        ,8) = {8};
-// Physical Volume("SCAT"          ,9) = {9};
-
-// 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 Volume("DEPOSIT"       ,10000) = {11};
-
-// pts_PMLBOT()   = PointsOf{ Physical Volume{1}; };
-// pts_LAYER_L6() = PointsOf{ Physical Volume{2}; };
-// pts_LAYER_L5() = PointsOf{ Physical Volume{3}; };
-// pts_LAYER_L4() = PointsOf{ Physical Volume{4}; };
-// pts_LAYER_L3() = PointsOf{ Physical Volume{5}; };
-// pts_LAYER_L2() = PointsOf{ Physical Volume{6}; };
-// pts_LAYER_L1() = PointsOf{ Physical Volume{7}; };
-// pts_PMLTOP()   = PointsOf{ Physical Volume{8}; };
-// pts_ROD()      = PointsOf{ Physical Volume{9}; };
-// // pts_DEPOSIT()  = PointsOf{ Physical Volume{10000}; };
-// // For k In {0:#pts_PMLBOT()-1}
-// //     Printf("blabal %g %f PML_lc;",pts_PMLBOT(k),PML_lc);
-// // EndFor
-
-// Characteristic Length{:}              = PML_lc;
-// Characteristic Length{pts_LAYER_L1()} =  L1_lc/refine_mesh_L1;
-// Characteristic Length{pts_LAYER_L2()} =  L2_lc/refine_mesh_L2;
-// Characteristic Length{pts_LAYER_L4()} =  L4_lc/refine_mesh_L4;
-// Characteristic Length{pts_LAYER_L5()} =  L5_lc/refine_mesh_L5;
-// Characteristic Length{pts_LAYER_L6()} =  L6_lc/refine_mesh_L6;
-// Characteristic Length{pts_LAYER_L3()} =  L3_lc/refine_mesh_L3;
-// Characteristic Length{pts_ROD()}      =  scat_lc;
-// If (tag_geom==3)
-//     Mesh.Algorithm = 6;
-// EndIf
-// Mesh.SurfaceEdges = 0;
-// Mesh.VolumeEdges = 0;
+Call SetPBCs;
+
+Physical Volume("PMLBOT"        ,1) = {1};
+Physical Volume("LAYER_L6_SUBS" ,2) = {2};
+Physical Volume("LAYER_L5_LOW"  ,3) = {3};
+Physical Volume("LAYER_L4_GUI"  ,4) = {4};
+Physical Volume("LAYER_L3_EMB"  ,5) = {10};
+Physical Volume("LAYER_L2_TOP"  ,6) = {6};
+Physical Volume("LAYER_L1_SUPER",7) = {7};
+Physical Volume("PMLTOP"        ,8) = {8};
+Physical Volume("SCAT"          ,9) = {9};
+
+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 Volume("DEPOSIT"       ,10000) = {11};
+
+pts_PMLBOT()   = PointsOf{ Physical Volume{1}; };
+pts_LAYER_L6() = PointsOf{ Physical Volume{2}; };
+pts_LAYER_L5() = PointsOf{ Physical Volume{3}; };
+pts_LAYER_L4() = PointsOf{ Physical Volume{4}; };
+pts_LAYER_L3() = PointsOf{ Physical Volume{5}; };
+pts_LAYER_L2() = PointsOf{ Physical Volume{6}; };
+pts_LAYER_L1() = PointsOf{ Physical Volume{7}; };
+pts_PMLTOP()   = PointsOf{ Physical Volume{8}; };
+pts_ROD()      = PointsOf{ Physical Volume{9}; };
+// pts_DEPOSIT()  = PointsOf{ Physical Volume{10000}; };
+// For k In {0:#pts_PMLBOT()-1}
+//     Printf("blabal %g %f PML_lc;",pts_PMLBOT(k),PML_lc);
+// EndFor
+
+Characteristic Length{:}              = PML_lc;
+Characteristic Length{pts_LAYER_L1()} =  L1_lc/refine_mesh_L1;
+Characteristic Length{pts_LAYER_L2()} =  L2_lc/refine_mesh_L2;
+Characteristic Length{pts_LAYER_L4()} =  L4_lc/refine_mesh_L4;
+Characteristic Length{pts_LAYER_L5()} =  L5_lc/refine_mesh_L5;
+Characteristic Length{pts_LAYER_L6()} =  L6_lc/refine_mesh_L6;
+Characteristic Length{pts_LAYER_L3()} =  L3_lc/refine_mesh_L3;
+Characteristic Length{pts_ROD()}      =  scat_lc;
+If (tag_geom==3)
+    Mesh.Algorithm = 6;
+EndIf
+Mesh.SurfaceEdges = 0;
+Mesh.VolumeEdges = 0;
+