// Number of band
N = 3 ;
// x and y size
H = 2.0 ;
dH = 0.5 ;
L = 1.0 ;
delta = L/N ;
// mesh element lenght
lc = delta ;
// mesh element lenght  at the surface
lcs = lc ;



Function bedrock
bed = L + dH*Sin(2*Pi*xx/L)*Sin(2*Pi*yy/L) ;
Return

Function surface 
Surf = H ; 
Return

// mesh the bed surface 6
y= 0 ;
x = 0 ;    
Nsurf0 = 15*N*N + 1 ;
Nsurf = Nsurf0 ;
For j In {1:N} 
For i In {1:N} 
//
//     p4 ---- l3 ---- p3
//     | \            / |
//     |  \          /  |
//     |   l8      l7   |
//     |     \    /     |
//     l4      p5       l2 
//     |     /    \     |
//     |   l5      l6   |
//     |  /          \  |
//     |/             \ |
//     p1 ---- l1 ---- p2
//
// Define 5 points 
xx = x ; yy = y ; Call bedrock ;
p1 = newp ; Point(p1) = {xx,yy,bed,lc} ;
xx = x + delta ; yy = y ; Call bedrock ;
p2 = newp ; Point(p2) = {xx,yy,bed,lc} ;
xx = x + delta ; yy = y + delta ; Call bedrock ;
p3 = newp ; Point(p3) = {xx,yy,bed,lc} ;
xx = x ; yy = y + delta ; Call bedrock ;
p4 = newp ; Point(p4) = {xx,yy,bed,lc} ;
xx = x + 0.5*delta ; yy = y + 0.5*delta ; Call bedrock ;
p5 = newp ; Point(p5) = {xx,yy,bed,lc} ;

// Define 8 lines 
l1 = newc ; Line(l1) = {p1,p2} ; 
l2 = newc ; Line(l2) = {p2,p3} ; 
l3 = newc ; Line(l3) = {p3,p4} ; 
l4 = newc ; Line(l4) = {p4,p1} ; 
l5 = newc ; Line(l5) = {p1,p5} ; 
l6 = newc ; Line(l6) = {p2,p5} ; 
l7 = newc ; Line(l7) = {p3,p5} ; 
l8 = newc ; Line(l8) = {p4,p5} ; 

// Define 4 Line Loop and 4 Ruled Surface 
ll1 = newc; Line Loop(ll1) = {l1,l6,-l5};   Ruled Surface(Nsurf) = {ll1};
Nsurf += 1 ;
ll2 = newc; Line Loop(ll2) = {l2,l7,-l6};   Ruled  Surface(Nsurf) = {ll2};
Nsurf += 1 ;
ll3 = newc; Line Loop(ll3) = {l3,l8,-l7};   Ruled  Surface(Nsurf) = {ll3};
Nsurf += 1 ;
ll4 = newc; Line Loop(ll4) = {l4,l5,-l8};   Ruled Surface(Nsurf) = {ll4};
Nsurf += 1 ;

x +=  delta ; 
EndFor
x=0.0 ;
y +=  delta ; 
EndFor

Coherence ;
Physical Surface(6) = {Nsurf0:Nsurf-1} ;

// mesh the 5 other surfaces 
xx = 0.0 ; yy = 0.0 ; Call bedrock ;
p1 = newp ; Point(p1) = {xx,yy,bed,lc} ;
xx = L ; yy = 0.0 ; Call bedrock ;
p2 = newp ; Point(p2) = {xx,yy,bed,lc} ;
xx = L ; yy = L ; Call bedrock ;
p3 = newp ; Point(p3) = {xx,yy,bed,lc} ;
xx = 0.0 ; yy = L ; Call bedrock ;
p4 = newp ; Point(p4) = {xx,yy,bed,lc} ;

xx = 0.0 ; yy = 0.0 ; Call surface ;
p5 = newp ; Point(p5) = {xx,yy,Surf,lcs} ;
xx = L ; yy = 0.0 ; Call surface ;
p6 = newp ; Point(p6) = {xx,yy,Surf,lcs} ;
xx = L ; yy = L ; Call surface ;
p7 = newp ; Point(p7) = {xx,yy,Surf,lcs} ;
xx = 0.0 ; yy = L ; Call surface ;
p8 = newp ; Point(p8) = {xx,yy,Surf,lcs} ;

Coherence;
Line(236) = {34,53};
Line(237) = {1,50};
Line(238) = {12,51};
Line(239) = {43,52};
Line(240) = {51,52};
Line(241) = {52,53};
Line(242) = {53,50};
Line(243) = {50,51};
Line Loop(244) = {203,167,4,237,-242,-236};
Plane Surface(245) = {244};
Line Loop(246) = {236,-241,-239,226,214,202};
Plane Surface(247) = {246};
Line Loop(248) = {239,-240,-238,153,189,225};
Plane Surface(249) = {248};
Line Loop(250) = {243,-238,-152,-140,-1,237};
Plane Surface(251) = {250};
Line Loop(252) = {242,243,240,241};
Plane Surface(253) = {252};
Surface Loop(254) = {253,245,-163,-160,-150,-151,-148,-138,-139,-136,-251,249,247,-170,-171,-165,-166,-167,-161,-162,-164,-154,-155,-149,-152,-142,-143,-137,-140,-141,-147,-144,-145,-146,-156,-157,-158,-168,-169,-159,-153};
Volume(255) = {254};
