Skip to content
Snippets Groups Projects
Select Git revision
  • 6f56a703591e4abb1d0b77f9d57e2330d2637664
  • master default protected
  • albertpiwonski-master-patch-57409
  • quadspheres
  • fix_Tmatrix_code_epsr_background
  • albertpiwonski-master-patch-12427
  • cavity
  • c1
8 results

grating3D.geo

Blame
  • grating3D.geo 8.43 KiB
    // test_case = "pyramid";
    // test_case = "hole";
    // test_case = "torus";
    // test_case = "halfellipsoid";
    // test_case = "checker";
    // test_case = "bisin";
    Include StrCat["grating3D_data_",test_case,".geo"];
    
    SetFactory("OpenCASCADE");
    
    Macro SetPBCs
      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)
      // 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
    
    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;
      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) = {97, 99, -46};
      Plane Surface(49) = {49};
      Curve Loop(50) = {99, 42, 100};
      Plane Surface(50) = {50};
      Curve Loop(51) = {98, -38, 97};
      Plane Surface(51) = {51};
      Curve Loop(52) = {98, 48, 100};
      Plane Surface(52) = {52};
      Curve Loop(53) = {38, 48, -42, -46};
      Plane Surface(53) = {53};
      Surface Loop(9) = {49, 51, 52, 50, 53};
      Volume(9) = {9};
    EndIf
    
    If (tag_geom==2)
      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-1*nm,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}; }
    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};}
    EndIf
    
    If (tag_geom==7)
      Box(9) = {-period_x/2,-ry/2, hh_L_3, period_x, ry, rz};
    EndIf
    
    
    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"      ,3) = {3};
    Physical Volume("LAYER_L4"      ,4) = {4};
    Physical Volume("LAYER_L3"      ,5) = {10};
    Physical Volume("LAYER_L2"      ,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};
    
    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}; };
    
    lc_PML = lambda_m/(paramaille*.5);
    list_lc(0) = lc_PML;
    For k In {1:6}
      n_L~{k}  = Sqrt(Abs(eps_re_L~{k}));
      lc_L~{k} = lambda_m/(paramaille*n_L~{k})/refine_mesh_L~{k};
      list_lc(7-k) = lc_L~{k};
    EndFor
    list_lc(7) = lc_PML;
    list_lc(8) = lc_scat;
    
    // This helps meshing: The default behavior of the PointsOf techinque
    // overides points belonging to several domains
    If (tag_geom==7)
      meshing_sequence() = {1,8,2,3,5,6,7,4,9};
    Else
        meshing_sequence() = {1,8,2,3,4,6,7,5,9};
    EndIf
    
    // Start with coarsest
    Characteristic Length{:} = lc_PML;
    
    For k In {0:8}
      which_dom = meshing_sequence(k);
      Characteristic Length{PointsOf{Physical Volume{which_dom};}} = list_lc(which_dom-1);
    EndFor
    
    
    If (tag_geom==3) // Split torus weird otherwise
      Mesh.Algorithm = 6;
    EndIf
    // Mesh.Optimimze3D=1;
    // Mesh.SurfaceEdges = 0;
    Mesh.VolumeEdges = 0;
    Mesh.ElementOrder = og;