mu = 1e-3;
rho_zeo = 1714;
rho_PU  = 1030;

p_vv=0.2;

// zeo_ww = p_ww;
// zeo_vv = zeo_ww * rho_PU / ((100 - zeo_ww) * rho_zeo);
zeo_vv = p_vv;
zeo_side = 1*mu;

short_side = 8 * mu;
long_side  = 6 * mu;
side_ratio = short_side / long_side;
mem_A = short_side * long_side;


//zeo_n      = mem_A * zeo_vv / (zeo_side^2 * (1+zeo_vv));
zeo_n      = mem_A * zeo_vv / (zeo_side^2);
zeo_x 	   = Ceil(Sqrt(zeo_n * side_ratio));
zeo_y	   = Ceil(zeo_n / zeo_x);

//real_load_vv  = (zeo_x * zeo_y)*zeo_side^2 / (mem_A - (zeo_x * zeo_y)*zeo_side^2);
real_load_vv  = (zeo_x * zeo_y)*zeo_side^2 / (mem_A);
real_load_ww  = real_load_vv*rho_zeo/rho_PU / (1 + real_load_vv*rho_zeo/rho_PU);

xint = (short_side - zeo_x * zeo_side) / (zeo_x + 1);
x0=xint;

yint = (long_side - zeo_y * zeo_side) / (zeo_y + 1);
y0=yint*1;

If (zeo_side < xint*10)
   lc = zeo_side;
EndIf
If (zeo_side >= xint*10)
   lc = zeo_side/4;
EndIf
If (lc >= 5*mu)
   lc = zeo_side;
EndIf
//lc = 1*mu;

zeo_lc=lc;
mem_lc = lc;


// Print info to file
Printf ("Zeolite side:   %5.1f mu
Membrane sides: %5.0f mu x %5.0f mu

Obtained loading(w/w):      %5.2f
Zeolite volume ratio (v/v): %5.2f
Obtained volume ratio (v/v):%5.2f
Number of Zeolites:       %5.0f
Zeolite in x and y:       %5.0f x %5.0f
Spacing in x and y:       %5.2f x %5.2f

", zeo_side/mu, short_side/mu, long_side/mu, real_load_ww*100, zeo_vv*100, real_load_vv*100, zeo_n, zeo_x, zeo_y, xint/mu, yint/mu) > "membrane.info";

prepost= 5 * lc;

Point(1) = {0,		-prepost,	0, mem_lc};
Point(2) = {short_side,	-prepost,	0, mem_lc};

Point(3) = {0, 		0, 		0, mem_lc};
Point(4) = {short_side, 0, 		0, mem_lc};
Point(5) = {short_side, long_side, 	0, mem_lc};
Point(6) = {0, 		long_side, 	0, mem_lc};

Point(7) = {short_side,	long_side+prepost,	0, mem_lc};
Point(8) = {0,		long_side+prepost,	0, mem_lc};

Line(1) = {1,2} ;
Line(2) = {2,4} ;
Line(3) = {3,4} ;
Line(4) = {3,1} ;
//Line Loop(5) = {1,2,-3,4} ;

Line(6) = {4,5} ;
Line(7) = {5,6} ;
Line(8) = {6,3} ;
//Line Loop(9) = {3,6,7,8} ;


Line(10) = {5,7} ;
Line(11) = {7,8} ;
Line(12) = {8,6} ;
Line Loop(13) = {-7,10,11,12} ;

Function CarvZeo
	p1 = newp;	Point(p1) = {x,			y,		0,  zeo_lc} ;
	p2 = newp;	Point(p2) = {x+zeo_side,	y,		0,  zeo_lc} ;
	p3 = newp;	Point(p3) = {x+zeo_side,	y+zeo_side,	0,  zeo_lc} ;
	p4 = newp;	Point(p4) = {x,  		y+zeo_side,	0,  zeo_lc} ;

	l1 = newl;	Line(l1) = {p1,p2} ;
	l2 = newl;	Line(l2) = {p2,p3} ;
	l3 = newl;	Line(l3) = {p3,p4} ;
	l4 = newl;	Line(l4) = {p4,p1} ;
	zeoliteline[t] = newll ;	Line Loop(zeoliteline[t]) = {l1,l2,l3,l4} ;
	thezeolites[t] = news ;		Plane Surface(thezeolites[t]) = {zeoliteline[t]} ;
Return

x = x0;
y = y0;
t=1;
trans=1;
yt[0]=t;
For i In {1:zeo_y}
	x = x0;
	For j In {1:zeo_x}
		Call CarvZeo ;
		x += xint+zeo_side;
		t += 1;

	EndFor
	yt[i]=t;
	yc[i] = y+zeo_side/2;
	y += yint+zeo_side;
	trans += 1;
EndFor
x=x0;
For j In {1:zeo_x}
	p1 = newp;	Point(p1) = {x,			0,		0,  zeo_lc} ;
	p2 = newp;	Point(p2) = {x+zeo_side,	0,		0,  zeo_lc} ;
	p3 = newp;	Point(p3) = {x+zeo_side+xint,	0,		0,  zeo_lc} ;
	If (j==1)
		l1 = newl;	Line(l1) = {3,p1} ;	pz1[0]=l1;
	EndIf
	l1 = newl;	Line(l1) = {p1,p2} ; 	l2 = newl;	Line(l2) = {p2,p3} ;
	pz1[j*2-1]=l1;	pz1[j*2]=l2;
	x += xint+zeo_side;
EndFor


small_lines=1;

If (small_lines==1)
	pz1[j*2+1]=-2;	pz1[j*2+2]=-1;	pz1[j*2+3]=-4;
	Line Loop(5) = {pz1[]} ;
	theinlet    = news ;
	Plane Surface(theinlet) = {5} ;
EndIf

If (small_lines==0)
	Line Loop(5) = {1,2,-3,4} ;
	theinlet    = news ;
	Plane Surface(theinlet) = {5} ;

	Line{pz1[1]} In Surface{theinlet};
EndIf



Line Loop(9) = {3,6,7,8} ;
//pz1[j*2+1]=6;	pz1[j*2+2]=7;	pz1[j*2+3]=8;
//Line Loop(9) = {pz1[]} ;

zeoliteline[0] = 9;
themembrane = news ;
Plane Surface(themembrane) = {zeoliteline[0]} ;

theoutlet   = news ;
Plane Surface(theoutlet) = {13} ;

Physical Line(1)    = {1} ;
Physical Surface(1) = {theinlet} ;
//Physical Line(2)    = {3} ;
Physical Surface(2) = {themembrane} ;
Printf ("Zeolite centers:") >> "membrane.info";

Physical Surface(3) = {thezeolites[yt[0]]:thezeolites[yt[1]-1]} ;
Physical Line(3)    = {13:zeoliteline[yt[1]-1]} ;
Printf ("%5.2f", yc[1]/mu ) >> "membrane.info";

For i In {2:zeo_y}
	Physical Surface(i+2) = {thezeolites[yt[i-1]]:thezeolites[yt[i]-1]} ;
	Physical Line(i+2)    = {zeoliteline[yt[i-1]-1]:zeoliteline[yt[i]-1]} ;
//	Physical Line(i+2)    = {zeoliteline[yt[i-1]-1]:zeoliteline[yt[i]-1]} ;
	Printf ("%5.2f", yc[i]/mu ) >> "membrane.info";
EndFor
Physical Line(zeo_y+3)    = {7} ;
Physical Surface(zeo_y+3) = {theoutlet} ;
Physical Line(zeo_y+4)    = {11} ;

Physical Line(zeo_y+5)    = {pz1[]} ;

// 1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad
Mesh.Algorithm = 6;
Mesh.CharacteristicLengthFactor = 1.0;
//Mesh.Optimize = 1;
//Mesh.OptimizeNetgen = 1;

//Recombine Surface {themembrane};


