// Keep the ratio of lc1/lc2=constant. Just for simple comparison base.

// Define characteristic length for rectangle
lc1=2;
// Define characteristic length for circle in 5 parts
lc2=0.1;
// Define characteristic length for outer fictious circle 11-12-13-14-15
lc3=0.6;

// Define the radius of the circle
r=1;

//first define the rectangle, including upstream (100r) and downstream(100r), width 100r
// 20 points
//  4------------19--18---17---------------3   --
//  |             |   |   |                |    |
//  |             |  13   |                |    |
//  |             |/  | \ |                |   100r
//  20-----------14   8   12---------------16   |
//  |           /   9   7   \              |    |
//  1----------15-10  5  6--11-------------2    --
//  
//  |<-----100r------>|<------100r-------->|

Point(1)={-100*r,0,-0.5,lc1};
Point(2)={100*r,0,-0.5,lc1};
Point(3)={100*r,100*r,-0.5,lc1};
Point(4)={-100*r,100*r,-0.5,lc1};

// Use a line loop to generate the isolated fiber in the rectangle:
// Diameter 1 μm

// fiber with its center after extrusion is (0 0 0)
p0=5;Point(p0) = {0,0,-0.5,lc2}; 
p1=6;Point(p1) = {r,0,-0.5,lc2};
p2=7;Point(p2) = {0.707107*r,0.707107*r,-0.5,lc2};
p3=8;Point(p3) = {0,r,-0.5,lc2};
p4=9;Point(p4) = {-0.707107*r,0.707107*r,-0.5,lc2};
p5=10;Point(p5) = {-r,0,-0.5,lc2};

// concentric circle
p1=11;Point(p1) = {6*r,0,-0.5,lc3};
p2=12;Point(p2) = {3*1.41421*r,3*1.41421*r,-0.5,lc3};
p3=13;Point(p3) = {0,6*r,-0.5,lc3};
p4=14;Point(p4) = {-3*1.41421*r,3*1.41421*r,-0.5,lc3};
p5=15;Point(p5) = {-6*r,0,-0.5,lc3};

// other points on the rectangular side, i.e. point 16, 17, 18, 19, 20;
Point(16)={100*r,3*1.41421*r,-0.5,lc1};
Point(17)={3*1.41421*r,100*r,-0.5,lc1};
Point(18)={0,100*r,-0.5,lc1};
Point(19)={-3*1.41421*r,100*r,-0.5,lc1};
Point(20)={-100*r,3*1.41421*r,-0.5,lc1};

// the lines between the inner fiber, outer circle, and rectangular points
// straight lines around the rectangle
Line(1) = {1,15};
Line(2) = {11,2};
Line(3) = {2,16};
Line(4) = {16,3};
Line(5) = {3,17};
Line(6) = {17,18};
Line(7) = {18,19};
Line(8) = {19,4};
Line(9) = {4,20};
Line(10) = {20,1};
// straight lines between outer circle and rectangular points
Line(11) = {12,16};
Line(12) = {12,17};
Line(13) = {13,18};
Line(14) = {14,19};
Line(15) = {14,20};
// the arc on outer circle
Circle(16) = {11,5,12};
Circle(17) = {12,5,13};
Circle(18) = {13,5,14};
Circle(19) = {14,5,15};
// straight lines between inner fiber and outer circle
Line(20) = {6,11};
Line(21) = {7,12};
Line(22) = {8,13};
Line(23) = {9,14};
Line(24) = {10,15};
// the arc on inner fiber
Circle(25) = {6,5,7};
Circle(26) = {7,5,8};
Circle(27) = {8,5,9};
Circle(28) = {9,5,10};

// the faces between the inner fiber, outer circle, and rectangular points
// the faces between inner fiber and outer circle
Line Loop(1) = {20,16,-21,-25};  //point 6-11-12-7; line 20,16,-21,-25
Plane Surface(1) = {1};
Line Loop(2) = {21,17,-22,-26};  //point 7-12-13-8; line 21,17,-22,-26
Plane Surface(2) = {2};
Line Loop(3) = {22,18,-23,-27};  //point 8-13-14-9; line 22,18,-23,-27
Plane Surface(3) = {3};
Line Loop(4) = {23,19,-24,-28};  //point 9-14-15-10;line 23,19,-24,-28
Plane Surface(4) = {4};
// the faces between outer circle and rectangule
Line Loop(5) = {2,3,-11,-16};     //point 11-2-16-12; line 2,3,-11,-16
Plane Surface(5) = {5};
Line Loop(6) = {11,4,5,-12};     //point 12-16-3-17; line 11,4,5,-12
Plane Surface(6) = {6};
Line Loop(7) = {12,6,-13,-17};   //point 12-17-18-13;line 12,6,-13,-17
Plane Surface(7) = {7};
Line Loop(8) = {13,7,-14,-18};   //point 13-18-19-14;line 13,7,-14,-18
Plane Surface(8) = {8};
Line Loop(9) = {14,8,9,-15};     //point 14-19-4-20; line 14,8,9,-15
Plane Surface(9) = {9};
Line Loop(10) = {15,10,1,-19};   //point 14-20-1-15; line 15,10,1,-19
Plane Surface(10) = {10};

// Extrude both the geometry and mesh
// Note:this will influence the Aref for the drag coefficients.
//h=1; //Max aspect ratio = 36.0988
h=0.1*r; //Max aspect ratio = 13.9788
//h=0.5*r; //Max aspect ratio = 18.0494

e1[] =Extrude{0,0,h}{Surface{1:10};Layers{1};Recombine;};
//e1[0-5]  --surface1
//e1[6-11] --surface2
//e1[12-17]--surface3
//e1[18-23]--surface4
//e1[24-29]--surface5
//e1[30-35]--surface6
//e1[36-41]--surface7
//e1[42-47]--surface8
//e1[48-53]--surface9
//e1[54-59]--surface10

//use surface1 as an example for illustration:
//e1[0]--the opposite face of Surface(1)
//e1[1]--the whole volume
//e1[2-5]--surfaces from bottom, right side, upper side, and left side
//according to the sequence of the nodes during the face construction
//where e1[5] is the 1/4 of fiber boundary


Physical Surface ("inlet") = {e1[52],e1[57]};

Physical Surface ("outlet") = {e1[27],e1[33]};

Physical Surface ("symmetry_up") = {e1[34],e1[39],e1[45],e1[51]};

Physical Surface ("symmetry_down") = {e1[58],e1[22],e1[2],e1[26]};

Physical Surface("fiber")={e1[5],e1[11],e1[17],e1[23]};

Physical Surface("front")={1:10};

Physical Surface("back")={e1[0],e1[6],e1[12],e1[18],e1[24],e1[30],e1[36],e1[42],e1[48],e1[54]};

Physical Volume ("rectangleFlowField") = {e1[1],e1[7],e1[13],e1[19],e1[25],e1[31],e1[37],e1[43],e1[49],e1[55]};

// At last, set the unit micrometer---Warn: this will make your geometry invisable in paraFoam
Mesh.ScalingFactor = 1e-06;


