//-*- C++ -*-

//!@file

Mesh.Algorithm3D=4; // frontal3D algorithm
Mesh.CharacteristicLengthFromPoints=0;
//Mesh.CharacteristicLengthExtendFromBoundary = 0;
//Mesh.CharacteristicLengthFromCurvature = 1;
//Mesh.CharacteristicLengthFactor=0.5;

factor=1;

//discretization definitions
lcOutlet = 0.2/factor;//level of refinement (associated to the poins)
lcSuddenExpansion= 0.04/factor;
lcThroatFine= 0.001/factor;
nb_points_on_the_arc=20*factor; // nb of points on the throat arc
nb_points_on_the_BL=5*factor;

// domain definitions
R1= 0.6;//radius of inlet/outlet
R2= 0.2;//radius of the throat
R22= 0.2-0.02/factor;//radius inside of the throat (R2-boundary layer))
LO = 10.0;//length of the outlet cylinder

//Center of the circles
Point(0) = {0, 0, 0, lcThroatFine};

//Points on the first (larger) circle with radius R1
Point(1) = {R1, 0, 0, lcSuddenExpansion};
Point(2) = {0, R1, 0, lcSuddenExpansion};
Point(3) = {-R1, 0, 0, lcSuddenExpansion};
Point(4) = {0, -R1, 0, lcSuddenExpansion};

//First larger circle (divided in 4 arcs)
Circle(1) = {1,0,2};
Circle(2) = {2,0,3};
Circle(3) = {3,0,4};
Circle(4) = {4,0,1};


// line loop generating the larger disk
Line Loop(5) = {1,2,3,4};

Plane Surface(1) = {5}; //on line loops

//Points on the large outlet circle
Point(699) = {0, 0, LO, lcSuddenExpansion};
Point(700) = {R1, 0, LO, lcOutlet};
Point(701) = {0, R1, LO, lcOutlet};
Point(702) = {-R1, 0, LO, lcOutlet};
Point(703) = {0, -R1, LO, lcOutlet};

// arcs or the circle
Circle(710) = {700,699,701};
Circle(711) = {701,699,702};
Circle(712) = {702,699,703};
Circle(713) = {703,699,700};
//outlet surface
Line Loop(714) = {710,711,712,713};
Plane Surface(715) = {714};

//straight lines
Line(716) = {700, 1};
Line(717) = {701, 2};
Line(718) = {702, 3};
Line(719) = {703, 4};

Line Loop(720) = {719, 4, -716, -713};
Ruled Surface(721) = {720};
Line Loop(722) = {716, 1, -717, -710};
Ruled Surface(723) = {722};
Line Loop(724) = {711, 718, -2, -717};
Ruled Surface(725) = {724};
Line Loop(726) = {712, 719, -3, -718};
Ruled Surface(727) = {726};
Surface Loop(728) = {715, 723, 721, 727, 1, 725};//137,133,135
Volume(729) = {728};

Field[1]=MathEval;
Field[1].F= Sprintf("(sqrt(z*z)+1)/80*%g", factor);
Background Field = 1;


