/* 
First try to represent a channel with solid microstructures (beams) on the lower-level plate
Author: Thibaut Francis Bourlet
*/

// 1. Geometry
// 1.1. Define the size of the channel
delta = 0.5; 		lcar_b = 0.1; 	//charac length of a beam
dx = 0.2; 		dz = 0.2;      //spacing between beams
Lx = 2*Pi*delta;	lx = 0.02;
Ly = 2*delta;		ly = 0.1;
Lz = 2*Pi*delta;	lz = 0.02;


// 1.2. Define the corners of the channel
Point(1) = {0, 0, 0};    Point(2) = {0, 0, Lz};
Point(3) = {Lx, 0, Lz};  Point(4) = {Lx, 0, 0};
Point(5) = {0, Ly, 0};   Point(6) = {0, Ly, Lz};
Point(7) = {Lx, Ly, Lz}; Point(8) = {Lx, Ly, 0};

// 1.3. Define the lines
Line(1) = {1,2};  Line(2) = {2,3};  Line(3) = {3,4};
Line(4) = {4,1};  Line(5) = {5,6};  Line(6) = {6,7};
Line(7) = {7,8};  Line(8) = {8,5};  Line(9) = {1,5};
Line(10) = {2,6}; Line(11) = {3,7}; Line(12) = {4,8};

// 1.4. Define the surfaces
Line Loop(13) = {1, 2, 3, 4};     Plane Surface(14) = {13};
Line Loop(15) = {5, 6, 7, 8};     Plane Surface(16) = {15};
Line Loop(17) = {2, 11, -6, -10}; Plane Surface(18) = {17};
Line Loop(19) = {-4, 12, 8, -9};  Plane Surface(20) = {19};
Line Loop(21) = {1, 10, -5, -9};  Plane Surface(22) = {21};
Line Loop(23) = {-3, 11, 7, -12}; Plane Surface(24) = {23};


// 2. Beam coating
// 2.1. Function 
Function Beam

	p1 = newp; Point(p1) = {x, y, z, lcar_b};
	p2 = newp; Point(p2) = {x, y, z+lz, lcar_b};
	p3 = newp; Point(p3) = {x+lx, y, z+lz, lcar_b};
	p4 = newp; Point(p4) = {x+lx, y, z, lcar_b};
	p5 = newp; Point(p5) = {x, y+ly, z, lcar_b};
	p6 = newp; Point(p6) = {x, y+ly, z+lz, lcar_b};
	p7 = newp; Point(p7) = {x+lx, y+ly, z+lz, lcar_b};
	p8 = newp; Point(p8) = {x+lx, y+ly, z, lcar_b};

	l1 =  newreg; Line(l1) =  {p1, p2};
	l2 =  newreg; Line(l2) =  {p2, p3};
	l3 =  newreg; Line(l3) =  {p3, p4};
	l4 =  newreg; Line(l4) =  {p4, p1};
	l5 =  newreg; Line(l5) =  {p5, p6};
	l6 =  newreg; Line(l6) =  {p6, p7};
	l7 =  newreg; Line(l7) =  {p7, p8};
	l8 =  newreg; Line(l8) =  {p8, p5};
	l9 =  newreg; Line(l9) =  {p1, p5};
	l10 = newreg; Line(l10) = {p2, p6};
	l11 = newreg; Line(l11) = {p3, p7};
	l12 = newreg; Line(l12) = {p4, p8};

	ll1 = newreg; Line Loop(ll1) = {l1, l2, l3, l4};
	ll2 = newreg; Line Loop(ll2) = {l5, l6, l7, l8};
	ll3 = newreg; Line Loop(ll3) = {l1, l10, -l5, -l9};
	ll4 = newreg; Line Loop(ll4) = {l2, l11, -l6, -l10};
	ll5 = newreg; Line Loop(ll5) = {-l3, l11, l7, -l12};
	ll6 = newreg; Line Loop(ll6) = {-l4, l12, l8, -l9};

	Plane Surface(news) =  {ll1}; //is saved later
	Plane Surface(news) =  {ll2};
	Plane Surface(news) =  {ll3};
	Plane Surface(news) =  {ll4};
	Plane Surface(news) =  {ll5};
	Plane Surface(news) =  {ll6};

// Store the surface loops identification numbers in a list for later reference
	theloops[t] = newreg;
	Surface Loop(theloops[t]) = {ll6+1, ll6+2,  ll6+3,  ll6+4,  ll6+5,  ll6+6}; 

	thebeam = newreg;
	Volume(thebeam) = theloops[t]; 

	cross_surface[t-1] = newreg;
	Plane Surface(cross_surface[t-1]) = ll1; //touches the plate

Return

// 2.2. Generation of the beam coating

x = 0.0; y = 0; z = 0; t=1;
max_x = Lx/dx; 	
max_z = Lz/dz;

For i In {1:max_x} //double For: row and columns

	x += dx;

	For j In {1:max_z}
 
		z += dz;
		Call Beam;
		Physical Volume(t) = thebeam;
		t += 1; 

	EndFor

	z = 0;
EndFor

//Merge cross_surface and #14
merged_s = newreg;
Physical Surface(merged_s) = {14, cross_surface[]};

theloops[0] = newreg;
Surface Loop(theloops[0]) = {22, 18, 16, 20, merged_s, 24};

Coherence;
channel_v = newreg;
Volume(channel_v) = {theloops[]};

channel_beamed = newreg;
Physical Volume(channel_beamed) = {channel_v}; 


