// Gmsh project created on Mon Sep 24 10:45:41 2012

d = 5;		// mesh 

//******************************************************************************
Function Quantum

// Trapezio da ruotare
p1 = newp;
Point(p1) = {x, y, z, d};
p2 = newp;
Point(p2) = {x+10, y, z, d};
p3 = newp;
Point(p3) = {x+5, y+10, z,d};
p4 = newp;
Point(p4) = {x, y+10, z, d};

l1 = newl;
Line(l1) = {p1,p2};
l2 = newl;
Line(l2) = {p2,p3};
l3 = newl;
Line(l3) = {p3,p4};
l4 = newl;
Line(l4) = {p4,p1};

// Costruiamo il tronco di cono
l5 = newl;
Line Loop(l5) = {l1,l2,l3,l4};
s1 = news;
Plane Surface(s1) = {l5};

t1[] = Extrude { {0,1,0}, {x,y,z}, 2*Pi/3 }	
	{ Surface{s1}; Recombine;  };
Vol_qt[s] = t1[1];
s += 1;	

lay[e] = t1[3];
e += 1;
lay[e] = t1[4];
e += 1;	
							
t2[] = Extrude { {0,1,0}, {x,y,z}, 2*Pi/3 }	
	{ Surface{t1[0]}; Recombine; };	
Vol_qt [s] = t2[1];
s += 1;	

lay[e] = t2[3];
e += 1;
lay[e] = t2[4];
e += 1;
						
t3[] = Extrude { {0,1,0}, {x,y,z}, 2*Pi/3 }	
	{ Surface{t2[0]}; Recombine; };	
Vol_qt[s] = t3[1];
s += 1;

lay[e] = t3[4];
e += 1;
lay[e] = t3[5];
e += 1;
	

Return		// Quantum


//***********************************************************************************************

Function Base

p5 = newp;
Point(p5) = {x-50,y,z+350,d};	// 17

p6 = newp;
Point(p6) = {x+350,y,z+350,d};	// 18

p7 = newp;
Point(p7) = {x+350,y,z-50,d};	// 19

p8 = newp;
Point(p8) = {x-50,y,z-50,d};	// 20

l5 = newl;
Line(l5) = {p5,p6};	// 56

l6 = newl;
Line(l6) = {p6,p7};	// 57

l7 = newl;
Line(l7) = {p7,p8};	// 58

l8 = newl;
Line(l8) = {p8,p5};	// 59

piano_app = newreg;
Line Loop(piano_app) = {l5,l6,l7,l8};

s1 = news;
Plane Surface(s1) = {piano_app};

b[] = Extrude {0,-4,0} 
	{ Surface{s1}; Recombine; };		// BASE ( b[1] )
Vol_qt[s] = b[1];
s += 1;

lay[e] = b[0];
e += 1;
lay[e] = b[2];
e += 1;
lay[e] = b[3];
e += 1;
lay[e] = b[4];
e += 1;
lay[e] = b[5];
e += 1;

Return


//**********************************************************************************************

Function Circonf

pc = newp;
Point(pc) = {x, y, z, d};
p10 = newp;
Point(p10) = {x+10, y, z, d};
p11 = newp;
Point(p11) = {x, y, z-10, d};
p12 = newp;
Point(p12) = {x-10, y, z, d};
p13 = newp;
Point(p13) = {x, y, z+10, d};
c1 = newl;
Circle(c1) = {p10, pc, p11};
c2 = newl;
Circle(c2) = {p11, pc, p12};
c3 = newl;
Circle(c3) = {p12, pc, p13};
c4 = newl;
Circle(c4) = {p13, pc, p10};
circ[c] = newl;
Line Loop(circ[c]) = {c1, c2, c3, c4};
c += 1;

Return

//***********************************************************************************************
Function HoleySurface


s10 = news;
Plane Surface(s10) = {piano_app, circ[]};
lay[e] = s10;


Return

//*********************************************************************************************


// Realizziamo i q.d. e le basi
x = 0; y = 0; z = 0; s = 1;
e = 0; // indice per l'array di superfici che realizzano il volume di uno strato


For v In {1:2}
	c = 0;	// indice dell'array delle circonferenze di ciascuna base d'appoggio
	For u In {1:4}
		For t In {1:4}
			Call Quantum;
			Call Circonf;	// Realizza la circonferenza di base del q.d. ...
					// ... e la mette nell' array circ[]
			x += 100;
		EndFor
		x = 0;
		z += 100;
	EndFor
	x = 0; z = 0;

	// Realizziamo la base corrispondente
	Call Base;
	Call HoleySurface;  	// realizza la superficie bucata
				// e la mette come ultimo elemento di lay[]
	theloops[v] = news;
	Surface Loop(theloops[v]) = {lay[]};
	y += 100; 
EndFor



// Volume fisico VOL_QT: volumi di tutti i q.d. + volumi di tutte le basi; 

// Sia i volumi dei q.t. che quelli delle basi sono ottenuti tramite estrusione; 
// bisogna ottenere un unico array di tutti questi volumi

VOL_QT = newv;
Physical Volume ("VOL_QT") = {Vol_qt[]};




// Calcoliamoci il volume della barriera
	
// Volume cubo esterno
p9 = newp;
Point(p9) = {-60, -10, 360, d};
p10 = newp;
Point(p10) = {360, -10, 360, d};
p11 = newp;
Point(p11) = {360, -10, -60, d};
p12 = newp;
Point(p12) = {-60, -10, -60, d};

l9 = newl;
Line(l9) = {p9,p10};
l10 = newl;
Line(l10) = {p10,p11};
l11 = newl;
Line(l11) = {p11,p12};
l12 = newl;
Line(l12) = {p12,p9};

l13 = newl;
Line Loop(l13) = {l9,l10,l11,l12};

s10 = news;
Plane Surface(s10) = {l13};

t4[] = Extrude {0,130,0} {Surface{s10}; Recombine;};	// t4[1]

theloops[0] = news;
Surface Loop(theloops[0]) = {s10, t4[0], t4[2], t4[3], t4[4], t4[5]};


/*
 theloops[1] e theloops[2] vengono calcolate all' interno del ciclo For.
 Adesso resta da aggiungere all' interno di theloops[1] e theloops[2] le 
 rispettive superfici bucate dove poggiano i q.d.
 Costruiamo questa superficie tramite le funzioni Circonf e HoleySurface
 che vengono chiamate all' interno del ciclo For
*/


barrier = newv;
Volume(barrier) = {theloops[]};

// We finally define a physical volume for barrier:

Physical Volume("BARRIER") = barrier;















