// naufrage

// paramètres
PARAMETER nu = 0.3 
PARAMETER pos = -0.1
PARAMETER Radius = 0.1
PARAMETER max_angle = 0.2;

//condition à la limite (pour ne pas briser la symetrie)
constraint 1
formula : x1=1

//condition de non-interpenetration suite à une contrainte appliquée selon l'axe de symétrie
constraint 2 nonnegative
formula: (x-1)^2+(y-pos)^2-Radius^2 

// Définitions energies
// 1-energie d'extension
quantity stretch energy modulus 10 method linear_elastic 
define facet attribute poisson_ratio real
define facet attribute form_factors real[3] 
// pour avoir un �tat de reference et faire des raffinement.
define vertex attribute ref_coord real[3]
define vertex attribute old_vid integer
define edge attribute	old_eid integer
// Pour info seulement, les composantes principales des d�formations
quantity stretch1 info_only method relaxed_elastic1_A global
quantity stretch2 info_only method relaxed_elastic2_A global

// 2-energie de courbure
quantity bending energy modulus 10 method star_perp_sq_mean_curvature 


//Géométrie
vertices
1	0	0	0  	fixed
2	0.49	0	0
3	0.5	0.5	0
4	0.51	0	0
5	1 	0	0	constraint 1
6	1	1	0	fixed constraint 1
7	0	1	0	fixed
8	0.5	1	0	fixed

edges
1 	1	2 color blue
2	2	3 color blue
3	3	4 color blue
4	4	5 color blue
5	5	6 color red constraint 1	
6	6	8 fixed color blue
7	7	1 fixed color blue
8	3	8
9	8	7 fixed color blue
	

faces
1	1 2 8 9 7 tension 0 
2	3 4 5 6 -8 tension 0


read
//********************************************************************
// Definir les form_fators � partir des positions de references connues.
set_form_factor:={
set vertex quantity bending;
foreach facet ff do 
	{ set ff.form_factors[1]  
       	   	(ff.vertex[2].ref_coord[1] - ff.vertex[1].ref_coord[1])^2
       	 	+ (ff.vertex[2].ref_coord[2] - ff.vertex[1].ref_coord[2])^2
       	 	+ (ff.vertex[2].ref_coord[3] - ff.vertex[1].ref_coord[3])^2;
  	set ff.form_factors[2]  
        	  (ff.vertex[2].ref_coord[1] - ff.vertex[1].ref_coord[1])
         	*(ff.vertex[3].ref_coord[1] - ff.vertex[1].ref_coord[1])
        	+ (ff.vertex[2].ref_coord[2] - ff.vertex[1].ref_coord[2])
         	*(ff.vertex[3].ref_coord[2] - ff.vertex[1].ref_coord[2])
        	+ (ff.vertex[2].ref_coord[3] - ff.vertex[1].ref_coord[3])
         	*(ff.vertex[3].ref_coord[3] - ff.vertex[1].ref_coord[3]);
  	set ff.form_factors[3]  
       	   (ff.vertex[3].ref_coord[1] - ff.vertex[1].ref_coord[1])^2
        	+ (ff.vertex[3].ref_coord[2] - ff.vertex[1].ref_coord[2])^2
        	+ (ff.vertex[3].ref_coord[3] - ff.vertex[1].ref_coord[3])^2;
	set ff.poisson_ratio nu;
   	set ff quantity stretch;
	};
	foreach edge ee where color = blue OR color = red do {unset ee.vertex quantity bending; };
};

//fixer la configuration de r�f�rence + d�finir �lasticit�
setformini := {
	set vertex ref_coord[1] x;
	set vertex ref_coord[2] y;
	set vertex ref_coord[3] z;
	set vertex old_vid id;
     	set edge old_eid id;
	set_form_factor;
};
setformini;

// Apr�s manipulation de la maille, red�finition de l'�lasticit� + contraint 1
setform := {
//Definir les coordonn�es initiales pour ces nouveaux vertex
     foreach vertex vv where old_vid == 0 do 
     { vv.ref_coord[1] := 
          avg(vv.edge ee where old_eid != 0,sum(ee.vertex where old_vid != 0,
                    ref_coord[1]));
       vv.ref_coord[2] := 
          avg(vv.edge ee where old_eid != 0,sum(ee.vertex where old_vid != 0,
                    ref_coord[2]));
       vv.ref_coord[3] := 
          avg(vv.edge ee where old_eid != 0,sum(ee.vertex where old_vid != 0,
                    ref_coord[3]));
     };
     set vertex old_vid id;
     set edge old_eid id;
//redefinir l'�lasticit� pour ce maillage en utilisant les positions initiales
     	set_form_factor;	
	recalc;
};

//rajoute des points l� o� la courbure est mal prise en compte (angle important entre les facettes)
rajoute:=
{
set edges no_refine;
foreach edge ee where color=black and old_eid!=0 do
	{
	if (1-abs((ee.facet[1].x * ee.facet[2].x + ee.facet[1].y * ee.facet[2].y + ee.facet[1].z * ee.facet[2].z)/ee.facet[1].area/ee.facet[2].area)) > max_angle then 
		{
		foreach ee.facet ff do
			{ 
			foreach ff.edge eeee where eeee.id != ee.id do {unset eeee no_refine};
			};
		};
	};
r;
unset edges no_refine;
setform;
};

// histogramme des "angles" entre facettes pour aider � choisir max_angle
hist_angle:=
{
histogram(edge ee where color=black, (1-abs((ee.facet[1].x * ee.facet[2].x + ee.facet[1].y * ee.facet[2].y + ee.facet[1].z * ee.facet[2].z)/ee.facet[1].area/ee.facet[2].area)))
};

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


// raffine orthoradialement autour du vertex 3
ortho_r:=
{
set edges no_refine;
foreach vertex[3].edge ee do { unset ee no_refine };
r;
unset edges no_refine;
setform;
};

//raffine radialement autour du vertex 3
radial_r:=
{
set edges no_refine;
foreach vertex[3].facet ff do 
	{ 
	foreach ff.edge ee do {unset ee no_refine};
	};
foreach vertex[3].edge ee do {set ee no_refine;};
r;
unset edges no_refine;
setform;
};


// Enregistrer les positions actuelles pour post-traitement.
Post :=
{foreach facet ff do
	{foreach ff.edge ee do 
		{foreach ee.vertex vv do 
			{print vv x >> "naufrage_fin";
			print vv y >> "naufrage_fin";
			print vv z >> "naufrage_fin";
			print 0.0 >> "naufrage_fin";
			print 0.0 >> "naufrage_fin";
               		print 0.0 >> "naufrage_fin";
			};
		};	
	};
foreach facet ff do
	{foreach ff.edge ee do 
		{foreach ee.vertex vv do 
			{print vv ref_coord[1] >> "naufrage_ini";
			print vv ref_coord[2] >> "naufrage_ini";
			print vv ref_coord[3] >> "naufrage_ini";
			print 0.0 >> "naufrage_ini";
			print 0.0 >> "naufrage_ini";
               		print 0.0 >> "naufrage_ini";
			};
		};	
	};
};

prepare := {
vertex[5].y:=0.2;
vertex[5].z:=0.2;
fix vertex[5];
};