/*This script calculates the interfacial energy of a lattice cell containing two H0 particles


*/
//********************************************************** Simulation Parameters **************************************

parameter d_sep = 2 + 0.1 //use this parameter to change distance between particles
parameter Amp = 0.1 //Amplitide of the undulations relative to the longest side length
parameter z_offset = 0
CONSTRAINT_TOLERANCE 1e-12

parameter eta = 20 //Sharpness of corners of in plane shape of polygonal plate
parameter delta = 0 //offset of the undulations wrt particle shape (in radians)

//******************************************* Particle configurations ***************************************

parameter rot1 = -pi/6
parameter tilt1 = 0
parameter p1x = sqrt(3)*d_sep
parameter p1y = d_sep

parameter rot2 = pi/6
parameter tilt2 = 0
parameter p2x = sqrt(3)*d_sep/2
parameter p2y = d_sep/2

//****************************************** Reference frame Transformations *********************************

#define x_prime1 (((y-p1y)*sin(rot1)*cos(tilt1))+(cos(rot1)*cos(tilt1)*(x-p1x))-(sin(tilt1)*(z-z_offset))) 	
#define y_prime1 (-(x-p1x)*sin(rot1) + (y-p1y)*cos(rot1))
#define z_prime1 (((y-p1y)*sin(rot1)*sin(tilt1))+(cos(rot1)*sin(tilt1)*(x-p1x))+(cos(tilt1)*(z-z_offset)))

#define x_prime2 (((y-p2y)*sin(rot2)*cos(tilt2))+(cos(rot2)*cos(tilt2)*(x-p2x))-(sin(tilt2)*(z-z_offset))) 	
#define y_prime2 (-(x-p2x)*sin(rot2) + (y-p2y)*cos(rot2))
#define z_prime2 (((y-p2y)*sin(rot2)*sin(tilt2))+(cos(rot2)*sin(tilt2)*(x-p2x))+(cos(tilt2)*(z-z_offset)))

//*********************************************** shape definition parameters **********************************

//Defining paramters for hexagon corner pointing in +ve x
parameter a1 = (3^0.5)/2
parameter a2 = (3^0.5)/2
parameter t1 = pi/2
parameter t2 = 7*pi/6
parameter t3 = 11*pi/6

parameter a_mean = (a1 + a2) / 2
parameter a_diff = (a1 - a2) / 2

//********************************************* barrier functions ***********************************************************************

constraint barrierfn1 //inplane shape for particle 1
formula: ((cos(t1)*x_prime1 + sin(t1)*y_prime1 - a_diff) / a_mean)^eta +
	((cos(t2)*x_prime1 + sin(t2)*y_prime1 - a_diff) / a_mean)^eta +
	((cos(t3)*x_prime1 + sin(t3)*y_prime1 - a_diff) / a_mean)^eta = 1

constraint barrierfn2 //inplane shape for particle 2
formula: ((cos(t1)*x_prime2 + sin(t1)*y_prime2 - a_diff) / a_mean)^eta +
	((cos(t2)*x_prime2 + sin(t2)*y_prime2 - a_diff) / a_mean)^eta +
	((cos(t3)*x_prime2 + sin(t3)*y_prime2 - a_diff) / a_mean)^eta = 1

//********************************************* Contact line undulations ***************************************************************

constraint z01
formula: z_prime1 = Amp * sqrt((x_prime1^2)+(y_prime1^2)) * cos(3 * (atan2(y_prime1,x_prime1)-delta)) //+ Amp * sqrt((x_prime1^2)+(y_prime1^2)) * cos(6 * (atan2(y_prime1,x_prime1)-delta))

constraint z02
formula: z_prime2 = Amp * sqrt((x_prime2^2)+(y_prime2^2)) * cos(3 * (atan2(y_prime2,x_prime2)-delta)) //+ Amp * sqrt((x_prime2^2)+(y_prime2^2)) * cos(6 * (atan2(y_prime2,x_prime2)-delta))

//********************************************* Torus model (translational symmetry) for periodic boundaries ************************************************************

TORUS_FILLED
periods
d_sep*sqrt(3) 0 0 //Lattice Vector a
d_sep*sqrt(3)/2 d_sep*3/2 0 //Lattice vector b
0 0 10*d_sep //Lattice vector c, used for 3D lattices, set to large value to model 2d lattices

display_origin 0 0 -1

//******************************************** Geometry ***************************************************************************************
/*The geometry for vertices are defined as usual, edges are defined as a start and end point followed by lattice vector a,b,c transformations
this allows our boundaries to share interfacial deformation by copying another edge which mimics a larger simulation, i.e., periodic boundarys
the notation goes as follows: * represents no translation along that lattice vector, + means one positive translation along that lattice vector
and - means one negative translation along that lattice vector. The order of these indicate which lattice vector you translate with resepct to
for example + - * means positive translation along lattice vector a, negative translation along lattice vector b and no translation along lattice vector c
*/

vertices
//Hexagon1 Parimeter
1 p1x+cos(1*pi/6)	p1y+sin(1*pi/6) 0 constraints barrierfn1,z01
2 p1x+cos(3*pi/6)	p1y+sin(3*pi/6) 0 constraints barrierfn1,z01
3 p1x+cos(5*pi/6)	p1y+sin(5*pi/6) 0 constraints barrierfn1,z01
4 p1x+cos(7*pi/6)	p1y+sin(7*pi/6) 0 constraints barrierfn1,z01
5 p1x+cos(9*pi/6)	p1y+sin(9*pi/6) 0 constraints barrierfn1,z01
6 p1x+cos(11*pi/6)	p1y+sin(11*pi/6) 0 constraints barrierfn1,z01

//Hexagon2 Parimeter
7 p2x+cos(1*pi/6)	p2y+sin(1*pi/6) 0 constraints barrierfn2,z02
8 p2x+cos(3*pi/6)	p2y+sin(3*pi/6) 0 constraints barrierfn2,z02
9 p2x+cos(5*pi/6)	p2y+sin(5*pi/6) 0 constraints barrierfn2,z02
10 p2x+cos(7*pi/6)	p2y+sin(7*pi/6) 0 constraints barrierfn2,z02
11 p2x+cos(9*pi/6)	p2y+sin(9*pi/6) 0 constraints barrierfn2,z02
12 p2x+cos(11*pi/6)	p2y+sin(11*pi/6) 0 constraints barrierfn2,z02

//Interface Parimiter
13 0 0 0

edges
//Hexagon1 Parimeter
1 1 2 * * * constraints barrierfn1,z01 
2 2 3 * * * constraints barrierfn1,z01
3 3 4 * * * constraints barrierfn1,z01
4 4 5 * * * constraints barrierfn1,z01
5 5 6 * * * constraints barrierfn1,z01
6 6 1 * * * constraints barrierfn1,z01

//Hexagon2 Parimeter
7 7 8 * * * constraints barrierfn2,z02 
8 8 9 * * * constraints barrierfn2,z02
9 9 10 * * * constraints barrierfn2,z02
10 10 11 * * * constraints barrierfn2,z02
11 11 12 * * * constraints barrierfn2,z02
12 12 7 * * * constraints barrierfn2,z02


//interface edge
13 13 13 + * *
14 13 13 * + *

//Interface Connecting lines
15 13 9 * * *
16 13 10 * * *
17 13 11 * * *
18 11 13 + * * 
19 6 13 + * *
20 6 13 + + *
21 1 13 + + *
22 2 13 + + *
23 2 13 * + *
24 9 13 * + *

//particle connections
25 8 3 * * *
26 7 4 * * *
27 12 5 * * *

faces
1 14 -20 19
2 20 -21 -6
3 21 -22 -1
4 22 -13 -23
5 23 -24 -8 25 -2
6 24 -14 15
7 16 -9 -15
8 17 -10 -16
9 13 -18 -17
10 18 -19 -5 -27 -11
11 26 -3 -25 -7
12 27 -4 -26 -12

//************************************************************ Code scripting / Execution ********************************************************

read

direct := ".txt" //File output directory / name

transform_expr "2a2b"; //sets number of lattice cells to show

//try to even out the mesh
groom_length := 0.1;
groom := {
   	refine edge where length > groom_length and not no_refine;
   	u;V;
   	delete edge where length < groom_length/5;
}

//Helps avoid getting stuck in local minimas
hess_loop := {
	k 0;
	fix vertex where on_constraint barrierfn1;
	fix vertex where on_constraint z01;
	fix vertex where on_constraint barrierfn2;
	fix vertex where on_constraint z02;
	fix z_offset;
	hessian_seek; hessian_seek;
	unfix vertex where on_constraint barrierfn1;
	unfix vertex where on_constraint z01;
	unfix vertex where on_constraint barrierfn2;
	unfix vertex where on_constraint z02;
	unfix z_offset;
	k 1;
}

//Allows the simulation to be reset easily
reset := {
	LOAD "D:/Evolver/Myfiles/Hex0 Honey.txt"
}

//Output configuration data and energy to text file
outp := { recalc;
	printf "%f | %f",total_energy,d_sep >> direct;
}

//Meshing, tessallation parameters and procedure
gogo := {unfix z_offset; 
   groom 8;
   {g 5; groom; } 20;
   hess_loop;
   { g 5; groom; } 5;
}

// standard tessallation and minimisation procuedure (step 1) re Ken Brakke
go_more := {
   { recalc; M 1;
     {g 5; groom } 10; 
      hess_loop;
     {g 5; groom } 10; 
     M 2; {g 5; u; V;} 10;
     refine edge where on_constraint barrierfn1;
     refine edge where on_constraint barrierfn2;
     u; V;  
     {g 5; u; V;} 50;
      }
}

//Command to run minimization procedure
run:= {gogo 2; go_more;outp}