
// Note : this file contains the definition of the problem geometry.
// To run this example you need to open the .pro file with GetDP.

// ----------------------------------------------------------------- //
//                ____           __          __                      //
//               /  _/___  _____/ /_  ______/ /__  _____             //
//               / // __ \/ ___/ / / / / __  / _ \/ ___/             //
//             _/ // / / / /__/ / /_/ / /_/ /  __(__  )              //
//            /___/_/ /_/\___/_/\__,_/\__,_/\___/____/               //
//                                                                   //
// ----------------------------------------------------------------- //

// Inclusion of the parameters already defined in the .par file ---- //

Include "2Daxi_single_pancake.par";

// ----------------------------------------------------------------- //
//       ____                                  __                    //
//      / __ \____  _________  ____ ___  ___  / /____  __________    //
//     / /_/ / __ `/ ___/ __ `/ __ `__ \/ _ \/ __/ _ \/ ___/ ___/    //
//    / ____/ /_/ / /  / /_/ / / / / / /  __/ /_/  __/ /  (__  )     //
//   /_/    \__,_/_/   \__,_/_/ /_/ /_/\___/\__/\___/_/  /____/      //
//                                                                   //
// ----------------------------------------------------------------- //

// Refer to 2Daxi_single_pancake.par for the definition of geometrical
// parameters ------------------------------------------------------ //

// Apply the scaling factor on the geometrical dimensions
TapeWidth     *= ScalingFactor ;
TapeThick     *= ScalingFactor ;
CoilIntRadius *= ScalingFactor ;
CoilStep      *= ScalingFactor ;
AirIntRadius  *= ScalingFactor ;
AirExtRadius  *= ScalingFactor ;

// Typical mesh sizes in the conductor and in the air
TapeMeshSize = TapeWidth / TapeWidthDivs;
AirMeshSize = ( AirExtRadius - AirIntRadius ) / AirDivs;

// Mesher options :
//  - We ask to save the mesh in a binary file (lightweight)
//  - We use ExpertMode so that Gmsh won't complain for mixed
//    structured-unstructured mesh
//  - We call the MeshAdapt algorithm and optimize its result
Mesh.Algorithm = 1;
Mesh.Binary = 1;
Mesh.Optimize = 1;
General.ExpertMode = 1;

// Some visualization options (in the interest of clarity)
Geometry.Color.Lines = Black;
Geometry.LineWidth = 1;
Geometry.PointSize = 1;

// ----------------------------------------------------------------- //
//            ______                          __                     //
//           / ____/__  ____  ____ ___  ___  / /________  __         //
//          / / __/ _ \/ __ \/ __ `__ \/ _ \/ __/ ___/ / / /         //
//         / /_/ /  __/ /_/ / / / / / /  __/ /_/ /  / /_/ /          //
//         \____/\___/\____/_/ /_/ /_/\___/\__/_/   \__, /           //
//                                                 /____/            //
//                                                                   //
// ----------------------------------------------------------------- //

p0 = newp;
Point(p0) = { 0 , 0 , 0 , AirMeshSize };

// Building the pancake structure ---------------------------------- //

// Prepare some lists to store elementary identifiers :
//   - sTapes[] will store the surface of the tapes
//   - lTapesSkin[] is for the edge contours of the tapes
//   - llTapesSkin[] contains edge loops around the tapes
sTapes[] = {} ; lTapesSkin[] = {} ; llTapesSkin = {} ;

For iTurn In { 1 : CoilTurns }
  
  // Computation of the radius of the current turn
  TurnRadius = CoilIntRadius + (iTurn - 1) * CoilStep ;

  // Create the inner side of the tape
  ps = newp;
  Point(ps)   = { TurnRadius ,  TapeWidth / 2 , 0 , TapeMeshSize } ;
  Point(ps+1) = { TurnRadius , -TapeWidth / 2 , 0 , TapeMeshSize } ;
  ls = newl ; Line(ls) = { ps , (ps+1) } ;
  
  // This extrusion generates the superconducting tape
  tmp[] = Extrude { TapeThick , 0 , 0 }
          { Line{ ls } ; Layers{ TapeThickDivs } ; Recombine ; } ;

  // We recover the information about the extruded geometry
  sTapes[] += tmp[1] ;
  lTapesSkin[] += { ls, tmp[3], -tmp[0], tmp[2] } ;
  lls = newll; Line Loop(lls) = { ls, tmp[3], -tmp[0], tmp[2] } ;
  llTapesSkin[] += { lls } ;

EndFor

// Building the first air enclosure -------------------------------- //

pa = newp ;
Point(pa)   = { 0             , -AirIntRadius , 0 , AirMeshSize } ;
Point(pa+1) = { AirIntRadius  , 0             , 0 , AirMeshSize } ;
Point(pa+2) = { 0             , AirIntRadius  , 0 , AirMeshSize } ;

la = newl ;
Circle(la)   = {  pa    , p0 , (pa+1) } ;
Circle(la+1) = { (pa+1) , p0 , (pa+2) } ;
Line(la+2)   = { (pa+2) , pa } ;

lla = newll ; Line Loop(lla) = { la, (la+1), (la+2) } ;
sa = news ; Plane Surface(sa) = { lla , llTapesSkin[] } ;

// Building the infinite domain ------------------------------------ //

pi = newp ;
Point(pi)   = { 0             , -AirExtRadius , 0 , AirMeshSize } ;
Point(pi+1) = { AirExtRadius  , 0             , 0 , AirMeshSize } ;
Point(pi+2) = { 0             , AirExtRadius  , 0 , AirMeshSize } ;

li = newl ;
Circle(li)   = {  pi    ,  p0 , (pi+1) } ;
Circle(li+1) = { (pi+1) ,  p0 , (pi+2) } ;
Line(li+2)   = { (pi+2) , (pa+2) } ;
Line(li+3)   = {  pa    ,  pi } ;

lli = newll ; Line Loop(lli) = {   li    , (li+1) , (li+2) ,
                                 -(la+1) , -la    , (li+3) } ;
si = news ; Plane Surface(si) = { lli } ;

// Final touches --------------------------------------------------- //

// Tag the physical regions with the IDs defined in the .par file 
Physical Surface(SUPRA)           = sTapes[] ;
For iTurn In { 1 : CoilTurns }
  Physical Surface(SUPRA+iTurn)   = sTapes[ (iTurn-1) ] ;
EndFor
Physical Surface(AIR)             = { sa } ;
Physical Surface(INFINITE_DOMAIN) = { si } ;
Physical Line(SUPRA_SKIN)         = lTapesSkin[] ;
Physical Line(FAR_LIMIT)          = { li , (li+1) } ;
Physical Line(SYMMETRY_AXIS)      = { (li+2) , (la+2) , (li+3) } ;

// Compute cuts in the air domain
For iTurn In { 1 : CoilTurns }
  TurnsToFill[] = {} ;
  If ( iTurn > 1 ) 
    TurnsToFill[] += { SUPRA+1:SUPRA+iTurn-1 } ;
  EndIf
  If ( iTurn < CoilTurns )
    TurnsToFill[] += { SUPRA+iTurn+1:SUPRA+CoilTurns } ;
  EndIf
  Cohomology(1) { { AIR , INFINITE_DOMAIN , TurnsToFill[] } , {} } ;
EndFor

// Some colors to help read the mesh
Color Orange  { Surface{ sTapes[] } ; }
Color SkyBlue { Surface{ sa , si } ; }
 