
lc = 1e2;

Point(1) = {76, 0, 0, lc};
Point(2) = {304, 560,  0, lc} ;
Point(3) = {304, -940+580, 0, lc} ;
Point(4) = {609, 810, 0, lc} ;
Point(5) = {609, -1330+790, 0, lc} ;
Point(6) = {914, 1040, 0, lc} ;
Point(7) = {914, -1630+1040, 0, lc} ;
Point(8)={1219,1200,0,lc};
Point(9)={1219,-1790+1210,0,lc};
Point(10)={1524,1250,0,lc};
Point(11)={1524,-1810+1250,0,lc};
Point(12)={1829,1280,0,lc};
Point(13)={1829,-1830+1300,0,lc};
Point(14)={2134,1225,0,lc};
Point(15)={2134,-1820+1225,0,lc};
Point(16)={2438,1150,0,lc};
Point(17)={2438,-1840+1130,0,lc};
Point(18)={2743,1060,0,lc};
Point(19)={2743,-1790+1050,0,lc};
Point(20)={3048,950,0,lc};
Point(21)={3048,-1700+950,0,lc};
Point(22)={3353,820,0,lc};
Point(23)={3353,-1410+820,0,lc};
// 11'5" Edge of weir
Point(24)={3479,770,0,lc};
Point(25)={3479,-1260+770,0,lc};
Point(26)={3657,810,0,lc};
Point(27)={3657,-1170+810,0,lc};
// 12'6" Other edge of weir
Point(28)={3810,0,0,lc};

// Back of weir.Pt24+654,248; & Pt28 Adjusted :-)
Point(29)={3479+604,770+100,0,lc};
Point(30)={3810+500,0+248,0,lc};

//Try a whopping Spline...
//Spline(30)={1,2,4,6,8,10,12,14,16,18,20,22,24,26,29,30,28,27,25,23,21,19,17,15,13,11,9,7,5,3,1};
Line(30)={1,2,4,6,8,10,12,14,16,18,20,22,24,26,29,30,28,27,25,23,21,19,17,15,13,11,9,7,5,3,1};

Line Loop(31) = {30};
Plane Surface(32) = {31};

// Duplicate every single point, at the correct depth.....
depth[]={340,430,460,500,520,540,550,600,630,620,540,200};
For i In {0:11}
  Translate{0,-200,-depth[i]}{Duplicata{Point{2*i+2};}}
  p=newp-1; nP[] += p;
  Translate{0,200,-depth[i]}{Duplicata{Point{2*i+3};}}
  p=newp-1; nP[] += p;
EndFor  

// For each PAIR, we define three planes.....
For i In {0:10}
  l1=newl; Line(l1)={2*i+2,2*(i+1)+2,nP[2*i+2],nP[2*i],2*i+2};
  ll1=newll; Line Loop(ll1)={l1};
  s1=news; Plane Surface(s1)={ll1};
  nS[]+=s1;
  l2=newl; Line(l2)={nP[2*i],nP[2*i+2],nP[2*i+3],nP[2*i+1],nP[2*i]};
  ll2=newll; Line Loop(ll2)={l2};
  s2=news; Plane Surface(s2)={ll2};
  nS[]+=s2;
  l3=newl; Line(l3)={2*i+3,2*(i+1)+3,nP[2*i+3],nP[2*i+1],2*i+3};
  ll3=newll; Line Loop(ll3)={l3};
  s3=news; Plane Surface(s3)={ll3};
  nS[]+=s3;
EndFor

// Triangular bit to start...
  Translate{20,0,-340}{Duplicata{Point{1};}}
  nP1=newp-1;
  ln=newl; Line(ln)={1,nP1,nP[0],2,1}; ll=newll; Line Loop(ll)={ln};
  sn=news; Plane Surface(sn)={ll}; nS[]+=sn;
  ln=newl; Line(ln)={1,nP1,nP[1],3,1}; ll=newll; Line Loop(ll)={ln};
  sn=news; Plane Surface(sn)={ll}; nS[]+=sn;
  ln=newl; Line(ln)={nP1,nP[0],nP[1],nP1}; ll=newll; Line Loop(ll)={ln};
  sn=news; Plane Surface(sn)={ll}; nS[]+=sn;

// Rectangular weir
  Translate{100,-200,-350}{Duplicata{Point{24};}}
  nW[]+=newp-1;
  Translate{-100,-200,-450}{Duplicata{Point{29};}}
  nW[]+=newp-1;
  Translate{-250,150,-450}{Duplicata{Point{30};}}
  nW[]+=newp-1;
  Translate{-100,200,-350}{Duplicata{Point{28};}}
  nW[]+=newp-1;
  ln=newl; Line(ln)={nW[0],nW[1],nW[2],nW[3],nW[0]}; ll=newll; Line Loop(ll)={ln};
  sn=news; Plane Surface(sn)={ll}; nS[]+=sn;
  ln=newl; Line(ln)={nW[0],nW[1],29,24,nW[0]}; ll=newll; Line Loop(ll)={ln};
  sn=news; Plane Surface(sn)={ll}; nS[]+=sn;
  ln=newl; Line(ln)={nW[1],nW[2],30,29,nW[1]}; ll=newll; Line Loop(ll)={ln};
  sn=news; Plane Surface(sn)={ll}; nS[]+=sn;
  ln=newl; Line(ln)={nW[2],nW[3],28,30,nW[2]}; ll=newll; Line Loop(ll)={ln};
  sn=news; Plane Surface(sn)={ll}; nS[]+=sn;

// Force a Patch. Tacky...... Numbers from GUI, Yuck...
  ln=newl; Line(ln)={53,56,24,53}; ll=newll; Line Loop(ll)={ln};
  sn=news; Plane Surface(sn)={ll}; nS[]+=sn;
  ln=newl; Line(ln)={59,56,53,54,59}; ll=newll; Line Loop(ll)={ln};
  sn=news; Plane Surface(sn)={ll}; nS[]+=sn;
  ln=newl; Line(ln)={25,54,59,28,27,25}; ll=newll; Line Loop(ll)={ln};
  sn=news; Plane Surface(sn)={ll}; nS[]+=sn;


Surface Loop(1)=nS;

Volume(1)={1};


Mesh 1 ;
Mesh 2 ;                                     // Mesh the surface of the object
Plugin(NewView).View = 0 ;                   // ID of the view = 0
Plugin(NewView).Dimension = 1 ;              // Dimension of the view = 1 (scalar values)
Plugin(NewView).Run ;                        // Create view 0, it has values on surface elements
Plugin(ModifyComponents).View = 0 ;          // Select view 0 for modification
Plugin(ModifyComponents).Expression0 = "0.000001" ; // Set the value "1" for each surface element 
Plugin(ModifyComponents).Run ;               // Apply the new value
Plugin(Integrate).View = 0 ;                 // Select view 0
Plugin(Integrate).Dimension = 2 ;            // Ask to integrate on surface elements
Plugin(Integrate).Run ;                      // Run the integration (it generates a new view with ID = 1)
Printf("Area = %g", View[1].Max) ;           // Print out the area


Mesh 3 ;                                     // Mesh the volume of the object
Plugin(NewView).View = 2 ;                   // ID of new view = 2
Plugin(NewView).Dimension = 1 ;              // Dimension of new view still 1 (scalar field)
Plugin(NewView).Run ;                        // Create view 2, this one has values on volume elements
Plugin(ModifyComponents).View = 2 ;          // Select view 2 for modification
Plugin(ModifyComponents).Expression0 = "1e-6" ; // Set the value "1" to all value elements
Plugin(ModifyComponents).Run ;               // Apply the new value
Plugin(Integrate).View = 2 ;                 // Select view 2
Plugin(Integrate).Dimension = 3 ;            // Set up a volume integration
Plugin(Integrate).Run ;                      // Run the integration (generating view 3)
Printf("Volume = %g", View[3].Max) ;         // Print out the volume

