// long microchannel between two reservoirs

ll_x = 1e3; // length of microchannel : Maximum dimension of system
delx_s = 1.e-5; // smallest delx within the stretch region : Minimum prescribed mesh element size
//ll_x = 5.e2; // length of microchannel - trialling shorter domain
ll_y = 2.; // width of microchannel
ll_r_x = 1.e1; // length of reservoir region
ll_r_y = ll_r_x; // width of reservoir region


delx_c = 1.; // target largest delx within the constant region
dely_s = delx_s; // smallest dely
deltax = 1.15e0; // stretch ratio in x direction along nanochannel
deltay = 1.15e0; // stretch ratio in y direction perpendicular to nanochannel walls
deltayr = 1.04e0; // stretch ratio in y direction in reservoirs

lc_small = delx_s;  // charateristic mesh length variable - small
lc_large = 0.25;  // charateristic mesh length variable - small
lc = lc_small; // temp

//---------------------------------------------------------
// dump data about user variables
Printf ("delx_s %f: delx_c = %f: deltax = %f",delx_s,delx_c,deltax);

// calculate some stretched mesh variables
// stretch in x direction along microchannel
n_s = Floor(Log(delx_c/delx_s)/Log(deltax)+1)+1;
ll_s = (1-deltax^(n_s))/(1-deltax)*delx_s;
ll_c = ll_x-ll_s/2;
n_c = Floor(ll_c/delx_c)+1;
delx_c = ll_c/n_c;
// dump data about generated variables
Printf ("n_s = %f: n_c = %f: ll_s = %f: ll_c = %f: delx_c = %f",n_s,n_c,ll_s,ll_c,delx_c);

// stretch in y direction along microchannel
// ll_s = (1-deltax^(n_s))/(1-deltax)*delx_s;
// ll_y/2*(1-deltay)/dely_s = 1-deltay^{n_y)
// deltay^{n_y} = 1-ll_y/2*(1-deltay)/dely_s
n_y = Floor(Log(1-ll_y/2*(1-deltay)/dely_s)/Log(deltay)) + 1;
dely_s = ll_y/2*(1-deltay)/(1-deltay^(n_y));
Printf ("n_y = %f: dely_s = %f",n_y,dely_s);

// stretch in y direction along reservoir wall
n_y2 = Floor(Log(1-(ll_r_y/2-ll_y/2)*(1-deltayr)/dely_s)/Log(deltayr)) + 1;
Printf ("n_y2 = %f",n_y2);

//---------------------------------------------------------

// setup top left quadrant of microchannel
Point(1) = {0, 0, 0, lc_large/2};
Point(2) = {ll_s,  0, 0, lc} ;
Point(3) = {ll_s,  ll_y/2, 0, lc} ;
Point(4) = {0,  ll_y/2, 0, lc} ;
Line(1) = {1,2};
Line(2) = {3,2};
Line(3) = {4,3};
Line(4) = {4,1};
Line Loop(1) = {1, -2, -3, 4};
Plane Surface(1) = {1};

// setup top middle section of microchannel
Point(5) = {ll_x-ll_s,  0, 0, lc} ;
Point(6) = {ll_x-ll_s,  ll_y/2, 0, lc} ;
Line(5) = {2,5};
Line(6) = {6,5};
Line(7) = {6,3};
Line Loop(2) = {5, -6, 7, 2};
Plane Surface(2) = {2};

// setup top right quadrant section
Point(7) = {ll_x,  0, 0, lc_large/2} ;
Point(8) = {ll_x,  ll_y/2, 0, lc} ;
Line(8) = {7,5};
Line(9) = {8,7};
Line(10) = {8,6};
Line Loop(3) = {-8, -9, 10, 6};
Plane Surface(3) = {3};

// mesh microchannel with structured grid
Transfinite Line {1} = n_s Using Progression deltax; // bottom
//Transfinite Line {1} = n_s ; // bottom
Transfinite Line {2} = n_y Using Progression deltay; // right
Transfinite Line {3} = n_s Using Progression deltax; // top
Transfinite Line {4} = n_y Using Progression deltay; // left
Transfinite Surface {1} = {1,2,3,4}; // these are the points
Recombine Surface{1};
Transfinite Line {5} = n_c; // bottom
Transfinite Line {6} = n_y Using Progression deltay; // right
Transfinite Line {7} = n_c; // top
Transfinite Surface {2} = {2,5,6,3}; // these are the points
Recombine Surface{2};
Transfinite Line {8} = n_s Using Progression deltax; // bottom
Transfinite Line {9} = n_y Using Progression deltay; // right
Transfinite Line {10} = n_s Using Progression deltax; // top
Transfinite Surface {3} = {5,7,8,6}; // these are the points
Recombine Surface{3};




// reflect geometry to create right half
Symmetry {1, 0, 0, -ll_x/2} {
  Duplicata { Surface{4}; }
}




