
// 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 "Bobinette_supra_2.par";

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

// Definition of the geometrical parameters of the winding

// Scaling all dimensions using the mm unit factor
TapeWidth = TapeWidth * mm;
SupraThick = SupraThick * mm;
MatrixThick = MatrixThick * mm;

// Typical mesh sizes in the conductor and in the air
TapeMeshSize = SupraThick / NbElemThick;
AirMeshSize = ( AirExtRadius - AirIntRadius ) / NbElemThick;

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

// Mesher setup
Mesh.Algorithm = 1;
Mesh.Binary = 1;
Mesh.Optimize = 1;
General.ExpertMode = 1;

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

// Construction of the coil turns ---------------------------------- //

SurfSupra[] = {};
SurfMatrix[] = {};
SurfAirInt[] = {};
LoopSupra[] = {};
LoopMatrix[] = {};
LoopAirInt[] = {};
ConductorSkin[] = {};

y1 = -TapeWidth/2 ; y2 = TapeWidth/2;
For iTurn In { 1 : NbTurns }
  x1 = CoilIntRadius + (iTurn - 1) * CoilStep;
  x2 = x1 + MatrixThick;
  x3 = x2 + SupraThick;

  pm1 = newp ; Point(pm1) = { x1 , y1 , 0 , TapeMeshSize };
  pm2 = newp ; Point(pm2) = { x2 , y1 , 0 , TapeMeshSize };
  pm3 = newp ; Point(pm3) = { x2 , y2 , 0 , TapeMeshSize };
  pm4 = newp ; Point(pm4) = { x1 , y2 , 0 , TapeMeshSize };
  ps1 = pm2;
  ps2 = newp ; Point(ps2) = { x3 , y1 , 0 , TapeMeshSize };
  ps3 = newp ; Point(ps3) = { x3 , y2 , 0 , TapeMeshSize };
  ps4 = pm3;

  lm1 = newl ; Line(lm1) = { pm1 , pm2 };
  lm2 = newl ; Line(lm2) = { pm2 , pm3 };
  lm3 = newl ; Line(lm3) = { pm3 , pm4 };
  lm4 = newl ; Line(lm4) = { pm4 , pm1 };
  llMatrix = newll ; Line Loop(llMatrix) = { lm1 , lm2 , lm3 , lm4 };
  Transfinite Line { -lm1 , lm3 } = ( NbElemThick + 1 );
  Transfinite Line { -lm2 , lm4 } = ( NbElemWidth + 1 )
				    Using Bump ( 1 / BumpAmount ); 
  ls1 = newl ; Line(ls1) = { ps1 , ps2 };
  ls2 = newl ; Line(ls2) = { ps2 , ps3 };
  ls3 = newl ; Line(ls3) = { ps3 , ps4 };
  ls4 = newl ; Line(ls4) = { ps4 , ps1 };
  llSupra = newll ; Line Loop(llSupra) = { ls1 , ls2 , ls3 , ls4 };
  Transfinite Line { -ls1 , ls3 } = ( NbElemThick + 1 );
  Transfinite Line { -ls2 , ls4 } = ( NbElemWidth + 1 )
				    Using Bump ( 1 / BumpAmount ) ;
  
  sMatrix = news ; Plane Surface(sMatrix) = { llMatrix };
  sSupra = news ; Plane Surface(sSupra) = { llSupra };

  Transfinite Surface { sSupra };
  Transfinite Surface { sMatrix };
  Recombine Surface { sMatrix , sSupra };

  Color Red     { Surface{sSupra}; }
  Color Yellow  { Surface{sMatrix}; }

  SurfMatrix[] += sMatrix;
  SurfSupra[] += sSupra;
  LoopMatrix[] += llMatrix;
  LoopSupra[] += llSupra;
  ConductorSkin[] += { lm1 , ls1 , ls2 , ls3 , lm3 , lm4 };

  // Add an air region if it is not the last turn ------------------ //

  If ( iTurn < NbTurns )
    x4 = CoilIntRadius + iTurn * CoilStep;
    
    pai1 = ps2;
    pai2 = newp ; Point(pai2) = { x4 , y1 , 0 , TapeMeshSize };
    pai3 = newp ; Point(pai3) = { x4 , y2 , 0 , TapeMeshSize };
    pai4 = ps3;

    lai1 = newl ; Line(lai1) = { pai1 , pai2 };
    lai2 = newl ; Line(lai2) = { pai2 , pai3 };
    lai3 = newl ; Line(lai3) = { pai3 , pai4 };
    lai4 = newl ; Line(lai4) = { pai4 , pai1 };
    llAirInt = newll ; Line Loop(llAirInt) = { lai1 , lai2 , lai3 , lai4 };
    Transfinite Line { -lai1 , lai3 } = ( NbElemAir + 1 )
                                        Using Bump ( 1 / BumpAmount );
    Transfinite Line { -lai2 , lai4 } = ( NbElemWidth + 1 )
				        Using Bump ( 1 / BumpAmount );

    sAirInt = news ; Plane Surface(sAirInt) = { llAirInt };
    Transfinite Surface { sAirInt };
    Recombine Surface { sAirInt };
    
    Hide { Line {lai1 , lai3}; }
    Color SkyBlue { Surface{sAirInt}; }

    SurfAirInt[] += sAirInt;
    LoopAirInt[] += llAirInt;
  EndIf
EndFor

// Enclosure points ------------------------------------------------ //

// Definition of the enclosure points 
pe0 = newp ; Point(pe0) = { 0            , 0             , 0 , AirMeshSize };
pe1 = newp ; Point(pe1) = { 0            , -AirIntRadius , 0 , AirMeshSize };
pe2 = newp ; Point(pe2) = { AirIntRadius , 0             , 0 , AirMeshSize };
pe3 = newp ; Point(pe3) = { 0            , AirIntRadius  , 0 , AirMeshSize };
pe4 = newp ; Point(pe4) = { 0            , -AirExtRadius , 0 , AirMeshSize };
pe5 = newp ; Point(pe5) = { AirExtRadius , 0             , 0 , AirMeshSize };
pe6 = newp ; Point(pe6) = { 0            , AirExtRadius  , 0 , AirMeshSize };

// Build the lines for the air regions
la1 = newl ; Circle(la1) = { pe1 , pe0 , pe2 };
la2 = newl ; Circle(la2) = { pe2 , pe0 , pe3 };
la3 = newl ; Line(la3)   = { pe3 , pe1 };
llAir = newll ; Line Loop(llAir) = { la1 , la2 , la3 };

// Construction of the pseudo-infinite zone
li1 = newl ; Circle(li1) = { pe4 , pe0 , pe5 };
li2 = newl ; Circle(li2) = { pe5 , pe0 , pe6 };
li3 = newl ; Line(li3)   = { pe6 , pe3 };
li4 = -la2 ; li5 = -la1;
li6 = newl ; Line(li6) = { pe1 , pe4 };
llInf = newll ; Line Loop(llInf) = { li1 , li2 , li3 , li4 , li5 , li6 };

sAir = news ; Plane Surface(sAir) = { llAir , LoopMatrix[] , LoopSupra[] , LoopAirInt[] };
sInf = news ; Plane Surface(sInf) = { llInf };

// Tag the physical regions with the IDs defined in the .par file
Physical Surface(SUPRA)           = SurfSupra[];
Physical Surface(MATRIX)          = SurfMatrix[];
Physical Surface(AIR)             = { sAir , SurfAirInt[] };
Physical Surface(INFINITE_DOMAIN) = { sInf };
Physical Line(CONDUCTOR_SKIN)     = ConductorSkin[];
Physical Line(FAR_LIMIT)          = { li1 , li2 };
Physical Line(SYMMETRY_AXIS)      = { li3 , la3 , li6 };

// Cohomology computation for the H-Phi formulation, this will be
// used in the .pro file to enforce constraints
Cohomology(1) { { AIR , INFINITE_DOMAIN },
                { } };

// Some colors to help read the mesh
Color SkyBlue { Surface{sAir}; }
Color RoyalBlue { Surface{sInf}; }

// Check the mesh coherence
Coherence;
