Geometrie Poly-Helices

trophime christophe trophime at labs.polycnrs-gre.fr
Mon Apr 2 17:03:24 CEST 2001


Salut,
come convenu je t'envoie en attachment le fichier de geometrie
d'un "insert polyhelice". Le maillage de l'ensemble des helices
et des bagues est fait en quadrangles. Le maillage du domaine
exterieur a ete decoupe en plusieurs parties (Cooling Channels, ...).
Ce sont ces surfaces que je n'arrive pas a mailler. Dans le fichier
toutes ces surfaces sont mises en commentaires sauf celle qui
se trouve entre l'axe Oz et la premiere helice (i.e 149).

Encore merci pour ton aide et a bientot.
Je t'envoie une liste des points que l'on aimerait aborder.

Christophe Trophime
-------------- next part --------------
/*******************************************************

Test de maillage d'un insert poly-helices
Parametrage a revoir pour eviter probleme....

Ch. Trophime 22/03/2001
*******************************************************/
weight = 0.01;

r0 = 2;       // rayon interne 1ere helice
alpha = 0;    // angle de la decoupe
glue = 0.2;   // epaisseur colle
h_ring = 1;
r_ext = 30;
r_inf = 35;

N_helices = 6;

N_turns[1] = 3;
N_turns[2] = 4;
N_turns[3] = 5;
N_turns[4] = 5;
N_turns[5] = 5;
N_turns[6] = 5;

epaisseur[1] = 1; // epaisseur helice
epaisseur[2] = 2;
epaisseur[3] = 3;
epaisseur[4] = 3;
epaisseur[5] = 3;
epaisseur[6] = 3;

hauteur[1] = 7; // hauteur helice
hauteur[2] = 8;
hauteur[3] = 9;
hauteur[4] = 10;
hauteur[5] = 11;
hauteur[6] = 12;

canal[1] = 1; // epaisseur canal hydraulique entre les helices i et i+1
canal[2] = 1;
canal[3] = 1;   
canal[4] = 1;   // a ajouter pour faire fonctionner correctement
canal[5] = 1;   
canal[6] = 0.;   // a ajouter pour faire fonctionner correctement

z0[1] = -(hauteur[1]/2.0+1); // ordonnees de la base de  helice
z0[2] = -(hauteur[2]/2.0+2);
z0[3] = z0[2];  
z0[4] = -(hauteur[4]/2.0+3);
z0[5] = z0[4];  
z0[6] = -(hauteur[6]/2.0+4);

z1[1] = (hauteur[1]/2.0+1); // ordonnees du sommet de  helice
z1[2] = z1[1];
z1[3] = (hauteur[3]/2.0+2);  
z1[4] = z1[3];
z1[5] = (hauteur[5]/2.0+3);  
z1[6] = z1[5];

/* Discretisation Parameters
  nr : nombre de noeuds dans epaisseur
  nh : ..................... hauteur par tour
 */

nr = 10;
nh = 10;
nh_glue = 2;
nh_end = 5;
nh_ring = 5;
nr_channel = 6;

n_circum = 40;
nr_inf =10;
n_axe = 40;
n_bore = 20;

Function Helix_Turn

  p1 = newp; Point(p1) = {x, y, 0, weight};  
  p2 = newp; Point(p2) = {x+r, y+r*Tan(alpha), 0, weight};
  p3 = newp; Point(p3) = {x+r, y+h+r*Tan(alpha), 0, weight};  
  p4 = newp; Point(p4) = {x, y+h, 0, weight};  

  c1 = newreg; Line(c1) = {p1,p2};  
  c2 = newreg; Line(c2) = {p2,p3};  
  c3 = newreg; Line(c3) = {p3,p4};  
  c4 = newreg; Line(c4) = {p4,p1};
  
  Transfinite Line{c1,-c3} = nr;  
  Transfinite Line{c2,-c4} = nh;

  l1 = newreg; Line Loop(l1) = {c1,c2,c3,c4};
    
  helix_turn = newreg; Plane Surface(helix_turn) = {l1}; 
  
  Transfinite Surface{helix_turn} = {p1,p2,p3,p4};
  Recombine Surface{helix_turn};
  
Return

Function Glue_Turn
  gl1 = newreg; Line(gl1) = {num_p1-1, num_p1+2};  
  gl2 = newreg; Line(gl2) = {num_p1, num_p1+1};  

  Transfinite Line{gl1,-gl2} = nh_glue;

  g1 = newreg; Line Loop(g1) = {-num_l1,gl1,-(num_l1+4),-gl2};
  //Printf("Lignes : %g %g %g %g",num_l1,gl1,-(num_l1+4),-gl2);
  glue_turn = newreg; Plane Surface(glue_turn) = {g1};
  
  Transfinite Surface{glue_turn} = {num_p1,num_p1-1,num_p1+2,num_p1+1};
  Recombine Surface{glue_turn};
  
Return

Function Helix_EndTurn
  p1 = newp; Point(p1) = {x, y, 0, weight};  
  p2 = newp; Point(p2) = {x+r, y, 0, weight};
  Helix_EndTurn_P1[num] = p1; 
  Helix_EndTurn_P2[num] = p2; 

  l1 = newreg; Line(l1) = {p1, p2};              //Printf("Ligne(%g) : %g %g",l1,p1,p2);  
  l2 = newreg; Line(l2) = {p2, num_p1+signe};    //Printf("Ligne(%g) : %g %g",l2,p2, num_p1+signe);  
  l3 = newreg; Line(l3) = {num_p1, p1};          //Printf("Ligne(%g) : %g %g",l3,num_p1, p1);  
  Helix_EndTurn_L1[num] = l1;
  
  Transfinite Line{signe*l1} = nr;
  Transfinite Line{signe*l2} = nh_end;
  Transfinite Line{-signe*l3} = nh_end;

  g1 = newreg; Line Loop(g1) = {l1,l2,-signe*num_l1,l3};
  //Printf("Lignes : %g %g %g %g",l1,l2,-signe*num_l1,l3);
  end_turn = newreg; Plane Surface(end_turn) = {g1};
  
  Transfinite Surface{end_turn} = {p1,p2,num_p1+signe,num_p1};
  //Printf("EndTurn : %g %g  %g %g",p1,p2,num_p1+signe,num_p1); 
  Recombine Surface{end_turn};
  
  num += 1;
Return

Function Ring
  p1 = newp; Point(p1) = {x, y+signe*h_ring, 0, weight};  
  p2 = newp; Point(p2) = {x+epaisseur[num], y+signe*h_ring, 0, weight};
  p3 = newp; Point(p3) = {x+epaisseur[num]+canal[num], y+signe*h_ring, 0, weight};
  p4 = newp; Point(p4) = {x+epaisseur[num]+canal[num]+epaisseur[num+1], y+signe*h_ring, 0, weight};

  l1 = newreg; Line(l1) = {p1, p2};  
  l2 = newreg; Line(l2) = {p2, p3};
  l3 = newreg; Line(l3) = {p3, p4}; 

  hl1 = newreg; Line(hl1) = {num_p1, p1};  
  hl2 = newreg; Line(hl2) = {num_p2, p2};
  hl3 = newreg; Line(hl3) = {num_p3, p3};
  hl4 = newreg; Line(hl4) = {p4, num_p4}; 
  hl5 = newreg; Line(hl5) = {num_p2, num_p3}; 
  
  Transfinite Line{signe*hl2} = nh_ring;
  Transfinite Line{-signe*l1} = nr;
  Transfinite Line{-signe*hl1} = nh_ring;
  
  Transfinite Line{signe*hl5} = nr_channel;
  Transfinite Line{signe*hl3} = nh_ring;
  Transfinite Line{-signe*l2} = nr_channel;
  Transfinite Line{-signe*hl2} = nh_ring;

  Transfinite Line{-signe*hl4} = nh_ring;
  Transfinite Line{-signe*l3} = nr;
  Transfinite Line{signe*hl3} = nh_ring;
  
  r1 = newreg; Line Loop(r1) = {signe*num_l1,signe*hl2,-signe*l1,-signe*hl1};
  r2 = newreg; Line Loop(r2) = {signe*hl5,signe*hl3,-signe*l2,-signe*hl2};
  r3 = newreg; Line Loop(r3) = {signe*num_l2,-signe*hl4,-signe*l3,-signe*hl3};
  
  //Printf("r1: %g(%g-%g), %g(%g-%g), %g(%g-%g), %g(%g-%g)",
  //   num_l1,num_p1, num_p2,
  //   signe*hl2,num_p2, p2,
  //   -signe*l1,p1, p2,
  //   -signe*hl1,num_p1, p1);
     
  r1_ring = newreg; Plane Surface(r1_ring) = {r1};
  Transfinite Surface{r1_ring} = {num_p1, num_p2, p2, p1};
  r2_ring = newreg; Plane Surface(r2_ring) = {r2};
  Transfinite Surface{r2_ring} = {num_p2, num_p3, p3, p2};
  r3_ring = newreg; Plane Surface(r3_ring) = {r3};
  Transfinite Surface{r3_ring} = {num_p3, num_p4, p4, p3};
  Recombine Surface{r1_ring, r2_ring, r3_ring};
  
Return

Function Cooling_Channel
  p1 = newp; Point(p1) = {x, y, 0, weight};
  
  l1 = newreg; Line(l1) = {num_p1, p1}; //Printf("L1: %g %g", num_p1, num_p2);
  l2 = newreg; Line(l2) = {p1, num_p2}; //Printf("L2: %g %g", num_p1, num_p2);


  Transfinite Line{signe*l1} = nh_ring; //Using Progression 2;
  Transfinite Line{signe*l2} = nr_channel; //Using Progression 2;
//   If (num % 2 != 0)
//     Transfinite Line{signe*l1} = nr_channel; //Using Progression 2;
//     Transfinite Line{signe*l2} = nh_ring; //Using Progression 2;
//   EndIf
  
Return

Function Exterior
   origin = newp; Point(origin) = {0,0,0, weight};
   p1 = newp; Point(p1) = {0, -r_ext, 0, weight};
   p2 = newp; Point(p2) = {0, -r_inf, 0, weight};
   //p3 = newp; Point(p3) = {r_ext, 0, 0, weight};
   //p4 = newp; Point(p4) = {r_inf, 0, 0, weight};
   p5 = newp; Point(p5) = {0, r_ext, 0, weight};
   p6 = newp; Point(p6) = {0, r_inf, 0, weight};
   p7 = newp; Point(p7) = {0, z0[1], 0, weight};
   p10 = newp; Point(p10) = {0, z1[1]+h_ring, 0, weight};
   
   x = r0;
   y = z1[N_helices] + h_ring;
   For i In{1:N_helices}
      r = epaisseur[i];
      x += r + canal[i];
   EndFor
   theta = Atan(y/x);
   
   p8 = newp; Point(p8) = {r_ext * Cos(theta), -r_ext * Sin(theta), 0, weight};
   p11 = newp; Point(p11) = {r_inf * Cos(theta), -r_inf * Sin(theta), 0, weight};


   y = z0[N_helices];
   p9 = newp; Point(p9) = {r_ext * Cos(theta), r_ext * Sin(theta), 0, weight};
   p12 = newp; Point(p12) = {r_inf * Cos(theta),  r_inf * Sin(theta), 0, weight};

   c1 = newreg; Circle(c1) = {p1, origin, p8};   
   c2 = newreg; Circle(c2) = {p2, origin, p11};   
   c3 = newreg; Circle(c3) = {p8, origin, p9};   
   c4 = newreg; Circle(c4) = {p11, origin, p12};   
   c5 = newreg; Circle(c5) = {p9, origin, p5};   
   c6 = newreg; Circle(c6) = {p12, origin, p6};   
   l1 = newreg; Line(l1) = {p2, p1};
   l2 = newreg; Line(l2) = {p11, p8};
   l3 = newreg; Line(l3) = {p9, p12};
   l4 = newreg; Line(l4) = {p5, p6};
   l5 = newreg; Line(l5) = {p10,p7};
   
   l6 = newreg; Line(l6) = {Helix_EndTurn_P1[1],p7};
   l7 = newreg; Line(l7) = {Helix_EndTurn_P2[2*N_helices-1],p8};
   l8 = newreg; Line(l8) = {p7,p1};
   
   l9 = newreg; Line(l9) = {Ring_P1[1],p10};
   l10 = newreg; Line(l10) = {Ring_P4[N_helices-1],p9};
   l11 = newreg; Line(l11) = {p5,p10};

   nbr_nodes_on_last = (nh_end-1) * 2 + 
                       N_turns[N_helices] * (nh-1) + (N_turns[N_helices] - 1) * (nh_glue-1) +
		       nh_ring;
		        
   nbr_nodes_on_first = (nh_end-1) * 2 + 
                       N_turns[1] * (nh-1) + (N_turns[1] - 1) * (nh_glue-1) +
		       nh_ring;
		       
   Transfinite Line{c1} = n_circum;
   Transfinite Line{c2} = n_circum;
   Transfinite Line{c3} = nbr_nodes_on_last;// nbre node on last helice
   Transfinite Line{c4} = nbr_nodes_on_last;// nbre node on last helice
   Transfinite Line{c5} = n_circum;
   Transfinite Line{c6} = n_circum;
   
   Transfinite Line{l1} = nr_inf;
   Transfinite Line{l2} = nr_inf;
   Transfinite Line{l3} = nr_inf;
   Transfinite Line{l4} = nr_inf;
   Transfinite Line{l5} = nbr_nodes_on_first;      // nbr de noeuds sur premiere helice
   
   Transfinite Line{l6} = n_bore;
   Transfinite Line{l7} = n_axe;
   Transfinite Line{l8} = n_axe;

   Transfinite Line{l9} = n_bore;
   Transfinite Line{l10} = n_axe;
   Transfinite Line{l11} = n_axe;

   infini_1 = newreg; Line Loop(infini_1) = {-l1, c2, l2, -c1}; //Printf("infini_1: %g %g %g %g", -l1, c2, l2, -c1);
   infini_2 = newreg; Line Loop(infini_2) = {-l2, c4, -l3, -c3}; //Printf("infini_2: %g %g %g %g",-l2, c4, -l3, -c3);
   infini_3 = newreg; Line Loop(infini_3) = {l3, c6, -l4, -c5}; //Printf("infini_3: %g %g %g %g",l3, c6, -l4, -c5);
      
   inf_1 =  newreg; Plane Surface(inf_1) = {infini_1};
   inf_2 =  newreg; Plane Surface(inf_2) = {infini_2};
   inf_3 =  newreg; Plane Surface(inf_3) = {infini_3};
   
   Transfinite Surface{inf_1} = {p1, p2, p11, p8};
   Transfinite Surface{inf_2} = {p8, p11, p12, p9};
   Transfinite Surface{inf_3} = {p9, p12, p6, p5};
   Recombine Surface{inf_1,inf_2,inf_3};
   
Return

x = r0;
num = 1;
num_int = 2;
num_ext = 4;
num_tot = 0;

Printf("Built Helices");
For i In{1:N_helices}
 y = -hauteur[i]/2.0;
 r = epaisseur[i];
 h = (hauteur[i]-glue*(N_turns[i]-1))/N_turns[i];
 
 Printf("Helix %g : ",i); // Printf("(%g)", h);
 For j In {1:N_turns[i]}

  Call Helix_Turn ;
  //Printf(" %g, ",j);
  Printf("internal[%g] : %g external[%g] : %g ", i, num_int, i, num_ext); 
  num_int += 6;
  num_ext += 6;
  
  y += h + glue;
  num += 1;
 EndFor
 x += r + canal[i];
 num_tot += 6*N_turns[i];
 //Printf("): %g  total : %g",6*N_turns[i],num_tot);
 //Physical Volume (i) = {num} ;
EndFor

num_p1 = 0;
num_l1 = 0;
Printf("Built Glue");
For i In{1:N_helices}
 num_p1 += 4;
 num_l1 += 3;
 //Printf("Built Glue %g : (Point P1:%g) (Ligne L1:%g)",i,num_p1,num_l1);
 
 For j In {1:N_turns[i]-1}

  Call Glue_Turn ;

  //Physical Volume (num) = helix_turn ;
  //Printf("[%g,%g] ",j,j+1);
  //Printf("internal : %g external : %g ",num_tot+1, num_tot+2); num_tot += 4;
 
  num_p1 += 4;
  num_l1 += 6;
 EndFor
 num_l1 += 3;//because of 3 calls to newreg after line creation in Helix_Turn
 
 //Printf(" ");
EndFor

x = r0;
num_first = 1;
num_last = 0;
num_l1 = 1;
num =1;

Printf("Build End Turns");
For i In{1:N_helices}
 r = epaisseur[i];
 num_points = 4*N_turns[i];
 num_last += num_points;
 
 y = z0[i];
 signe = 1;
 num_p1 = num_first;
 l2 = 0;
 Printf("Built Helix Lower End Turns [%g] : ",i);
 //Printf("num_first %g %g: ",num_first, num_last);
 Call Helix_EndTurn ;
 //Printf("internal: %g external: %g bottom : %g", l2, l3, l1);
 Helix_Lower_EndTurn_Int[i] = signe*l2;
 Helix_Lower_EndTurn_Ext[i] = signe*l3;
 Helix_Bottom[i] = signe*l1;
 
 y = z1[i];
 signe = -1;
 num_p1 = num_last;
 num_l1 += 4*N_turns[i]+2*(N_turns[i]-1)-2;
 l2 = 0;
 Printf("Built Helix Upper End Turns [%g] : ",i);
 Call Helix_EndTurn ;
 //Printf("internal: %g external: %g", l2, l3);
 Helix_Upper_EndTurn_Int[i] = signe*l2;
 Helix_Upper_EndTurn_Ext[i] = signe*l3;
 
 x += r + canal[i];
 num_l1 += 4;
 
 num_first = num_last+1; 
EndFor

x = r0;

num_helix = 2;
signe = 1;

Printf("Built Upper Rings");
For i In{1:N_helices/2}
 num = 1+2*(i-1); 
 y = z1[num];
 
 num_p1 = Helix_EndTurn_P1[num_helix];
 num_p2 = Helix_EndTurn_P2[num_helix];
 num_p3 = Helix_EndTurn_P1[num_helix+2];
 num_p4 = Helix_EndTurn_P2[num_helix+2];
 
 num_l1 = Helix_EndTurn_L1[num_helix];
 num_l2 = Helix_EndTurn_L1[num_helix+2];

 Printf("Built Ring [%g-%g] : ",num,num+1);
 Call Ring ;
 //Printf("internal: %g external: %g %g %g %g %g", hl5, hl1,l1,l2,l3,hl4);
 Ring_P1[num] = p1;
 Ring_P4[num] = p4;
 //Printf("Ring[%g] P1=%g P4=%g", num, p1, p4);

 Ring_Int[num] = hl5;
 Ring_Ext_L1[num] = hl1;
 Ring_Ext_L2[num] = l1;
 Ring_Ext_L3[num] = l2;
 Ring_Ext_L4[num] = l3;
 Ring_Ext_L5[num] = hl4;
 
 x += epaisseur[num] + canal[num];
 x += epaisseur[num+1] + canal[num+1];
 num_helix += 4;
EndFor

x = r0 + epaisseur[1] + canal[1];

signe = -1;
num_helix = 3;
Printf("Built Lower Rings");
For i In{1:N_helices/2-1}
 num = 2*i;
 y = z0[num];
 
 num_p1 = Helix_EndTurn_P1[num_helix];
 num_p2 = Helix_EndTurn_P2[num_helix];
 num_p3 = Helix_EndTurn_P1[num_helix+2];
 num_p4 = Helix_EndTurn_P2[num_helix+2];

 num_l1 = Helix_EndTurn_L1[num_helix];
 num_l2 = Helix_EndTurn_L1[num_helix+2];

 Printf("Built Ring [%g-%g] : ",num,num+1);
 Call Ring ;
 //Printf("internal: %g external: %g %g %g %g %g", hl5, hl1,l1,l2,l3,hl4);
 Ring_P1[num] = p1;
 Ring_P4[num] = p4;
 //Printf("Ring[%g] P1=%g P4=%g", num, p1, p4);
 
 Ring_Int[num] = hl5;
 Ring_Ext_L1[num] = hl1;
 Ring_Ext_L2[num] = l1;
 Ring_Ext_L3[num] = l2;
 Ring_Ext_L4[num] = l3;
 Ring_Ext_L5[num] = hl4;

 x += epaisseur[num] + canal[num];
 x += epaisseur[num+1] + canal[num+1];
 num_helix += 4;
EndFor

Printf("Building Cooling Channels");

signe = 1;
cpte = 1;
x_num = r0;

Printf("Built Lines");
num_ring = 1;
num_endturn = 1; 
For num In{1:N_helices-1}

 x_num += epaisseur[num];
 x_next = x_num + canal[num];
 
 //Printf("Built Cooling Channel [%g-%g] Ring[%g]: %g",num,num+1, num_ring, l1);
 If (num % 2 != 0) 
   y = (z0[num]-h_ring <= z0[num+1]-h_ring) ? z0[num]-h_ring : z0[num+1]-h_ring;
   x = (z0[num]-h_ring <= z0[num+1]-h_ring) ? x_next  : x_num;
   If (num == 1)
     y = (z0[num] <= z0[num+1]-h_ring) ? z0[num] : z0[num+1]-h_ring;
     x = (z0[num] <= z0[num+1]-h_ring) ? x_next  : x_num;
   EndIf
   If (num == N_helices-1)
   y = (z0[num]-h_ring <= z0[num+1]) ? z0[num]-h_ring : z0[num+1];
   x = (z0[num]-h_ring <= z0[num+1]) ? x_next  : x_num;
   EndIf
 EndIf
 If (num % 2 == 0) 
   y = (z1[num]+h_ring <= z1[num+1]+h_ring) ? z1[num+1]+h_ring : z1[num]+h_ring;
   x = (z1[num]+h_ring <= z1[num+1]+h_ring) ? x_num  : x_next;
 EndIf
 
 If ( num_ring > 1 && num_ring <= N_helices-2)
   num_p1 = Ring_P4[num_ring-1];
   num_p2 = Ring_P1[num_ring+1];
 EndIf
 
 If ( num == N_helices-1)
   num_p2 = Helix_EndTurn_P1[num_endturn+2];
   num_p1 = Ring_P4[num_ring-1];
 EndIf
 If ( num == 1)
   num_p1 = Helix_EndTurn_P2[num_endturn];
   num_p2 = Ring_P1[num_ring+1];
 EndIf
 
 Call Cooling_Channel;
 //Printf("Line(%g): %g %g", l1, num_p1, num_p2);
 Channel_L1[num_ring] = l1;
 Channel_L2[num_ring] = l2;
 If (num == N_helices -1)
   Channel_L1[num_ring] = l2;
   Channel_L2[num_ring] = l1;
 EndIf
 
 num_endturn += 2;
 cpte += 1;
 If (cpte ==2)
   num_ring += 1;
   cpte = 1;
 EndIf
 x_num = x_next;
EndFor

Printf("Create Surfaces");
 
num_tot = 0;
For i In{1:N_helices}
  num_tot += 6*N_turns[i];
EndFor
//Printf("num_tot= %g",num_tot);

orient = 1;
num_helix = 0;
num_channel = 1;
n = 0;
For i In{1:N_helices}
  If ( i % 2 != 0) 
    //Printf("Odd Helix (%g)", i);
    //Printf("Channel_L1=%g", -Channel_L1[num_channel]); 
    Channel_Loop[n] = -Channel_L1[num_channel]; n+=1;
    //Printf("Channel_L2=%g", -Channel_L2[num_channel]); 
    Channel_Loop[n] = -Channel_L2[num_channel]; n+=1;
    If ( i > 1)
       //Printf(" %g", Ring_Ext_L5[i-1]);
       Channel_Loop[n] = Ring_Ext_L5[i-1]; n += 1;
    EndIf
    //Printf(" %g ",Helix_Lower_EndTurn_Int[i]);
    Channel_Loop[n] = Helix_Lower_EndTurn_Int[i]; n+=1;
    num_int = num_helix + 2;
    For j In {1:N_turns[i]-1}
       //Printf(" %g ",num_int);
       Channel_Loop[n] = num_int; n+=1;
       num_int += 6;
       //Printf(" %g",num_tot+1); 
       Channel_Loop[n] = num_tot+1; n+=1;
       num_tot += 4;
    EndFor
    //Printf(" %g ",num_int);
    Channel_Loop[n] = num_int; n+=1;
    num_int += 6;
    //Printf(" %g",Helix_Upper_EndTurn_Int[i]);
    Channel_Loop[n] = Helix_Upper_EndTurn_Int[i]; n+=1;
    //Printf(" %g",Ring_Int[i]);
    Channel_Loop[n] = Ring_Int[i]; n+=1;
    num_helix += 6*N_turns[i];

    num_channel += 1;
    
  EndIf
  If ( i % 2 == 0) 
    
    //
    // Odd Chanel Loop
    
    num_ext = num_helix + 6*N_turns[i] - 2;
    num_tot += 4*(N_turns[i]-2);
    //Printf("Even Helix(%g)", i);
    //Printf(" %g ",Helix_Upper_EndTurn_Ext[i]);
    Channel_Loop[n] = Helix_Upper_EndTurn_Ext[i]; n+=1;
    For j In {1:N_turns[i]-1}
       //Printf(" %g",num_ext); 
       Channel_Loop[n] = num_ext; n+=1;
       num_ext -= 6;
       //Printf(" %g",-(num_tot+2));
       Channel_Loop[n] = -(num_tot+2); n+=1;
       num_tot -= 4;
    EndFor
    //Printf(" %g ",num_ext);
    Channel_Loop[n] = num_ext; n+=1;
    num_ext -= 6;
    //Printf("Lower_EndTurn=%g",Helix_Lower_EndTurn_Ext[i]);
    Channel_Loop[n] = Helix_Lower_EndTurn_Ext[i]; n+=1;
    If ( i <= N_helices-1)
      //Printf("Lower_Ring=%g",Ring_Ext_L1[i]);//??
      Channel_Loop[n] = Ring_Ext_L1[i]; n+=1;
    EndIf
    num_helix += 6*N_turns[i];
    num_tot += 4*N_turns[i];

    num_channel += 1;
    // reverse Order of Channel_Loop
    For m In{0:n-1}
      Reverse_Channel_Loop[m] = -Channel_Loop[n-1-m];
    EndFor

    Channel_Boundary = newreg; Line Loop(Channel_Boundary) = {Reverse_Channel_Loop[]};
    Channel = newreg; //Plane Surface(Channel) = {Channel_Boundary};
    Printf("Create Channel[%g,%g] #%g  Boundary[%g] Surface[%g]", i-1, i, num_channel -1, Channel_Boundary, Channel);
    n = 0;
    
    //
    // Even Channel Loop
  EndIf
EndFor

Printf(" ");
num_tot = 0;
num_helix = 0;
For i In{1:N_helices}
  num_helix += 6*N_turns[i];
  For j In{1:N_turns[i]-1}
     num_tot += 4;
  EndFor
EndFor

num_tot += num_helix - 4*(N_turns[i]-1);
n = 0;

num_channel = N_helices -1;
//Printf("num_channel = %g", num_channel);

For i In{1:N_helices}
  num_i = N_helices+1 -i;
  If ( num_i % 2 != 0 && num_i != 1) 
    //n = 0;
    num_helix -= 6*N_turns[num_i];
    //Printf("Odd Helix (%g) %g", num_tot, num_i);
    //Printf(" %g",-Ring_Ext_L1[num_i]);
    Channel_Loop[n] = -Ring_Ext_L1[num_i]; n += 1;
    
   
    //Printf(" %g",Helix_Upper_EndTurn_Ext[num_i]);  
    Channel_Loop[n] = Helix_Upper_EndTurn_Ext[num_i]; n += 1;

    num_ext = num_helix+6*N_turns[num_i]-2;
    For j In {1:N_turns[num_i]-1}
       //Printf(" %g",num_ext);   
       Channel_Loop[n] = num_ext; n += 1;
       num_ext -= 6;
       //Printf(" %g",-(num_tot+2));
       Channel_Loop[n] = -(num_tot+2); n += 1;
       num_tot -= 4;
    EndFor
    //Printf(" %g ", num_ext);
    Channel_Loop[n] = num_ext; n += 1;
    num_ext -= 6;
    //Printf(" %g",Helix_Lower_EndTurn_Ext[num_i]);
    Channel_Loop[n] = Helix_Lower_EndTurn_Ext[num_i]; n+=1;
    //Printf(" %g",-Ring_Int[num_i-1]);
    Channel_Loop[n] = -Ring_Int[num_i-1]; n += 1;
    num_tot -= 4*(N_turns[num_i-1]-1)-4;
    
    num_channel -= 1;
  EndIf
  
  If ( num_i % 2 == 0 ) 
    num_helix -= 6*N_turns[num_i];
    num_int = num_helix + 2;
    //Printf("Even Helix(%g) %g", num_tot, num_i);
    If ( num_i != N_helices) 
      //Printf(" %g ",Helix_Lower_EndTurn_Int[num_i]);
      Channel_Loop[n] = Helix_Lower_EndTurn_Int[num_i]; n += 1;
    EndIf   
    For j In {1:N_turns[num_i]-1}
       If ( num_i != N_helices) 
         //Printf(" %g ",num_int);
         Channel_Loop[n] = num_int; n += 1;
       EndIf  
       num_int += 6;
       If ( num_i != N_helices) 
          //Printf(" %g",num_tot+1);   
          Channel_Loop[n] = num_tot+1; n += 1;
       EndIf
       num_tot += 4;
    EndFor
    If ( num_i != N_helices) 
      //Printf(" %g ",num_int); 
      Channel_Loop[n] = num_int; n += 1;
    EndIf  
    num_int += 6;
    num_tot -= 4*N_turns[num_i];
    If ( num_i != N_helices) 
      //Printf(" %g",Helix_Upper_EndTurn_Int[num_i]);
      Channel_Loop[n] = Helix_Upper_EndTurn_Int[num_i]; n += 1;
      //Printf(" %g",-Ring_Ext_L5[num_i-1]);
      Channel_Loop[n] = -Ring_Ext_L5[num_i-1]; n += 1;
      //Printf(" %g",Channel_L1[num_channel]);
      Channel_Loop[n] = Channel_L1[num_channel];
      //Printf(" %g",Channel_L2[num_channel]);
      Channel_Loop[n] = Channel_L2[num_channel];
      // reverse Order of Channel_Loop
      For m In{0:n}
        Reverse_Channel_Loop[m] = -Channel_Loop[n-m];
	//Printf("Channel_Loop[%g] = %g ", m, Channel_Loop[m]);
      EndFor
      Channel_Boundary = newreg; Line Loop(Channel_Boundary) = {Reverse_Channel_Loop[]};
      Channel = newreg; //Plane Surface(Channel) = {Channel_Boundary};
      Printf("Create Channel[%g,%g] #%g  Boundary[%g] Surface[%g]", i-1, i, num_channel, Channel_Boundary, Channel);
      num_channel -= 1;
    EndIf
    //num_channel -= 1;
    n = 0;
  EndIf
EndFor

Printf("Create Exterior Boundary");
Call Exterior;

Printf("Create insert boundary Loop");
n = 0;

num_tot = 0;
For i In{1:N_helices}
  num_tot += 6*N_turns[i];
EndFor
//Printf("num_tot= %g",num_tot);

orient = 1;
num_helix = 0;
num_channel = 1;
n = 0;

//Printf(" %g",Helix_Bottom[1]);
External_Loop[n] = Helix_Bottom[1]; n += 1;
//Printf("Channel_L1=%g", Channel_L1[num_channel]);
External_Loop[n] = Channel_L1[num_channel]; n += 1;
//Printf("Channel_L2=%g", Channel_L2[num_channel]); 
External_Loop[n] = Channel_L2[num_channel]; n += 1;
For i In{1:N_helices}
  If ( i % 2 != 0) 
    //Printf("Odd Helix (%g)", i);
    num_int = num_helix + 2;
    For j In {1:N_turns[i]-1}
       num_int += 6;
       num_tot += 4;
    EndFor
    num_int += 6;
    num_helix += 6*N_turns[i];
  EndIf
  If ( i % 2 == 0) 
    num_ext = num_helix + 6*N_turns[i] - 2;
    num_tot += 4*(N_turns[i]-2);
    //Printf("Even Helix(%g)", i);
    For j In {1:N_turns[i]-1}
       num_ext -= 6;
       num_tot -= 4;
    EndFor
    num_ext -= 6;
    If ( i <= N_helices-1)
      //Printf(" %g",Ring_Ext_L2[i]);
      External_Loop[n] = Ring_Ext_L2[i]; n += 1;
      //Printf(" %g",Ring_Ext_L3[i]);
      External_Loop[n] = Ring_Ext_L3[i]; n += 1;
      //Printf(" %g",Ring_Ext_L4[i]);
      External_Loop[n] = Ring_Ext_L4[i]; n += 1;
    EndIf
    num_helix += 6*N_turns[i];
    num_tot += 4*N_turns[i];
    If (num_channel+2 == N_helices-1 && i != N_helices)
      //Printf("Channel_L1=%g", Channel_L2[num_channel+2]);
      External_Loop[n] = Channel_L2[num_channel+2]; n += 1;
      //Printf("Channel_L2=%g", Channel_L1[num_channel+2]); 
      External_Loop[n] = Channel_L1[num_channel+2]; n += 1;
    EndIf
    If (num_channel+2 != N_helices-1 && i != N_helices)
      //Printf("Channel_L1=%g", Channel_L1[num_channel+2]);
      External_Loop[n] = Channel_L1[num_channel+2]; n += 1;
      //Printf("Channel_L2=%g", Channel_L2[num_channel+2]); 
      External_Loop[n] = Channel_L2[num_channel+2]; n += 1;
      num_channel += 2;
    EndIf
  EndIf
EndFor
//Printf(" %g",Helix_Bottom[N_helices]);
External_Loop[n] = Helix_Bottom[N_helices];

// reverse Order of Channel_Loop
For m In{0:n}
   Reverse_External_Loop[m] = -External_Loop[n-m];
EndFor

Ext = newreg; Line Loop(Ext) = {c1, -l7 ,Reverse_External_Loop[], -l6, l8};
Ext_Mesh = newreg; //Plane Surface(Ext_Mesh) = {Ext};

Printf(" ");
n = 0;
num_tot = 0;
num_helix = 0;
For i In{1:N_helices}
  num_helix += 6*N_turns[i];
  For j In{1:N_turns[i]-1}
     num_tot += 4;
  EndFor
EndFor

num_tot += num_helix - 4*(N_turns[i]-1);

Printf("New External Contour");

num_channel = N_helices -2;

For i In{1:N_helices}
  num_i = N_helices+1 -i;
  If ( num_i % 2 != 0) 
    num_helix -= 6*N_turns[num_i];
    //Printf("Odd Helix (%g) %g", num_tot, num_i);
    If ( num_i == N_helices - 1)
      //Printf(" %g",-Ring_Ext_L5[num_i]);
      External_Loop[n] = -Ring_Ext_L5[num_i];
      
      //Printf("loop:");
      //Printf(" %g %g %g", l7 ,c3, -l10);
      // reverse Order of Channel_Loop
      For m In{0:n}
         Reverse_External_Loop[m] = -External_Loop[n-m];
	 //Printf(" %g", Reverse_External_Loop[m]);
      EndFor
      
      Ext = newreg; Line Loop(Ext) = {l7 ,c3, -l10, Reverse_External_Loop[]};
      Ext_Mesh = newreg; 
      //Plane Surface(Ext_Mesh) = {Ext};
      //Transfinite Surface{Ext_Mesh} = {Helix_EndTurn_P2[2*N_helices-1], p8, p9, Ring_P4[N_helices-1]};
      //Recombine Surface{Ext_Mesh};
      
      n = 0;
      Printf("New External Contour");
    EndIf
    If ( num_i >= 1)
      //Printf(" %g",-Ring_Ext_L4[num_i]);
      External_Loop[n] = -Ring_Ext_L4[num_i]; n += 1;
      //Printf(" %g",-Ring_Ext_L3[num_i]);
      External_Loop[n] = -Ring_Ext_L3[num_i]; n += 1;
      //Printf(" %g",-Ring_Ext_L2[num_i]);
      External_Loop[n] = -Ring_Ext_L2[num_i]; n += 1;
      If (num_i != 1) 
        //Printf(" %g",-Channel_L2[num_channel]);
        External_Loop[n] = -Channel_L2[num_channel]; n += 1;
        //Printf(" %g",-Channel_L1[num_channel]);
        External_Loop[n] = -Channel_L1[num_channel]; n += 1;
        num_channel -= 2;
      EndIf
    EndIf
   
    num_ext = num_helix+6*N_turns[num_i]-2;
    For j In {1:N_turns[num_i]-1}
       num_ext -= 6;
       num_tot -= 4;
    EndFor
    num_ext -= 6;
    If ( num_i >= 2)
      num_tot -= 4*(N_turns[num_i-1]-1)-4;
    EndIf
  EndIf
  If ( num_i % 2 == 0) 
    num_helix -= 6*N_turns[num_i];
    num_int = num_helix + 2;
    If ( num_i == N_helices) 
      //Printf(" %g ",Helix_Lower_EndTurn_Int[num_i]);
      External_Loop[n] = Helix_Lower_EndTurn_Int[num_i]; n += 1;
    EndIf
    For j In {1:N_turns[num_i]-1}
       If ( num_i == N_helices)
         //Printf(" %g ",num_int); 
         External_Loop[n] = num_int; n += 1;
       EndIf
       num_int += 6;
       If ( num_i == N_helices)
         //Printf(" %g",num_tot+1); 
         External_Loop[n] = num_tot+1; n += 1;
       EndIf
       num_tot += 4;
    EndFor
    If ( num_i == N_helices)
      //Printf(" %g ",num_int);
      External_Loop[n] = num_int; n += 1;
    EndIf
    num_int += 6;
    If ( num_i == N_helices)
      //Printf(" %g",Helix_Upper_EndTurn_Int[num_i]);
      External_Loop[n] = Helix_Upper_EndTurn_Int[num_i]; n += 1;
    EndIf
    num_tot -= 4*N_turns[num_i];
  EndIf
EndFor
//Printf(" %g", l9);
External_Loop[n] = l9;

// Reverse Loop
For m In{0:n}
   Reverse_External_Loop[m] = -External_Loop[n-m];
EndFor

Ext = newreg; Line Loop(Ext) = {l10, c5, l11, Reverse_External_Loop[]};
Ext_Mesh = newreg; //Plane Surface(Ext_Mesh) = {Ext};
Printf("New External Contour");

n = 0;
Printf(" %g", -(Helix_EndTurn_L1[1]+2));
External_Loop[n] = -(Helix_EndTurn_L1[1]+2); n += 1;
num = 4;
num_tot = 0;
For i In{1:N_helices}
  num_tot += 6*N_turns[i];
EndFor

//Printf("num_tot= %g",num_tot);
For i In{1:N_turns[1]-1}
   Printf(" %g", -num);
   External_Loop[n] = -num; n += 1;
   num += 6;
   Printf(" %g", (num_tot+2));
   External_Loop[n] = (num_tot+2); n += 1;
EndFor   
Printf(" %g", -num);
External_Loop[n] = -num; n += 1;
Printf(" %g", Helix_EndTurn_L1[2]+2);
External_Loop[n] = Helix_EndTurn_L1[2]+2; n += 1;
Printf(" %g", Ring_Ext_L1[1]);
External_Loop[n] = Ring_Ext_L1[1]; n += 1;
Printf(" %g", l9);
External_Loop[n] = l9; n += 1;
Printf(" %g", l5);
External_Loop[n] = l5; n += 1;
Printf(" %g", -l6);
External_Loop[n] = -l6;

Ext = newreg; Line Loop(Ext) = {External_Loop[]};
Ext_Mesh = newreg; Plane Surface(Ext_Mesh) = {Ext};
Transfinite Surface{Ext_Mesh} = {p7, Helix_EndTurn_P1[1], Ring_P1[1], p10};
Recombine Surface{Ext_Mesh};