// Model of ... epistructure for simulations
// all length units are in 1.e-6 m




// --- GEOMETRICAL DATA (fixed) ---

// Sizes in y-direction (thickness of layers)
h_sub 	= 100.;
h_epi	= 7.209;
h_n		= 0.5 * h_epi;
h_act	= 0.1 * h_epi;
h_p		= 0.4 * h_epi;

// Sizes in x-direction
d_mesa 	= 150.;
d_gap	= 15.;
d_dcontact1 = 85.;
d_dcontact3 = 75.;
d_tot	= d_mesa + 2 * d_gap + 2 * d_dcontact1 + 2 * d_dcontact3;

// Shift in positive x-direction 
// (such that the left lower corner of the substrate is at ( 0, y_0, 0 ) 
d_shift = d_tot / 2;

// Lowest y coordinate
y_0 = 0.;
// -------------------------




// -------- MESH RESOLUTION (variable) --------

// Nomenclature
//		n:	number of layers
//		f:	mesh resolution scale factor for domains (~ # vertical points per sublayer)
//		p is for the p-DBR (in and out), n is for the n-DBR (in and out), act is the active layer


// Number of layers in each structure and their resolution (characteristic length, smaller -> finer resolution)

n_sub		= 1;
f_sub		= 0.5;

n_n_in 		= 1;
f_n_in		= 1;

n_act_in	= 2;
f_act_in	= 1;

n_p_in 		= 1;
f_p_in		= 1;

// define values like n_p_out, f_act_out here if needed 
// (i.e. if the resolution should be different for in and out)

// -------------------------

// Global: Scale  mesh resolution (characteristic length, smaller -> finer resolution)
Mesh.ScalingFactor = 1;


// Anisotropy factor (resolution in y-direction aniso times finer than in x-direction)
aniso = 10;

// -------------------------


// Vectors for iteration 
y[] = { h_sub, h_n, h_act, h_p , h_n, h_act, h_p};
x[] = { d_tot, d_mesa, d_mesa, d_mesa, d_dcontact3, d_dcontact3, d_dcontact3 };
n[] = { n_sub, n_n_in, n_act_in, n_p_in, n_n_in, n_act_in, n_p_in };	//n[] = { n_sub, n_n_in, n_act_in, n_p_in, n_n_out, n_act_out, n_p_out };
f[] = { f_sub, f_n_in, f_act_in, f_p_in, f_n_in, f_act_in, f_p_in };	//f[] = { f_sub, f_n_in, f_act_in, f_p_in, f_n_out, f_act_out, f_p_out };
// -------------------------





// -------- MISC OPTIONS --------

// Mesh algorithm
Mesh.Algorithm = 7; //2D mesh algorithm (1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=BAMG, 8=DelQuad) Default value: 2, BAMG for anisotropic meshes
Mesh.Binary = 1;
//Mesh.HighOrderOptimize = 1;
Mesh.SaveParametric = 1;
Mesh.CharacteristicLengthExtendFromBoundary = 1;

// Should all duplicate entities be automatically removed?
Geometry.AutoCoherence = 1; 

// Mesh random factor to allow finer meshes
Mesh.RandomFactor = 1e-7;			// default: 1e-9;

// (Don't) extend the elements sizes from the boundary inside the domain
Mesh.CharacteristicLengthExtendFromBoundary = 1;

// The default characteristic size of cells
cell_size = 250.0;

// -------------------------




// Def. function rect
// input: 2 x values, 2 y values, e.g. rectangle
Function rect
	p = newp;
	Point(p) = {r_x_0,r_y_0,0}; Point(p+1) = {r_x_1,r_y_0,0}; Point(p+2) = {r_x_1,r_y_1,0}; Point(p+3) = {r_x_0,r_y_1,0};
	points[] = { p, p+1, p+2, p+3 };
	l = newl;
	Line(l) = {p, p+1}; Line(l+1) = {p+1, p+2}; Line(l+2) = {p+2, p+3}; Line(l+3) = {p+3, p};
	//Printf("rect: line =%g ", l);
	ll = newll; Line Loop(ll) = {l,l+1,l+2,l+3};
	//Printf("rect: line loop =%g ", ll);
	surf = news; Plane Surface(surf) = {ll}; // surf = newreg
	//Printf("rect: surface = %g ", surf);
Return




// Def. function mesh
// input: 2 x values, 2 y values, e.g. rectangle + mesh_fac + k
Function mesh
	Field[k] = AttractorAnisoCurve ;
	Field[k].NNodesByEdge 	= 800; // # attractors on the edges
	Field[k].EdgesList 		= { l - 6 * 2 / 2 + 2 };				// ...
	Field[k].dMin 			= y[i] / 2;
	Field[k].dMax 			= y[i] / 2 * 1.5;
	Field[k].lMinNormal 	= y[i] / n[i] / f[i];					// Resolution normal to the edge (y), inside dMin: ~ f[i] nodes in y-direction per layer
	Field[k].lMinTangent 	= y[i] / n[i] / f[i] * aniso;			// Resolution tangential to the edge (x), inside dMin: 10	 times finer in y-direction
	Field[k].lMaxNormal 	= y[i] / n[i] / f[i] * 1000;			// Resolution normal to the edge (y), outside dMin
	Field[k].lMaxTangent 	= y[i] / n[i] / f[i] * 1000 * aniso;	// Resolution tangential to the edge (x), outside dMin
	Printf('Line for AttractorAnisoCurve: %g, field: %g, normal resolution: %g',l - 6 * 2 / 2 + 2,k,y[i] / n[i] / f[i]);
Return 




// Use functions rect and mesh to iterate through physical surfaces
// First substrate and inner n-DBR, p-DBR
ps 		= 1;
k 		= 0;
k_ 		= 0;
y_cum 	= y_0;

For i In { 0 : 3 }	 //3	//#y[]-1
	For j In { 0 : 1 }
		r_y_0 = y_cum;				//+ 0.00001; // 10pm
		r_y_1 = y_cum + y[i] / 2;
		y_cum += y[i] / 2;
		r_x_0 = - x[i] / 2.;
		r_x_1 = x[i] / 2.;
		r_x_0 += d_shift;
		r_x_1 += d_shift;
		Call rect;
		k_+= 1;
	EndFor
	Call mesh;
	k += 1;	

	Physical Surface(ps) = {surf}; ps += 1;
EndFor


Physical Line(0) = {1:l};



// Background field

Field[k] = MinAniso;
Field[k].FieldsList = {0,1,2,3}; // highest number in array must be geq k 
Background Field = k;