
/********************************************************************************/
/* ** Choose general settings ** */

Mesh.CharacteristicLengthMin = 0.01;
Mesh.CharacteristicLengthMax = 0.01;
//Mesh.Algorithm3D = 4;
Mesh.Smoothing = 100;
Mesh.Optimize=1;
Mesh.OptimizeNetgen=1;
//Geometry.Tolerance = 1e-3; // adjust value here for correct merge result

Geometry.NumSubEdges = 100; // nicer display of curve
General.Color.Text = Black;
Geometry.Color.Points = Purple;
General.Color.Background = White;
General.Color.Text = Black;
Mesh.Color.Points = {255,0,0};


/* ** Define constants ** */

DefineConstant[
r = {0.38, Name "Bend radius"}
rc = {0.125, Name "Radius of the cross-section"}
N = {20, Name "Number of viscous layers"}
Nv1 = {28, Name "Number of layers in extruded section v1"}
Nv2 = {80, Name "Number of layers in extruded section v2"}
Ne = {40, Name "Number of layers in the elbows"}
Nh = {60, Name "Number of layers in horizontal section"}
lw ={0.01, Name "Thickness of the boundary layer"}
mes = {1e-3, Name "Mesh element size"}
lh = {1.49, Name "Length of the horizontal section"}
lv1 = {0.74, Name "Length of the 1st vertical section"}
lv2 = {2.07, Name "Length of the 2nd vertical section"}];

/********************************************************************************/

/* ** Create the 1st elbow ** */

// Points in the elbow
rd = rc - lw; // Radius of the inner mesh = radius - boundary layer

np1 = newp; Point(np1) = {r, 0, 0, mes}; //middle
np2 = newp; Point(np2) = {(r - rd), 0, 0, mes};
np3 = newp; Point(np3) = {r, 0, rd, mes};
np4 = newp; Point(np4) = {(r+rd), 0, 0, mes};
np5 = newp; Point(np5) = {r, 0, -rd , mes};

// Circles in the elbow
// elbow 1
nl1 = newl; Circle(nl1) = {2, 1, 3};
nl2 = newl; Circle(nl2) = {3, 1, 4};

nl3 = newl; Circle(nl3) = {4, 1, 5};
nl4 = newl; Circle(nl4) = {5, 1, 2};

nll1 = newll; Line Loop(nll1) = {nl1,nl2,nl3,nl4};

ns1 = news; Surface(ns1) = {nll1};


elbow[] = Extrude { {0,0,1}, // direction of rotation axis
                     {0,0,0}, // a point on the rotation axis
                     -Pi/2 } { // the rotation angle
    Surface{ns1}; Layers{Ne}; Recombine;
};


// Boundary layer

boundary_layer_e1[] = Extrude {
Surface{-elbow[2]:-elbow[5]}; Layers{N, lw}; Recombine;};

/* ** Create the horizontal section ** */

// Create the element by extrusion

horizontal[] = Extrude {-lh, 0, 0} {Surface{elbow[0]}; Layers{Nh}; Recombine;};

// Boundary layer

boundary_layer_h[] = Extrude {Surface{-horizontal[2]:-horizontal[5]}; Layers{N, lw}; Recombine;};

/* ** Create the 2nd elbow ** */

// Create the element by extrusion

elbow2[] = Extrude { {0,0,1}, // direction of rotation axis
                     {-lh,0,0}, // a point on the rotation axis
                     -Pi/2 } { // the rotation angle
    Surface{horizontal[0]}; Layers{Ne}; Recombine;
};

// Boundary layer

boundary_layer_e2[] = Extrude {
Surface{-elbow2[2]:-elbow2[5]}; Layers{N, lw}; Recombine;};

/* ** Create the outflow section ** */

// Create the element by extrusion

outflow[] = Extrude {0, lv2, 0} {Surface{ns1}; Layers{Nv2}; Recombine;};

// Boundary layer

boundary_layer_o[] = Extrude {
Surface{outflow[2]:outflow[5]}; Layers{N, lw}; Recombine;};

/* ** Create the inflow section ** */

// Create the element by extrusion

inflow[] = Extrude {0, lv1, 0} {Surface{elbow2[0]}; Layers{Nv1}; Recombine;};

// Boundary layer

boundary_layer_i[] = Extrude {
Surface{-inflow[2]:-inflow[5]}; Layers{N, lw}; Recombine;};


/* ** Assign Physical entities ** */
Physical Volume("PIPE") = {inflow[1], horizontal[1], outflow[1], 
elbow[1], elbow2[1]}; 
Physical Surface("INLET") = {inflow[0]};
Physical Surface("OUTLET") = {outflow[0]};
Physical Surface("WALLS") = {inflow[2], inflow[3], inflow[4], inflow[5], 
outflow[2], outflow[3], outflow[4], outflow[5],
horizontal[2], horizontal[3], horizontal[4], horizontal[5],
elbow[2], elbow[3], elbow2[4], elbow2[5],
elbow2[2], elbow2[3], elbow2[4], elbow2[5]
};
Physical Surface("BOUNDARY LAYER") = {boundary_layer_i[], boundary_layer_o[], boundary_layer_e1[], boundary_layer_e2[], boundary_layer_h[]};
/********************************************************************************/
General.ExpertMode = 1; // disable warning about Delaunay

