######################## Mesh Generation function with GMSH ####################
function MeshGenerator(L,h1,h2,h3,hd,dpml,l0,ld,lpml)
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 1)
    gmsh.option.setNumber("Mesh.Algorithm", 6)
    gmsh.clear()
    gmsh.model.add("geometry") # name it whatever you want

    # Add points
    gmsh.model.geo.addPoint(-L/2-dpml, -hd-h2-h3-dpml, 0, lpml,  1)
    gmsh.model.geo.addPoint(-L/2     , -hd-h2-h3-dpml, 0, lpml,  2)
    gmsh.model.geo.addPoint(-L/2-dpml, -hd-h2-h3     , 0, lpml,  3)
    gmsh.model.geo.addPoint(-L/2     , -hd-h2-h3     , 0,   l0,  4)
    gmsh.model.geo.addPoint( L/2     , -hd-h2-h3-dpml, 0, lpml,  5)
    gmsh.model.geo.addPoint( L/2     , -hd-h2-h3     , 0,   l0,  6)
    gmsh.model.geo.addPoint( L/2+dpml, -hd-h2-h3-dpml, 0, lpml,  7)
    gmsh.model.geo.addPoint( L/2+dpml, -hd-h2-h3     , 0, lpml,  8)
    gmsh.model.geo.addPoint(-L/2-dpml, -hd-h2        , 0,   l0,  9)
    gmsh.model.geo.addPoint(-L/2     , -hd-h2        , 0,   l0, 10)
    gmsh.model.geo.addPoint( L/2     , -hd-h2        , 0,   l0, 11)
    gmsh.model.geo.addPoint( L/2+dpml, -hd-h2        , 0,   l0, 12)
    gmsh.model.geo.addPoint(-L/2-dpml, -hd           , 0,   ld, 13)
    gmsh.model.geo.addPoint(-L/2     , -hd           , 0,   ld, 14)
    gmsh.model.geo.addPoint( L/2     , -hd           , 0,   ld, 15)
    gmsh.model.geo.addPoint( L/2+dpml, -hd           , 0,   ld, 16)
    gmsh.model.geo.addPoint(-L/2-dpml, 0.            , 0,   ld, 17)
    gmsh.model.geo.addPoint(-L/2     , 0.            , 0,   ld, 18)
    gmsh.model.geo.addPoint( L/2     , 0.            , 0,   ld, 19)
    gmsh.model.geo.addPoint( L/2+dpml, 0.            , 0,   ld, 20)
    gmsh.model.geo.addPoint(-L/2-dpml, h1            , 0, lpml, 21)
    gmsh.model.geo.addPoint(-L/2     , h1            , 0,   l0, 22)
    gmsh.model.geo.addPoint( L/2     , h1            , 0,   l0, 23)
    gmsh.model.geo.addPoint( L/2+dpml, h1            , 0, lpml, 24)
    gmsh.model.geo.addPoint(-L/2-dpml, h1+dpml       , 0, lpml, 25)
    gmsh.model.geo.addPoint(-L/2     , h1+dpml       , 0, lpml, 26)
    gmsh.model.geo.addPoint( L/2     , h1+dpml       , 0, lpml, 27)
    gmsh.model.geo.addPoint( L/2+dpml, h1+dpml       , 0, lpml, 28)
    #gmsh.model.geo.addPoint(-Lt/2, h1/2-ht/2, 0, l0, 29)
    #gmsh.model.geo.addPoint( Lt/2, h1/2-ht/2, 0, l0, 30)
    #gmsh.model.geo.addPoint( Lt/2, h1/2+ht/2, 0, l0, 31)
    #gmsh.model.geo.addPoint(-Lt/2, h1/2+ht/2, 0, l0, 32)
    #gmsh.model.geo.addPoint(x_s,y_s, 0, l_s, 17)
    # Add lines
    gmsh.model.geo.addLine(  1,  2,  1)
    gmsh.model.geo.addLine(  2,  4,  2)
    gmsh.model.geo.addLine(  4,  3,  3)
    gmsh.model.geo.addLine(  3,  1,  4)
    gmsh.model.geo.addLine(  2,  5,  5)
    gmsh.model.geo.addLine(  5,  6,  6)
    gmsh.model.geo.addLine(  6,  4,  7)
    gmsh.model.geo.addLine(  5,  7,  8)
    gmsh.model.geo.addLine(  7,  8,  9)
    gmsh.model.geo.addLine(  8,  6, 10)
    gmsh.model.geo.addLine(  4, 10, 11)
    gmsh.model.geo.addLine( 10,  9, 12)
    gmsh.model.geo.addLine(  9,  3, 13)
    gmsh.model.geo.addLine(  6, 11, 14)
    gmsh.model.geo.addLine( 11, 10, 15)
    gmsh.model.geo.addLine(  8, 12, 16)
    gmsh.model.geo.addLine( 12, 11, 17)
    gmsh.model.geo.addLine( 10, 14, 18)
    gmsh.model.geo.addLine( 14, 13, 19)
    gmsh.model.geo.addLine( 13,  9, 20)
    gmsh.model.geo.addLine( 11, 15, 21)
    gmsh.model.geo.addLine( 15, 14, 22)
    gmsh.model.geo.addLine( 12, 16, 23)
    gmsh.model.geo.addLine( 16, 15, 24)
    gmsh.model.geo.addLine( 14, 18, 25)
    gmsh.model.geo.addLine( 18, 17, 26)
    gmsh.model.geo.addLine( 17, 13, 27)
    gmsh.model.geo.addLine( 15, 19, 28)
    gmsh.model.geo.addLine( 19, 18, 29)
    gmsh.model.geo.addLine( 16, 20, 30)
    gmsh.model.geo.addLine( 20, 19, 31)
    gmsh.model.geo.addLine( 18, 22, 32)
    gmsh.model.geo.addLine( 22, 21, 33)
    gmsh.model.geo.addLine( 21, 17, 34)
    gmsh.model.geo.addLine( 19, 23, 35)
    gmsh.model.geo.addLine( 23, 22, 36)
    gmsh.model.geo.addLine( 20, 24, 37)
    gmsh.model.geo.addLine( 24, 23, 38)
    gmsh.model.geo.addLine( 22, 26, 39)
    gmsh.model.geo.addLine( 26, 25, 40)
    gmsh.model.geo.addLine( 25, 21, 41)
    gmsh.model.geo.addLine( 23, 27, 42)
    gmsh.model.geo.addLine( 27, 26, 43)
    gmsh.model.geo.addLine( 24, 28, 44)
    gmsh.model.geo.addLine( 28, 27, 45)
    #gmsh.model.geo.addLine( 29, 30, 46)
    #gmsh.model.geo.addLine( 30, 31, 47)
    #gmsh.model.geo.addLine( 31, 32, 48)
    #gmsh.model.geo.addLine( 32, 29, 49)
    # Construct curve loops and surfaces 
    gmsh.model.geo.addCurveLoop([  1,  2,  3,  4], 1)
    gmsh.model.geo.addCurveLoop([  5,  6,  7, -2], 2)
    gmsh.model.geo.addCurveLoop([  8,  9, 10, -6], 3)
    gmsh.model.geo.addCurveLoop([ 11, 12, 13, -3], 4)
    gmsh.model.geo.addCurveLoop([ -7, 14, 15,-11], 5)
    gmsh.model.geo.addCurveLoop([ 16, 17,-14,-10], 6)
    gmsh.model.geo.addCurveLoop([ 18, 19, 20,-12], 7)
    gmsh.model.geo.addCurveLoop([-15, 21, 22,-18], 8)
    gmsh.model.geo.addCurveLoop([-17, 23, 24,-21], 9)
    gmsh.model.geo.addCurveLoop([-19, 25, 26, 27],10)
    gmsh.model.geo.addCurveLoop([-22, 28, 29,-25],11)
    gmsh.model.geo.addCurveLoop([-24, 30, 31,-28],12)
    gmsh.model.geo.addCurveLoop([-26, 32, 33, 34],13)
    gmsh.model.geo.addCurveLoop([-29, 35, 36,-32],14)
    gmsh.model.geo.addCurveLoop([-31, 37, 38,-35],15)
    gmsh.model.geo.addCurveLoop([-33, 39, 40, 41],16)
    gmsh.model.geo.addCurveLoop([-36, 42, 43,-39],17)
    gmsh.model.geo.addCurveLoop([-38, 44, 45,-42],18)
    #gmsh.model.geo.addCurveLoop([ 46, 47, 48, 49],19)
    #gmsh.model.geo.addCurveLoop([29, 30], 11)
    gmsh.model.geo.addPlaneSurface([ 1],  1)
    gmsh.model.geo.addPlaneSurface([ 2],  2)
    gmsh.model.geo.addPlaneSurface([ 3],  3)
    gmsh.model.geo.addPlaneSurface([ 4],  4)
    gmsh.model.geo.addPlaneSurface([ 5],  5)
    gmsh.model.geo.addPlaneSurface([ 6],  6)
    gmsh.model.geo.addPlaneSurface([ 7],  7)
    gmsh.model.geo.addPlaneSurface([ 8],  8)
    gmsh.model.geo.addPlaneSurface([ 9],  9)
    gmsh.model.geo.addPlaneSurface([10], 10)
    gmsh.model.geo.addPlaneSurface([11], 11)
    gmsh.model.geo.addPlaneSurface([12], 12)
    gmsh.model.geo.addPlaneSurface([13], 13)
    gmsh.model.geo.addPlaneSurface([14], 14)
    gmsh.model.geo.addPlaneSurface([15], 15)
    gmsh.model.geo.addPlaneSurface([16], 16)
    gmsh.model.geo.addPlaneSurface([17], 17)
    gmsh.model.geo.addPlaneSurface([18], 18)
    #gmsh.model.geo.addPlaneSurface([19], 19)
    # Physical groups
    gmsh.model.addPhysicalGroup(0, [1,3,21,25], 1)
    gmsh.model.setPhysicalName(0, 1, "LeftNodes")
    gmsh.model.addPhysicalGroup(0, [7,8,24,28], 2)
    gmsh.model.setPhysicalName(0, 2, "RightNodes")
    gmsh.model.addPhysicalGroup(0, [2,5], 3)
    gmsh.model.setPhysicalName(0, 3, "BottomNodes")
    gmsh.model.addPhysicalGroup(0, [26,27], 4)
    gmsh.model.setPhysicalName(0, 4, "TopNodes")
    gmsh.model.addPhysicalGroup(0, [13,17], 5)
    gmsh.model.setPhysicalName(0, 5, "LeftDNodes")
    gmsh.model.addPhysicalGroup(0, [16,20], 6)
    gmsh.model.setPhysicalName(0, 6, "RightDNodes")
    gmsh.model.addPhysicalGroup(0, [14,15,18,19], 7)
    gmsh.model.setPhysicalName(0, 7, "TBDNodes")
    gmsh.model.addPhysicalGroup(0, [4,6,22,23], 8)
    gmsh.model.setPhysicalName(0, 8, "InnerNodes")
    gmsh.model.addPhysicalGroup(0, [9,10,11,12], 9)
    gmsh.model.setPhysicalName(0, 9, "SourceNodes")
    
    gmsh.model.addPhysicalGroup(1, [4,13,20,34,41], 10)
    gmsh.model.setPhysicalName(1, 10, "LeftEdge")
    gmsh.model.addPhysicalGroup(1, [9,16,23,37,44], 11)
    gmsh.model.setPhysicalName(1, 11, "RightEdge")
    gmsh.model.addPhysicalGroup(1, [1,5,8], 12)
    gmsh.model.setPhysicalName(1, 12, "BottomEdge")
    gmsh.model.addPhysicalGroup(1, [40,43,45], 13)
    gmsh.model.setPhysicalName(1, 13, "TopEdge")
    gmsh.model.addPhysicalGroup(1, [27], 14)
    gmsh.model.setPhysicalName(1, 14, "LeftDEdge")
    gmsh.model.addPhysicalGroup(1, [30], 15)
    gmsh.model.setPhysicalName(1, 15, "RightDEdge")
    gmsh.model.addPhysicalGroup(1, [19,22,24,26,29,31], 16)
    gmsh.model.setPhysicalName(1, 16, "TBDEdge")
    gmsh.model.addPhysicalGroup(1, [2,3,6,7,10,11,14,18,21,25,28,32,33,35,36,38,39,42], 17)
    gmsh.model.setPhysicalName(1, 17, "InnerEdge")
    gmsh.model.addPhysicalGroup(1, [46,47,48,49], 18)
    gmsh.model.setPhysicalName(1, 18, "TargetEdge")
    gmsh.model.addPhysicalGroup(1, [12,15,17], 19)
    gmsh.model.setPhysicalName(1, 19, "SourceEdge")
    
    gmsh.model.addPhysicalGroup(2, [10,11,12], 20)
    gmsh.model.setPhysicalName(2, 20, "Design")
    gmsh.model.addPhysicalGroup(2, [16,17,18], 21)
    gmsh.model.setPhysicalName(2, 21, "TopPML")
    gmsh.model.addPhysicalGroup(2, [1,2,3], 22)
    gmsh.model.setPhysicalName(2, 22, "BottomPML")
    gmsh.model.addPhysicalGroup(2, [14], 23)
    gmsh.model.setPhysicalName(2, 23, "Target")
    gmsh.model.addPhysicalGroup(2, [4,5,6,7,8,9], 24)
    gmsh.model.setPhysicalName(2, 24, "Incident")
    gmsh.model.addPhysicalGroup(2, [13,15], 25)
    gmsh.model.setPhysicalName(2, 25, "LRTarget")
    gmsh.model.geo.synchronize()
    gmsh.model.mesh.generate(2)
    
    # ... and save it to disk
    gmsh.write("geometry.msh")
    gmsh.finalize()
end