// Essai d'ecrasement
// realise sur un tuyau 72'' (design 8A) sans renforcement TFC

// The anisotropic 2D mesh generator can be selected with:
//Mesh.Algorithm = 1 ;

// One can force a 4 step Laplacian smoothing of the mesh with:
//Mesh.Smoothing = 4 ;


// Centre du tuyau :
Point(1) = {0,0,0,1};

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

* Rayons des differentes couches

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

// nb de couches de TFC
nb_couche = 1;

// on definit les epaisseurs des differentes couches :
Eliner = 0.00152;
Ecore = 7.5*0.0254;
Ecoating = 0.05;
Ewire = 0.00156;
Etfc  = nb_couche * 0.00042;

// d'ou les rayons
Rexttfc = 72*0.0254/2;
Rint = Rexttfc - Etfc;
Rextliner = 76*0.0254/2;
Rintliner = Rextliner - Eliner;
Rextcore = Rexttfc + Ecore;
Rextwire = Rextcore + Ewire;
Rext = Rextcore + Ecoating;


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

* Parametres de l'essai

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

// largeur de l'appui en clef (appuis supposes centres sur l'axe (Oy), permettant d'imposer la force)
l_clef = 4*0.0254;

// angle forme par les axes medians des appuis (supposes admettre (Oy) pour axe de symetrie)
theta = 20.2674 * Pi / 180;
// largeur des appuis (theta), suppose tangents a l'exterieur du pipe
l_appui = l_clef;

// d'ou, secteurs angulaires concernes 
alpha_l_clef = Atan(l_clef / (2*Rext));
alpha_l_appui = Atan(l_appui / Rext);
// angle entre (Oy) et un appui
alpha_appui = Pi - theta/2;


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

* Nombre de mailles par couche et redecoupage

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

// nb secteurs angulaires pour le maillage de la section du tuyau :
// !!!!!!!!! nbmailles < 1000 !!!!!!!!!!
// !!!!!!!!! nbmailles pair !!!!!!!!!!
nbmailles = 144;

// pour les mailles physicals
nb_coating = 1;
nb_wire = 1;


// On d�ermine les indices des points en contacts avec les appuis (clef et "reins")
n_l_clef = Ceil(alpha_l_clef/(2*Pi/nbmailles));
n_appui = Floor(alpha_appui/(2*Pi/nbmailles));
n_l_appui = Floor(alpha_l_appui/(2*Pi/nbmailles));
Printf("nb de mailles a l'appui de clef : 2 x %g ; nb de mailles a l'appui : %g", n_l_clef, n_l_appui);

// parametres utilises pour definir le nb de mailles par secteurs angulaires
// b : selon l'arc de cercle
// c : selon le rayon
b = 3;
c = 1;

// longueur caracteristique pour la creation du maillage
lc = 1;
// coefficient de raffinage (raison d'une suite geometrique) pour les mailles, radialement
refine = 0.5;
refine_core = 0.1;
// nb de mailles raffinees, orthoradialement, au voisinage des mailles d'appuis
nb_maill_raff = 1;

// angle des coupures pour la suppression de la precontrainte
nb_coup = 4;
angl_coup[1] = Pi / 4;
angl_coup[2] = 3*Pi / 4;
angl_coup[3] = 5*Pi / 4;
angl_coup[4] = 7*Pi / 4;
//secteur angulaire concerne par la coupure
// !!!!! < 2Pi / nbmailles !!!!!
// hyp : 3mm sur la surface exterieure
alph_coup = 2.6e-3;

For i In {1 : nb_coup}
	a_coup[i] = Ceil(angl_coup[i]*nbmailles / (2*Pi));
	If ( angl_coup[i]*nbmailles / (2*Pi) - Floor(angl_coup[i]*nbmailles / (2*Pi)) < 0.5 )
		a_coup[i] = Floor(angl_coup[i]*nbmailles / (2*Pi));
	EndIf
EndFor


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

* Points de base

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

// NB : 4 "coordonnees" pour les points : la 4eme est la "longueur caracteristique" de la maille que l'on veut en ce point
// On commence sur l'axe (0y)
Point(10000) = {0,Rext,0,lc};
Point(20000) = {0,Rextwire,0,lc};
Point(30000) = {0,Rextcore,0,lc};
Point(40000) = {0,Rextliner,0,lc};
Point(50000) = {0,Rintliner,0,lc};
Point(60000) = {0,Rexttfc,0,lc};
Point(70000) = {0,Rint,0,lc};

	Line(15000) = {10000,20000};
	Line(25000) = {20000,30000};
	Line(35000) = {30000,40000};
	Line(45000) = {40000,50000};
	Line(55000) = {50000,60000};
	Line(65000) = {60000,70000};

		Transfinite Line (15000) = 4*c Using Progression refine;
		Transfinite Line (25000) = 4*c;
		Transfinite Line (35000) = 9*c Using Bump refine_core;
		Transfinite Line (45000) = 4*c;
		Transfinite Line (55000) = 5*c Using Bump refine_core;
		Transfinite Line (65000) = 3*c;

//indice de parcours du tableau contenant les coupures
i_coup = 1;
nb_coup_done = 0;
next_of_coup = 0;

For a In {1 : nbmailles}

i = 10000 + a + nb_coup_done;
j = 20000 + a + nb_coup_done;
k = 30000 + a + nb_coup_done;
l = 40000 + a;
m = 50000 + a;
n = 60000 + a;
o = 70000 + a;

alpha = -2*a*Pi/(nbmailles);
Printf("maille %g", a);

p = i + 5000;
q = j + 5000;
r = k + 5000 - nb_coup_done;
s = l + 5000;
t = m + 5000;
u = n + 5000;

// sens trigonometrique
	Point(i) = {Rext*Sin(alpha),Rext*Cos(alpha),0,lc};
	Point(j) = {Rextwire*Sin(alpha),Rextwire*Cos(alpha),0,lc};
	If (a == a_coup[i_coup] )
		Point(i+1) = {Rext*Sin(alpha-alph_coup),Rext*Cos(alpha-alph_coup),0,lc};
		Point(j+1) = {Rextwire*Sin(alpha-alph_coup),Rextwire*Cos(alpha-alph_coup),0,lc};
		Point(k+1) = {Rextcore*Sin(alpha-alph_coup),Rextcore*Cos(alpha-alph_coup),0,lc};
	EndIf
	Point(k) = {Rextcore*Sin(alpha),Rextcore*Cos(alpha),0,lc};
	Point(l) = {Rextliner*Sin(alpha),Rextliner*Cos(alpha),0,lc};
	Point(m) = {Rintliner*Sin(alpha),Rintliner*Cos(alpha),0,lc};
	Point(n) = {Rexttfc*Sin(alpha),Rexttfc*Cos(alpha),0,lc};
	Point(o) = {Rint*Sin(alpha),Rint*Cos(alpha),0,lc};

		If (a == a_coup[i_coup] )
			Printf("Line Loop %g", p+1);
			Line(p+1) = {i+1,j+1};
			Line(q+1) = {j+1,k+1};
		EndIf
		Line(p) = {i,j};
		Line(q) = {j,k};
		Line(r) = {k,l};
		Line(s) = {l,m};
		Line(t) = {m,n};
		Line(u) = {n,o};

			If (a == a_coup[i_coup] )
				Transfinite Line (p+1) = 4*c Using Progression refine;
				Transfinite Line (q+1) = 4*c;
			EndIf
			Transfinite Line (p) = 4*c Using Progression refine;
			Transfinite Line (q) = 4*c;
			Transfinite Line (r) = 9*c Using Bump refine_core;
			Transfinite Line (s) = 4*c;
			Transfinite Line (t) = 5*c Using Bump refine_core;
			Transfinite Line (u) = 3*c;

		If (a == a_coup[i_coup] )
			Circle(i+1) = {i,1,i+1};
			Circle(j+1) = {j,1,j+1};
			Circle(k+1) = {k,1,k+1};
		EndIf
		Circle(i) = {i-1,1,i};
		Circle(j) = {j-1,1,j};
		Circle(k) = {k-1,1,k};
		Circle(l) = {l-1,1,l};
		Circle(m) = {m-1,1,m};
		Circle(n) = {n-1,1,n};
		Circle(o) = {o-1,1,o};

			If (a == a_coup[i_coup] )
				Transfinite Line (i+1) = 1;
				Transfinite Line (j+1) = 1;
				Transfinite Line (k+1) = 1;
			EndIf
			Transfinite Line (i) = b;
			Transfinite Line (j) = b;
			Transfinite Line (k) = b;
/*			If (next_of_coup)
				Transfinite Line (l) = 2*b;
			EndIf
			If (!next_of_coup)
*/				Transfinite Line (l) = b;
//			EndIf
			Transfinite Line (m) = b;
			Transfinite Line (n) = b;
			Transfinite Line (o) = b;

// !!! le coup du "+1" est très crade... !!! - on compense les erreurs d'arrondi pour CE cas, à la clé et à l'appui droit...
// !!! on suppose qu'on n'a pas de coupure superposée aux appuis... !!!!
			If ( (a <= n_l_clef + nb_maill_raff) || (a >= nbmailles-n_l_clef+1 - nb_maill_raff) || ( (a >= n_appui - nb_maill_raff) && (a <= n_appui+n_l_appui+nb_maill_raff) ) || ( (a >= nbmailles-n_appui-n_l_appui+1 - nb_maill_raff) && (a <= nbmailles-n_appui+1 +nb_maill_raff) ) )
				Transfinite Line (i) = 3*b;
				Transfinite Line (j) = 3*b;
				Transfinite Line (k) = 3*b;
				Transfinite Line (l) = 3*b;
				Printf("appui");
			EndIf			


		If (a == a_coup[i_coup] )
			Line Loop (p+1) = {p,j+1,-(p+1),-(i+1)};
			Line Loop (q+1) = {q,k+1,-(q+1),-(j+1)};
		EndIf
		Line Loop (p) = {p-1,j,-p,-i};
		Line Loop (q) = {q-1,k,-q,-j};
		If (next_of_coup)
			Line Loop (r) = {r-1,l,-r,-k,-(k-1)};
		EndIf
		If (!next_of_coup)
			Line Loop (r) = {r-1,l,-r,-k};
		EndIf
		Line Loop (s) = {s-1,m,-s,-l};
		Line Loop (t) = {t-1,n,-t,-m};
		Line Loop (u) = {u-1,o,-u,-n};

	If (a == a_coup[i_coup] )
		Plane Surface (p+1) = {p+1};
		Plane Surface (q+1) = {q+1};
	EndIf
	Plane Surface (p) = {p};
	Plane Surface (q) = {q};
	Plane Surface (r) = {r};
	Plane Surface (s) = {s};
	Plane Surface (t) = {t};
	Plane Surface (u) = {u};

		If (a == a_coup[i_coup] )
			Transfinite Surface {p+1} = {i+1,j+1,j,i};
			Transfinite Surface {q+1} = {j+1,k+1,k,j};
			
			// listes pour les physicals
			list_coup_coating[nb_coating] = p+1;
			list_coup_wire[nb_wire] = q+1;
			nb_coating ++;
			nb_wire ++;
		EndIf
		Transfinite Surface {p} = {i,j,j-1,i-1};
		Transfinite Surface {q} = {j,k,k-1,j-1};
		Transfinite Surface {r} = {k,l,l-1,k-1};
		Transfinite Surface {s} = {l,m,m-1,l-1};
		Transfinite Surface {t} = {m,n,n-1,m-1};
		Transfinite Surface {u} = {n,o,o-1,n-1};
		list_coating[a] = p;
		list_wire[a] = q;

		If (a == a_coup[i_coup] )
			Recombine Surface {p+1};
			Recombine Surface {q+1};
		EndIf
		Recombine Surface {p};
		Recombine Surface {q};
		If (!next_of_coup)
			Recombine Surface {r};
		EndIf
		Recombine Surface {s};
		Recombine Surface {t};
		Recombine Surface {u};

	// on etait sur une coupure : on actualise les differents indices pour la prochaine etape...
	If (next_of_coup)
		next_of_coup = 0;
	EndIf
	If (a == a_coup[i_coup] )
		Printf("coupure nb %g done", i_coup);
		If (i_coup < nb_coup)
			i_coup ++;
		EndIf
		nb_coup_done ++;
		next_of_coup = 1;
	EndIf

EndFor

// on supprime automatiquement les points en doublons dans le maillage 
// (d'o, pas de conditions sur les points alpha = {0, 2Pi}
Coherence;

 
	// clef
	Physical Point (100001) = {10000};

// !!!!! attention �l'orientation des Physical pour l'application ult�ieure, dans Aster, de la pression
// Soit : pour l' ordre de parcours noeuds n (1, 2) -> vecteur tangent t : 
// vecteur normal n, tq (n, t) rep�e direct.
// Puis, pour une pression p > 0, on a une force S x -p.n

	// la paroi exterieure (Rext)
	Physical Line (100000) = {10001 : 10000+nbmailles+nb_coup};
	// la ligne entre le coating et le core exterieur
	Physical Line (300000) = {30001 : 30000+nbmailles+nb_coup};


	// l'appui a la clef
	Physical Line (100011) = {10001 : 10000+n_l_clef, 10000+(nbmailles+nb_coup)-n_l_clef+1 : 10000+nbmailles+nb_coup};

// !!!!! On suppose ici qu'on a 2 coupures avant l'appui gauche, et aucune entre les 2 appuis !!!!!
	nb_coup_avt_g = 2;
	// appui gauche
	Physical Line (100012) = {10000+nb_coup_avt_g+n_appui:10000+nb_coup_avt_g+n_appui+n_l_appui};
		//et pour les conditions d'appui : point sur la face exterieure...
		// !! ici, on ne prend pas en compte, dans le modele aster, la couche de TFC => indice de la face interne = 50000
		Physical Point (100002) = {10000+nb_coup_avt_g+n_appui+n_l_appui};

	// appui droit
	Physical Line (100013) = {10000+nb_coup_avt_g+nbmailles-n_appui-n_l_appui +1 : 10000+nb_coup_avt_g+nbmailles-n_appui +1};
		// et pour les conditions d'appui
		Physical Point (100003) = {10000+nb_coup_avt_g+nbmailles-n_appui-n_l_appui};

	Physical Surface (150001) = {list_coating[]};
	Physical Surface (250001) = {list_wire[]};
	Physical Surface (150002) = {list_coup_coating[]};
	Physical Surface (250002) = {list_coup_wire[]};
	Physical Surface (350001) = {35001:35000+nbmailles};
	Physical Surface (450001) = {45001:45000+nbmailles};
	Physical Surface (550001) = {55001:55000+nbmailles};
	Physical Surface (650001) = {65001:65000+nbmailles};
