diff --git a/Conductors3D/w2d.geo b/Conductors3D/w2d.geo
new file mode 100644
index 0000000000000000000000000000000000000000..daa57e2d257728ae71e93e5a43d2469eebc6d4c6
--- /dev/null
+++ b/Conductors3D/w2d.geo
@@ -0,0 +1,125 @@
+Include "wires_common.pro";
+
+// global Gmsh options
+Mesh.Optimize = 1; // optimize quality of tetrahedra
+Mesh.VolumeEdges = 0; // hide volume edges
+Geometry.ExactExtrusion = 0; // to allow rotation of extruded shapes
+Solver.AutoMesh = 2; // always remesh if necessary (don't reuse mesh on disk)
+
+lc1 = (Flag_Thin==0) ? MeshSizeWire*mm : rs ;
+lc2 = box/10; // mesh size at outer surface
+
+llWires[] = {};
+centerWires[] = {};
+surfWires[] = {};
+
+If( Flag_Thin != 1)
+// then circles are in the geometry
+// their radius is R= WX~{i} in the real case (Flag_Thin == 0)
+// or R=rs in the "regular sleeve case" (Flag_Thin == 2)
+ For i In {1:NumWires}
+ R = (Flag_Thin == 0 ) ? rw : rs;
+ pC = newp; Point(pC) = {WX~{i} , WY~{i} , 0, lc1};
+ p1 = newp; Point(p1) = {WX~{i}+R, WY~{i} , 0, lc1};
+ p2 = newp; Point(p2) = {WX~{i} , WY~{i}+R, 0, lc1};
+ p3 = newp; Point(p3) = {WX~{i}-R, WY~{i} , 0, lc1};
+ p4 = newp; Point(p4) = {WX~{i} , WY~{i}-R, 0, lc1};
+ c1 = newl; Circle(c1) = {p1,pC,p2};
+ c2 = newl; Circle(c2) = {p2,pC,p3};
+ c3 = newl; Circle(c3) = {p3,pC,p4};
+ c4 = newl; Circle(c4) = {p4,pC,p1};
+ l1 = newl; Line(l1) = {pC,p1};
+ l2 = newl; Line(l2) = {pC,p2};
+ l3 = newl; Line(l3) = {pC,p3};
+ l4 = newl; Line(l4) = {pC,p4};
+ ll1 = newll; Line Loop(ll1) = {l1,c1,-l2};
+ s1 = news; Plane Surface(s1) = {ll1};
+ ll2 = newll; Line Loop(ll2) = {l2,c2,-l3};
+ s2 = news; Plane Surface(s2) = {ll2};
+ ll3 = newll; Line Loop(ll3) = {l3,c3,-l4};
+ s3 = news; Plane Surface(s3) = {ll3};
+ ll4 = newll; Line Loop(ll4) = {l4,c4,-l1};
+ s4 = news; Plane Surface(s4) = {ll4};
+
+ llWires[] += {ll1,ll2,ll3,ll4};
+ surfWires~{i-1}[] = {s1,s2,s3,s4};
+ surfWires[] += surfWires~{i-1}[];
+ centerWires[] += {pC};
+ EndFor
+Else
+ For i In {1:NumWires}
+ pC = newp; Point(pC) = {WX~{i} , WY~{i} , 0, lc1};
+ centerWires[] += {pC};
+ EndFor
+EndIf
+
+// create air box around wires
+BoundingBox; // recompute model bounding box
+cx = (General.MinX + General.MaxX) / 2;
+cy = (General.MinY + General.MaxY) / 2;
+cz = (General.MinZ + General.MaxZ) / 2;
+
+/* lx = 2*inf + General.MaxX - General.MinX; */
+/* ly = 2*inf + General.MaxY - General.MinZ; */
+lx = 2*box;
+ly = 2*box;
+
+p1 = newp; Point (p1) = {cx-lx/2, cy-ly/2, 0, lc2};
+p2 = newp; Point (p2) = {cx+lx/2, cy-ly/2, 0, lc2};
+p3 = newp; Point (p3) = {cx+lx/2, cy+ly/2, 0, lc2};
+p4 = newp; Point (p4) = {cx-lx/2, cy+ly/2, 0, lc2};
+l1 = newl; Line(l1) = {p1, p2};
+l2 = newl; Line(l2) = {p2, p3};
+l3 = newl; Line(l3) = {p3, p4};
+l4 = newl; Line(l4) = {p4, p1};
+ll1 = newll; Line Loop(ll1) = {l1,l2,l3,l4};
+Physical Line("INF", 2) = { l1, l2, l3, l4 };
+
+If( !Flag_3D )
+
+ If( Flag_Thin == 0 ) // Real
+ s1 = news; Plane Surface(s1) = {ll1,-llWires[]};
+ Physical Surface("AIR", 1) = { s1 };
+ For i In {1:NumWires}
+ Physical Surface(Sprintf("VWIRE_%g",i), 10+i) = { surfWires~{i-1}[] };
+ Physical Point (Sprintf("LWIRE_%g",i), 50+i) = { centerWires[i-1] };
+ EndFor
+ ElseIf( Flag_Thin == 1 ) // 1D raw sleeve
+ s1 = news; Plane Surface(s1) = {ll1};
+ Physical Surface("AIR", 1) = { s1 };
+ For i In {1:NumWires}
+ Physical Point (Sprintf("LWIRE_%g",i), 50+i) = { centerWires[i-1] };
+ Point { centerWires[i-1] } In Surface { s1 };
+ EndFor
+ Else // 1D regular sleeve
+ s1 = news; Plane Surface(s1) = {ll1,-llWires[]};
+ Physical Surface("AIR", 1) = { s1, surfWires[] };
+ For i In {1:NumWires}
+ Physical Point (Sprintf("LWIRE_%g",i), 50+i) = { centerWires[i-1] };
+ Point { centerWires[i-1] } In Surface { s1 };
+ EndFor
+ EndIf
+
+Else
+
+ For i In {1:NumWires}
+ e[] = Extrude {0,0,Lz} { Surface{ surfWires[i-1] } ; Layers{1}; Recombine; } ;
+ Printf("%g", #e[]);
+ Printf("%g %g %g %g %g %g", e[0], e[1], e[2], e[3], e[4], e[5]);
+ Physical Volume (Sprintf("VWIRE_%g",i), 10+i) = { e[1] };
+ Physical Surface(Sprintf("SANODE_%g",i), 20+i) = { surfWires[i-1] };
+ Physical Surface(Sprintf("SCATHODE_%g",i), 30+i) = { e[0] };
+ e[] = Extrude {0,0,Lz} { Point{ centerWires[i-1] } ; Layers{1}; Recombine; } ;
+ Physical Line (Sprintf("LWIRE_%g",i), 50+i) = { e[1] };
+ Physical Point (Sprintf("PANODE_%g",i), 60+i) = { centerWires[i-1] };
+ Physical Point (Sprintf("PCATHODE_%g",i), 70+i) = { e[0] };
+ EndFor
+
+ e[] = Extrude {0,0,Lz} { Surface{ s1 } ; Layers{1}; Recombine; } ;
+ Physical Volume ("AIR", 1) = { e[1] };
+ Physical Surface("INF", 2) = { s1, e[0], e[2], e[3], e[4], e[5] };
+ //Physical Surface("INF", 1) = { CombinedBoundary{ Volume{ e[1]}; } };
+
+EndIf
+
+
diff --git a/Conductors3D/w2d.pro b/Conductors3D/w2d.pro
new file mode 100644
index 0000000000000000000000000000000000000000..3205664ae70ff57fc67ecd97ab1240b54dc32347
--- /dev/null
+++ b/Conductors3D/w2d.pro
@@ -0,0 +1,707 @@
+Include "wires_common.pro";
+
+Struct S::SOLVER [ Enum, mumps, gmres, fgmres,bcgs ];
+DefineConstant[
+ FileId = ""
+ DataDir = StrCat["data/" , FileId]
+ ElekinDir = StrCat["Elekin_Line_Const/" , FileId]
+ MagDynDir = StrCat["MagDyn_Line_Const/" , FileId]
+ FullWaveDir = StrCat["FullWave_Line_Const/", FileId]
+ type_solver = {
+ S::SOLVER.mumps,
+ Choices {
+ S::SOLVER.mumps ="mumps",
+ S::SOLVER.gmres="gmres",
+ S::SOLVER.fgmres= "fgmres",
+ S::SOLVER.bcgs="bcgs"
+ },
+ Name "Solver/type_solver"
+ }
+];
+
+DefineConstant[
+ R_ = {"Dynamic", Name "GetDP/1ResolutionChoices", Visible 0},
+ C_ = {"-sol -pos", Name "GetDP/9ComputeCommand", Visible 0},
+ P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}
+ ResDir = "Res/"
+];
+
+
+Group{
+ // Geometrical regions
+ AIR = Region[ 1 ];
+ INF = Region[ 2 ];
+
+ VWIRES = Region[ {} ];
+ LWIRES = Region[ {} ];
+ For i In {1:NumWires}
+ LWIRE~{i} = Region[ (50+i) ];
+ LWIRES += Region[ LWIRE~{i} ];
+ If( Flag_Thin == 0 )
+ VWIRE~{i} = Region[ (10+i) ];
+ VWIRES += Region[ VWIRE~{i} ];
+ SWIRE~{i} = Region[ {} ];
+ Else
+ VWIRE~{i} = Region[ {} ];
+ SWIRE~{i} = ElementsOf[ AIR, OnOneSideOf LWIRE~{i} ];
+ EndIf
+ EndFor
+
+ // Abstract regions
+ // Cette formulation peut contenir des conducteurs volumiques Vol_C
+ // et lin�iques lVol_C.
+
+ If( Flag_Thin == 0 )
+ Vol_C = Region[ VWIRES ];
+ Vol_CC = Region[ AIR ];
+ lVol_C = Region[ {} ];
+ Else
+ Vol_C = Region[ {} ];
+ Vol_CC = Region[ {VWIRES, AIR} ];
+ lVol_C = Region[ LWIRES ];
+ EndIf
+
+ Vol_nu = Region[ {Vol_C, Vol_CC, lVol_C} ];
+ Sur_Dirichlet_a = Region[ INF ];
+
+ // 3D
+ If( Flag_3D )
+ Sur_Dirichlet_a += Region[ { BOT, TOP} ];
+ EndIf
+ Sur_Complete_Tree = Region[ { lVol_C, Sur_Dirichlet_a } ];
+ Electrodes = Region[ {} ];
+
+ Dom_Hcurl_a = Region[ Vol_nu ];
+ Dom_Hthin_a = ElementsOf[ Vol_nu, OnOneSideOf lVol_C ]; // Sleeve
+ Dom_Hgrad_u = Region[ Vol_C ];
+ Dom_Hregion_i = Region[ lVol_C ];
+
+ // For integration on the volume elements (excluding line elements) of the sleeve
+ Vol_Sleeve = ElementsOf[ {Vol_CC}, OnOneSideOf lVol_C ];
+}
+
+Function{
+ omega = 2*Pi*Freq;
+ mu0 = Pi*4e-7;
+ mur = 1.;
+ sigma = 5.96e7;
+
+ i[] = Complex[0,1];
+
+ mu[] = mu0;
+ nu[] = 1/mu[];
+ sigma[] = sigma;
+ skin_depth = Sqrt[ 2/(omega*sigma*mu0*mur) ];
+ tau[] = Complex[-1,1]/skin_depth*rw;
+ J0[] = JnComplex[0,$1];
+ J1[] = JnComplex[1,$1];
+ dJ1[] = JnComplex[0,$1]-JnComplex[1,$1]/$1;
+
+ For i In {1:NumWires}
+ // Current density
+ J[ Region[ {VWIRE~{i}, LWIRE~{i}} ] ] =
+ Vector[ 0, 0, WI~{i}/A_c * Exp[ i[]*WP~{i}*deg ] ];
+
+ // local radial coordinate with origin on wire i
+ R~{i}[] = Hypot[ X[]-WX~{i} , Y[]-WY~{i} ];
+ EndFor
+
+ // shape function of the current, wI[r]
+ wI[] = tau[] * J0[tau[]*$1/rw] / J1[tau[]] /(2*A_c);
+
+ //radial analytical solutions with a zero flux condition imposed at $1(=r)=rs
+ Analytic_A[] = // per Amp
+ (($1>rw) ? mu0/(2*Pi) * Log[rs/$1] :
+ mu0/(2*Pi) * Log[rs/rw] + i[]*( wI[$1]-wI[rw] )/(sigma*omega) );
+ AnalyticStatic_A[] = mu0/(2*Pi) * // per Amp
+ (($1>rw) ? Log[rs/$1] : Log[rs/rw] + mur*(1-($1/rw)^2)/2 );
+
+ // Impedance of thin wire p.u. length
+ R_DC = 1./(sigma*A_c);
+ Analytic_R[] = Re[ wI[rw]/sigma ];
+ Analytic_L[] = mu0/(2*Pi) * Log[rs/rw] + Re[ wI[rw]/(i[]*omega*sigma) ];
+
+ If( NumWires == 1 )
+ Exact_A[] = WI_1*(Analytic_A[R_1[]] + mu0/(2*Pi) * Log[rout/rs]) ;
+ Correction_A[] = ((R_1[]<rs) ? WI_1*Analytic_A[R_1[]] : 0 ) ;
+ EndIf
+ If ( NumWires == 2 )
+ Exact_A[] = WI_1*Analytic_A[R_1[]] + WI_2*Analytic_A[R_2[]]
+ + (WI_1+WI_2)*mu0/(2*Pi) * NumWires * Log[rout/rs] ;
+ Correction_A[] = ((R_1[]<rs) ? WI_1*Analytic_A[R_1[]] :
+ ((R_2[]<rs) ? WI_2*Analytic_A[R_2[]] : 0 ));
+ EndIf
+ If ( NumWires == 3 )
+ Exact_a[] = WI_1*Analytic_A[R_1[]] + WI_2*Analytic_A[R_2[]]
+ + WI_3*Analytic_A[R_3[]] + (WI1_+WI_2+WI_3)*mu0/(2*Pi) * NumWires * Log[rout/rs] ;
+ Correction_A[] = ((R_1[]<rs) ? WI_1*Analytic_A[R_1[]] :
+ ((R_2[]<rs) ? WI_2*Analytic_A[R_2[]] :
+ ((R_3[]<rs) ? WI_3*Analytic_A[R_3[]] : 0 )));
+ EndIf
+}
+
+
+Jacobian {
+ { Name Vol ;
+ Case {
+ { Region lVol_C ; Jacobian Lin ; Coefficient A_c; }
+ { Region All ; Jacobian Vol ; }
+ }
+ }
+}
+
+Integration {
+ { Name I1 ;
+ Case {
+ { Type Gauss ;
+ Case {
+ { GeoElement Point ; NumberOfPoints 1 ; }
+ { GeoElement Line ; NumberOfPoints 3 ; }
+ { GeoElement Triangle ; NumberOfPoints 3 ; }
+ { GeoElement Quadrangle ; NumberOfPoints 4 ; }
+ { GeoElement Prism ; NumberOfPoints 21 ; }
+ { GeoElement Tetrahedron ; NumberOfPoints 4 ; }
+ { GeoElement Hexahedron ; NumberOfPoints 6 ; }
+ { GeoElement Pyramid ; NumberOfPoints 8 ; }
+ }
+ }
+ }
+ }
+}
+
+Constraint{
+ // 2D
+
+ { Name Hcurl_a_2D_ac; Type Assign;
+ Case {
+ { Region Region[ Sur_Dirichlet_a ]; Value 0; }
+ }
+ }
+
+ { Name Impose_I ;
+ Case {
+ For i In {1:NumWires}
+ { Region VWIRE~{i} ; Value WI~{i} ; }
+ { Region LWIRE~{i} ; Value WI~{i} ; }
+ EndFor
+ }
+ }
+
+ // 3D
+ { Name Hcurl_a_3D_ac; Type Assign;
+ Case {
+ { Region Region[ Sur_Dirichlet_a ]; Value 0.; }
+ }
+ }
+ { Name GaugeCondition_a; Type Assign;
+ Case {
+ { Region Vol_CC; SubRegion Sur_Complete_Tree; Value 0.; }
+ }
+ }
+}
+
+FunctionSpace {
+ // 2D
+
+ { Name Hcurl_a_2D; Type Form1P;
+ BasisFunction {
+ { Name sc; NameOfCoef ac; Function BF_PerpendicularEdge;
+ Support Dom_Hcurl_a; Entity NodesOf[All]; }
+ { Name sw; NameOfCoef aw; Function BF_PerpendicularEdge;
+ Support Dom_Hthin_a; Entity NodesOf[lVol_C]; }
+ }
+ SubSpace {
+ { Name ac ; NameOfBasisFunction { sc } ; }
+ { Name aw ; NameOfBasisFunction { sw } ; }
+ }
+ Constraint {
+ { NameOfCoef ac;
+ EntityType Auto; NameOfConstraint Hcurl_a_2D_ac; }
+ }
+ }
+
+ { Name Hgrad_u_2D ; Type Form1P ;
+ BasisFunction {
+ { Name sr ; NameOfCoef ur ; Function BF_RegionZ ;
+ Support Dom_Hgrad_u ; Entity Vol_C ; }
+ }
+ GlobalQuantity {
+ { Name U ; Type AliasOf ; NameOfCoef ur ; }
+ { Name I ; Type AssociatedWith ; NameOfCoef ur ; }
+ }
+ Constraint {
+ { NameOfCoef I ; EntityType GroupsOfNodesOf ; NameOfConstraint Impose_I ; }
+ }
+ }
+
+ { Name Hregion_i_2D; Type Vector;
+ BasisFunction {
+ { Name sr; NameOfCoef ir; Function BF_RegionZ;
+ Support Dom_Hregion_i; Entity lVol_C; }
+ }
+ GlobalQuantity {
+ { Name Is; Type AliasOf; NameOfCoef ir; }
+ { Name Us; Type AssociatedWith; NameOfCoef ir; }
+ }
+ Constraint {
+ { NameOfCoef Is;
+ EntityType Region; NameOfConstraint Impose_I; }
+ }
+ }
+
+ // 3D
+ { Name Hcurl_a_3D; Type Form1;
+ BasisFunction {
+ { Name sc; NameOfCoef ac; Function BF_Edge;
+ Support Dom_Hcurl_a; Entity EdgesOf[All]; }
+ { Name sw; NameOfCoef aw; Function BF_Edge;
+ Support Dom_Hthin_a; Entity EdgesOf[lVol_C]; }
+ }
+ SubSpace {
+ { Name ac ; NameOfBasisFunction { sc } ; }
+ { Name aw ; NameOfBasisFunction { sw } ; }
+ }
+ Constraint {
+ { NameOfCoef ac;
+ EntityType Auto; NameOfConstraint Hcurl_a_3D_ac; }
+ }
+ }
+
+ { Name Hgrad_u_3D; Type Form0;
+ BasisFunction {
+ { Name sn; NameOfCoef un; Function BF_Node;
+ Support Dom_Hgrad_u; Entity NodesOf[ Vol_C, Not Electrodes]; }
+ { Name sn2; NameOfCoef un2; Function BF_GroupOfNodes;
+ Support Dom_Hgrad_u; Entity GroupsOfNodesOf[ Electrodes ]; }
+ }
+ GlobalQuantity {
+ { Name U; Type AliasOf ; NameOfCoef un2; }
+ { Name I; Type AssociatedWith; NameOfCoef un2; }
+ }
+ Constraint {
+ //{ NameOfCoef U; EntityType Auto; NameOfConstraint Impose_U; }
+ { NameOfCoef I; EntityType Auto; NameOfConstraint Impose_I; }
+ }
+ }
+}
+
+Formulation {
+ { Name MagnetoDynamics; Type FemEquation;
+ If( !Flag_3D )
+ Quantity {
+ { Name a; Type Local; NameOfSpace Hcurl_a_2D; }
+ { Name ac; Type Local; NameOfSpace Hcurl_a_2D[ac]; }
+ { Name aw; Type Local; NameOfSpace Hcurl_a_2D[aw]; }
+
+ { Name dv; Type Local; NameOfSpace Hgrad_u_2D; }
+ { Name U; Type Global; NameOfSpace Hgrad_u_2D [U]; }
+ { Name I; Type Global; NameOfSpace Hgrad_u_2D [I]; }
+
+ { Name i; Type Local; NameOfSpace Hregion_i_2D; }
+ { Name Is; Type Global; NameOfSpace Hregion_i_2D [Is]; }
+ { Name Us; Type Global; NameOfSpace Hregion_i_2D [Us]; }
+ }
+
+ Equation {
+
+ Integral { [ nu[] * Dof{d ac} , {d ac} ];
+ In Vol_nu; Jacobian Vol; Integration I1; }
+ Integral { [ -Dof{i}/A_c , {ac} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+
+ Integral { [ nu[] * Dof{d aw} , {d aw} ];
+ In Vol_nu; Jacobian Vol; Integration I1; }
+ Integral { [-Dof{i}/A_c , {aw} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+
+ /*
+ Integral { [ -J[] , {a} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ */
+
+ Integral { DtDof[ sigma[] * Dof{a} , {a} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ Integral { [ sigma[] * Dof{dv} , {a} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+
+ Integral { DtDof[ sigma[] * Dof{a} , {dv} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ Integral { [ sigma[] * Dof{dv} , {dv} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ GlobalTerm { [ Dof{I} , {U} ]; In Vol_C; }
+
+ GlobalTerm { [ Dof{Us} , {Is} ];
+ In lVol_C; }
+ GlobalTerm { [ Analytic_R[] * Dof{Is} , {Is} ];
+ In lVol_C;}
+ Integral { DtDof [ Dof{ac}/A_c , {i} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+ If( Flag_Thin != 3)
+ Integral { DtDof [ -Dof{aw}/A_c , {i} ];
+ In lVol_C; Jacobian Vol; Integration I1;}
+ GlobalTerm { DtDof [ Analytic_L[] * Dof{Is} , {Is} ];
+ In lVol_C;}
+ EndIf
+ }
+ Else // ############### 3D
+ Quantity {
+ { Name a; Type Local; NameOfSpace Hcurl_a_3D; }
+ { Name ac; Type Local; NameOfSpace Hcurl_a_3D[ac]; }
+ { Name aw; Type Local; NameOfSpace Hcurl_a_3D[aw]; }
+ // { Name v; Type Local; NameOfSpace Hgrad_u_3D; }
+ // { Name U; Type Global; NameOfSpace Hgrad_u_3D [U]; }
+ // { Name I; Type Global; NameOfSpace Hgrad_u_3D [I]; }
+ }
+
+ Equation {
+ /*
+ Integral { [ nu[] * Dof{d ac} , {d ac} ];
+ In Vol_nu; Jacobian Vol; Integration I1; }
+ Integral { [ nu[] * Dof{d aw} , {d aw} ];
+ In Vol_nu; Jacobian Vol; Integration I1; }
+
+ Integral { [ -Vector[ 0, 0, I0/A_c] , {ac} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+ Integral { [ -Vector[ 0, 0, I0/A_c] , {aw} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+
+ Integral { [ -Vector[ 0, 0, I0/A_c] , {a} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+
+ Integral { DtDof[ sigma[] * Dof{a} , {a} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ Integral { [ sigma[] * Dof{d v} , {a} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ Integral { DtDof[ sigma[] * Dof{a} , {d v} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ Integral { [ sigma[] * Dof{d v} , {d v} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ GlobalTerm { [Dof{U} , {I} ]; In Electrodes; }
+ */
+
+ }
+ EndIf
+ }
+}
+
+Resolution {
+ { Name Dynamic;
+ System {
+ { Name S; NameOfFormulation MagnetoDynamics;
+ Type ComplexValue; Frequency Freq; }
+ }
+ Operation {
+ CreateDir[DataDir];
+
+ // SOLVER OPTION
+ If(type_solver == S::SOLVER.gmres)
+ SetGlobalSolverOptions[
+ "-pc_type ilu -ksp_type gmres -ksp_max_it 1000 -ksp_rtol 1e-6 -ksp_atol 1e-6"];
+ ElseIf(type_solver == S::SOLVER.fgmres)
+ SetGlobalSolverOptions[
+ "-pc_type ilu -ksp_type fgmres -ksp_max_it 1000 -ksp_rtol 1e-6 -ksp_atol 1e-6"];
+ ElseIf(type_solver == S::SOLVER.bcgs)
+ SetGlobalSolverOptions[
+ "-pc_type ilu -ksp_type bcgs -ksp_max_it 1000 -ksp_rtol 1e-6 -ksp_atol 1e-6"];
+ EndIf
+ SetGlobalSolverOptions["-petsc_prealloc 900"];
+
+ Generate[S]; Solve[S]; SaveSolution[S];
+
+ PostOperation[integaz];
+
+ If( Flag_Maps == 1 )
+ PostOperation[mapaz];
+ PostOperation[cutaz];
+
+ DeleteFile["U_full.dat"];
+ DeleteFile["U_raw.dat"];
+ DeleteFile["U_reg.dat"];
+ DeleteFile["U_naive.dat"];
+ DeleteFile["joule_full.dat"];
+ DeleteFile["joule_raw.dat"];
+ DeleteFile["joule_reg.dat"];
+ DeleteFile["joule_naive.dat"];
+ EndIf
+
+ If( Flag_Thin != 0 )
+ For i In {1:NumWires}
+ Print[ {i, rw, Sqrt[ $sarea~{i}/Pi], rs, $sarea~{i}, Pi*rs^2, skin_depth, $intby~{i}/$sarea~{i}},
+ Format "WIRE %2g: rw = %7.3e rs = %7.3e (%7.3e) as = %7.3e (%7.3e) delta = %7.3e bave = %7.3e" ] ;
+ EndFor
+ EndIf
+ }
+ }
+}
+
+PostProcessing {
+ { Name PostProcessings; NameOfFormulation MagnetoDynamics;
+ PostQuantity {
+ { Name az;
+ Value { Local { [ CompZ[{ac}] - (Flag_Thin!=0)*(CompZ[{aw}]-Correction_A[])];
+ In Dom_Hcurl_a; Jacobian Vol; }}}
+ { Name acz;
+ Value { Local { [ CompZ[ {ac}] ];
+ In Dom_Hcurl_a; Jacobian Vol; }}}
+ { Name awz;
+ Value { Local { [ CompZ[ {aw}] ];
+ In Dom_Hthin_a; Jacobian Vol; }}}
+ { Name acwz;
+ Value { Local { [ CompZ[ {ac} - {aw} ] ];
+ In Dom_Hcurl_a; Jacobian Vol; }}}
+ { Name exact;
+ Value { Local { [ Exact_A[] ] ;
+ In Dom_Hcurl_a; Jacobian Vol; }}}
+ { Name corr;
+ Value { Local { [ Correction_A[] ] ;
+ In Dom_Hcurl_a; Jacobian Vol; }}}
+ { Name bcw;
+ Value { Local { [ {d ac} - {d aw} ];
+ In Dom_Hcurl_a; Jacobian Vol; }}}
+
+ { Name j;
+ Value { Local { [ -sigma[] * ( Dt[{a}] + {dv} ) ];
+ In Vol_C; Jacobian Vol; }}}
+
+ { Name un;
+ Value{ Local { [ 1 ] ;
+ In Vol_nu; Integration I1 ; Jacobian Vol; }}}
+ { Name surf;
+ Value{ Integral{ [ 1 ] ;
+ In Vol_nu; Integration I1 ; Jacobian Vol; }}}
+
+
+ If( Flag_Thin == 0 )
+ { Name flux;
+ Value { Integral { [ CompZ[ {a} ] / A_c ];
+ In Vol_C; Integration I1 ; Jacobian Vol; }}}
+ { Name JouleLosses;
+ Value { Integral { [ sigma[]*SquNorm[ Dt[{a}] + {dv} ] ] ;
+ In Vol_C ; Integration I1 ; Jacobian Vol ; }}}
+ { Name U;
+ Value { Term { [ Cart2Pol [ {U} ] ]; In Vol_C; }}}
+ { Name L;
+ Value { Term { [ Im [ {U}/{I}/omega ] ]; In Vol_C; }}}
+
+ Else
+ { Name flux; // naive flux
+ Value { Integral { [ CompZ[ {ac} ] ];
+ In lVol_C; Integration I1 ; Jacobian Vol; }}}
+
+ { Name JouleLosses;
+ Value { Term { [ Analytic_R[]*{Is}*{Is} ] ; In lVol_C; }}}
+ { Name U;
+ Value { Term { [ Cart2Pol [ {Us} ] ] ; In lVol_C; }}}
+ { Name L; // only valid with one current
+ Value { Term { [ Im [ {Us}/{Is}/omega ] ] ; In lVol_C; }}}
+
+ For i In {1:NumWires}
+ { Name area~{i};
+ Value{ Integral{ [ 1 ] ;
+ In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
+ { Name intbx~{i}; // integral of b, 1st step for computing b average
+ Value { Integral { [ Norm [ CompX[ {d ac} - {d aw} ] ] ] ;
+ In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
+ { Name intby~{i}; // integral of b, 1st step for computing b average
+ Value { Integral { [ Norm [ CompY[ {d ac} - {d aw} ] ] ] ;
+ In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
+ EndFor
+
+ EndIf
+ }
+ }
+}
+
+PostOperation mapaz UsingPost PostProcessings {
+ Print [az, OnElementsOf Vol_nu, File "az.pos"];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].IntervalsType = 3;",
+ "View[l].NbIso = 30;",
+ "View[l].RaiseZ = 15000;",
+ "View[l].NormalRaise = 0;"],
+ File "tmp.geo", LastTimeStepOnly];
+
+ If(Flag_Thin != 0)
+ Print [bcw, OnElementsOf Vol_Sleeve, File "b.pos"];
+
+ Print [acwz, OnElementsOf Vol_nu, File "acz.pos"];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].IntervalsType = 3;",
+ "View[l].NbIso = 30;",
+ "View[l].RaiseZ = 15000;",
+ "View[l].NormalRaise = 0;"],
+ File "tmp.geo", LastTimeStepOnly] ;
+ EndIf
+}
+
+
+NbPoints = 1000;
+Xmax = 0.15*box;
+
+PostOperation cutaz UsingPost PostProcessings {
+ // cut of az displayed in GUI
+ If( Flag_Thin == 0 )
+ Print [ az, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
+ Format Gmsh, File "Cut_az.pos" ];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].Name = 'cut az full';",
+ "View[l].Axes = 3;",
+ "View[l].LineWidth = 3;",
+ "View[l].Type = 2;"],
+ File "tmp.geo", LastTimeStepOnly];
+
+ Else
+
+ Print [ acz, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
+ Format Gmsh, File "Cut_az.pos" ];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].Name = 'cut az not corrected';",
+ "View[l].Axes = 3;",
+ "View[l].LineWidth = 3;",
+ "View[l].Type = 2;"],
+ File "tmp.geo", LastTimeStepOnly];
+
+ Print [ az, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
+ Format Gmsh, File "Cut_az.pos" ];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].Name = 'cut az corrected';",
+ "View[l].Axes = 3;",
+ "View[l].LineWidth = 3;",
+ "View[l].Type = 2;"],
+ File "tmp.geo", LastTimeStepOnly];
+ EndIf
+
+ // cuts in txt format for Gnuplot
+ Print [ exact, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
+ Format SimpleTable, File "Cut_analytic.txt" ];
+
+ If( Flag_Thin == 0 )
+ Print [ az, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
+ Format SimpleTable, File "Cut_fullfem.txt" ];
+ Else
+ Print [ az, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
+ Format SimpleTable, File "Cut_thinfem.txt" ];
+ Print [ acz, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
+ Format SimpleTable, File "Cut_ac.txt" ];
+ Print [ awz, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
+ Format SimpleTable, File "Cut_aw.txt" ];
+ Print [ acwz, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
+ Format SimpleTable, File "Cut_acw.txt" ];
+ Print [ corr, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
+ Format SimpleTable, File "Cut_correction.txt" ];
+ EndIf
+}
+
+
+PostOperation integaz UsingPost PostProcessings {
+ For i In {1:NumWires}
+
+ If( Flag_Thin == 0 )
+ Print[ U, OnRegion Vol_C, File > "U_full.dat", Format Table,
+ SendToServer Sprintf["}Voltage/Wire %g/1Full",i]];
+ Print[ L, OnRegion Vol_C, File > "L_full.dat", Format Table,
+ SendToServer Sprintf["}Inductance/Wire %g/1Full",i]];
+ Print[ JouleLosses[ VWIRE~{i}], OnGlobal,
+ Format Table, File > "joule_full.dat",
+ SendToServer Sprintf["}Joule Losses/Wire %g/1Full",i]];
+ EndIf
+ If( Flag_Thin == 1 )
+ Print[ U, OnRegion lVol_C, File > "U_raw.dat", Format Table,
+ SendToServer Sprintf["}Voltage/Wire %g/2Raw",i]];
+ Print[ L, OnRegion lVol_C, File > "L_raw.dat", Format Table,
+ SendToServer Sprintf["}Inductance/Wire %g/2Raw",i]];
+ Print[ JouleLosses, OnRegion LWIRE~{i},
+ Format Table, File > "joule_raw.dat",
+ SendToServer Sprintf["}Joule Losses/Wire %g/2Raw",i]];
+ EndIf
+ If( Flag_Thin == 2 )
+ Print[ U, OnRegion lVol_C, File > "U_reg.dat", Format Table,
+ SendToServer Sprintf["}Voltage/Wire %g/3Reg",i]];
+ Print[ L, OnRegion lVol_C, File > "L_reg.dat", Format Table,
+ SendToServer Sprintf["}Inductance/Wire %g/3Reg",i]];
+ Print[ JouleLosses, OnRegion LWIRE~{i},
+ Format Table, File > "joule_reg.dat",
+ SendToServer Sprintf["}Joule Losses/Wire %g/3Reg",i]];
+ EndIf
+ If( Flag_Thin == 3 )
+ Print[ U, OnRegion lVol_C, File > "U_naive.dat", Format Table,
+ SendToServer Sprintf["}Voltage/Wire %g/4Naive",i]];
+ Print[ L, OnRegion lVol_C, File > "L_naive.dat", Format Table,
+ SendToServer Sprintf["}Inductance/Wire %g/4naive",i]];
+ Print[ JouleLosses, OnRegion LWIRE~{i},
+ Format Table, File > "joule_naive.dat",
+ SendToServer Sprintf["}Joule Losses/Wire %g/4Naiv",i]];
+ EndIf
+
+ If( Flag_Thin != 0 )
+ Print[ area~{i}[SWIRE~{i}], OnGlobal,
+ Format Table, File >"zzz",
+ StoreInVariable $sarea~{i}];
+ Print[ intby~{i}[SWIRE~{i}], OnGlobal,
+ Format Table, File > "zzz",
+ StoreInVariable $intby~{i}];
+ EndIf
+
+ EndFor
+}
+
+
+
+
+ /*{ Name flux_corr; // corrected flux
+ Value {
+ Integral { [ CompZ[ {ac} - {aw} ] ];
+ In lVol_C; Integration I1 ; Jacobian Vol; }
+ Integral { [ Analytic_L[]*{Is} ];
+ In lVol_C; Integration I1 ; Jacobian Vol; }
+ }
+ }*/
+
+ /* Print [awz, OnElementsOf Vol_nu, File "awz.pos"];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].IntervalsType = 3;",
+ "View[l].NbIso = 30;",
+ "View[l].RaiseZ = 15000;",
+ "View[l].NormalRaise = 0;"],
+ File "tmp.geo", LastTimeStepOnly] ;
+
+ Print [acz, OnElementsOf Vol_nu, File "acz.pos"];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].IntervalsType = 3;",
+ "View[l].NbIso = 30;",
+ "View[l].RaiseZ = 15000;",
+ "View[l].NormalRaise = 0;"],
+ File "tmp.geo", LastTimeStepOnly] ;
+
+
+
+ Print[ flux[ VWIRE~{i} ], OnGlobal,
+ Format Table, File > "flux_full.dat",
+ SendToServer Sprintf["}Flux/Wire %g/1Full",i]];
+ Print[ flux[ LWIRE~{i} ], OnGlobal,
+ Format Table, File > "flux_raw.dat",
+ SendToServer Sprintf["}Flux/Wire %g/2Raw",i]];
+ Print[ flux[ LWIRE~{i} ], OnGlobal,
+ Format Table, File > "flux_reg.dat",
+ SendToServer Sprintf["}Flux/Wire %g/4Reg",i]];
+*/
+ /*
+ X = rw/skin_depth;
+ pI[] = (X/2)*Re[(1+i[])*(J_0[(1+i[])*X]/J_1[(1+i[])*X])];
+ Resist[] = Re[pI[]*R_dc];
+ // radial analytical solution with a zero flux condition imposed at r=rs
+ Analytic_a~{i}[] = mu0/(2*Pi) * WI~{i} *
+ ( (R~{i}[]>WR~{i}) ? Log[rs/R~{i}[]] :
+ Log[rs/WR~{i}] + (J0[tau[]*R~{i}[]/rw] - J0[tau[]]) / J1[tau[]] / tau[] );
+
+ Analytic_L~{i}[] = mu0 / (2*Pi) *
+ ( 2/tau[]/tau[]- J0[tau[]] / J1[tau[]] / tau[] - Log[rs/WR~{i}] ) ;
+
+ // Idem in static case (thus particular case of the above)
+ AnalyticStatic_a~{i}[] = mu0 * WI~{i} / (2*Pi) *
+ ( (R~{i}[]>WR~{i}) ? Log[rs/R~{i}[]] : Log[rs/WR~{i}] + (1-(R~{i}[]/WR~{i})^2)/2. );
+ Analytic_e~{i}[] = mu0 * WI~{i} / (2*Pi) * i[] * omega
+ * J0[tau[]*R~{i}[]/WR~{i}] / J1[tau[]] / tau[] ;
+ Analytic_flux~{i}[] = mu0 * WI~{i} / (2*Pi) *
+ ( i[] * (skin_depth/WR~{i})^2 + Log[rs/WR~{i}] - J0[tau[]] / J1[tau[]] / tau[] ) ;
+ */
diff --git a/Conductors3D/w2s_common.pro b/Conductors3D/w2s_common.pro
new file mode 100644
index 0000000000000000000000000000000000000000..69bd48127bb139c534405924453a72a9ef129171
--- /dev/null
+++ b/Conductors3D/w2s_common.pro
@@ -0,0 +1,47 @@
+
+mm = 1e-3;
+deg = 180/Pi;
+
+DefineConstant[
+ NumWires = {1, Name "Parameters/0Number of wires",
+ Choices {1="1", 2="2", 3="3"}, Visible 1}
+
+ Flag_Thin = {1, Name "Parameters/1Wire representation",
+ Choices {0="Real", 1="1D (raw sleeve)", 2="1D (regular sleeve)", 3="1D (naive)"}}
+ Flag_3D = {0, Name "Parameters/2Extruded model", Choices {0,1}, Visible 0}
+ Flag_Maps = {1, Name "Input/2Display maps", Choices {0,1}, Visible 1}
+ LogFreq = {4, Min 0, Max 8, Step 1, Name "Input/4Log of frequency"}
+ Freq = 10^LogFreq
+
+ WireRadius = {0.564189, Name "Parameters/5Wire radius [mm]"}
+ MeshSizeWire = {0.2, Name "Parameters/6Mesh size on wire [mm]"}
+ SleeveRadius = {2, Name "Parameters/7Sleeve radius [mm]"}
+ box = {100*mm, Name "Parameters/7box half-width [m]"}
+
+ rw = WireRadius*mm
+ rs = SleeveRadius*mm
+ A_c = Pi*rw^2 // wire cross section
+ rout = box
+ Lz = 1*mm
+
+ /* 'WireRadius' is the actual radius of the wire.
+ 'MeshSizeWire' is the imposed mesh size at the nodes of the wire,
+ and thus also the practical approximative value of the radius of the sleeve.
+ The 'rs' and 'rw' parameters are assumed identical for all wires.
+ */
+];
+
+
+
+Xlocs() = { 0, 5*mm, -5*mm};
+Ylocs() = { 0, 0, 0};
+For i In {1:NumWires}
+ DefineConstant[
+ WX~{i} = { Xlocs(i-1), // place wires symmetrically around x=0
+ Name Sprintf("Parameters/Wire %g/0X position [m]", i), Closed }
+ WY~{i} = { Ylocs(i-1), Name Sprintf("Parameters/Wire %g/0Y position [m]", i) }
+ WI~{i} = { 1, Name Sprintf("Parameters/Wire %g/2Current [A]", i) }
+ WP~{i} = { 0*NumWires*(i-1),
+ Name Sprintf("Parameters/Wire %g/3Phase [deg]", i) }
+ ];
+EndFor
diff --git a/Conductors3D/w3d.geo b/Conductors3D/w3d.geo
new file mode 100644
index 0000000000000000000000000000000000000000..d103979ff78346c5d89fb5c552ecd49531e9ce0a
--- /dev/null
+++ b/Conductors3D/w3d.geo
@@ -0,0 +1,149 @@
+Include "w3d_common.pro";
+
+// global Gmsh options
+Mesh.Optimize = 1; // optimize quality of tetrahedra
+Mesh.VolumeEdges = 0; // hide volume edges
+Geometry.ExactExtrusion = 0; // to allow rotation of extruded shapes
+Solver.AutoMesh = 2; // always remesh if necessary (don't reuse mesh on disk)
+
+lc1 = (Flag_Thin) ? MeshSizeThinWire*mm : MeshSizeThinWire*mm ;
+lc2 = box/10; // mesh size at outer surface
+
+llWires[] = {};
+centerWires[] = {};
+surfWires[] = {};
+
+If( !Flag_Thin || (Flag_Thin && Flag_Corr==2))
+ // then circles are in the geometry
+ // their radius is R=rw if Flag_Thin != 0
+ // or R=rs in the "regular sleeve case" (Flag_Corr == 2)
+
+ R = ( !Flag_Thin ) ? rw : rs;
+
+ For i In {1:NumWires}
+ pC = newp; Point(pC) = {WX~{i} , WY~{i} , 0, lc1};
+ p1 = newp; Point(p1) = {WX~{i}+R, WY~{i} , 0, lc1};
+ p2 = newp; Point(p2) = {WX~{i} , WY~{i}+R, 0, lc1};
+ p3 = newp; Point(p3) = {WX~{i}-R, WY~{i} , 0, lc1};
+ p4 = newp; Point(p4) = {WX~{i} , WY~{i}-R, 0, lc1};
+ c1 = newl; Circle(c1) = {p1,pC,p2};
+ c2 = newl; Circle(c2) = {p2,pC,p3};
+ c3 = newl; Circle(c3) = {p3,pC,p4};
+ c4 = newl; Circle(c4) = {p4,pC,p1};
+ l1 = newl; Line(l1) = {pC,p1};
+ l2 = newl; Line(l2) = {pC,p2};
+ l3 = newl; Line(l3) = {pC,p3};
+ l4 = newl; Line(l4) = {pC,p4};
+ ll1 = newll; Line Loop(ll1) = {l1,c1,-l2};
+ s1 = news; Plane Surface(s1) = {ll1};
+ ll2 = newll; Line Loop(ll2) = {l2,c2,-l3};
+ s2 = news; Plane Surface(s2) = {ll2};
+ ll3 = newll; Line Loop(ll3) = {l3,c3,-l4};
+ s3 = news; Plane Surface(s3) = {ll3};
+ ll4 = newll; Line Loop(ll4) = {l4,c4,-l1};
+ s4 = news; Plane Surface(s4) = {ll4};
+
+ ll5 = newll;
+ Line Loop(ll5) = {c1,c2,c3,c4};
+ llWires[] += { ll5 };
+ surfWires~{i-1}[] = {s1,s2,s3,s4};
+ surfWires[] += surfWires~{i-1}[];
+ centerWires[] += {pC};
+ EndFor
+Else
+ // then geometry only contains line conductors
+ For i In {1:NumWires}
+ pC = newp; Point(pC) = {WX~{i}, WY~{i}, 0, lc1};
+ centerWires[] += {pC};
+ EndFor
+
+EndIf
+
+// create air box around wires
+BoundingBox; // recompute model bounding box
+cx = (General.MinX + General.MaxX) / 2;
+cy = (General.MinY + General.MaxY) / 2;
+cz = (General.MinZ + General.MaxZ) / 2;
+
+// lx = 2*inf + General.MaxX - General.MinX;
+// ly = 2*inf + General.MaxY - General.MinZ;
+lx = 2*box;
+ly = 2*box;
+
+p1 = newp; Point (p1) = {cx-lx/2, cy-ly/2, 0, lc2};
+p2 = newp; Point (p2) = {cx+lx/2, cy-ly/2, 0, lc2};
+p3 = newp; Point (p3) = {cx+lx/2, cy+ly/2, 0, lc2};
+p4 = newp; Point (p4) = {cx-lx/2, cy+ly/2, 0, lc2};
+l1 = newl; Line(l1) = {p1, p2};
+l2 = newl; Line(l2) = {p2, p3};
+l3 = newl; Line(l3) = {p3, p4};
+l4 = newl; Line(l4) = {p4, p1};
+ll1 = newll; Line Loop(ll1) = {l1,l2,l3,l4};
+
+Skin[] = {};
+
+If( !Flag_Thin ) // Wires with finite radius
+
+ For i In {1:NumWires}
+ e[] = Extrude {0,0,Lz}
+ { Surface{ surfWires~{i-1}[] } ; /*Layers{1}; Recombine;*/ } ;
+ Physical Volume (Sprintf("VWIRE_%g",i), 10+i) = { e[ {1,6,11,16} ] };
+ Physical Surface(Sprintf("SANODE_%g",i), 20+i) = { surfWires~{i-1}[] };
+ Physical Surface(Sprintf("SCATHODE_%g",i), 30+i) = { e[ {0,5,10,15} ] };
+ Skin[] += e[ {3,8,13,18} ];
+ EndFor
+ s1 = news; Plane Surface(s1) = {ll1,-llWires[]};
+ e[] = Extrude {0,0,Lz} { Surface{ s1 } ; /*Layers{1}; Recombine;*/ } ;
+ Physical Volume ("AIR", 1) = { e[1] };
+
+ElseIf( Flag_Corr == 2 ) // regular sleeve
+
+ Wires[]= {};
+ For i In {1:NumWires}
+ e[] = Extrude {0,0,Lz} { Point{ centerWires[i-1] } ; /*Layers{1};*/ } ;
+ Physical Line (Sprintf("LWIRE_%g",i), 50+i) = { e[1] };
+ Physical Point (Sprintf("PANODE_%g",i), 60+i) = { centerWires[i-1] };
+ Physical Point (Sprintf("PCATHODE_%g",i), 70+i) = { e[0] };
+
+ e[] = Extrude {0,0,Lz} { Surface{ surfWires~{i-1}[] } ;
+ /*Layers{1}; Recombine;*/ } ;
+ Wires() += e[ {1,6,11,16} ];
+ Physical Volume (Sprintf("VWIRE_%g",i), 10+i) = { e[ {1,6,11,16} ] };
+ Physical Surface(Sprintf("SANODE_%g",i), 20+i) = { surfWires~{i-1}[] };
+ Physical Surface(Sprintf("SCATHODE_%g",i), 30+i) = { e[ {0,5,10,15} ] };
+ EndFor
+ s1 = news; Plane Surface(s1) = {ll1,-llWires[]};
+ e[] = Extrude {0,0,Lz} { Surface{ s1 } ; /*Layers{1}; Recombine;*/ } ;
+ Physical Volume ("AIR", 1) = { e[1], Wires[] };
+
+Else
+
+ s1 = news; Plane Surface(s1) = {ll1};
+ e[] = Extrude {0,0,Lz} { Surface{ s1 } ; /*Layers{1}; Recombine;*/ } ;
+ s2 = e[0];
+ v1 = e[1];
+ For i In {1:NumWires}
+ e[] = Extrude {0,0,Lz} { Point{ centerWires[i-1] } ; /*Layers{1};*/ } ;
+ //Printf("e=", e[]);
+ Physical Line (Sprintf("LWIRE_%g",i), 50+i) = { e[1] };
+ Physical Point (Sprintf("PANODE_%g",i), 60+i) = { centerWires[i-1] };
+ Physical Point (Sprintf("PCATHODE_%g",i), 70+i) = { e[0] };
+ Physical Volume ("AIR", 1) = { v1 };
+ // embed line conductors in mesh
+ Point { centerWires[i-1] } In Surface { s1 };
+ Point { e[0] } In Surface { s2 };
+ Curve { e[1] } In Volume { v1 };
+ EndFor
+
+EndIf
+
+Physical Surface("INF", 2) = { CombinedBoundary{ Volume{ : }; } };
+Physical Surface("SKIN", 3) = Skin[];
+
+Electrodes[] = {};
+If( !Flag_Thin ) // Wires with finite radius
+ For i In {1:NumWires}
+ Electrodes[] += {20+i, 30+i};
+ EndFor
+EndIf
+Physical Line("LINTREE", 4) = { Boundary { Physical Surface { Electrodes[] }; } };
diff --git a/Conductors3D/w3d.pro b/Conductors3D/w3d.pro
new file mode 100644
index 0000000000000000000000000000000000000000..d2fc8df9b4173794a84dd5b169ef7bdf2f190de9
--- /dev/null
+++ b/Conductors3D/w3d.pro
@@ -0,0 +1,1015 @@
+Include "w3d_common.pro";
+
+Struct S::SOLVER [ Enum, mumps, gmres, fgmres, bcgs ];
+
+DefineConstant[
+ FileId = ""
+ DataDir = StrCat["data/" , FileId]
+ ElekinDir = StrCat["Elekin_Line_Const/" , FileId]
+ MagDynDir = StrCat["MagDyn_Line_Const/" , FileId]
+ FullWaveDir = StrCat["FullWave_Line_Const/", FileId]
+ type_solver = {
+ S::SOLVER.mumps,
+ Choices {
+ S::SOLVER.mumps ="mumps",
+ S::SOLVER.gmres="gmres",
+ S::SOLVER.fgmres= "fgmres",
+ S::SOLVER.bcgs="bcgs"
+ },
+ Name "Solver/type_solver"
+ }
+];
+
+DefineConstant[
+ R_ = {"Dynamic", Name "GetDP/1ResolutionChoices", Visible 0},
+ C_ = {"-sol -pos", Name "GetDP/9ComputeCommand", Visible 0},
+ P_ = {"", Name "GetDP/2PostOperationChoices", Visible 0}
+ ResDir = "Res/"
+];
+
+
+Group{
+ // Geometrical regions
+ AIR = Region[ 1 ];
+ INF = Region[ 2 ];
+ SKIN = Region[ 3 ];
+ LINTREE = Region[ 4 ]; // not used
+
+ VWIRES = Region[ {} ];
+ LWIRES = Region[ {} ];
+ ANODES = Region[ {} ];
+ CATHODES = Region[ {} ];
+
+ For i In {1:NumWires}
+ LWIRE~{i} = Region[ (50+i) ];
+ LWIRES += Region[ LWIRE~{i} ];
+
+ If( !Flag_Thin )
+ VWIRE~{i} = Region[ (10+i) ];
+ VWIRES += Region[ VWIRE~{i} ];
+ SWIRE~{i} = Region[ {} ];
+ ANODE~{i} = Region[ (20+i) ];
+ CATHODE~{i} = Region[ (30+i) ];
+ Else
+ VWIRE~{i} = Region[ {} ];
+ ANODE~{i} = Region[ (60+i) ];
+ CATHODE~{i} = Region[ (70+i) ];
+ SWIRE~{i} = ElementsOf[ AIR, OnOneSideOf LWIRE~{i} ];
+ EndIf
+
+ ANODES += Region[ ANODE~{i} ];
+ CATHODES += Region[ CATHODE~{i} ];
+ EndFor
+
+
+ // Abstract regions
+ // Cette formulation peut contenir des conducteurs volumiques Vol_C
+ // et lin�iques lVol_C.
+
+ If( !Flag_Thin )
+ Vol_C = Region[ VWIRES ];
+ Vol_CC = Region[ AIR ];
+ //lVol_C = Region[ {} ];
+ Else
+ Vol_C = Region[ {LWIRES} ];
+ Vol_CC = Region[ {VWIRES, AIR} ];
+ //lVol_C = Region[ LWIRES ];
+ EndIf
+
+ Vol_nu = Region[ {Vol_C, Vol_CC} ];
+ Sur_Dirichlet_a = Region[ INF ];
+ Electrodes = Region[ {ANODES, CATHODES} ];
+
+ Vol_Tree = Region[ { Vol_nu } ];
+ Sur_Tree = Region[ { Sur_Dirichlet_a /*, SKIN*/ } ];
+ // If( Flag_Corr!= 0 )
+ // Lin_Tree = Region[ { LWIRES /*, LINTREE*/ } ];
+ // Else
+ Lin_Tree = Region[ {} ];
+
+
+ Dom_Hcurl_a = Region[ Vol_nu ];
+ Dom_Hgrad_u = Region[ Vol_C ];
+ //Dom_Hregion_i = Region[ Vol_nu ];
+
+ // additional Groups for the local correction method
+ //Dom_Hthin_a = ElementsOf[ Vol_nu, OnOneSideOf LWIRES ];
+
+ Dom_Hthin_a = Region[ Vol_nu ];
+
+ // For integration on the volume elements (excluding line elements)
+ // of the sleeve
+ Vol_Sleeve = ElementsOf[ {Vol_CC}, OnOneSideOf LWIRES ];
+}
+
+
+Function{
+ omega = 2*Pi*Freq;
+ mu0 = Pi*4e-7;
+ mur = 1.;
+ sigma = 5.96e7;
+
+ i[] = Complex[0,1];
+
+ mu[] = mu0;
+ nu[] = 1/mu[];
+ sigma[] = sigma;
+ skin_depth = Sqrt[ 2/(omega*sigma*mu0*mur) ];
+
+ tau[] = Complex[-1,1]/skin_depth*rw;
+ J0[] = JnComplex[0,$1];
+ J1[] = JnComplex[1,$1];
+ dJ1[] = JnComplex[0,$1]-JnComplex[1,$1]/$1;
+
+ For i In {1:NumWires}
+ // Current density
+ J[ Region[ {VWIRE~{i}, LWIRE~{i}} ] ] =
+ Vector[ 0, 0, WI~{i}/A_c * Exp[ i[]*WP~{i}*deg ] ];
+
+ // local radial coordinate with origin on wire i
+ R~{i}[] = Hypot[ X[]-WX~{i} , Y[]-WY~{i} ];
+ EndFor
+
+ // shape function of the current, wI[r]
+ wI[] = tau[] * J0[tau[]*$1/rw] / J1[tau[]] /(2*A_c);
+
+ //radial analytical solutions with a zero flux condition imposed at $1(=r)=rs
+ Analytic_A[] = // per Amp
+ (($1>rw) ? mu0/(2*Pi) * Log[rs/$1] :
+ mu0/(2*Pi) * Log[rs/rw] + i[]*( wI[$1]-wI[rw] )/(sigma*omega) );
+ AnalyticStatic_A[] = mu0/(2*Pi) * // per Amp
+ (($1>rw) ? Log[rs/$1] : Log[rs/rw] + mur*(1-($1/rw)^2)/2 );
+ Analytic_B[] = // per Amp
+ (($1>rw) ? mu0/(2*Pi*$1) :
+ mu0*mur/(2*Pi) * J1[tau[]*$1/rw] / J1[tau[]] );
+
+ // Impedance of thin wire p.u. length
+ R_DC = 1./(sigma*A_c);
+ Analytic_R[] = Re[ wI[rw]/sigma ];
+ Analytic_L[] = mu0/(2*Pi) * Log[rs/rw] + Re[ wI[rw]/(i[]*omega*sigma) ];
+
+ If( NumWires == 1 )
+ Exact_B[] = Analytic_B[R_1[]] ;
+ Correction_B[] = ((R_1[]<rs) ? Analytic_B[R_1[]] : 0 ) ;
+ EndIf
+ If ( NumWires == 2 )
+ Exact_B[] = Analytic_B[R_1[]] + Analytic_B[R_2[]] ;
+ Correction_B[] = ((R_1[]<rs) ? Analytic_B[R_1[]] :
+ ((R_2[]<rs) ? Analytic_B[R_2[]] : 0 ));
+ EndIf
+ If ( NumWires == 3 )
+ Exact_B[] = Analytic_B[R_1[]] + Analytic_B[R_2[]] + Analytic_B[R_3[]] ;
+ Correction_A[] = ((R_1[]<rs) ? Analytic_B[R_1[]] :
+ ((R_2[]<rs) ? Analytic_B[R_2[]] :
+ ((R_3[]<rs) ? Analytic_B[R_3[]] : 0 )));
+ EndIf
+}
+
+
+Jacobian {
+ { Name Vol ;
+ Case {
+ { Region LWIRES ; Jacobian Lin ; Coefficient A_c; }
+ { Region All ; Jacobian Vol ; }
+ }
+ }
+}
+
+Integration {
+ { Name I1 ;
+ Case {
+ { Type Gauss ;
+ Case {
+ { GeoElement Point ; NumberOfPoints 1 ; }
+ { GeoElement Line ; NumberOfPoints 3 ; }
+ { GeoElement Triangle ; NumberOfPoints 3 ; }
+ { GeoElement Quadrangle ; NumberOfPoints 4 ; }
+ { GeoElement Prism ; NumberOfPoints 21 ; }
+ { GeoElement Tetrahedron ; NumberOfPoints 4 ; }
+ { GeoElement Hexahedron ; NumberOfPoints 6 ; }
+ { GeoElement Pyramid ; NumberOfPoints 8 ; }
+ }
+ }
+ }
+ }
+}
+
+Constraint{
+ { Name Impose_U ;
+ Case {
+ For i In {1:NumWires}
+ If( Flag_Form == 1) // massive
+ { Region ANODE~{i} ; Value 0 ; }
+ EndIf
+ If( Flag_U )
+
+ If( Flag_Form == 1) // massive
+ { Region CATHODE~{i} ; Value -R_DC*Lz*WI~{i} ; }
+ Else // stranded
+ { Region VWIRE~{i} ; Value -R_DC*Lz*WI~{i} ; }
+ EndIf
+
+ EndIf
+ EndFor
+ }
+ }
+ { Name Impose_I ;
+ Case {
+ For i In {1:NumWires}
+ If( !Flag_U )
+
+ If( Flag_Form == 1) // massive
+ { Region ANODE~{i} ; Value WI~{i} ; }
+ { Region CATHODE~{i} ; Value WI~{i} ; }
+ Else // stranded
+ { Region VWIRE~{i} ; Value WI~{i} ; }
+ { Region LWIRE~{i} ; Value WI~{i} ; }
+ EndIf
+
+ EndIf
+ EndFor
+ }
+ }
+
+ { Name Hcurl_a_3D_ac; Type Assign;
+ Case {
+ { Region Region[ Sur_Dirichlet_a ]; Value 0.; }
+ }
+ }
+
+ { Name GaugeCondition_a; Type Assign;
+ Case {
+ { Region Vol_Tree; SubRegion Region[ { Sur_Tree, Lin_Tree } ];
+ Value 0.; }
+ }
+ }
+}
+
+
+
+FunctionSpace {
+ // Magnetic vector potential
+ { Name Hcurl_a_3D; Type Form1;
+ BasisFunction {
+ { Name sc; NameOfCoef ac; Function BF_Edge;
+ Support Dom_Hcurl_a; Entity EdgesOf[All]; }
+ }
+ Constraint {
+ { NameOfCoef ac;
+ EntityType Auto; NameOfConstraint Hcurl_a_3D_ac; }
+ { NameOfCoef ac; EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
+ NameOfConstraint GaugeCondition_a ; }
+ }
+ }
+
+ // Electric scalar potential for massive conductors
+ { Name Hgrad_u_3D; Type Form0;
+ BasisFunction {
+ { Name sn; NameOfCoef un; Function BF_Node;
+ Support Dom_Hgrad_u; Entity NodesOf[ Dom_Hgrad_u, Not Electrodes]; }
+ { Name sn2; NameOfCoef un2; Function BF_GroupOfNodes;
+ Support Dom_Hgrad_u; Entity GroupsOfNodesOf[ Electrodes ]; }
+ }
+ GlobalQuantity {
+ { Name U; Type AliasOf ; NameOfCoef un2; }
+ { Name I; Type AssociatedWith; NameOfCoef un2; }
+ }
+ Constraint {
+ { NameOfCoef U; EntityType GroupsOfNodesOf; NameOfConstraint Impose_U; }
+ { NameOfCoef I; EntityType GroupsOfNodesOf; NameOfConstraint Impose_I; }
+ }
+ }
+
+ { Name Hregion_i_3D; Type Vector;
+ BasisFunction {
+ { Name sr; NameOfCoef ir; Function BF_RegionZ;
+ Support Dom_Hthin_a; Entity Vol_C; }
+ }
+ GlobalQuantity {
+ { Name Is; Type AliasOf ; NameOfCoef ir; }
+ { Name Us; Type AssociatedWith; NameOfCoef ir; }
+ }
+ Constraint {
+ { NameOfCoef Us; EntityType Region; NameOfConstraint Impose_U; }
+ { NameOfCoef Is; EntityType Region; NameOfConstraint Impose_I; }
+ }
+ }
+
+ // For the local correction method
+ { Name Hcurl_acorr_3D; Type Form1;
+ BasisFunction {
+ { Name sw; NameOfCoef aw; Function BF_Edge;
+ Support Vol_nu; Entity EdgesOf[ Vol_nu, ConnectedTo LWIRES ]; }
+ }
+ Constraint {
+ // { NameOfCoef aw; EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
+ // NameOfConstraint GaugeCondition_a ; }
+ }
+ }
+}
+
+Formulation {
+ { Name MagnetoDynamics; Type FemEquation;
+
+ If( Flag_Form == 1) // Massive conductor
+ Quantity {
+ { Name a; Type Local; NameOfSpace Hcurl_a_3D; }
+
+ { Name v; Type Local; NameOfSpace Hgrad_u_3D; }
+ { Name U; Type Global; NameOfSpace Hgrad_u_3D [U]; }
+ { Name I; Type Global; NameOfSpace Hgrad_u_3D [I]; }
+ }
+
+ Equation {
+ Integral { [ nu[] * Dof{d a} , {d a} ];
+ In Vol_nu; Jacobian Vol; Integration I1; }
+
+ Integral { DtDof[ sigma[] * Dof{a} , {a} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ Integral { [ sigma[] * Dof{d v} , {a} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ Integral { DtDof[ sigma[] * Dof{a} , {d v} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ Integral { [ sigma[] * Dof{d v} , {d v} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+
+ GlobalTerm { [Dof{I} , {U} ]; In CATHODES; }
+ }
+
+ Else // stranded conductor
+ If( Flag_Corr==0 ) // naive thin wire model
+ Quantity {
+ { Name a; Type Local; NameOfSpace Hcurl_a_3D; }
+
+ { Name i; Type Local; NameOfSpace Hregion_i_3D; }
+ { Name Is; Type Global; NameOfSpace Hregion_i_3D [Is]; }
+ { Name Us; Type Global; NameOfSpace Hregion_i_3D [Us]; }
+ }
+
+ Equation {
+ Integral { [ nu[] * Dof{d a} , {d a} ];
+ In Vol_nu; Jacobian Vol; Integration I1; }
+ Integral { [ -Dof{i}/A_c , {a} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+
+ // Integral { [ -J[] , {a} ];
+ // In Vol_C; Jacobian Vol; Integration I1; }
+
+ // Us is the voltage drop
+ // grad v = -Dt a - j/sigma
+ Integral { DtDof [ Dof{a} , {i} ];
+ In Vol_C; Jacobian Vol; Integration I1;}
+ Integral { [ Dof{i}/(sigma[]*A_c) , {i} ];
+ In Vol_C; Jacobian Vol; Integration I1;}
+ GlobalTerm { [ Dof{Us}*A_c , {Is} ]; In Vol_C; }
+ }
+ Else // use local correction method
+ Quantity {
+ { Name a; Type Local; NameOfSpace Hcurl_a_3D; }
+ { Name as; Type Local; NameOfSpace Hcurl_acorr_3D; }
+
+ { Name i; Type Local; NameOfSpace Hregion_i_3D; }
+ { Name Is; Type Global; NameOfSpace Hregion_i_3D [Is]; }
+ { Name Us; Type Global; NameOfSpace Hregion_i_3D [Us]; }
+ }
+
+ Equation {
+ /*
+ Integral { [ nu[] * Dof{d a} , {d a} ];
+ In Vol_nu; Jacobian Vol; Integration I1; }
+ Integral { [ -Dof{i}/A_c , {a} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ */
+ Integral { [ nu[] * Dof{d as} , {d as} ];
+ In Vol_nu; Jacobian Vol; Integration I1; }
+ Integral { [ -Dof{i}/A_c , {as} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+
+ /*
+ GlobalTerm { [ Analytic_R[] * Dof{Is} , {Is} ];
+ In Vol_C;}
+ Integral { DtDof [ Dof{a}/A_c , {i} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ Integral { DtDof [ -Dof{as}/A_c , {i} ];
+ In Vol_C; Jacobian Vol; Integration I1;}
+ GlobalTerm { DtDof [ Analytic_L[] * Dof{Is} , {Is} ];
+ In Vol_C;}
+ GlobalTerm { [ Dof{Us} , {Is} ];
+ In Vol_C; }
+ */
+ }
+ EndIf
+ EndIf
+ }
+}
+
+
+Resolution {
+ { Name Dynamic;
+ System {
+ { Name S; NameOfFormulation MagnetoDynamics;
+ Type ComplexValue; Frequency Freq; }
+ }
+ Operation {
+ CreateDir[DataDir];
+
+ Generate[S]; Solve[S]; SaveSolution[S];
+
+ If( Flag_Maps == 1 )
+ // DeleteFile["U.dat"];
+ // DeleteFile["I.dat"];
+ // DeleteFile["Z.dat"];
+
+ PostOperation[map];
+ EndIf
+ PostOperation[integaz];
+ PostOperation[cut];
+
+ For i In {1:NumWires}
+ Print[ {i, rw, rs, A_c, skin_depth, R_DC, WI~{i}/A_c, mu0*WI~{i}/(2*Pi*rw)},
+ Format "WIRE %2g: rw = %7.3e rs = %7.3e A_c = %7.3e delta = %7.3e R_DC = %7.3e J = %7.3e Bmax = %7.3e"] ;
+ // Print[ {i, rw, Sqrt[ $sarea~{i}/Pi], rs, $sarea~{i}, Pi*rs^2, skin_depth, $intby~{i}/$sarea~{i}},
+ // Format "WIRE %2g: rw = %7.3e rs = %7.3e (%7.3e) as = %7.3e (%7.3e) delta = %7.3e bave = %7.3e" ] ;
+ EndFor
+ }
+ }
+}
+
+PostProcessing {
+ { Name MagnetoDynamics; NameOfFormulation MagnetoDynamics;
+ PostQuantity {
+ { Name a;
+ Value { Local { [ {a}]; In Dom_Hcurl_a; Jacobian Vol; }}}
+ { Name b;
+ Value { Local { [ {d a}]; In Dom_Hcurl_a; Jacobian Vol; }}}
+ { Name by;
+ Value { Local { [ Re[Cart2Pol[CompY [ {d a}]]]];
+ In Dom_Hcurl_a; Jacobian Vol; }}}
+
+ If( Flag_Form == 1 ) // massive
+ { Name J;
+ Value { Local { [ -sigma[] * ( Dt[{a}] + {d v} ) ];
+ In Vol_C; Jacobian Vol; }}}
+ { Name U;
+ Value { Term { [ {U} ]; In Electrodes; }}}
+ { Name I;
+ Value { Term { [ {I} ]; In Electrodes; }}}
+ { Name R;
+ Value { Term { [ Re[ -{U}/{I}/Lz ] ]; In Electrodes; }}}
+ { Name L; // actually total flux/I
+ Value { Term { [ Im[ -{U}/{I}/omega/Lz ] ]; In Electrodes; }}}
+ Else // stranded
+ { Name J;
+ Value { Local { [ {i}/A_c ];
+ In Vol_C; Jacobian Vol; }}}
+ { Name U;
+ Value { Term { [ {Us} ]; In Vol_C; }}}
+ { Name I;
+ Value { Term { [ {Is} ]; In Vol_C; }}}
+ { Name R;
+ Value { Term { [ Re[ -{Us}/{Is}/Lz ] ]; In Vol_C; }}}
+ { Name L; // actually total flux/I
+ Value { Term { [ Im[ -{Us}/{Is}/omega/Lz ] ]; In Vol_C; }}}
+ EndIf
+
+ If( Flag_Corr )
+ { Name bs;
+ Value { Local { [ {d as}]; In Dom_Hthin_a; Jacobian Vol; }}}
+ { Name bcorry;
+ Value { Local { [ Re[Cart2Pol[CompY [{d a}-{d as}]]] + Correction_B[] ];
+ In Dom_Hthin_a; Jacobian Vol; }}}
+ EndIf
+
+ }
+ }
+}
+
+PostOperation map UsingPost MagnetoDynamics {
+ Print[ J, OnElementsOf Vol_C, File "j.pos"];
+ Print[ b, OnElementsOf Vol_nu, File "b.pos"];
+ If( Flag_Corr )
+ Print[ bs, OnElementsOf Vol_nu, File "bs.pos"];
+ PrintGroup[ EdgesOf[ Vol_nu, ConnectedTo LWIRES ],
+ In Vol_nu, File "group.pos"];
+ EndIf
+}
+
+PostOperation integaz UsingPost MagnetoDynamics {
+ For i In {1:NumWires}
+ If( Flag_Form == 1 ) // massive
+
+ Print[ U, OnRegion CATHODES, File > "U.dat", Format Table,
+ SendToServer Sprintf["Results/1Voltage/Wire %g",i]];
+ Print[ I, OnRegion CATHODES, File > "I.dat", Format Table,
+ SendToServer Sprintf["Results/2Current/Wire %g",i]];
+ Print[ R, OnRegion CATHODES, File > "Z.dat", Format Table,
+ SendToServer Sprintf["Results/3Resistance p.u.l./Wire %g",i]];
+ Print[ L, OnRegion CATHODES, File > "Z.dat", Format Table,
+ SendToServer Sprintf["Results/4Inductance p.u.l./Wire %g",i]];
+ Else
+ Print[ U, OnRegion Vol_C, File > "U.dat", Format Table,
+ SendToServer Sprintf["Results/1Voltage/Wire %g",i]];
+ Print[ I, OnRegion Vol_C, File > "I.dat", Format Table,
+ SendToServer Sprintf["Results/2Current/Wire %g",i]];
+ Print[ R, OnRegion Vol_C, File > "Z.dat", Format Table,
+ SendToServer Sprintf["Results/3Resistance p.u.l./Wire %g",i]];
+ Print[ L, OnRegion Vol_C, File > "Z.dat", Format Table,
+ SendToServer Sprintf["Results/4Inductance p.u.l./Wire %g",i]];
+ EndIf
+ EndFor
+}
+
+NbPoints = 400;
+Xmax = 0.15*box;
+
+PostOperation cut UsingPost MagnetoDynamics {
+ Print [ by, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
+ Format Gmsh, File "by.pos" ];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].Name = 'cut by';",
+ "View[l].Axes = 3;",
+ "View[l].LineWidth = 3;",
+ "View[l].Type = 2;"],
+ File "tmp.geo", LastTimeStepOnly];
+ If( Flag_Corr )
+ Print [ bcorry, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
+ Format Gmsh, File "by.pos" ];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].Name = 'cut by';",
+ "View[l].Axes = 3;",
+ "View[l].LineWidth = 3;",
+ "View[l].Type = 2;"],
+ File "tmp.geo", LastTimeStepOnly];
+ EndIf
+}
+
+/*
+
+PostProcessing {
+ { Name MagnetoDynamics_corr; NameOfFormulation MagnetoDynamics;
+ PostQuantity {
+ { Name a;
+ Value { Local { [ {a}]; In Dom_Hcurl_a; Jacobian Vol; }}}
+ { Name b;
+ Value { Local { [ {d a}]; In Dom_Hcurl_a; Jacobian Vol; }}}
+ { Name by;
+ Value { Local { [ Cart2Pol[CompY [ {d a}]]]; In Dom_Hcurl_a; Jacobian Vol; }}}
+
+ If( Flag_Form <= 1 )
+ { Name U;
+ Value { Term { [ {U} ]; In CATHODES; }}}
+ { Name flux;
+ Value { Term { [ Im[ {U}/omega ] ]; In CATHODES; }}}
+
+ Else
+
+ { Name bc;
+ Value { Local { [ {d ac}]; In Dom_Hcurl_a; Jacobian Vol; }}}
+ { Name bw;
+ Value { Local { [ {d aw}]; In Dom_Hthin_a; Jacobian Vol; }}}
+ { Name byw;
+ Value { Local { [ Cart2Pol[CompY [ {d aw}]]]; In Dom_Hthin_a; Jacobian Vol; }}}
+
+ { Name U;
+ Value { Term { [ {Us} ]; In CATHODES; }}}
+ { Name flux;
+ Value { Term { [ Im[ {Us}/omega ] ]; In CATHODES; }}}
+
+ EndIf
+
+ }
+ }
+}
+
+PostOperation map_corr UsingPost MagnetoDynamics {
+
+ Print[ b, OnElementsOf Vol_nu, File "b.pos"];
+ Print[ U, OnRegion CATHODES, File > "U.dat", Format Table,
+ SendToServer Sprintf["}Voltage/Wire %g",i]];
+ Print[ flux, OnRegion CATHODES, File > "L.dat", Format Table,
+ SendToServer Sprintf["}Inductance/Wire %g",i]];
+
+ If( Flag_Form >= 2 )
+ PrintGroup[ EdgesOf[ Vol_nu, ConnectedTo lVol_C ],
+ In Vol_nu, File "group.pos"];
+ Print [bw, OnElementsOf Vol_nu, File "b.pos"];
+ EndIf
+
+}
+
+
+
+FunctionSpace {
+ { Name Hcurl_ai_3D; Type Form1;
+ BasisFunction {
+ { Name sc; NameOfCoef ac; Function BF_Edge;
+ Support Dom_Hcurl_a; Entity EdgesOf[All]; }
+ { Name sw; NameOfCoef aw; Function BF_Edge;
+ Support Dom_Hthin_a; Entity EdgesOf[ Vol_nu, ConnectedTo lVol_C ]; }
+ }
+ SubSpace {
+ { Name ac ; NameOfBasisFunction { sc } ; }
+ { Name aw ; NameOfBasisFunction { sw } ; }
+ }
+ Constraint {
+ { NameOfCoef ac;
+ EntityType Auto; NameOfConstraint Hcurl_a_3D_ac; }
+ { NameOfCoef ac; EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
+ NameOfConstraint GaugeCondition_a ; }
+ { NameOfCoef aw; EntityType EdgesOfTreeIn ; EntitySubType StartingOn ;
+ NameOfConstraint GaugeCondition_a ; }
+ }
+ }
+}
+
+Formulation {
+ { Name Corrected; Type FemEquation;
+ If( Flag_Form == 1) // A-v formulation
+ // Conventional A-v formulation
+ // Works with FE discretized true wires of finite radius (Flag_Thin==0)
+ // and idealized thin wires (Flag_Thin!=0)
+ Quantity {
+ { Name a; Type Local; NameOfSpace Hcurl_a_3D; }
+
+ { Name v; Type Local; NameOfSpace Hgrad_u_3D; }
+ { Name U; Type Global; NameOfSpace Hgrad_u_3D [U]; }
+ { Name I; Type Global; NameOfSpace Hgrad_u_3D [I]; }
+ }
+
+ Equation {
+ Integral { [ nu[] * Dof{d a} , {d a} ];
+ In Vol_nu; Jacobian Vol; Integration I1; }
+
+ Integral { DtDof[ sigma[] * Dof{a} , {a} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ Integral { [ sigma[] * Dof{d v} , {a} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ Integral { DtDof[ sigma[] * Dof{a} , {d v} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ Integral { [ sigma[] * Dof{d v} , {d v} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ //GlobalTerm { [Dof{I} , {U} ]; In CATHODES; }
+
+ Integral { DtDof[ sigma[] * Dof{a} , {a} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+ Integral { [ sigma[] * Dof{d v} , {a} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+ Integral { DtDof[ sigma[] * Dof{a} , {d v} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+ Integral { [ sigma[] * Dof{d v} , {d v} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+
+ GlobalTerm { [Dof{I} , {U} ]; In CATHODES; }
+ }
+ Else // A-I formulation
+ Quantity {
+ { Name a; Type Local; NameOfSpace Hcurl_ai_3D; }
+ { Name ac; Type Local; NameOfSpace Hcurl_ai_3D[ac]; }
+ { Name aw; Type Local; NameOfSpace Hcurl_ai_3D[aw]; }
+
+ { Name i; Type Local; NameOfSpace Hregion_i_3D; }
+ { Name Is; Type Global; NameOfSpace Hregion_i_3D [Is]; }
+ { Name Us; Type Global; NameOfSpace Hregion_i_3D [Us]; }
+ }
+
+ Equation {
+ Integral { [ nu[] * Dof{d ac} , {d ac} ];
+ In Vol_nu; Jacobian Vol; Integration I1; }
+ Integral { [ nu[] * Dof{d aw} , {d aw} ];
+ In Vol_nu; Jacobian Vol; Integration I1; }
+
+ If( 0 ) // impose current density
+ Integral { [ -J[], {ac} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+ Integral { [ -J[] , {aw} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+ Integral { [ -J[] , {a} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+
+ Else
+
+ Integral { [ -Dof{i}/A_c , {ac} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+ Integral { [ -Dof{i}/A_c , {aw} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+
+ GlobalTerm { [ Dof{Us} , {Is} ]; In lVol_C; }
+ GlobalTerm { [ Analytic_R[] * Dof{Is} , {Is} ]; In lVol_C;}
+ Integral { DtDof [ Dof{ac}/A_c , {i} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+ If( Flag_Form != 2)
+ Integral { DtDof [ -Dof{aw}/A_c , {i} ];
+ In lVol_C; Jacobian Vol; Integration I1;}
+ GlobalTerm { DtDof [ Analytic_L[] * Dof{Is} , {Is} ];
+ In lVol_C;}
+ EndIf
+ EndIf
+ }
+ EndIf
+ }
+}
+
+*/
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+// ################################
+
+/*
+
+NbPoints = 1000;
+Xmax = 0.15*box;
+PostOperation cut UsingPost PostProcessings {
+ // cut of by displayed in GUI
+ If( Flag_Thin == 0 )
+ Print [ by, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
+ Format Gmsh, File "Cut_az.pos" ];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].Name = 'cut |by|';",
+ "View[l].Axes = 3;",
+ "View[l].LineWidth = 3;",
+ "View[l].Type = 2;"],
+ File "tmp.geo", LastTimeStepOnly];
+
+ Else
+
+ Print [ acz, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
+ Format Gmsh, File "Cut_az.pos" ];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].Name = 'cut az not corrected';",
+ "View[l].Axes = 3;",
+ "View[l].LineWidth = 3;",
+ "View[l].Type = 2;"],
+ File "tmp.geo", LastTimeStepOnly];
+
+ Print [ az, OnLine { {-Xmax,0,0} {Xmax,0,0} }{NbPoints},
+ Format Gmsh, File "Cut_az.pos" ];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].Name = 'cut az corrected';",
+ "View[l].Axes = 3;",
+ "View[l].LineWidth = 3;",
+ "View[l].Type = 2;"],
+ File "tmp.geo", LastTimeStepOnly];
+ EndIf
+
+ // cuts in txt format for Gnuplot
+ Print [ exact, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
+ Format SimpleTable, File "Cut_analytic.txt" ];
+
+ If( Flag_Thin == 0 )
+ Print [ az, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
+ Format SimpleTable, File "Cut_fullfem.txt" ];
+ Else
+ Print [ az, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
+ Format SimpleTable, File "Cut_thinfem.txt" ];
+ Print [ acz, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
+ Format SimpleTable, File "Cut_ac.txt" ];
+ Print [ awz, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
+ Format SimpleTable, File "Cut_aw.txt" ];
+ Print [ acwz, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
+ Format SimpleTable, File "Cut_acw.txt" ];
+ Print [ corr, OnLine { {0,0,0} {Xmax,0,0} } {NbPoints},
+ Format SimpleTable, File "Cut_correction.txt" ];
+ EndIf
+}
+
+PostProcessing {
+ { Name PostProcessings; NameOfFormulation MagnetoDynamics;
+ PostQuantity {
+ If( Flag_3D)
+ { Name a;
+ Value { Local { [ {a}];
+ In Dom_Hcurl_a; Jacobian Vol; }}}
+ // { Name ac;
+ // Value { Local { [ {ac} ];
+ // In Dom_Hcurl_a; Jacobian Vol; }}}
+ // { Name aw;
+ // Value { Local { [ {aw} ];
+ // In Dom_Hthin_a; Jacobian Vol; }}}
+ { Name b;
+ Value { Local { [ {d a}];
+ In Dom_Hcurl_a; Jacobian Vol; }}}
+ { Name U;
+ Value { Term { [ {U} ]; In CATHODES; }}}
+ { Name L;
+ Value { Term { [ Im[ {U}/{I}/omega ] ]; In CATHODES; }}}
+
+
+ Else ## 2D
+ { Name az;
+ Value { Local { [ CompZ[{ac}] - (Flag_Thin!=0)*(CompZ[{aw}]-Correction_A[])];
+ In Dom_Hcurl_a; Jacobian Vol; }}}
+ { Name acz;
+ Value { Local { [ CompZ[ {ac}] ];
+ In Dom_Hcurl_a; Jacobian Vol; }}}
+ { Name awz;
+ Value { Local { [ CompZ[ {aw}] ];
+ In Dom_Hthin_a; Jacobian Vol; }}}
+ { Name acwz;
+ Value { Local { [ CompZ[ {ac} - {aw} ] ];
+ In Dom_Hcurl_a; Jacobian Vol; }}}
+ { Name exact;
+ Value { Local { [ Exact_A[] ] ;
+ In Dom_Hcurl_a; Jacobian Vol; }}}
+ { Name corr;
+ Value { Local { [ Correction_A[] ] ;
+ In Dom_Hcurl_a; Jacobian Vol; }}}
+ { Name bcw;
+ Value { Local { [ {d ac} - {d aw} ];
+ In Dom_Hcurl_a; Jacobian Vol; }}}
+
+ { Name j;
+ Value { Local { [ -sigma[] * ( Dt[{a}] + {dv} ) ];
+ In Vol_C; Jacobian Vol; }}}
+
+ { Name un;
+ Value{ Local { [ 1 ] ;
+ In Vol_nu; Integration I1 ; Jacobian Vol; }}}
+ { Name surf;
+ Value{ Integral{ [ 1 ] ;
+ In Vol_nu; Integration I1 ; Jacobian Vol; }}}
+
+
+ If( Flag_Thin == 0 )
+ { Name flux;
+ Value { Integral { [ CompZ[ {a} ] / A_c ];
+ In Vol_C; Integration I1 ; Jacobian Vol; }}}
+ { Name JouleLosses;
+ Value { Integral { [ sigma[]*SquNorm[ Dt[{a}] + {dv} ] ] ;
+ In Vol_C ; Integration I1 ; Jacobian Vol ; }}}
+ // { Name U;
+ // Value { Term { [ Cart2Pol [ {U} ] ]; In Vol_C; }}}
+
+ { Name U;
+ Value { Term { [ Im [ {U}/omega ] ]; In Vol_C; }}}
+
+ Else
+ { Name flux; // naive flux
+ Value { Integral { [ CompZ[ {ac} ] ];
+ In lVol_C; Integration I1 ; Jacobian Vol; }}}
+
+ { Name JouleLosses;
+ Value { Term { [ Analytic_R[]*{Is}*{Is} ] ; In lVol_C; }}}
+
+ // { Name U;
+ // Value { Term { [ Cart2Pol [ {Us} ] ] ; In lVol_C; }}}
+
+ { Name U;
+ Value { Term { [ Im [ {Us}/omega ] ] ; In lVol_C; }}}
+
+ For i In {1:NumWires}
+ { Name area~{i};
+ Value{ Integral{ [ 1 ] ;
+ In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
+ { Name intbx~{i}; // integral of b, 1st step for computing b average
+ Value { Integral { [ Norm [ CompX[ {d ac} - {d aw} ] ] ] ;
+ In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
+ { Name intby~{i}; // integral of b, 1st step for computing b average
+ Value { Integral { [ Norm [ CompY[ {d ac} - {d aw} ] ] ] ;
+ In SWIRE~{i}; Integration I1 ; Jacobian Vol; }}}
+ EndFor
+
+ EndIf
+ EndIf
+ }
+ }
+}
+
+PostOperation map UsingPost PostProcessings {
+ Print [az, OnElementsOf Vol_nu, File "az.pos"];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].IntervalsType = 3;",
+ "View[l].NbIso = 30;",
+ "View[l].RaiseZ = 15000;",
+ "View[l].NormalRaise = 0;"],
+ File "tmp.geo", LastTimeStepOnly];
+
+ If(Flag_Thin != 0)
+ Print [bcw, OnElementsOf Vol_Sleeve, File "b.pos"];
+
+ Print [acwz, OnElementsOf Vol_nu, File "acz.pos"];
+ Echo[ StrCat["l=PostProcessing.NbViews-1;",
+ "View[l].IntervalsType = 3;",
+ "View[l].NbIso = 30;",
+ "View[l].RaiseZ = 15000;",
+ "View[l].NormalRaise = 0;"],
+ File "tmp.geo", LastTimeStepOnly] ;
+ EndIf
+}
+
+PostOperation integaz UsingPost PostProcessings {
+ For i In {1:NumWires}
+
+ If( Flag_Thin == 0 )
+ Print[ U, OnRegion Vol_C, File > "U_full.dat", Format Table,
+ SendToServer Sprintf["}Voltage/Wire %g/1Full",i]];
+ Print[ JouleLosses[ VWIRE~{i}], OnGlobal,
+ Format Table, File > "joule_full.dat",
+ SendToServer Sprintf["}Joule Losses/Wire %g/1Full",i]];
+ EndIf
+ If( Flag_Thin == 1 )
+ Print[ U, OnRegion lVol_C, File > "U_raw.dat", Format Table,
+ SendToServer Sprintf["}Voltage/Wire %g/2Raw",i]];
+ Print[ JouleLosses, OnRegion LWIRE~{i},
+ Format Table, File > "joule_raw.dat",
+ SendToServer Sprintf["}Joule Losses/Wire %g/2Raw",i]];
+ EndIf
+ If( Flag_Thin == 2 )
+ Print[ U, OnRegion lVol_C, File > "U_reg.dat", Format Table,
+ SendToServer Sprintf["}Voltage/Wire %g/3Reg",i]];
+ Print[ JouleLosses, OnRegion LWIRE~{i},
+ Format Table, File > "joule_reg.dat",
+ SendToServer Sprintf["}Joule Losses/Wire %g/3Reg",i]];
+ EndIf
+ If( Flag_Thin == 3 )
+ Print[ U, OnRegion lVol_C, File > "U_naive.dat", Format Table,
+ SendToServer Sprintf["}Voltage/Wire %g/4Naive",i]];
+ Print[ JouleLosses, OnRegion LWIRE~{i},
+ Format Table, File > "joule_naive.dat",
+ SendToServer Sprintf["}Joule Losses/Wire %g/4Naiv",i]];
+ EndIf
+
+ If( Flag_Thin != 0 )
+ Print[ area~{i}[SWIRE~{i}], OnGlobal,
+ Format Table, File >"zzz",
+ StoreInVariable $sarea~{i}];
+ Print[ intby~{i}[SWIRE~{i}], OnGlobal,
+ Format Table, File > "zzz",
+ StoreInVariable $intby~{i}];
+ EndIf
+
+ EndFor
+}
+
+
+*/
+
+ /* 2D case
+ Quantity {
+ { Name a; Type Local; NameOfSpace Hcurl_a_2D; }
+ { Name ac; Type Local; NameOfSpace Hcurl_a_2D[ac]; }
+ { Name aw; Type Local; NameOfSpace Hcurl_a_2D[aw]; }
+
+ { Name dv; Type Local; NameOfSpace Hgrad_u_2D; }
+ { Name U; Type Global; NameOfSpace Hgrad_u_2D [U]; }
+ { Name I; Type Global; NameOfSpace Hgrad_u_2D [I]; }
+
+ { Name i; Type Local; NameOfSpace Hregion_i_2D; }
+ { Name Is; Type Global; NameOfSpace Hregion_i_2D [Is]; }
+ { Name Us; Type Global; NameOfSpace Hregion_i_2D [Us]; }
+ }
+
+ Equation {
+
+ Integral { [ nu[] * Dof{d ac} , {d ac} ];
+ In Vol_nu; Jacobian Vol; Integration I1; }
+ Integral { [ -Dof{i}/A_c , {ac} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+
+ Integral { [ nu[] * Dof{d aw} , {d aw} ];
+ In Vol_nu; Jacobian Vol; Integration I1; }
+ Integral { [-Dof{i}/A_c , {aw} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+
+
+ Integral { [ -J[] , {a} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+
+ Integral { DtDof[ sigma[] * Dof{a} , {a} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ Integral { [ sigma[] * Dof{dv} , {a} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+
+ Integral { DtDof[ sigma[] * Dof{a} , {dv} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ Integral { [ sigma[] * Dof{dv} , {dv} ];
+ In Vol_C; Jacobian Vol; Integration I1; }
+ GlobalTerm { [ Dof{I} , {U} ]; In Vol_C; }
+
+ GlobalTerm { [ Dof{Us} , {Is} ];
+ In lVol_C; }
+ GlobalTerm { [ Analytic_R[] * Dof{Is} , {Is} ];
+ In lVol_C;}
+ Integral { DtDof [ Dof{ac}/A_c , {i} ];
+ In lVol_C; Jacobian Vol; Integration I1; }
+ If( Flag_Thin != 3)
+ Integral { DtDof [ -Dof{aw}/A_c , {i} ];
+ In lVol_C; Jacobian Vol; Integration I1;}
+ GlobalTerm { DtDof [ Analytic_L[] * Dof{Is} , {Is} ];
+ In lVol_C;}
+ EndIf
+ }
+ */
+
diff --git a/Conductors3D/w3d_common.pro b/Conductors3D/w3d_common.pro
new file mode 100644
index 0000000000000000000000000000000000000000..e8dc467f0dfc6edd0cfbabd4403ae550b31be86b
--- /dev/null
+++ b/Conductors3D/w3d_common.pro
@@ -0,0 +1,57 @@
+
+mm = 1e-3;
+deg = 180/Pi;
+
+/*
+ Expected inductance of a rectilinear 1mm radius wire: 4 nH/cm
+ */
+
+DefineConstant[
+ NumWires = {1, Name "Parameters/0Number of wires",
+ Choices {1="1", 2="2", 3="3"}, Visible 1}
+ Flag_Form = {2, Name "Parameters/1Conductor type", Visible 1,
+ Choices {1="Massive", 2="Stranded"}}
+ Flag_Thin = {1, Name "Parameters/2Use thin wires", Choices {0,1}, Visible 1}
+ Flag_Corr = {1, Name "Parameters/3Correction method",
+ Visible (Flag_Form==2) && Flag_Thin,
+ Choices {0="None", 1="Raw sleeve", 2="regular sleeve"}}
+ Flag_U = {0, Name "Parameters/4Impose U", Choices {0,1}, Visible 1}
+ Flag_Maps = {1, Name "Input/2Display maps", Choices {0,1}, Visible 1}
+
+ LogFreq = {4, Min 1, Max 8, Step 0.25, Name "Input/4Log of frequency"}
+ Freq = 10^LogFreq
+ WireRadius = {1, Name "Parameters/5Wire radius [mm]"}
+ MeshSizeWire = {1, Name "Parameters/6Mesh size in wire [mm]",
+ Visible !Flag_Thin}
+ MeshSizeThinWire = {1, Name "Parameters/7Mesh size on thin wire [mm]",
+ Visible Flag_Thin}
+ SleeveRadius = {2, Name "Parameters/8Sleeve radius [mm]"}
+ box = {100*mm, Name "Parameters/9box half-width [m]"}
+
+ rw = WireRadius*mm
+ rs = SleeveRadius*mm
+ A_c = Pi*rw^2 // wire cross section
+ rout = box
+ Lz = 10*mm
+
+ /* 'WireRadius' is the actual radius of the wire.
+ 'MeshSizeWire' is the imposed mesh size at the nodes of the wire,
+ and thus also the practical approximative value of the radius of the sleeve.
+ The 'rs' and 'rw' parameters are assumed identical for all wires.
+ */
+];
+
+
+
+Xlocs() = { 0, 5*mm, -5*mm};
+Ylocs() = { 0, 0, 0};
+For i In {1:NumWires}
+ DefineConstant[
+ WX~{i} = { Xlocs(i-1), // place wires symmetrically around x=0
+ Name Sprintf("Parameters/Wire %g/0X position [m]", i), Closed }
+ WY~{i} = { Ylocs(i-1), Name Sprintf("Parameters/Wire %g/0Y position [m]", i) }
+ WI~{i} = { 1.23456, Name Sprintf("Parameters/Wire %g/2Current [A]", i) }
+ WP~{i} = { 0*NumWires*(i-1),
+ Name Sprintf("Parameters/Wire %g/3Phase [deg]", i) }
+ ];
+EndFor