Include "3D_single_pancake_emag_parameters.geo" ;

r_int = SP_INT_RADIUS ;
r_ext = r_int + SP_NB_OF_TURNS * SP_COND_THICK ;
r_air = AIR_RADIUS ;
r_inf = INF_RADIUS ;
z_top = SP_COND_WIDTH / 2. ;
t_int = SP_COND_WIDTH / ( 2. * r_int ) ;
t_ext = SP_COND_WIDTH / ( 2. * r_ext ) ;
l_adi = ADI_LENGTH ;
e_adi = ADI_THICK ;
lc_sp = ( r_ext - r_int ) / 3. ;
lc_hdl = e_adi ;
lc_air = ( r_inf - r_air ) / 2. ;

p0 = newp ;
Point( p0   ) = {  0., 0., -z_top } ;

p_int = newp ;
Point( p_int   ) = {  r_int * Cos( -t_int ),  r_int * Sin( -t_int ), -z_top, lc_sp } ;
Point( p_int+1 ) = {  r_int * Cos(  t_int ),  r_int * Sin(  t_int ), -z_top, lc_sp } ;
Point( p_int+2 ) = {  0.                   ,  r_int                , -z_top, lc_sp } ;
Point( p_int+3 ) = { -r_int                ,  0.                   , -z_top, lc_sp } ;
Point( p_int+4 ) = {  0.                   , -r_int                , -z_top, lc_sp } ;
p_adi_int = newp ;
Point( p_adi_int   ) = { r_int * Cos( -t_int ) - e_adi, r_int * Sin( -t_int ), -z_top, lc_sp } ;
Point( p_adi_int+1 ) = { r_int * Cos(  t_int ) - e_adi, r_int * Sin(  t_int ), -z_top, lc_sp } ;
p_ext = newp ;
Point( p_ext   ) = { r_ext                    ,  0.                       , -z_top, lc_sp } ;
Point( p_ext+1 ) = { 0.                       ,  r_ext                    , -z_top, lc_sp } ;
Point( p_ext+2 ) = { r_ext * Cos( Pi - t_ext ),  r_ext * Sin( Pi - t_ext ), -z_top, lc_sp } ;
Point( p_ext+3 ) = { r_ext * Cos( Pi + t_ext ),  r_ext * Sin( Pi + t_ext ), -z_top, lc_sp } ;
Point( p_ext+4 ) = { 0.                       , -r_ext                    , -z_top, lc_sp } ;
p_adi_ext = newp ;
Point( p_adi_ext   ) = { -r_ext - e_adi, r_ext * Sin( Pi - t_ext ), -z_top, lc_sp } ;
Point( p_adi_ext+1 ) = { -r_ext - e_adi, r_ext * Sin( Pi + t_ext ), -z_top, lc_sp } ;

l_int = newl ;
Circle( l_int   ) = { p_int  , p0, p_int+1 } ;
Circle( l_int+1 ) = { p_int+1, p0, p_int+2 } ;
Circle( l_int+2 ) = { p_int+2, p0, p_int+3 } ;
Circle( l_int+3 ) = { p_int+3, p0, p_int+4 } ;
Circle( l_int+4 ) = { p_int+4, p0, p_int   } ;
l_adi_int = newl ;
Line( l_adi_int   ) = { p_int      , p_adi_int   } ;
Line( l_adi_int+1 ) = { p_adi_int  , p_adi_int+1 } ;
Line( l_adi_int+2 ) = { p_adi_int+1, p_int+1     } ;
l_ext = newl ;
Circle( l_ext   ) = { p_ext  , p0, p_ext+1 } ;
Circle( l_ext+1 ) = { p_ext+1, p0, p_ext+2 } ;
Circle( l_ext+2 ) = { p_ext+2, p0, p_ext+3 } ;
Circle( l_ext+3 ) = { p_ext+3, p0, p_ext+4 } ;
Circle( l_ext+4 ) = { p_ext+4, p0, p_ext   } ;
l_adi_ext = newl ;
Line( l_adi_ext   ) = { p_ext+2    , p_adi_ext   } ;
Line( l_adi_ext+1 ) = { p_adi_ext  , p_adi_ext+1 } ;
Line( l_adi_ext+2 ) = { p_adi_ext+1, p_ext+3     } ;

ll_int = newll ;
Line Loop( ll_int ) = { l_int:l_int+4 } ;
ll_adi_int = newll ;
Line Loop( ll_adi_int ) = { l_adi_int:l_adi_int+2, -l_int } ;
ll_ext = newll ;
Line Loop( ll_ext ) = { l_ext:l_ext+4 } ;
ll_adi_ext = newll ;
Line Loop( ll_adi_ext ) = { l_adi_ext:l_adi_ext+2, -(l_ext+2) } ;

s_sp = news ;
Plane Surface( s_sp ) = { ll_ext, ll_int } ;
s_adi_int = news ;
Plane Surface( s_adi_int ) = { ll_adi_int } ;
s_adi_ext = news ;
Plane Surface( s_adi_ext ) = { ll_adi_ext } ;
s_cond_domain_bnd[] = { s_sp, s_adi_int, s_adi_ext } ;

tmp[] = Extrude { 0., 0., 2. * z_top } { Surface{ s_sp, s_adi_int, s_adi_ext } ; } ;
v_sp = tmp[1] ;
v_inner_cl[] = { tmp[13] } ;
v_outer_cl[] = { tmp[19] } ;
s_junct_inner_cl = tmp[11] ;
s_junct_outer_cl = tmp[4] ;
s_cond_domain_bnd[] += { tmp[0] , tmp[2] , tmp[3] , tmp[5] , tmp[6]  , 
                         tmp[7] , tmp[8] , tmp[9] , tmp[10], tmp[14] ,
			 tmp[15], tmp[16], tmp[20], tmp[21], tmp[22] } ;
tmp[] = Extrude { 0., 0., l_adi - 2. * e_adi } { Surface{ tmp[12], tmp[18] } ; } ;
v_inner_cl[] += { tmp[ 1 ] } ;
v_outer_cl[] += { tmp[ 7 ] } ;
s_cond_domain_bnd[] += { tmp[2], tmp[3], tmp[4] , tmp[5]  ,
                         tmp[8], tmp[9], tmp[10], tmp[11] } ;
tmp[] = Extrude { 0., 0., e_adi } { Surface{ tmp[0], tmp[6] } ; } ;
v_inner_cl[] += { tmp[ 1 ] } ;
v_outer_cl[] += { tmp[ 7 ] } ;
s_cond_domain_bnd[] += { tmp[0], tmp[2], tmp[4] , tmp[5]  ,
                         tmp[6], tmp[8], tmp[9], tmp[10] } ;

l_br = newl ;
Line( l_br   ) = { 94 , 72  } ;
Line( l_br+1 ) = { 76 , 85  } ;
Line( l_br+2 ) = { 118, 96  } ;
Line( l_br+3 ) = { 100, 109 } ;
ll_br = newll ;
Line Loop( ll_br   ) = { l_br, 171, -(l_br+2), -201 } ;
Line Loop( ll_br+1 ) = { l_br+2, 166, l_br+3, -190 } ;
Line Loop( ll_br+2 ) = { l_br+1, 192, -(l_br+3), -175 } ;
Line Loop( ll_br+3 ) = { -l_br, 146, -(l_br+1), -122 } ;
s_br = news ;
Plane Surface( s_br   ) = { ll_br   } ;
Plane Surface( s_br+1 ) = { ll_br+1 } ;
Plane Surface( s_br+2 ) = { ll_br+2 } ;
Plane Surface( s_br+3 ) = { ll_br+3 } ;
sl_br = news ;
Surface Loop( sl_br ) = { s_br:s_br+3, 176, -206 } ;
v_br = newv ;
Volume( v_br ) = { sl_br } ;
s_cond_domain_bnd[] += { s_br:s_br+3 } ;

p_air = newp ;
Point( p_air   ) = { 0.   , 0.,  0.   , lc_air } ;
Point( p_air+1 ) = { 0.   , 0.,  r_air, lc_air } ;
Point( p_air+2 ) = { r_air, 0.,  0.   , lc_air } ;
Point( p_air+3 ) = { 0.   , 0., -r_air, lc_air } ;
p_inf = newp ;
Point( p_inf   ) = { 0.   , 0.,  r_inf, lc_air } ;
Point( p_inf+1 ) = { r_inf, 0.,  0.   , lc_air } ;
Point( p_inf+2 ) = { 0.   , 0., -r_inf, lc_air } ;
l_air = newl ;
Circle( l_air   ) = { p_air+1, p_air, p_air+2 } ;
Circle( l_air+1 ) = { p_air+2, p_air, p_air+3 } ;
l_inf = newl ;
Circle( l_inf   ) = { p_inf  , p_air, p_inf+1 } ;
Circle( l_inf+1 ) = { p_inf+1, p_air, p_inf+2 } ;
tmp[] = Extrude { { 0., 0., 1. }, { 0., 0., 0. }, Pi / 2. } { Line{ l_air, l_air+1, l_inf, l_inf+1 } ; } ;
s_air[] = { tmp[1], tmp[ 4] } ;
s_inf[] = { tmp[7], tmp[10] } ;
tmp[] = Extrude { { 0., 0., 1. }, { 0., 0., 0. }, Pi / 2. } { Line{ tmp[0], tmp[3], tmp[6], tmp[9] } ; } ;
s_air[] += { tmp[1], tmp[ 4] } ;
s_inf[] += { tmp[7], tmp[10] } ;
tmp[] = Extrude { { 0., 0., 1. }, { 0., 0., 0. }, Pi / 2. } { Line{ tmp[0], tmp[3], tmp[6], tmp[9] } ; } ;
s_air[] += { tmp[1], tmp[ 4] } ;
s_inf[] += { tmp[7], tmp[10] } ;
tmp[] = Extrude { { 0., 0., 1. }, { 0., 0., 0. }, Pi / 2. } { Line{ tmp[0], tmp[3], tmp[6], tmp[9] } ; } ;
s_air[] += { tmp[1], tmp[ 4] } ;
s_inf[] += { tmp[7], tmp[10] } ;
ss_air = news ;
Surface Loop( ss_air ) = { s_air[], s_cond_domain_bnd[] } ;
ss_inf = news ;
Surface Loop( ss_inf ) = { s_inf[], s_air[] } ;
v_air = newv ;
Volume( v_air ) = { ss_air } ;
v_inf = newv ;
Volume( v_inf ) = { ss_inf } ;

// Assign physical tags to important regions of the geometrical model
Physical Point( "reference Phi point", TAG_REF_PHI_PT ) = { p_inf+2 } ;
Physical Surface( "inner current lead junction", TAG_JUNCT_INNER_CL ) = { s_junct_inner_cl } ;
Physical Surface( "outer current lead junction", TAG_JUNCT_OUTER_CL ) = { s_junct_outer_cl } ;
Physical Surface( "conducting domain boundary", TAG_COND_DOMAIN_BND ) = { s_cond_domain_bnd[] } ;
Physical Surface( "domain boundary", TAG_DOMAIN_BND ) = { s_inf[] } ;
Physical Volume( "superconducting coil", TAG_SC_COIL ) = { v_sp } ;
Physical Volume( "inner current lead", TAG_INNER_CL ) = { v_inner_cl[] } ;
Physical Volume( "outer current lead", TAG_OUTER_CL ) = { v_outer_cl[] } ;
Physical Volume( "bridge", TAG_BRIDGE ) = { v_br } ;
Physical Volume( "air", TAG_AIR ) = { v_air } ;
Physical Volume( "pseudo-infinite air", TAG_INFINITE_AIR ) = { v_inf } ;

// Assign some colors to the mesh elements to make the model more readable
Recursive Color SkyBlue { Volume{ v_air, v_inf } ; }
Recursive Color Orange { Volume{ v_br, v_inner_cl[], v_outer_cl[] } ; }
Recursive Color Red { Volume{ v_sp } ; }
Recursive Color Grey { Surface{ s_junct_inner_cl, s_junct_outer_cl } ; }

Recursive Hide { Volume{ v_inf } ; }
Hide { Volume{ v_air } ; }

// Compute cohomology of the air domain.
Cohomology(1) { { TAG_AIR, TAG_INFINITE_AIR }, { } } ;
