Hi,
I have some trouble generating the mesh of a "slice" of our magnets
especially for the region "Air' defined in HR_New_Virtual_Sphere.geo.
The geometry is defined in HR_New_GMSH.geo which includes several files.
The construction is parametrized (see data in HR_New-DAT.dat) : first we
build the magnet (see functions defined in HR_New-GEO.geo) then we
define the air around (see HR_New_Virtual_Sphere.geo which defines Air
and AirInf ) to compute the magnetic field generated by the magnet.
The trouble is that I cannot get a volumic mesh of the Air region...
I only have a mesh of the "Air" boundaries.
I must do something wrong in the definition of the Air region?
I would be grateful for any help to fix this problem
Best
C
/// 2D mesh algorithm (1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
Mesh.Algorithm=2;
// 3D mesh algorithm (1=Delaunay, 4=Frontal, 5=Frontal Delaunay, 6=Frontal Hex, 7=MMG3D)
Mesh.Algorithm3D=6;
// Remeshing algorithm (0=no split, 1=automatic, 2=automatic only with metis)
Mesh.RemeshAlgorithm=0;
// Remesh Parametrization (0=harmonic, 1=conformal, 2=rbf harmonic)
Mesh.RemeshParametrization=0;
// Mesh recombination algorithm (0=standard, 1=blossom)
Mesh.RecombinationAlgorithm=0;
// Mesh coloring (0=by element type, 1=by elementary entity, 2=by physical entity, 3=by partition)
Mesh.ColorCarousel=2;
//Mesh.NbPartitions = 1;
//Mesh.Smoothing = 5;
//Mesh.Optimize = 1;
/*-----------------------------------------------
-------------PARAMETRES UTILISATEURS-------------
------------------------------------------------*/
Include "HR_New-DAT.dat";
/*-----------------------------------------------
-------CONSTRUCTION SECTEUR AVEC KAPTON----------
------------------------------------------------*/
Include "HR_New-GEO.geo";
/*-----------------------------------------------
------------DEFINITIONS--------------------
------------------------------------------------*/
num_bord_r = 0;
Bord_RI[] = {};
Bord_RE[] = {};
num_n = 0;
Bord_Neumann[] = {};
num_h = 0;
Bord_RChannel[]={};
num_r = 0;
Bord_RChannel[]={};
num_i = 0;
Bord_Isolant[]={};
/*-----------------------------------------------
------------BOUCLE PRINCIPALE--------------------
------------------------------------------------*/
z=z0-Ep_Isolant;
t=0;
EpaisseurSpire=h[0];
Printf("-------------------------------");
Printf("spire=%g",t);Printf("h[%g]=%g (m)",t,h[t]);
t_end=Nb_spires-1;
If (ends)
Call GeomSpire;
t_end +=1;
EndIf
z=z0+EpaisseurSpire;
For t In {1:(t_end)} //Nb_spires
EpaisseurSpire=h[t];
If ( t >= Start_spire_number && t <= End_spire_number )
Call GeomSpire;
If(isolant)
Call GeomIsol;
EndIf
Coherence;
EndIf
z +=EpaisseurSpire+Ep_Isolant; // Printf("z=%g",z);
If ( ends==0 && t == End_spire_number )
If (isolant)
Call GeomIsol;
EndIf
Coherence;
EndIf
EndFor
Call Physical_Region;
/*-----------------------------------------------
------------Surrounding Air --------------------
------------------------------------------------*/
If(sphere)
LcAxis = 0.1*Unit ;
LcBox = 5*Unit ;
LcInf = 10*Unit ;
R=0.8*Height; //200*Unit;
Delta_R=2*Height; //10*Unit;
Include "HR_New_Virtual_Sphere.geo";
Call Virtual_Sphere ;
EndIf
/******************************************************************************************
PARAMETRES UTILISATEURS
*******************************************************************************************/
hexas=0; // Maillage Hexas -- to be implemented using blossom and quad. support
isolant=0; // Presence de l'isolant
sphere=1; // Precense de l'infini
ends=1; // Presence des extremites helices
/******************************************************************************************
Definition unite de longueur
*******************************************************************************************/
Unit = 1.e-3; //(conversion mm --> m)
/******************************************************************************************
Definition de la geometrie
*******************************************************************************************/
Rint = 61.2*0.5*Unit ; // Valeur du Rayon interieur (mm)
Rext = 106.4*0.5*Unit ; // Valeur du Rayon exterieur (mm)
Nb_spires = 34; //34 ; // Nombres de spires
Nb_isolants = 8 ; // Nombre d'isolants par spire
Ep_Isolant = 0.2*Unit; // Valeur de l'epaisseur de l'isolant (mm)
/* Eventually choose a specific spire from 0 to Nb_spires+1 */
Start_spire_number = 0;
End_spire_number = Nb_spires+1;
/******************************************************************************************
Definition des hauteurs des spires
*******************************************************************************************/
h[1] = 4.61*Unit;h[2] = 4.61*Unit;h[3] = 4.96*Unit;h[4] = 5.33*Unit;h[5] = 5.11*Unit;
h[6] = 4.61*Unit;h[7] = 4.61*Unit;h[8] = 4.61*Unit;h[9] = 4.61*Unit;h[10] = 4.61*Unit;
h[11] = 4.61*Unit;h[12] = 4.61*Unit;h[13] = 4.61*Unit;h[14] = 4.61*Unit;h[15] = 4.61*Unit;
h[16] = 4.61*Unit;h[17] = 4.61*Unit;h[18] = 4.61*Unit;h[19] = 4.61*Unit;h[20] = 4.61*Unit;
h[21] = 4.61*Unit;h[22] = 4.61*Unit;h[23] = 4.61*Unit;h[24] = 4.61*Unit;h[25] = 4.61*Unit;
h[26] = 4.61*Unit;h[27] = 4.61*Unit;h[28] = 4.61*Unit;h[29] = 4.61*Unit;h[30] = 4.75*Unit;
h[31] = 5.33*Unit;h[32] = 5.32*Unit;h[33] = 4.61*Unit;h[34] = 4.61*Unit;
// 1st Nth spire
Nb_spires++;
// Extremites
If (ends)
h[0]=42.54*Unit;
h[Nb_spires]=52.76*Unit;
EndIf
Spire_Milieu = (Nb_spires-1)/2+1;//18
Height = 0;
For i In {0:#h[]-1}
Height += h[i];
EndFor
Height += (#h[]-1) * Ep_Isolant;
Printf("Height=%g", Height);
Active_Height = Height;
If (!ends)
Active_Height -= h[0]+h[Nb_spires]-2*Ep_Isolant);
EndIf
Printf("Active_Height = %g m", Active_Height);
Z_Spire_Milieu = Active_Height /2.;
Printf("Z_Spire_Milieu = %g m", Z_Spire_Milieu);
Printf("%g %g %g", End_spire_number, Start_spire_number, Nb_spires);
If ( End_spire_number-Start_spire_number != Nb_spires )
ends = 0;
sphere = 0;
EndIf
/*******************************************************************************************
Definition de la "TUILE"
Forme generique :
AngleIsolant_Int
<---->
|____|..................
___/ \___...............|H_TUILE
| | | |
| <------> |
| AngleIsolant_Mil
| |
<------------->
AngleIsolant_Ext=Longueur angulaire de l'isolant en degre
*******************************************************************************************/
AngleIsolant_Ext= 10; // Angle Exterieur de l'isolant (deg�)
AngleIsolant_Mil= 7.5; // Angle Mileu de l'isolant (deg�)
AngleIsolant_Int= 4; // Angle Interieur de l'isolant (deg�)
TUILE=1;
H_TUILE=0.9*Unit; //Valeur hauteur tuile (mm)
/*******************************************************************************************
Definition du canal
******************************************************************************************/
Canal = 1;
H_canal= 0.2*Unit; // Hauteur du canal de refroidissement
lg_canal= 1.5*Unit; // Valeur Largeur du canal (mm)
/* Position en z des spires*/
z0 =0;
/*****************************************************************************************
Caracteristiques du Maillage
*******************************************************************************************/
LcSpire = 1*Unit ;
LcCanal = 0.5*Unit;
/*****************************************************************************************
Affichage caracteristiques Helices
*******************************************************************************************/
Printf("Rayon Interieur (m) = %g",Rint);
Printf("Rayon Exterieur(m) = %g",Rext);
Printf("Nombre de spires = %g",Nb_spires);
Printf("Nombre d'isolants par spire = %g",Nb_isolants);
Printf("Epaisseur isolant (m) = %g",Ep_Isolant);
/***
For t In {0:#h[]-1}
Printf("Spire[%g] h = %g m",t, h[t]);
EndFor
****/
Printf("Angle Exterieur de l'isolant (�) = %g",AngleIsolant_Ext);
Printf("Angle Mileu de l'isolant (�) = %g",AngleIsolant_Mil);
Printf("Angle Interieur de l'isolant (�) = %g",AngleIsolant_Int);
/*-------------------------------------------------
-----------------FONCTIONS------------------------
------------------------------------------------*/
in_isolant = 0;
/* Position en z des spires*/
z0 = -Z_Spire_Milieu;
TetaKapton_ext = (Pi*AngleIsolant_Ext)/180. ; // Angle Exterieur de l'isolant (rad)
TetaKapton_mid = (Pi*AngleIsolant_Mil)/180. ; // Angle Mileu de l'isolant (rad)
TetaKapton_int = (Pi*AngleIsolant_Int)/180. ; // Angle Mileu de l'isolant (rad)
lg_canal2 = lg_canal+(0.1*Unit); // largeur du canal2----> points necessaires pour un maillage Transfinite
/* Store angles into a table */
alpha[] = {};
is_tuile[] = {};
is_canal[] = {};
is_line[] = {};
alpha[0] = -(2*Pi)/Nb_isolants /2.; is_tuile[0] = 0; is_canal[0] = 0; is_line[0] = 0; //P0
alpha[1] = -(TetaKapton_ext/2); is_tuile[1] = 0; is_canal[1] = 0; is_line[1] = 0; //P1
i=1;
If ( TUILE )
i++; alpha[i] = -(TetaKapton_mid/2); is_tuile[i] = 0; is_canal[i] = 0; is_line[i] = 0; //P2
i++; alpha[i] = -(TetaKapton_int/2); is_tuile[i] = 1; is_canal[i] = 0; is_line[i] = 1; //P3
If ( Canal )
i++; alpha[i] = NaN; is_tuile[i] = 1; is_canal[i] = 1; is_line[i] = 0; //P4
i++; alpha[i] = NaN; is_tuile[i] = 0; is_canal[i] = 1; is_line[i] = 1; //P5
EndIf
EndIf
i++; alpha[i] = 0; mid_point = i; // P6
If ( TUILE )
is_tuile[i] = 1; is_canal[i] = 0; is_line[i] = 0;
If ( Canal )
is_canal[i] = 1;
i++; alpha[i] = NaN; is_tuile[i] = 0; is_canal[i] = 1; is_line[i] = 0; // P7
i++; alpha[i] = NaN; is_tuile[i] = 1; is_canal[i] = 1; is_line[i] = 1; // P8
EndIf
i++; alpha[i] = (TetaKapton_int/2); is_tuile[i] = 1; is_canal[i] = 0; is_line[i] = 0; // P9
i++; alpha[i] = (TetaKapton_mid/2); is_tuile[i] = 0; is_canal[i] = 0; is_line[i] = 1; // P10
EndIf
i++; alpha[i] = (TetaKapton_ext/2); is_tuile[i] = 0; is_canal[i] = 0; is_line[i] = 0; //P11
i++; alpha[i] = -alpha[0]; is_tuile[i] = 0; is_canal[i] = 0; is_line[i] = 0; //P12
/* Store dz for corresponding nodes */
dz0[] = {};
For j In {0 : #alpha[]-1}
dz0[j] = 0;
If ( is_tuile[j] || is_canal[j] )
dz0[j] = H_TUILE;
EndIf
EndFor
/* define Characteristic Length on Nodes */
Lc[] = {};
For j In {0 : #alpha[]-1}
Lc[j] = LcSpire;
If ( is_canal[j] )
Lc[j] = LcCanal ;
EndIf
EndFor
/* Create Quadrangle mesh*/
If (hexas)
Mesh.RecombinationAlgorithm=1; // Activate Blossom support
Mesh.RecombineAll=1;//recombine all defined surfaces
Mesh.Algorithm=8;//delquad mesher
EndIf
If (!hexas)
Mesh.RecombinationAlgorithm=0; // Deactivate Blossom support
Mesh.RemeshAlgorithm=1; //(0=no split, 1=automatic, 2=automatic only with metis)
EndIf
/*----------------------------------------------
Update Dz depending on Spire postion
----------------------------------------------*/
Function UpdateDz
If(t < Spire_Milieu)
For j In {0 : #dz0[]-1}
dz[j] = -dz0[j];
EndFor
EndIf
If (t==Spire_Milieu)
For j In {0 : #dz0[]-1}
dz[j] = bas ? -dz0[j] : dz0[j];
EndFor
EndIf
If(t > Spire_Milieu)
For j In {0 : #dz0[]-1}
dz[j] = dz0[j];
EndFor
EndIf
If (t==0)
For j In {0 : #dz0[]-1}
dz[j] = bas ? 0: -dz0[j];
EndFor
EndIf
If (t==Nb_spires)
For j In {0 : #dz0[]-1}
dz[j] = bas ? dz0[j] : 0;
EndFor
EndIf
For j In {0 : #dz0[]-1}
If ( !in_isolant && (is_canal[j] && !is_tuile[j]) )
dz[j] -= bas ? -H_canal : H_canal;
EndIf
If ( !in_isolant && (Canal && j==mid_point) )
dz[j] -= bas ? -H_canal : H_canal;
EndIf
EndFor
Return
/*----------------------------------------------
Create Points and Lines on Secteur Edge R,Z
----------------------------------------------*/
Function Lines
step = 1;
j0 = 0;
jend = #alpha[]-1;
If (in_isolant )
j0 = 1;
jend = #alpha[]-2;
EndIf
For j In {j0 : jend: step}
If ( !is_canal[j] || j==mid_point )
P[j] = newp; Point(P[j]) = {R*Cos[alpha[j]], R*Sin[alpha[j]], Z + dz[j], Lc[j]};
//Printf("P[%g]=[%g]=[%g,%g,%g]", j, P[j], R*Cos[alpha[j]], R*Sin[alpha[j]], Z + dz[j]);
EndIf
If ( is_canal[j] && j!=mid_point)
If ( j < mid_point )
P[j] = newp; Point(P[j]) = {Sqrt[R*R-lg_canal/2.*lg_canal/2.], -lg_canal/2., Z + dz[j], Lc[j]};
//Printf("P[%g]=[%g]=[%g,%g,%g] *", j, P[j], Sqrt[R*R-lg_canal/2.*lg_canal/2.], -lg_canal/2., Z + dz[j]);
EndIf
If ( j >= mid_point )
P[j] = newp; Point(P[j]) = {Sqrt[R*R-lg_canal/2.*lg_canal/2.], lg_canal/2., Z + dz[j], Lc[j]};
//Printf("P[%g]=[%g]=[%g,%g,%g] x", j, P[j], Sqrt[R*R-lg_canal/2.*lg_canal/2.], lg_canal/2., Z + dz[j]);
EndIf
EndIf
EndFor
/* Point on Axis in plane z=Z */
O = newp; Point(O) = {0, 0, Z, 1*Unit};
Nb_Line = #P[]-1;
For j In {j0 : Nb_Line-1: step}
L[j] = 0;
EndFor
i = -1;
For j In {j0 : #P[]-2}
If ( !is_line[j+1] )
//Printf("P[%g]P[%g]=[%g, O, %g]", j, j+1, P[j], P[j+1]);
i++; L[i] = newl; Circle(L[i]) = {P[j], O, P[j+1]};
EndIf
If ( is_line[j+1] )
//Printf("P[%g]P[%g]=[%g, %g]", j, j+1, P[j], P[j+1]);
i++; L[i] = newl; Line(L[i]) = {P[j], P[j+1]};
EndIf
EndFor
//Printf("N_P=%g", #P[]);
//Printf("N_L=%g", #L[]);
/***********
If ( in_isolant)
Printf("in isolant : Reset Points %g, %g, ", P[mid_point-1], P[mid_point+1]);
Printf("in isolant : Reset Lines %g %g", L[mid_point-2], L[mid_point+1]);
P[mid_point-1] = 0;
P[mid_point+1] = 0;
L[mid_point-2] = 0;
L[mid_point+1] = 0;
EndIf
*****/
Return
/*------------------------------------------------
-------------------GEOMETRIE CUIVRE---------------
------------------------------------------------*/
Function GeomSpire //FONCTION DE DUPLICATION
/*----------------------------------------------
- SCHEMA "HR_NEW" -
- -
- Identifications des elements: -
- -
- Points: pXX -
- -
- Lignes: -
- - suivant x : r_pXX_pXXX -
- - suivant y : c_pXX_pXXX -
- - suivant z : h_pXX_pXXX -
exemple: ligne entres points p20 & p30: c_p20_p30
- -
- Faces: -
- - Normale x: RXX -
- - Normale y: LXX -
- - Normale z: HXX ou BXX -
- -
- Volumes: VXX -
------------------------------------------------
IH7_______________IH6 IH4_____________IH3
/ | | \
/ | | \
/ IH12_____IH5_____IH11 \
/ \
/ \
IH10________IH9______________IH8 \IH2_________IH1__________IH0
| |
| |
| |
| |
INLET | | OUTLET
| |
| |
| |
| |
| |
|___________ _____________ ______________ ___________|
IB10 IB9 \IB8 /IB2 IB1 IB0
\ /
\ IB12____IB5_____ IB11 /
z| \ | | /
| \ | | /
|____theta \ ________________| |________________ /
IB7 IB6 IB4 IB3
EH7_______________EH6 EH4_____________EH3
/ | | \
/ | | \
/ EH12_____EH5_____EH11 \
/ \
/ \
EH10________EH9______________EH8 \EH2_________EH1__________EH0
| |
| |
| |
| |
INLET | | OUTLET
| |
| |
| |
| |
| |
|___________ _____________ ______________ ___________|
EB10 EB9 \EB8 /EB2 EB1 EB0
\ /
\ EB12____EB5_____ EB11 /
z| \ | | /
| \ | | /
|____theta \ ________________| |________________ /
EB7 EB6 EB4 EB3
*/
/*---------GEOMETRIE SPIRE--------------------*/
// Face du haut
Axe_Z = z+EpaisseurSpire+Ep_Isolant;
Printf("Axe_Z (m) = %g",Axe_Z);
bas = 0;
dz[] = {};
Call UpdateDz;
sens = 1;
R = Rint;
Z = Axe_Z;
IH[] = {};
LIH[] = {};
P[] = {};
L[] = {};
Call Lines;
For j In {0 : #P[]-1}
IH[j] = P[j];
EndFor
For j In {0 : #L[]-1}
LIH[j] = L[j];
EndFor
N_IH=#IH[]-1; Printf("N_IH=%g", N_IH);
If ( Canal )
N_IH -=2;
EndIf
sens = -1;
R = Rext;
Z = Axe_Z;
EH[] = {};
LEH[] = {};
P[] = {};
L[] = {};
Call Lines;
For j In {0 : #P[]-1}
EH[j] = P[j];
EndFor
For j In {0 : #L[]-1}
LEH[j] = L[j];
EndFor
N_EH=#EH[]-1;
If ( Canal )
N_EH -=2;
EndIf
LIEH[] = {};
For j In {0 : #IH[]-1}
LIEH[j] = newl; Line(LIEH[j]) = {IH[j], EH[j]};
EndFor
// Face du bas
Axe_Z = z+Ep_Isolant;
Printf("Axe_Z (m) = %g",Axe_Z);
bas = 1;
dz[] = {};
Call UpdateDz;
sens = -1;
R = Rint;
Z = Axe_Z;
IB[] = {};
LIB[] = {};
P[] = {};
L[] = {};
Call Lines;
For j In {0 : #P[]-1}
IB[j] = P[j];
EndFor
For j In {0 : #L[]-1}
LIB[j] = L[j];
EndFor
N_IB=#IB[]-1;
If ( Canal )
N_IB -=2;
EndIf
sens = 1;
R = Rext;
Z = Axe_Z;
EB[] = {};
LEB[] = {};
P[] = {};
L[] = {};
Call Lines;
For j In {0 : #P[]-1}
EB[j] = P[j];
EndFor
For j In {0 : #L[]-1}
LEB[j] = L[j];
EndFor
N_EB=#EB[]-1;
If ( Canal )
N_EB -=2;
EndIf
LIEB[] = {};
For j In {0 : #IB[]-1}
LIEB[j] = newl; Line(LIEB[j]) = {IB[j], EB[j]};
EndFor
LIHIB[] = {};
LEHEB[] = {};
For j In {0 : #IB[]-1}
LIHIB[j] = newl; Line(LIHIB[j]) = {IH[j], IB[j]};
LEHEB[j] = newl; Line(LEHEB[j]) = {EH[j], EB[j]};
EndFor
/* Surfaces */
N_Surf = #IH[] -1;
N_R_Surf = N_Surf;
If ( Canal )
N_R_Surf -=2;
N_R_Surf += 3 * 2;
EndIf
Printf("N_Surf = %g", N_Surf);
Printf("N_R_Surf = %g", N_R_Surf);
/*
For j In {0 : #LIH[]-1}
Printf("LIH[%g]=%g", j, LIH[j]);
EndFor
For j In {0 : #LIB[]-1}
Printf("LIB[%g]=%g", j, LIB[j]);
EndFor
For j In {0 : #LIH[]-1}
Printf("LIH[%g]=%g", j, LIH[j]);
Printf("LEH[%g]=%g", j, LEH[j]);
EndFor
*/
SH[] = {};
For j In {0 : #IH[]-2}
//Printf("Line Loop = %g, %g, %g, %g", LIEH[j], LEH[j], -LIEH[j+1], -LIH[j]);
IEH[j] = newl; Line Loop(IEH[j]) = {LIEH[j], LEH[j], -LIEH[j+1], -LIH[j]};
SH[j]=newreg; Plane Surface(SH[j]) = {IEH[j]};
EndFor
//H = newreg; Compound Surface(H) = {SH[]};
SB[] = {};
For j In {0 : #IB[]-2}
//Printf("Line Loop = %g, %g, %g, %g", LIEB[j], LEB[j], -LIEB[j+1], -LIB[j]);
IEB[j] = newl; Line Loop(IEB[j]) = {LIEB[j], LEB[j], -LIEB[j+1], -LIB[j]};
SB[j]=newreg; Plane Surface(SB[j]) = {IEB[j]};
EndFor
//B = newreg; Compound Surface(B) = {SB[]};
SRI[] = {}; // Create Coumpound Surface
i = 0;
For j In {0 : #IB[]-2}
//Printf("Line Loop = %g, %g, %g, %g", LIH[j], LIHIB[j+1], -LIB[j], -LIHIB[j]);
dSRI[j] = newl; Line Loop(dSRI[j]) = {LIH[j], LIHIB[j+1], -LIB[j], -LIHIB[j]};
If ( !Canal )
SRI[i]=newreg; Ruled Surface(SRI[i]) = {dSRI[j]}; i++;
EndIf
If ( Canal )
/* Remove surface with 0 area */
If ( j != 7 && j != 4)
If ( j != 3 && j != 8)
SRI[i]=newreg; Ruled Surface(SRI[i]) = {dSRI[j]}; i++;
EndIf
/* Split the surface into 3 parts */
If ( j == 3 )
L1 = newl; Line(L1) = {IB[3], IB[5]}; //Printf("L1=%g", L1);
L2 = newl; Line(L2) = {IH[3], IH[5]}; //Printf("L2=%g", L2);
dSRI_reg = newl; Line Loop(dSRI_reg) = {L2, LIHIB[5], -L1, -LIHIB[3]};
//Printf("Line Loop = %g, %g, %g, %g", L2, LIHIB[5], -L1, -LIHIB[3]);
SRI[i]=newreg; Ruled Surface(SRI[i]) = {dSRI_reg}; i++;
dSRI_reg = newl; Line Loop(dSRI_reg) = {L2, -LIH[4], -LIH[3]};
//Printf("Line Loop = %g, %g, %g", L2, -LIH[4], -LIH[3]);
SRI[i]=newreg; Ruled Surface(SRI[i]) = {dSRI_reg}; i++;
dSRI_reg = newl; Line Loop(dSRI_reg) = {L1, -LIB[4], -LIB[3]};
//Printf("Line Loop = %g, %g, %g", L1, -LIB[4], -LIB[3]);
SRI[i]=newreg; Ruled Surface(SRI[i]) = {dSRI_reg}; i++;
EndIf
If ( j == 8 )
L1 = newl; Line(L1) = {IB[9], IB[7]}; //Printf("L1=%g", L1);
L2 = newl; Line(L2) = {IH[9], IH[7]}; //Printf("L2=%g", L2);
dSRI_reg = newl; Line Loop(dSRI_reg) = {-L2, LIHIB[9], L1, -LIHIB[7]};
//Printf("Line Loop = %g, %g, %g, %g", -L2, LIHIB[9], L1, -LIHIB[7]);
SRI[i]=newreg; Ruled Surface(SRI[i]) = {dSRI_reg}; i++;
dSRI_reg = newl; Line Loop(dSRI_reg) = {L2, LIH[7], LIH[8]};
//Printf("Line Loop = %g, %g, %g", L2, LIH[7], LIH[8]);
SRI[i]=newreg; Ruled Surface(SRI[i]) = {dSRI_reg}; i++;
dSRI_reg = newl; Line Loop(dSRI_reg) = {L1, LIB[7], LIB[8]};
//Printf("Line Loop = %g, %g, %g", L1, LIB[7], LIB[8]);
SRI[i]=newreg; Ruled Surface(SRI[i]) = {dSRI_reg}; i++;
EndIf
EndIf
EndIf
EndFor
RI[] = {};
If ( !Canal )
RI[0] = newreg; Compound Surface(RI[0]) = {SRI[]};
EndIf
If ( Canal )
RI[0] = newreg; Compound Surface(RI[0]) = {SRI[0], SRI[1], SRI[2], SRI[3], SRI[4], SRI[5]};
RI[1] = newreg; Compound Surface(RI[1]) = {SRI[6]};
RI[2] = newreg; Compound Surface(RI[2]) = {SRI[7]};
RI[3] = newreg; Compound Surface(RI[3]) = {SRI[8], SRI[9], SRI[10], SRI[11], SRI[12], SRI[13]};
EndIf
SRE[] = {}; // Create Coumpound Surface
i = 0;
For j In {0 : #EB[]-2}
//Printf("Line Loop = %g, %g, %g, %g", LEH[j], LEHEB[j+1], -LEB[j], -LEHEB[j]);
dSRE[j] = newl; Line Loop(dSRE[j]) = {LEH[j], LEHEB[j+1], -LEB[j], -LEHEB[j]};
If ( !Canal )
SRE[i]=newreg; Ruled Surface(SRE[i]) = {dSRE[j]}; i++;
EndIf
If ( Canal )
/* Remove surface with 0 area */
If ( j != 7 && j != 4 )
If ( j != 3 && j != 8)
SRE[i]=newreg; Ruled Surface(SRE[i]) = {dSRE[j]}; i++;
EndIf
/* Split the surface into 3 parts */
If ( j == 3 )
L1 = newl; Line(L1) = {EB[3], EB[5]}; //Printf("L1=%g", L1);
L2 = newl; Line(L2) = {EH[3], EH[5]}; //Printf("L2=%g", L2);
dSRE_reg = newl; Line Loop(dSRE_reg) = {L2, LEHEB[5], -L1, -LEHEB[3]};
//Printf("Line Loop = %g, %g, %g, %g", L2, LEHEB[5], -L1, -LEHEB[3]);
SRE[i]=newreg; Ruled Surface(SRE[i]) = {dSRE_reg}; i++;
dSRE_reg = newl; Line Loop(dSRE_reg) = {L2, -LEH[4], -LEH[3]};
//Printf("Line Loop = %g, %g, %g", L2, -LEH[4], -LEH[3]);
SRE[i]=newreg; Ruled Surface(SRE[i]) = {dSRE_reg}; i++;
dSRE_reg = newl; Line Loop(dSRE_reg) = {L1, -LEB[4], -LEB[3]};
//Printf("Line Loop = %g, %g, %g", L1, -LEB[4], -LEB[3]);
SRE[i]=newreg; Ruled Surface(SRE[i]) = {dSRE_reg}; i++;
EndIf
If ( j == 8 )
L1 = newl; Line(L1) = {EB[9], EB[7]}; //Printf("L1=%g", L1);
L2 = newl; Line(L2) = {EH[9], EH[7]}; //Printf("L2=%g", L2);
dSRE_reg = newl; Line Loop(dSRE_reg) = {-L2, LEHEB[9], L1, -LEHEB[7]};
//Printf("Line Loop = %g, %g, %g, %g", -L2, LEHEB[9], L1, -LEHEB[7]);
SRE[i]=newreg; Ruled Surface(SRE[i]) = {dSRE_reg}; i++;
dSRE_reg = newl; Line Loop(dSRE_reg) = {L2, LEH[7], LEH[8]};
//Printf("Line Loop = %g, %g, %g", L2, LEH[7], LEH[8]);
SRE[i]=newreg; Ruled Surface(SRE[i]) = {dSRE_reg}; i++;
dSRE_reg = newl; Line Loop(dSRE_reg) = {L1, LEB[7], LEB[8]};
//Printf("Line Loop = %g, %g, %g", L1, LEB[7], LEB[8]);
SRE[i]=newreg; Ruled Surface(SRE[i]) = {dSRE_reg}; i++;
EndIf
EndIf
EndIf
EndFor
Printf("SRI[] : %g", #SRI[]);
Printf("SRE[] : %g", #SRE[]);
RE[] = {};
If ( !Canal )
RE[0] = newreg; Compound Surface(RE[0]) = {SRE[]};
EndIf
If ( Canal )
RE[0] = newreg; Compound Surface(RE[0]) = {SRE[0], SRE[1], SRE[2], SRE[3], SRE[4], SRE[5]};
RE[1] = newreg; Compound Surface(RE[1]) = {SRE[6]};
RE[2] = newreg; Compound Surface(RE[2]) = {SRE[7]};
RE[3] = newreg; Compound Surface(RE[3]) = {SRE[8], SRE[9], SRE[10], SRE[11], SRE[12], SRE[13]};
EndIf
//Printf("Line Loop = %g, %g, %g, %g", LIEB[0], -LEHEB[0], -LIEH[0], LIHIB[0]);
V0[t] = newl; Line Loop(V0[t]) = {LIEB[0], -LEHEB[0], -LIEH[0], LIHIB[0]};
TOP[t] = newreg; Plane Surface(TOP[t]) = {V0[t]};
//Printf("Line Loop = %g, %g, %g, %g", LIEB[#LIEB[]-1], -LEHEB[#LEHEB[]-1], -LIEH[#LIEH[]-1], LIHIB[#LIHIB[]-1]);
V1[t] = newl; Line Loop(V1[t]) = {LIEB[#LIEB[]-1], -LEHEB[#LEHEB[]-1], -LIEH[#LIEH[]-1], LIHIB[#LIHIB[]-1]};
BOTTOM[t] = newreg; Plane Surface(BOTTOM[t]) = {V1[t]};
/* Create Volumes */
Bord = newreg; Surface Loop(Bord) = {RI[], RE[], SB[], SH[], TOP[t], BOTTOM[t]};
//Bord = newreg; Surface Loop(Bord) = {SRI[], SRE[], SB[], SH[], TOP[t], BOTTOM[t]};
Copper[t] = newreg; Volume(Copper[t]) = {Bord};
Printf("Defining volume");
/* Update RI and RE Boundaries */
//For j In {0:#SRI[]-1}
// Bord_RI[num_bord_r] = SRI[j];
// Bord_RE[num_bord_r] = SRE[j];
// num_bord_r++;
//EndFor
For j In {0:#RI[]-1}
Bord_RI[num_bord_r] = RI[j];
Bord_RE[num_bord_r] = RE[j];
num_bord_r++;
EndFor
/* Update HChannels */
If ( t!=0 && t!=Nb_spires)
Bord_HChannel[num_h] = SB[0]; num_h++;
Bord_HChannel[num_h] = SB[#SB[]-1]; num_h++;
Bord_HChannel[num_h] = SH[0]; num_h++;
Bord_HChannel[num_h] = SH[#SH[]-1]; num_h++;
EndIf
If ( t==0)
For j In {0:#SB[]-1}
Bord_Neumann[num_n] = SB[j]; num_n++;
EndFor
Bord_HChannel[num_h] = SH[0]; num_h++;
Bord_HChannel[num_h] = SH[#SH[]-1]; num_h++;
EndIf
If ( t==Nb_spires)
For j In {0:#SB[]-1}
Bord_Neumann[num_n] = SH[j]; num_n++;
EndFor
Bord_HChannel[num_h] = SB[0]; num_h++;
Bord_HChannel[num_h] = SB[#SB[]-1]; num_h++;
EndIf
If ( !Canal )
For j In {1 : #SB[]-2}
Bord_Isolant[num_i] = SB[j]; num_i++;
Bord_Isolant[num_i] = SH[j]; num_i++;
EndFor
EndIf
If ( Canal )
For j In {1 : 3}
Bord_Isolant[num_i] = SB[j]; num_i++;
Bord_Isolant[num_i] = SH[j]; num_i++;
EndFor
For j In {4 : 7}
Bord_RChannel[num_r] = SB[j]; num_r++;
Bord_RChannel[num_r] = SH[j]; num_r++;
EndFor
For j In {8 : #SB[]-2}
Bord_Isolant[num_i] = SB[j]; num_i++;
Bord_Isolant[num_i] = SH[j]; num_i++;
EndFor
EndIf
Printf("GeomSpire done");
Return
/*------------------------------------------------
-------------------GEOMETRIE ISOLANT--------------
------------------------------------------------*/
/*----------------------------------------------
- SCHEMA "HR_NEW" -
- -
- Identifications des elements: -
- -
- Points: piXX -
- -
- Lignes: -
- - suivant x : ri_piXX_piXXX -
- - suivant y : ci_piXX_piXXX -
- - suivant z : hi_piXX_piXXX -
exemple: ligne entres points pi20 & pi30: ci_pi20_pi30
- -
- Faces: -
- - Normale x: RiXX -
- - Normale y: LiXX -
- - Normale z: HiXX ou BiXX -
- -
- Volumes: ViXX -
------------------------------------------------
IH7____________IH6___________________IH4_____________IH3
/| | | |\
/ | | | | \
/ | | | | \
IH9___________IH8 | | | | IH2_____________IH1
| | | | | | | |
| | | | | | | |
| | | | | | | |
| | | | | | | |
INLET | | IB7____________IB6___________________IB4_____________IB3 | | OUTLET
| | / \ | |
| | / \ | |
| | / \ | |
IB9___________IB8 IB2_____________IB1
|z
|
|_____theta
EH7____________EH6___________________EH4_____________EH3
/| | | |\
/ | | | | \
/ | | | | \
EH9___________IH8 | | | | EH2_____________EH1
| | | | | | | |
| | | | | | | |
| | | | | | | |
| | | | | | | |
INLET | | EB7____________EB6___________________EB4_____________EB3 | | OUTLET
| | / \ | |
| | / \ | |
| | / \ | |
EB9___________EB8 EB2_____________EB1
|z
|
|_____theta
*/
Function GeomIsol //FONCTION DE DUPLICATION
in_isolant = 1;
//face du haut
Axe_Z = z;
dz[] = {};
Call UpdateDz;
sens = 1;
R = Rint;
Z = Axe_Z;
IH[] = {};
LIH[] = {};
P[] = {};
L[] = {};
Call Lines;
//Printf("N_P=%g (Isolant)", #P[]-1);
//Printf("N_L=%g (Isolant)", #L[]-1);
i = 0;
For j In {0 : #P[]-1}
If ( P[j] != 0 )
IH[i] = P[j]; i++;
EndIf
EndFor
i = 0;
For j In {0 : #L[]-1}
If ( L[j] != 0 )
LIH[i] = L[j]; i++;
EndIf
EndFor
N_IH=#IH[]-1;
//Printf("N_IH=%g (Isolant)", #IH[]);
//Printf("N_LIH=%g (Isolant)", #LIH[]);
sens = -1;
R = Rext;
EH[] = {};
LEH[] = {};
P[] = {};
L[] = {};
Call Lines;
i = 0;
For j In {0 : #P[]-1}
If ( P[j] != 0 )
EH[i] = P[j]; i++;
EndIf
EndFor
i = 0;
For j In {0 : #L[]-1}
If ( L[j] != 0 )
LEH[i] = L[j]; i++;
EndIf
EndFor
N_EH=#EH[]-1;
//Printf("N_IE=%g (Isolant)", #EH[]);
//Printf("N_LEH=%g (Isolant)", #LEH[]);
LIEH[] = {};
For j In {0 : #IH[]-1}
LIEH[j] = newl; Line(LIEH[j]) = {IH[j], EH[j]};
EndFor
//Printf("N_LIEH=%g (Isolant)", #LIEH[]);
// Face du bas
Axe_Z = z+Ep_Isolant;
bas = 1;
dz[] = {};
Call UpdateDz;
sens = -1;
R = Rint;
Z = Axe_Z;
IB[] = {};
LIB[] = {};
P[] = {};
L[] = {};
Call Lines;
i = 0;
For j In {0 : #P[]-1}
If ( P[j] != 0 )
IB[i] = P[j]; i++;
EndIf
EndFor
i = 0;
For j In {0 : #L[]-1}
If ( L[j] != 0 )
LIB[i] = L[j]; i++;
EndIf
EndFor
N_IB=#IB[]-1;
//Printf("N_IB=%g (Isolant)", #IB[]);
//Printf("N_LIB=%g (Isolant)", #LIB[]);
sens = 1;
R = Rext;
Z = Axe_Z;
EB[] = {};
LEB[] = {};
P[] = {};
L[] = {};
Call Lines;
i = 0;
For j In {0 : #P[]-1}
If ( P[j] != 0 )
EB[i] = P[j]; i++;
EndIf
EndFor
i = 0;
For j In {0 : #L[]-1}
If ( L[j] != 0 )
LEB[i] = L[j]; //Printf("%g LEB[%g]=%g", j, i, LEB[i]);
/*******************
If ( Canal )
If ( j == 3 || j == 6)
LEB[i] = 0;
EndIf
EndIf
********************/
i++;
EndIf
EndFor
N_EB=#EB[]-1;
//Printf("N_EB=%g (Isolant)", #EB[]);
//Printf("N_LEB=%g (Isolant)", #LEB[]);
LIEB[] = {};
For j In {0 : #IB[]-1}
LIEB[j] = newl; Line(LIEB[j]) = {IB[j], EB[j]};
EndFor
//Printf("N_LIEB=%g (Isolant)", #LIEB[]);
LIHIB[] = {};
LEHEB[] = {};
For j In {0 : #IB[]-1}
LIHIB[j] = newl; Line(LIHIB[j]) = {IH[j], IB[j]};
LEHEB[j] = newl; Line(LEHEB[j]) = {EH[j], EB[j]};
EndFor
/* Surfaces */
iSH[] = {};
i = 0;
For j In {0 : #IH[]-2}
If ( !Canal )
//Printf("SH[%g]=%g, %g, %g, %g", j, LIEH[j], LEH[j], -LIEH[j+1], -LIH[j]);
IEH[i] = newl; Line Loop(IEH[i]) = {LIEH[j], LEH[j], -LIEH[j+1], -LIH[j]};
iSH[i]=newreg; Plane Surface(iSH[i]) = {IEH[i]}; i++;
EndIf
If ( Canal )
If ( j != 3 && j != 6)
//Printf("SH[%g]=%g, %g, %g, %g", j, LIEH[j], LEH[j], -LIEH[j+1], -LIH[j]);
IEH[i] = newl; Line Loop(IEH[i]) = {LIEH[j], LEH[j], -LIEH[j+1], -LIH[j]};
iSH[i]=newreg; Plane Surface(iSH[i]) = {IEH[i]}; i++;
EndIf
EndIf
EndFor
//Printf("SH Isolant=%g ", #iSH[]);
iSB[] = {};
i = 0;
For j In {0 : #IB[]-2}
If ( !Canal )
//Printf("SB[%g]=%g, %g, %g, %g", j, LIEB[j], LEB[j], -LIEB[j+1], -LIB[j]);
IEB[i] = newl; Line Loop(IEB[i]) = {LIEB[j], LEB[j], -LIEB[j+1], -LIB[j]};
iSB[i]=newreg; Plane Surface(iSB[i]) = {IEB[i]}; i++;
EndIf
If ( Canal )
If ( j != 3 && j != 6)
//Printf("SB[%g]=%g, %g, %g, %g", j, LIEB[j], LEB[j], -LIEB[j+1], -LIB[j]);
IEB[i] = newl; Line Loop(IEB[i]) = {LIEB[j], LEB[j], -LIEB[j+1], -LIB[j]};
iSB[i]=newreg; Plane Surface(iSB[i]) = {IEB[i]}; i++;
EndIf
EndIf
EndFor
//Printf("SB Isolant=%g ", #iSB[]);
//Printf("IB Isolant=%g ", #IB[]);
SRI[] = {}; // Create Coumpound Surface
i = 0;
For j In {0 : #IB[]-2}
dSRI[j] = newl; Line Loop(dSRI[j]) = {LIH[j], LIHIB[j+1], -LIB[j], -LIHIB[j]};
If ( Canal )
If ( j!= 3 && j!=6 )
SRI[i]=newreg; Ruled Surface(SRI[i]) = {dSRI[j]}; //Printf("%g SRI[%g]=%g", j, i, SRI[i]);
i++;
EndIf
EndIf
If (!Canal)
SRI[i]=newreg; Ruled Surface(SRI[i]) = {dSRI[j]}; //Printf("%g SRI[%g]=%g", j, i, SRI[i]);
i++;
EndIf
EndFor
Printf("SRI Isolant=%g ", #SRI[]);
SRE[] = {}; // Create Coumpound Surface
i = 0;
For j In {0 : #IB[]-2}
dSRE[j] = newl; Line Loop(dSRE[j]) = {LEH[j], LEHEB[j+1], -LEB[j], -LEHEB[j]};
If ( Canal )
If ( j!= 3 && j!=6 )
SRE[i]=newreg; Ruled Surface(SRE[i]) = {dSRE[j]}; i++;
EndIf
EndIf
If (!Canal)
SRE[i]=newreg; Ruled Surface(SRE[i]) = {dSRE[j]}; i++;
EndIf
EndFor
Printf("SRE Isolant=%g ", #SRE[]);
//RI = newreg; Compound Surface(RI) = {SRI[]};
//Printf("Line Loop = %g, %g, %g, %g", LIEB[0], -LEHEB[0], -LIEH[0], LIHIB[0]);
iV0[t] = newl; Line Loop(iV0[t]) = {LIEB[0], -LEHEB[0], -LIEH[0], LIHIB[0]};
ITOP[t] = newreg; Plane Surface(ITOP[t]) = {iV0[t]};
Printf("ITOP Isolant=%g ", ITOP[t]);
//Printf("Line Loop = %g, %g, %g, %g", LIEB[#LIEB[]-1], -LEHEB[#LEHEB[]-1], -LIEH[#LIEH[]-1], LIHIB[#LIHIB[]-1]);
iV1[t] = newl; Line Loop(iV1[t]) = {LIEB[#LIEB[]-1], -LEHEB[#LEHEB[]-1], -LIEH[#LIEH[]-1], LIHIB[#LIHIB[]-1]};
IBOTTOM[t] = newreg; Plane Surface(IBOTTOM[t]) = {iV1[t]};
Printf("IBOTTOM Isolant=%g ", IBOTTOM[t]);
/* Create Volumes */
Bord = newreg; Surface Loop(Bord) = {SRI[], SRE[], iSB[], iSH[], ITOP[t], IBOTTOM[t]};
Kapton[t] = newreg; Volume(Kapton[t]) = {Bord};
/* Update HChannels */
Bord_HChannel[num_h] = ITOP[t]; num_h++;
Bord_HChannel[num_h] = IBOTTOM[t]; num_h++;
If ( Canal )
For j In {3 : 4}
Bord_RChannel[num_r] = iSB[j]; num_r++; //Printf("Bord_RChannel : %g", SB[j]);
Bord_RChannel[num_r] = iSH[j]; num_r++; //Printf("Bord_RChannel : %g", SH[j]);
EndFor
EndIf
For j In {0:#SRI[]-1}
Bord_RI[num_bord_r] = SRI[j];
Bord_RE[num_bord_r] = SRE[j];
num_bord_r++;
EndFor
// reset isolant indicator
in_isolant = 0;
Printf("GeomIsolant done");
Return
/*-----------------------------------------------
DEFINITION DES "PHYSICALS"
-------------------------------------------------*/
Function Physical_Region
Printf("Defining Physical");
// Rint and R_ext
Physical Surface("Rint") = {Bord_RI[]} ; //Rint_Cuivre
Physical Surface("Rext") = {Bord_RE[]} ; //Rext_Cuivre
// HChannels
Physical Surface("Hchannels")={Bord_HChannel[]};
If ( Canal )
Physical Surface("Rchannels")={Bord_RChannel[]};
EndIf
If (ends)
Physical Surface("Neumann")={Bord_Neumann[]};
EndIf
// Interface Isolants
If(!isolant)
Physical Surface("Isolants")={Bord_Isolant[]};
EndIf
Physical Volume("Copper") = {Copper[]}; // SPIRE
If ( isolant )
Physical Volume("Kapton") = {Kapton[]}; // KAPTON
EndIf
// Top_Cu
TopCu=100;
For t In {Start_spire_number:#TOP[]-1}
Physical Surface(TopCu)={TOP[t]}; TopCu++;
EndFor
// Bottom_Cu
BottomCu=200;
For t In {Start_spire_number:#BOTTOM[]-1}
Physical Surface(BottomCu)={BOTTOM[t]};BottomCu++;
EndFor
Return
Function Virtual_Sphere
Z_C = (z+z0)/2.;
Printf("z0=%g", z0);
Printf("z=%g", z);
Printf("z+z0=%g", Fabs[z-z0]);
Printf("Z_C=%g", (z+z0)/2.);
// Define Points
P0=newp; Point(P0) = {0, 0, Z_C, LcAxis};
P1=newp; Point(P1) = {R * Cos[alpha[0]], R * Sin[alpha[0]], Z_C, LcBox};
P2=newp; Point(P2) = {R * Cos[alpha[#alpha[]-1]], R *
Sin[alpha[#alpha[]-1]], Z_C, LcBox};
P3=newp; Point(P3) = {0, 0, Z_C+R, LcBox};
P4=newp; Point(P4) = {0, 0, Z_C-R, LcBox};
C12=newl; Circle(C12) = {P1, P0, P2};
C23=newl; Circle(C23) = {P2, P0, P3};
C31=newl; Circle(C31) = {P3, P0, P1};
C14=newl; Circle(C14) = {P1, P0, P4};
C24=newl; Circle(C24) = {P2, P0, P4};
//L01=newl; Line(L01) = {P0, P1};
//L02=newl; Line(L02) = {P0, P2};
L03=newl; Line(L03) = {P0, P3};
L04=newl; Line(L04) = {P0, P4};
Hemi1=newl; Line Loop(Hemi1) = {C12, C23, C31};
HemiS1=newl; Ruled Surface(HemiS1) = {Hemi1};
Printf("HemiS1=%g", HemiS1);
lb230=newreg; Line Loop(lb230)={-C31, -L03, L04, -C14};
Lb230=newreg; Plane Surface(Lb230)={lb230, -V0[]};
Printf("Lb230=%g", Lb230);
Hemi2=newl; Line Loop(Hemi2) = {C12, C24, -C14};
HemiS2=newl; Ruled Surface(HemiS2) = {Hemi2};
Printf("HemiS2=%g", HemiS2);
lb380=newreg; Line Loop(lb380)={-C23, C24, -L04, L03};
Lb380=newreg; Plane Surface(Lb380)={lb380, -V1[]};
Printf("Lb380=%g", Lb380);
box = newreg ;
If(isolant)
Surface Loop(box) = {HemiS1, HemiS2, Lb230, Lb380, Bord_RI[], Bord_RE[],
Bord_HChannel[], Bord_RChannel[], Bord_Neumann[], Bord_Isolant[]};
EndIf
If (!isolant)
Surface Loop(box) = {HemiS1, HemiS2, Lb230, Lb380, Bord_RI[], Bord_RE[],
Bord_HChannel[], Bord_RChannel[], Bord_Neumann[]};
EndIf
Printf("box=%g", box);
Box = newreg ; Volume(Box) = {box} ;
Printf("Box=%g", Box);
Physical Volume("Air") = {Box}; //Volume Box1
// Define Shell Points
X_R_S = (R+Delta_R) * Cos[alpha[0]];
Y_R_S = (R+Delta_R) * Sin[alpha[0]];
P10=newp; Point(P10) = {(R+Delta_R) * Cos[alpha[0]], (R+Delta_R) *
Sin[alpha[0]], Z_C, LcInf};
P20=newp; Point(P20) = {(R+Delta_R) * Cos[alpha[#alpha[]-1]], (R+Delta_R)
* Sin[alpha[#alpha[]-1]], Z_C, LcInf};
P30=newp; Point(P30) = {0, 0, Z_C+(R+Delta_R), 5*LcInf};
P40=newp; Point(P40) = {0, 0, Z_C-(R+Delta_R), 5*LcInf};
C1020=newl; Circle(C1020) = {P10, P0, P20};
C2030=newl; Circle(C2030) = {P20, P0, P30};
C3010=newl; Circle(C3010) = {P30, P0, P10};
C1040=newl; Circle(C1040) = {P10, P0, P40};
C2040=newl; Circle(C2040) = {P20, P0, P40};
L110=newl; Line(L110) = {P1, P10};
L220=newl; Line(L220) = {P2, P20};
L330=newl; Line(L330) = {P3, P30};
L440=newl; Line(L440) = {P4, P40};
HemiR1=newl; Line Loop(HemiR1) = {C1020, C2030, C3010};
HemiRS1=newl; Ruled Surface(HemiRS1) = {HemiR1};
HemiR2=newl; Line Loop(HemiR2) = {C1020, C2040, -C1040};
HemiRS2=newl; Ruled Surface(HemiRS2) = {HemiR2};
OxOz1=newl; Line Loop(OxOz1) = {C14, L440, -C1040, -L110};
Sym_OxOz1=newl; Plane Surface(Sym_OxOz1) = {OxOz1};
OxOz2=newl; Line Loop(OxOz2) = {C31, L110, -C3010, -L330};
Sym_OxOz2=newl; Plane Surface(Sym_OxOz2) = {OxOz2};
OyOz1=newl; Line Loop(OyOz1) = {-C24, L220, C2040, -L440};
Sym_OyOz1=newl; Plane Surface(Sym_OyOz1) = {OyOz1};
OyOz2=newl; Line Loop(OyOz2) = {C23, L330, -C2030, -L220};
Sym_OyOz2=newl; Plane Surface(Sym_OyOz2) = {OyOz2};
shell=newreg; Surface Loop(shell) = {Sym_OxOz1, Sym_OxOz2, HemiRS1,
HemiRS2, HemiS1, HemiS2, Sym_OyOz1, Sym_OyOz2};
Shell = newreg ; Volume(Shell) = {shell} ;
////////PHYSICAL REGION /////////
Physical Surface("Infinity") = {HemiRS1, HemiRS2} ; // Infini
Physical Surface("Symetrie1") = {Lb230, Sym_OxOz1, Sym_OxOz2} ; // Sym OxOz
Physical Surface("Symetrie2") = {Lb380, Sym_OyOz1, Sym_OyOz2} ; // Sym
OyOzFace Bas Box 2
Physical Volume("AirInf") = {Shell}; //Volume Box2
Return
_______________________________________________
gmsh mailing list
[email protected]
http://www.geuz.org/mailman/listinfo/gmsh
