diff --git a/DiffractionGratings/grating3D.geo b/DiffractionGratings/grating3D.geo
index 50dfc42dfc38d34a188725db830bf0e204bb9367..e028fcbbb3a4ed64b7c0d06fcfc435105ddf3549 100644
--- a/DiffractionGratings/grating3D.geo
+++ b/DiffractionGratings/grating3D.geo
@@ -1,97 +1,96 @@
-// Include "grating3D_data_torus.geo";
-// Include "grating3D_data_hole.geo";
-// Include "grating3D_data_pyramid.geo";
-// Include "grating3D_data_halfellipsoid.geo";
Include StrCat["grating3D_data_",test_case,".geo"];
SetFactory("OpenCASCADE");
-// Geometry.LineWidth = 5;
-// Geometry.Points = 0;
-// Geometry.Color.Lines = {255,66,55};
-
-// General.BackgroundImageFileName="struct2.png";
Macro SetPBCs
- e = 1e-6;
- 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};
- Periodic Surface{slaveX()} = {masterX()} Translate{period_x, 0, 0};
- Periodic Surface{slaveY()} = {masterY()} Translate{ 0, period_y, 0};
+ e = 10*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};
+ If (tag_geom==6) // bi-sin : BoundingBox does not catch BSpline Surfaces?
+ masterX()+={48,52};
+ slaveX()+={50,54};
+ masterY()+={51,55};
+ slaveY()+={49,53};
+ EndIf
+ Periodic Surface{slaveX()} = {masterX()} Translate{period_x, 0, 0};
+ 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};
+ // Bi-sinusoidal profile by Spline surface filling
+ // not the most general way to implement a freesurface height=f(x,y) and relies on multiple Coherence calls
+ 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);
@@ -115,56 +114,50 @@ 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};
+ 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 };
-w_arm1 = 35*nm;
-w_arm2 = 35*nm;
-w_arm3 = 35*nm;
-
-// tag_geom = 1; // half ellipsoid
-// tag_geom = 1; // pyramid
If (tag_geom==1)
- p=newp;
- Point(p) = {0,0,hh_L_3+rz};
- Line(97) = {25, 65};
- Line(98) = {65, 27};
- Line(99) = {65, 29};
- Line(100) = {31, 65};
- Curve Loop(49) = {38, -98, -97};
- Plane Surface(49) = {49};
- Curve Loop(50) = {98, 48, 100};
- Plane Surface(50) = {50};
- Curve Loop(51) = {99, 42, 100};
- Plane Surface(51) = {51};
- Curve Loop(52) = {99, -46, 97};
- Plane Surface(52) = {52};
- Surface Loop(9) = {52, 51, 50, 49, 24};
- Volume(9) = {9};
+ p=newp;
+ Point(p) = {0,0,hh_L_3+rz};
+ Line(97) = {25, 65};
+ Line(98) = {65, 27};
+ Line(99) = {65, 29};
+ Line(100) = {31, 65};
+ Curve Loop(49) = {38, -98, -97};
+ Plane Surface(49) = {49};
+ Curve Loop(50) = {98, 48, 100};
+ Plane Surface(50) = {50};
+ Curve Loop(51) = {99, 42, 100};
+ Plane Surface(51) = {51};
+ Curve Loop(52) = {99, -46, 97};
+ Plane Surface(52) = {52};
+ Surface Loop(9) = {52, 51, 50, 49, 24};
+ Volume(9) = {9};
EndIf
If (tag_geom==2)
- Cylinder(9) = {0,0,hh_L_3,0,0,thick_L_3,rx};
+ Cylinder(9) = {0,0,hh_L_3,0,0,thick_L_3,rx};
EndIf
If (tag_geom==3)
- Torus(9) = {0,0,hh_L_3+ry-2,rx,ry};
- Dilate { { 0,0,hh_L_3 }, { 1, 1, rz/ry } } { Volume{9}; }
- BooleanDifference{ Volume{9}; Delete; }{ Volume{4}; }
+ Torus(9) = {0,0,hh_L_3+ry-2,rx,ry};
+ Dilate { { 0,0,hh_L_3 }, { 1, 1, rz/ry } } { Volume{9}; }
+ BooleanDifference{ Volume{9}; Delete; }{ Volume{4}; }
EndIf
If (tag_geom==4)
- Sphere(9) = {0,0,hh_L_3,rx};
- Dilate { { 0,0,hh_L_3 }, { 1, ry/rx, rz/rx } } { Volume{9}; }
- BooleanDifference{ Volume{9}; Delete; }{ Volume{4}; }
+ Sphere(9) = {0,0,hh_L_3,rx};
+ Dilate { { 0,0,hh_L_3 }, { 1, ry/rx, rz/rx } } { Volume{9}; }
+ BooleanDifference{ Volume{9}; Delete; }{ Volume{4}; }
EndIf
If (tag_geom==5)
- Box(9) = {-rx/2,-ry/2, hh_L_2-rz, rx, ry, rz};
- Rotate {{0, 0, 1}, {0,0,0}, Pi/4} {Volume{9};}
+ Box(9) = {-rx/2,-ry/2, hh_L_2-rz, rx, ry, rz};
+ Rotate {{0, 0, 1}, {0,0,0}, Pi/4} {Volume{9};}
EndIf
@@ -194,9 +187,6 @@ Physical Surface("SBOT",302 ) = Surface In BoundingBox{-period_x/2-e,-period_y/2
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}; };
@@ -206,22 +196,17 @@ 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;
+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) // Split torus weird otherwise
+ Mesh.Algorithm = 6;
EndIf
Mesh.SurfaceEdges = 0;
Mesh.VolumeEdges = 0;
-