// Define geometry of the plane
L = 20;
W = 10;
Lg = 10;
Ls = 5;
Ld = 5;

// Define the unit cell coordiante for Graphene
acc = 0.142; // nm
// define the coordiantes of the basis points
xpb1 = 0.0;
ypb1 = acc*Sin(Pi/3);
zpb1 = 0.0;

xpb2 = acc*Cos(Pi/3.0);
ypb2 = 0.0;
zpb2 = 0.0;

xpb3 = acc + acc*Cos(Pi/3.0);
ypb3 = 0.0;
zpb3 = 0.0;

xpb4 = acc+acc*Cos(Pi/3.0)+acc*Cos(Pi/3.0);
ypb4 = acc*Sin(Pi/3.0);
zpb4 = 0.0;

// size of the unit cell
delX = 2.0*(acc+acc*Cos(Pi/3.0));
delY = 2.0*acc*Sin(Pi/3.0);

// Get the Number of unit cells in x and Y direction
Nx = Ceil((L-0.5)/delX);
Ny = Ceil((W-0.5)/delY);
// 
originx = (0.5)/2.0;
originy =  (0.5)/2.0;
// Define the boundary
lcg = acc;
Point(1) = {0, 0, 0,lcg};
Point(2) = {Ls,  0, 0, lcg};
Point(3) = {L-Ld, 0, 0, lcg};
Point(4) = {L,  0, 0,lcg};

Point(5) = {L,  W, 0,lcg};
Point(6) = {L-Ld,  W, 0,lcg};
Point(7) = {Ls,  W, 0, lcg};
Point(8) = {0,  W, 0,lcg};

Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 5};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 1};

Line Loop(24)={1,2,3,4,5,6,7,8};
Plane Surface(34) = {24};
// Embed attom positions in the mesh
For ic In {1:Nx}
For jc In {1:Ny}
 xp1 = xpb1 + (ic-1)*delX + originx;
 yp1 = ypb1 + (jc-1)*delY + originy;
 zp1 = zpb1;
 p1=newp;Point(p1)={xp1,yp1,0,lcg};
 
 xp2 = xpb2 + (ic-1)*delX + originx;
 yp2 = ypb2 + (jc-1)*delY + originy;
 zp2 = zpb2;
 p2=newp;Point(p2)={xp2,yp2,0,lcg};
 
 xp3 = xpb3 + (ic-1)*delX + originx;
 yp3 = ypb3 + (jc-1)*delY + originy;
 zp3 = zpb3;
 p3=newp;Point(p3)={xp3,yp3,0,lcg};
 
 xp4 = xpb4 + (ic-1)*delX + originx;
 yp4 = ypb4 + (jc-1)*delY + originy;
 zp4 = zpb4;
 p4=newp;Point(p4)={xp4,yp4,0,lcg};
 
Point {p1,p1,p3,p4} In Surface{34}; 
EndFor
EndFor

