// This script calculates the interfacial energy of two interacting H0 particles with centre-centre distance of 2*p1x

//**************************************** Simulation Parameters ***********************************************************************
parameter x_len = 8 //half length of simulation box x axis (set to 8*rc where rc is the radius of the circumscribed circle)
parameter y_len = 4 //half length of simulation box y axis (set to 4*rc where rc is the radius of the circumscribed circle)

parameter eta = 20 // Sharpness of corners of in plane shape of polygonal plate
parameter Amp = 0.1 //Amplitide of the undulations relative to the longest side length

parameter delta = 0 //offset of the undulations wrt particle shape (in radians)

CONSTRAINT_TOLERANCE 1e-6

//************************************** Particle 1 shape definition (left) *************************************************************

//Defining paramters for Hexagon corner pointing in +ve x
parameter p1a1 = (3^0.5)/2
parameter p1a2 = p1a1
parameter p1t1 = pi/6
parameter p1t2 = 3*pi/6
parameter p1t3 = 5*pi/6

parameter p1a_mean = (p1a1 + p1a2) / 2
parameter p1a_diff = (p1a1 - p1a2) / 2

//********************************************************** Particle 2 shape definition (Right) *********************************************************

//Defining paramters for Hexagon corner pointing in +ve x
parameter p2a1 = (3^0.5)/2
parameter p2a2 = p2a1
parameter p2t1 = pi/6
parameter p2t2 = 3*pi/6
parameter p2t3 = 5*pi/6

parameter p2a_mean = (p2a1 + p2a2) / 2
parameter p2a_diff = (p2a1 - p2a2) / 2

//************************************* Initial starting conditions *****************************************************************

//p1 (left) parameters
parameter p1x = -1.1 //Change to move particles appart
parameter p1y = 0
parameter p1z = 0 
parameter p1rot = 0 //Change to rotate particle 1 wrt lab frame
parameter p1tilt = 0

//p2 (right) parameters
parameter p2x = -p1x
parameter p2y = 0
parameter p2z = 0
parameter p2rot = 0 //Change to rotate particle 1 wrt lab frame
parameter p2tilt = 0

//************************************* Co-ordinate Transforms *****************************************************************

#define p1_x_prime (((y-p1y)*sin(p1rot)*cos(p1tilt))+(cos(p1rot)*cos(p1tilt)*(x-p1x))-(sin(p1tilt)*(z-p1z))) 	
#define p1_y_prime (-(x-p1x)*sin(p1rot) + (y-p1y)*cos(p1rot))
#define p1_z_prime (((y-p1y)*sin(p1rot)*sin(p1tilt))+(cos(p1rot)*sin(p1tilt)*(x-p1x))+(cos(p1tilt)*(z-p1z)))

#define p2_x_prime (((y-p2y)*sin(p2rot)*cos(p2tilt))+(cos(p2rot)*cos(p2tilt)*(x-p2x))-(sin(p2tilt)*(z-p2z))) 	
#define p2_y_prime (-(x-p2x)*sin(p2rot) + (y-p2y)*cos(p2rot))
#define p2_z_prime (((y-p2y)*sin(p2rot)*sin(p2tilt))+(cos(p2rot)*sin(p2tilt)*(x-p2x))+(cos(p2tilt)*(z-p2z)))

//************************************* Particle 1 barrier finction ************************************************

constraint barrierfn1 //for particles with 6 fold symetry
formula: ((cos(p1t1)*p1_x_prime + sin(p1t1)*p1_y_prime - p1a_diff) / p1a_mean)^eta +
	((cos(p1t2)*p1_x_prime + sin(p1t2)*p1_y_prime - p1a_diff) / p1a_mean)^eta +
	((cos(p1t3)*p1_x_prime + sin(p1t3)*p1_y_prime - p1a_diff) / p1a_mean)^eta = 1

//*********************************************** Particle 2 barrier function ***************************************************

constraint barrierfn2 //for particles with 6 fold symetry
formula: ((cos(p2t1)*p2_x_prime + sin(p2t1)*p2_y_prime - p2a_diff) / p2a_mean)^eta +
	((cos(p2t2)*p2_x_prime + sin(p2t2)*p2_y_prime - p2a_diff) / p2a_mean)^eta +
	((cos(p2t3)*p2_x_prime + sin(p2t3)*p2_y_prime - p2a_diff) / p2a_mean)^eta = 1

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

constraint z1
formula: p1_z_prime = Amp * ((p1_x_prime^2 + p1_y_prime^2)^0.5) * cos( 3 * (atan2(p1_y_prime,p1_x_prime)-delta))

constraint z2
formula: p2_z_prime = Amp * ((p2_x_prime^2 + p2_y_prime^2)^0.5) * cos( 3 * (atan2(p2_y_prime,p2_x_prime)-delta)) 

//********************************************** Interface pinning constraints *****************************************************

constraint x_min
formula: x = -x_len

constraint x_max
formula: x = x_len

constraint y_min
formula: y = -y_len

constraint y_max
formula: y = y_len

//********************************************** Geometry **************************************************************************

vertices
//Droplet perimiter
1 x_len y_len 0 constraints x_max, y_max
2 -x_len y_len 0 constraints x_min, y_max
3 -x_len 0 0 constraints x_min
4 -x_len -y_len 0 constraints x_min, y_min
5 x_len -y_len 0 constraints x_max, y_min
6 x_len 0 0 constraints x_max

//p 2 (Right)
7 cos(0*pi/3)+p2x	sin(0*pi/3)+p2y p2z constraints z2,barrierfn2
8 cos(1*pi/3)+p2x	sin(1*pi/3)+p2y p2z constraints z2,barrierfn2
9 cos(2*pi/3)+p2x	sin(2*pi/3)+p2y p2z constraints z2,barrierfn2
10 cos(3*pi/3)+p2x	sin(3*pi/3)+p2y p2z constraints z2,barrierfn2
11 cos(4*pi/3)+p2x	sin(4*pi/3)+p2y p2z constraints z2,barrierfn2
12 cos(5*pi/3)+p2x	sin(5*pi/3)+p2y p2z constraints z2,barrierfn2

//p 1 (left)
13 cos(0*pi/3)+p1x	sin(0*pi/3)+p1y p1z constraints z1,barrierfn1
14 cos(1*pi/3)+p1x	sin(1*pi/3)+p1y p1z constraints z1,barrierfn1
15 cos(2*pi/3)+p1x	sin(2*pi/3)+p1y p1z constraints z1,barrierfn1
16 cos(3*pi/3)+p1x	sin(3*pi/3)+p1y p1z constraints z1,barrierfn1
17 cos(4*pi/3)+p1x	sin(4*pi/3)+p1y p1z constraints z1,barrierfn1
18 cos(5*pi/3)+p1x	sin(5*pi/3)+p1y p1z constraints z1,barrierfn1

edges
//Droplet Perimiter
1 1 2 constraints y_max
2 2 3 constraints x_min
3 3 4 constraints x_min
4 4 5 constraints y_min
5 5 6 constraints x_max
6 6 1 constraints x_max

//Connecting lines
7 6 7
11 10 13
15 3 16

//p 2 (right)
8 7 8 constraints z2, barrierfn2
9 8 9 constraints z2, barrierfn2
10 9 10 constraints z2, barrierfn2
19 10 11 constraints z2, barrierfn2
20 11 12 constraints z2, barrierfn2
21 12 7 constraints z2, barrierfn2

//p 1 (left)
12 13 14 constraints z1, barrierfn1
13 14 15 constraints z1, barrierfn1
14 15 16 constraints z1, barrierfn1
16 16 17 constraints z1, barrierfn1
17 17 18 constraints z1, barrierfn1
18 18 13 constraints z1, barrierfn1


//Additional connecting lines
22 2 15
23 14 9
24 1 8
25 4 17
26 18 11
27 5 12

faces

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

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

read

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

//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 z1;
	fix vertex where on_constraint barrierfn2;
	fix vertex where on_constraint z2;
	hessian_seek; hessian_seek;
	unfix vertex where on_constraint barrierfn1;
	unfix vertex where on_constraint z1;
	unfix vertex where on_constraint barrierfn2;
	unfix vertex where on_constraint z2;
	k 1;
}

//Writes all vertices of surface around particle to a CSV
all_data :={recalc;
	Ad := sprintf "D:/Evolver/Myfiles/All Hex, Amp=%5.2f, Eta=%d, delta1=%d.csv", Amp, eta, (delta*180/pi);
	groom;
	printf "x-posn,y-posn,z-posn,\n" >> Ad;
	foreach vertex do {printf "%f,%f,%f,\n",x,y,z} >> Ad
}

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

//Allows the simulation to be reset easily
reset := {
	LOAD "D:/Evolver/Myfiles/Polygonal Plates interaction.txt"
}

//Meshing, tessallation parameters and procedure
gogo := { 
   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 8; go_more; outp}