/////////////
/// GMSH ///
/////////////

// CAD and meshing with cthe creation of cuts


/////////////////////
/// PARAMETERS ///
////////////////////

General.ExpertMode = 1;

/// Geometry options:
Geometry.Color.Lines = Black;
Geometry.LineWidth = 1;
Geometry.PointSize = 1;

/// Remove duplicate geometrical elementary entities:
Coherence;

/// Remove duplicate nodes of the mesh:
Coherence Mesh;

/// Optimize mesh:
Mesh.Optimize = 1;
Mesh.RandomFactor = 1e-9;
//Mesh.RemeshAlgorithm = 1;
Mesh.Algorithm = 6; // Use Frontal 2D algorithm

// Save mesh as binary file (to save storage):
Mesh.Binary = 1;

Include "hts_bench_3.par";

// Quadrangle mesh size:
If (isQuadMeshWanted == 1)
	nQuad_wire = 6;
	nQuadWidth = 180;
	bumpFactor = 1; // to refine progrssively the mesh at the extremities of the supercodncuting, silver and substrate layers
EndIf

// Local mesh size:
nn_wire = 0.008*wireWidth; // quadrangles = 0.004, triangles = 0.008
nn_iinf = 0.06*infiniteBoundaryInnerRadius; // quadrangles = 0.06, triangles = 0.08 (to check)
nn_oinf = 0.1*infiniteBoundaryOuterRadius;


////////////////////
/// GEOMETRIES ///
////////////////////

boundary[] = {}; // creation of an empty list

/// BEGIN LOOP ///
For i In {0:numberOfWires-1}

/*
Beware: various schemes have been tried out to obtained the thick cuts as expected. 
The creation of the wires as given below proved to be the most reliable without using extrusion.

Construction of the silver, SC and substrate layers following a counterclock wise rotation:
	     
  	   L3
  P4_______P3
  |	       |
L4|        | L2     
	|	Sc     |    
	|	       |
	P1_______P2 
	    L1       	
*/

// Wire:
  x_1 =  distanceToAxis+i*distanceBetweenWires; y_1 = (-1)*wireWidth/2; 	n_1 = nn_wire;
  x_2 = x_1+wireThickness;                      y_2 = y_1;                n_2 = nn_wire; 
  x_3 = x_2;		                                y_3 = (-1)*y_1;		        n_3 = nn_wire;
  x_4 = x_1;                                    y_4 = y_3;			          n_4 = nn_wire;	

	p~{1+i*4} = newp;	Point(p~{1+i*4}) = {x_1, y_1, 0, n_1}; 
	p~{2+i*4} = newp;	Point(p~{2+i*4}) = {x_2, y_2, 0, n_2};
	p~{3+i*4} = newp;	Point(p~{3+i*4}) = {x_3, y_3, 0, n_3}; 
	p~{4+i*4} = newp;	Point(p~{4+i*4}) = {x_4, y_4, 0, n_4};
	
	l~{1+i*4} = newl; Line(l~{1+i*4}) = {p~{1+i*4}, p~{2+i*4}};
	l~{2+i*4} = newl; Line(l~{2+i*4}) = {p~{2+i*4}, p~{3+i*4}};
	l~{3+i*4} = newl; Line(l~{3+i*4}) = {p~{3+i*4}, p~{4+i*4}};
	l~{4+i*4} = newl; Line(l~{4+i*4}) = {p~{4+i*4}, p~{1+i*4}};
		
	ll~{i} = newll; Line Loop(ll~{i}) = {l~{1+i*4}, l~{2+i*4}, l~{3+i*4}, l~{4+i*4}};
	boundary[] += ll~{i}; // tool for the hole of the insulation
	ps~{i} = news; Plane Surface(ps~{i}) = {ll~{i}};
	
  If (isQuadMeshWanted == 1)
	  Transfinite Line{l~{1+i*4}, l~{3+i*4}} = nQuad_wire;
    Transfinite Line{l~{2+i*4}, l~{4+i*4}} = nQuadWidth Using Bump bumpFactor; // Refine mesh at the extremities (progression from sides to inner tape width)
    Transfinite Surface{ps~{i}} = {p~{1+4*i}, p~{2+4*i}, p~{3+4*i}, p~{4+4*i}}; // Left: surface id , right: list of points at corners of the surface
    Recombine Surface {ps~{i}}; // just quad, surface id
  EndIf
	
  Physical Surface(wireID~{i}) = {ps~{i}};
	
	Color Black {Surface {ps~{i}};}
	
	// Definition of thermal and mechanical boundaries:
	Physical Line(nodeEdgeBoundaryID~{i}) = {l~{1+i*4}, l~{2+i*4}, l~{3+i*4}, l~{4+i*4}}; // boundary of wire
	
EndFor
/// END LOOP ///

/// Surrounding:
x_1 = 0; 						                y_1 =(-1)*infiniteBoundaryInnerRadius;  n_1 = nn_iinf;
x_2 = infiniteBoundaryInnerRadius;  y_2 = 0;                  							n_2 = nn_iinf; 
x_3 = 0;						                y_3 = infiniteBoundaryInnerRadius;	    n_3 = nn_iinf;

p_0 = newp; Point(p_0) = {0, 0, 0, n_1};
p_1 = newp;	Point(p_1) = {x_1, y_1, 0, n_1}; 
p_2 = newp;	Point(p_2) = {x_2, y_2, 0, n_2};
p_3 = newp;	Point(p_3) = {x_3, y_3, 0, n_3}; 
	
l_1 = newl; Circle(l_1) = {p_1, p_0, p_2};
l_2 = newl; Circle(l_2) = {p_2, p_0, p_3};
l_3 = newl; Line(l_3) = {p_3, p_1};

ll_1 = newll; Line Loop(ll_1) = {l_1, l_2, l_3};
ps~{1} = news; Plane Surface(ps~{1}) = {ll_1, boundary[]}; 
Physical Surface("Surroundings", surroundingID) = {ps~{1}};

/// Outer boundary (inifinite space) for shell transformation:
x_1 = 0; 						               y_1 =(-1)*infiniteBoundaryOuterRadius; n_1 = nn_oinf;
x_2 = infiniteBoundaryOuterRadius; y_2 = 0;                               n_2 = nn_oinf; 
x_3 = 0;						               y_3 = infiniteBoundaryOuterRadius;	    n_3 = nn_oinf;

p_4 = newp;	Point(p_4) = {x_1, y_1, 0, n_1}; 
p_5 = newp;	Point(p_5) = {x_2, y_2, 0, n_2};
p_6 = newp;	Point(p_6) = {x_3, y_3, 0, n_3}; 
	
l_4 = newl; Line(l_4) = {p_1, p_4};
l_5 = newl; Circle(l_5) = {p_4, p_0, p_5};
l_6 = newl; Circle(l_6) = {p_5, p_0, p_6};
l_7 = newl; Line(l_7) = {p_6, p_3};

ll_2 = newll; Line Loop(ll_2) = {l_4, l_5, l_6, l_7, -l_2, -l_1};
ps~{2} = news; Plane Surface(ps~{2}) = {ll_2};
Physical Surface("Infinite space", shellID) = {ps~{2}};

Color Blue {Surface {ps~{1}}; Surface {ps~{2}};} // surrounding and infinite space


///////////////////////////////////////
/// COHOMOLOGY AND HOMOLOGY ///
//////////////////////////////////////

/// Cohomology and Homology computation:
/*
Since we have a hole in a hole, the insulation beeing a multipe-connected domain as well as the surrounding, we need one cut for the insulation and one cut for the surrounding that we have to pair them in the solver to define the impressed current per tape (| : edge, \ : cut).
*/

For i In {0:numberOfWires-1}
listOfGeometries[] = {surroundingID, shellID};
	For j In {0:numberOfWires-1}
		If (j != i)
			listOfGeometries[] += {wireID~{j}}; // adding the wires to the group expect for the one for which we need to generate the class of homology group H_1
		EndIf
	EndFor
Printf("Cohomology for wire %g", wireID~{i});
Cohomology(1){ {listOfGeometries[]}, {}};
//Homology(1){ {listOfGeometries[]}, {}}; // Beware: uncomment to see the class of H_1. However, comment out to generate the cuts for the solver otherwise the labelling is wrong
EndFor

