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};