
Mesh.RecombinationAlgorithm = 0; // 0 standard - 1 blossom
Mesh.RecombineAll = 1;

symm = DefineNumber[1, Name "Foam/Geometry/symm", 
		    Label "Symmetry (planar/axial)",
		    Choices{1="Rectangular", 2="Tubular"}];
clearance = DefineNumber[2, Min 2, Max 10, Step 1, 
			 Name "Foam/Geometry/clearance",
			 Label "Clearance [mm]"];
diam = DefineNumber[40, Name "Foam/Geometry/diam", 
		    Label "Pipe diameter [mm]"];
width = DefineNumber[5, Name "Foam/Geometry/width",
		     Label "Model thickness [mm]"];
angle = DefineNumber[5, Name "Foam/Geometry/angle",
		     Label "Model aperture angle [deg]"];

clearance = clearance*1e-3; // convert to [m]
diam = diam*1e-3;
width = width*1e-3; 
angle = angle*Pi/180.; // convert to [rad]

d = clearance;
r1 = diam/2;  // inlet/outlet pipe radius & valve inner radius
r2 = r1+0.005;  // valve outer radius
r3 = r1+0.020;  // cavity pipe radius 
l1 = 0.15;  // inlet pipe length
l2 = 0.075; // cavity pipe length
l3 = 0.15;   // outlet pipe length
l4 = 0.005; // valve thickness

// Meshsizes
lc = r1/4;
lc2 = d/3.;

// small eccentricity to allow azimutal extrusion
eps=0.0;
If(symm == 2)
eps= 1e-3; // avoid singular geometry on axis
EndIf

Point(1) =  {0.0,           eps,        0.0, lc}; 
Point(2) =  {l1+d,          eps,        0.0 ,lc}; 
Point(3) =  {l1+d,          r1,         0.0, lc2}; 
Point(4) =  {l1+d+l4,       r2,         0.0, lc2}; 
Point(5) =  {l1+d+l4,       eps,        0.0, lc}; 
Point(6) =  {l1+d+l4+l2+l3, eps,        0.0, lc}; 
Point(7) =  {l1+d+l4+l2+l3, r1,         0.0, lc}; 
Point(8) =  {l1+d+l4+l2,    r1,         0.0, lc}; 
Point(9) =  {l1+(r3-r1),    r3,         0.0, lc}; 
Point(10) = {l1,            r1,         0.0, lc2}; 
Point(11) = {0.0,           r1,         0.0, lc}; 

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

Line Loop(12) = {1, 2, 3, -4, 5, 6, -7, 8, 9, -10, 11};
Plane Surface(13) = {12};

If(symm == 1)
  Extrude {0, 0, width} {
    Surface{13};
    Layers {{1},{1}};
    Recombine;        
  }
EndIf

If(symm == 2)
  Extrude { {1,0,0}, {0,0,0} , angle } {
    Surface{13};
    Layers {{1},{1}};
    Recombine;        
  }
EndIf

// Boundary condition
Physical Surface("frontAndBack") = {70, 13};
Physical Surface("inlet") = {69};
Physical Surface("outerWall") = {65, 61, 57, 53};
Physical Surface("outlet") = {49};
Physical Surface("valve") = {37, 33, 41};
Physical Surface("symetry") = {29, 45};
Physical Volume("fluid") = {1};

