diff --git a/ECE_full_wave/FullWave_E_ece_SISO.pro b/ECE_full_wave/FullWave_E_ece_SISO.pro
index 356d651972a2c28a4a26e7fe05c255179246477f..a8eb1b687ef79f1c61e450ff38678532843fc03c 100644
--- a/ECE_full_wave/FullWave_E_ece_SISO.pro
+++ b/ECE_full_wave/FullWave_E_ece_SISO.pro
@@ -13,10 +13,6 @@
Frequency domain simulations.
*/
-DefineConstant[
- h2Ddepth= 1 // Default value in 3D, to be adapted if axisymmetry or 2D planar
-];
-
Jacobian {
{ Name Vol;
@@ -167,9 +163,9 @@ PostProcessing {
{ Name I;
Value {
- Term { [ -1*{I}*h2Ddepth ]; In Sur_Terminals_FWece; }
+ Term { [ -1*{I} ]; In Sur_Terminals_FWece; }
}
- } // h2Ddepth is the depth for 2D problems, and it has to be 1 for 3D problems
+ }
{ Name VnotTerminals;
Value {
diff --git a/ECE_full_wave/Ishape3d.geo b/ECE_full_wave/Ishape3d.geo
index 0bd01a9620a02b2c3300c23b944723aaed49a802..8a95d1191c1d6bbaa628d00515adcd85484ac6dc 100644
--- a/ECE_full_wave/Ishape3d.geo
+++ b/ECE_full_wave/Ishape3d.geo
@@ -1,81 +1,123 @@
/* Ishape3d.geo
Geometrical description (for gmsh) of Ishape3D test for ECE
For details, see comments in the Ishape3d_data.pro file
-Meshing information is also defined here.
+Meshing information is also defined here.
+
+This uses OpenCascade
*/
Include "Ishape3d_data.pro";
+SetFactory("OpenCASCADE");
-/* Definition of parameters for local mesh dimensions */
-lcar = s*l/10; // characteristic length of mesh element
-
-
-p0 = newp; Point(p0) = {0,0,0,lcar};
-p1 = newp; Point(p1) = {a,0,0,lcar};
-p2 = newp; Point(p2) = {0,a,0,lcar};
-p3 = newp; Point(p3) = {-a,0,0,lcar};
-p4 = newp; Point(p4) = {0,-a,0,lcar};
-
-p0l = newp; Point(p0l) = {0,0,l,lcar};
-p1l = newp; Point(p1l) = {a,0,l,lcar};
-p2l = newp; Point(p2l) = {0,a,l,lcar};
-p3l = newp; Point(p3l) = {-a,0,l,lcar};
-p4l = newp; Point(p4l) = {0,-a,l,lcar};
-
-
-c1 = newc; Circle(c1) = {p1,p0,p2};
-c2 = newc; Circle(c2) = {p2,p0,p3};
-c3 = newc; Circle(c3) = {p3,p0,p4};
-c4 = newc; Circle(c4) = {p4,p0,p1};
-
-c1l = newc; Circle(c1l) = {p1l,p0l,p2l};
-c2l = newc; Circle(c2l) = {p2l,p0l,p3l};
-c3l = newc; Circle(c3l) = {p3l,p0l,p4l};
-c4l = newc; Circle(c4l) = {p4l,p0l,p1l};
-
-l1 = newl; Line(l1) = {p1, p1l};
-l2 = newl; Line(l2) = {p2, p2l};
-l3 = newl; Line(l3) = {p3, p3l};
-l4 = newl; Line(l4) = {p4, p4l};
-
-cL = newll; Curve Loop(cL) = {c1,c2,c3,c4}; scL = news; Surface(scL) = {cL};
-cLl = newll; Curve Loop(cLl) = {c1l,c2l,c3l,c4l}; scLl = news; Surface(scLl) = {cLl};
-cL12 = newll; Curve Loop(cL12) = {c1,l2,-c1l,-l1}; scL12 = news; Surface(scL12) = {cL12};
-cL23 = newll; Curve Loop(cL23) = {c2,l3,-c2l,-l2}; scL23 = news; Surface(scL23) = {cL23};
-cL34 = newll; Curve Loop(cL34) = {c3,l4,-c3l,-l3}; scL34 = news; Surface(scL34) = {cL34};
-cL41 = newll; Curve Loop(cL41) = {c4,l1,-c4l,-l4}; scL41 = news; Surface(scL41) = {cL41};
-
-closedS = newsl; Surface Loop(closedS) = {scL,scLl,scL12,scL23,scL34,scL41};
-volDom = newv; Volume(volDom) = {closedS};
-
-/*
-Geometry.PointNumbers = 1;
-Geometry.Color.Points = Orange;
-General.Color.Text = White;
-Mesh.Color.Points = {255, 0, 0}; //red
-*/
+xc = 0; yc = 0; zc = 0;
+vx = 0; vy = 0; vz = l;
+radius = a;
-/*
-If(_use_transfinite)
- Transfinite Line {p1, p1l} = 3*s;
- Transfinite Line {p2, p2l} = 3*s;
- Transfinite Line {p3, p3l} = 3*s;
- Transfinite Line {p4, p4l} = 3*s;
- Transfinite Surface {scL12};
- Transfinite Surface {scL23};
- Transfinite Surface {scL34};
- Transfinite Surface {scL41};
- Transfinite Volume{volDom};
-EndIf
-*/
+volDom = newv; Cylinder(newv) = {xc, yc, zc, vx, vy, vz, radius};
-Physical Volume ("MaterialX", 100) = {volDom} ; /* MaterialX */
+b() = Boundary{ Volume{volDom}; };
+Printf(" surfaces", b());
+
+// In order to identify which surface is which I played with the GUI!!!
+
+lateralCyl = b(0);
+bottomSurf = b(1);
+topSurf = b(2);
-Physical Surface ("Ground", 120) = {scLl} ; /* Ground */
-Physical Surface ("Terminal", 121) = {scL} ; /* Terminal */
-Physical Surface ("BoundaryNotTerminal", 131) = {scL12,scL23,scL34,scL41} ;
+//bs() = Boundary{ Surface{bottomSurf}; };
+//Printf(" curves", bs());
+//circle1 = bs(0);
+lc() = Boundary{ Surface{lateralCyl}; };
+Printf(" curves", lc());
+circle1 = lc(0);
+lineCyl = lc(1);
+circle2 = lc(2);
+//Transfinite Line {lineCyl} = 20;
+//Transfinite Curve{circle1} = 50;
+// Physical regions
+// --------------------------
+Physical Volume ("MaterialX", 100) = {volDom} ; /* MaterialX */
+
+Physical Surface ("Ground", 120) = {bottomSurf} ; /* Ground */
+Physical Surface ("Terminal", 121) = {topSurf} ; /* Terminal */
+Physical Surface ("BoundaryNotTerminal", 131) = {lateralCyl} ;
+
+
+/* Definition of parameters for local mesh dimensions */
+//lcar = s*l/10; // characteristic length of mesh element
+
+// lcar dep on Freq, wavelenth - 10 elements per wavelength and dimension
+//lambda = c0/Freq; // Label "Wavelength [m]"}
+//lcar = lambda/nbla
+// nbDelta from _data.geo, number of elements per skin depth
+
+delta = Sqrt(2.0/(2*Pi*Freq*mu*sigma));
+lcarDelta = delta/nbDelta;
+lcar = s*l/10;
+
+
+
+ //Field[1] = Distance;
+ //Field[1].CurvesList = {circle1};
+ //Field[1].SurfacesList = {lateralCyl};
+ //Field[1].Sampling = 200;
+
+ //Field[2] = Threshold;
+ //Field[2].InField = 1;
+ //Field[2].SizeMin = a/10; //lcar;
+ //Field[2].SizeMax = a/2;
+ //Field[2].DistMin = a/8; //delta;
+ //Field[2].DistMax = a/4; //delta;
+
+ //Field[2] = Box;
+ //Field[2].VIn = a/2;
+ //Field[2].VOut = a/10; //lcar;
+ //Field[2].XMin = -2*a/3;
+ //Field[2].XMax = 2*a/3;
+ //Field[2].YMin = -2*a/3;
+ //Field[2].YMax = 2*a/3;
+ //Field[2].ZMin = 0;
+ //Field[2].ZMax = l;
+ //Field[2].Thickness = 0.3;
+
+ Field[1] = Cylinder;
+ Field[1].VIn = lcar;
+ Field[1].VOut = lcar;
+ Field[1].Radius = a;
+ Field[1].XCenter = 0;
+ Field[1].YCenter = 0;
+ Field[1].ZCenter = l/2;
+ Field[1].XAxis = 0;
+ Field[1].YAxis = 0;
+ Field[1].ZAxis = l;
+
+ Field[2] = Cylinder;
+ Field[2].VIn = a/3;
+ Field[2].VOut = lcarDelta;
+ Field[2].Radius = a - delta;
+ Field[2].XCenter = 0;
+ Field[2].YCenter = 0;
+ Field[2].ZCenter = l/2;
+ Field[2].XAxis = 0;
+ Field[2].YAxis = 0;
+ Field[2].ZAxis = l;
+
+
+ Field[3] = Min;
+ Field[3].FieldsList = {1,2};
+ Background Field = 3;
+ Mesh.MeshSizeExtendFromBoundary = 0;
+ Mesh.MeshSizeFromPoints = 0;
+ Mesh.MeshSizeFromCurvature = 0;
+ Mesh.Algorithm = 5;
+
+
+/* Definition of parameters for local mesh dimensions */
+//lcar = s*l/10; // characteristic length of mesh element
+// Assign a mesh size to all the points of all the volumes:
+//MeshSize{ PointsOf{ Volume{:}; } } = lcar;
diff --git a/ECE_full_wave/Ishape3d.pro b/ECE_full_wave/Ishape3d.pro
index 62a2d0e824b6c71ce443194e6672b89ea479b2fd..8ec428bc6d23bc18e2a8faf93b6c0fae874c89c3 100644
--- a/ECE_full_wave/Ishape3d.pro
+++ b/ECE_full_wave/Ishape3d.pro
@@ -58,7 +58,7 @@ Group{
{ Name SetTerminalCurrent; Type Assign; // current excited terminals
Case {
If(Flag_AnalysisType==1)
- { Region Terminal; Value ITerminal1[]/h2Ddepth; } // here the depth is needed
+ { Region Terminal; Value ITerminal1[]; }
EndIf
}
}
@@ -90,34 +90,51 @@ Echo[Str[ "nv = PostProcessing.NbViews;",
], File "post.opt"];
Return
- PostOperation {
- { Name Maps; NameOfPostProcessing FW_E_ece ;
- Operation {
- Print [ gradV, OnElementsOf ElementsOf[Vol_FW,OnOneSideOf Sur_FW], File StrCat[Dir,"map_gradV.pos" ]];
+Macro WriteZ_header
+ Echo[ Str[Sprintf("# Hz Z RI R %g", 50)], Format Table, File > StrCat[Dir, "test_Z_RI.s1p"] ];
+Return
+Macro WriteY_header
+ Echo[ Str[Sprintf("# Hz Y RI R %g", 50)], Format Table, File > StrCat[Dir, "test_Y_RI.s1p"] ];
+Return
- Print [ E, OnElementsOf Vol_FW, File StrCat[Dir,"map_E.pos" ]];
- // Print [ J, OnElementsOf Vol_FW, File StrCat[Dir,"map_J.pos" ]];
- // Print [ rmsJ, OnElementsOf Vol_FW, File StrCat[Dir,"map_absJ.pos" ]];
- Print [ H, OnElementsOf Vol_FW, File StrCat[Dir,"map_H.pos" ]];
- Print [ Vterminals, OnElementsOf ElementsOf[Sur_Terminals_FWece], File StrCat[Dir,"map_Vterminals.pos" ]];
- Print [ VnotTerminals, OnElementsOf Vol_FW, File StrCat[Dir,"map_VnotTerminals.pos" ]];
- Call Change_post_options;
- }
+PostOperation {
+ { Name Maps; NameOfPostProcessing FW_E_ece ;
+ Operation {
+ Print [ gradV, OnElementsOf ElementsOf[Vol_FW,OnOneSideOf Sur_FW], File StrCat[Dir,"map_gradV.pos" ]];
+
+ Print [ E, OnElementsOf Vol_FW, File StrCat[Dir,"map_E.pos" ]];
+ // Print [ J, OnElementsOf Vol_FW, File StrCat[Dir,"map_J.pos" ]];
+ // Print [ rmsJ, OnElementsOf Vol_FW, File StrCat[Dir,"map_absJ.pos" ]];
+
+ Print [ H, OnElementsOf Vol_FW, File StrCat[Dir,"map_H.pos" ]];
+ Print [ Vterminals, OnElementsOf ElementsOf[Sur_Terminals_FWece], File StrCat[Dir,"map_Vterminals.pos" ]];
+ Print [ VnotTerminals, OnElementsOf Vol_FW, File StrCat[Dir,"map_VnotTerminals.pos" ]];
+
+ Call Change_post_options;
+ }
}
- { Name TransferMatrix; NameOfPostProcessing FW_E_ece ;
- Operation {
- If(Flag_AnalysisType==1) // current excitation
- Print [ Vterminals, OnRegion Terminal,
- Format Table , File > StrCat[Dir, "test_Z_RI.s1p"]] ;
- Else
- Print [ I, OnRegion Terminal,
- Format Table , File > StrCat[Dir, "test_Y_RI.s1p"]] ;
- EndIf
- }
- }
+ { Name TransferMatrix; NameOfPostProcessing FW_E_ece ;
+ Operation {
+ If(Flag_AnalysisType==1) // current excitation
+ If (Freq==freqs(0))
+ Call WriteZ_header;
+ Printf("======================> Writting header with first Freq=%g",Freq);
+ EndIf
+ Print [ Vterminals, OnRegion Terminal,
+ Format Table , File > StrCat[Dir, "test_Z_RI.s1p"]] ;
+ Else
+ If (Freq==freqs(0))
+ Call WriteY_header;
+ Printf("======================> Writting header with first Freq=%g",Freq);
+ EndIf
+ Print [ I, OnRegion Terminal,
+ Format Table , File > StrCat[Dir, "test_Y_RI.s1p"]] ;
+ EndIf
+ }
+ }
}
diff --git a/ECE_full_wave/Ishape3d_data.pro b/ECE_full_wave/Ishape3d_data.pro
index 074bf739a7bdf7376ad9458697fd29b5e185e756..f5588e5082b2346df48de05dc3ffe8994ce7d33a 100644
--- a/ECE_full_wave/Ishape3d_data.pro
+++ b/ECE_full_wave/Ishape3d_data.pro
@@ -64,7 +64,7 @@ mu = mu0*mur; // F/m - absolute permeability
nu = 1/mu; // m/H - absolute reluctivity
-/* new menu for the interface - type of BC */
+// new menu for the interface - type of BC
mTypeBC = StrCat[cGUI,"Type of BC/0"];
colorMTypeBC = "Red";
DefineConstant[
@@ -77,59 +77,46 @@ DefineConstant[
];
-
-/* new menu for the interface - values of excitations - significance depends on the type of BC */
+// new menu for the interface - values of excitations - significance depends on the type of BC
mValuesBC = StrCat[cGUI,"Values of BC (frequency analysis)/0"];
colorMValuesBC = "Ivory";
fmin = 1e7; // Hz
-//fmed = 3e9; // Hz
-//fmax = 100e9; // Hz
+fmax = 100e9; // Hz
+nop = 10;
+//freqs()=LogSpace[Log10[fmin],Log10[fmax],nop];
+freqs()=LinSpace[fmin,fmax,nop];
DefineConstant[
- Freq = {fmin, Name StrCat[mValuesBC, "0Working Frequency"], Units "Hz", Highlight Str[colorMValuesBC], Closed !close_menu }
-];
-// You may also use UndefineConstant[]
+ Freq = {freqs(0), Choices{freqs()}, Loop 0, Name StrCat[mValuesBC, "0Working Frequency"],
+ Units "Hz", Highlight Str[colorMValuesBC], Closed !close_menu }
-DefineConstant[
// Flag_AnalysisType
// ==0 bottom ev
// ==1 bottom ec
- VGroundRMS = {0, Name StrCat[mValuesBC, "2Top/1Potential on bottom (rms)"], Units "V" , ReadOnly 1, Highlight Str[colorro],
+ VGroundRMS = {0, Name StrCat[mValuesBC, "2Top/1Potential on bottom (rms)"],
+ Units "V", ReadOnly 1, Highlight Str[colorro],
Visible (Flag_AnalysisType==0) || (Flag_AnalysisType==1)}
- VGroundPhase = {0, Name StrCat[mValuesBC, "2Top/1Potential on bottom, phase"], Units "rad", ReadOnly 1, Highlight Str[colorro],
+ VGroundPhase = {0, Name StrCat[mValuesBC, "2Top/1Potential on bottom, phase"],
+ Units "rad", ReadOnly 1, Highlight Str[colorro],
Visible (Flag_AnalysisType==0) || (Flag_AnalysisType==1)}
- VTermRMS = {1, Name StrCat[mValuesBC, "1Bottom/1Potential on Term1 (rms)"], Units "V" , ReadOnly 0, Highlight Str[colorMValuesBC],
+ VTermRMS = {1, Name StrCat[mValuesBC, "1Bottom/1Potential on Term1 (rms)"],
+ Units "V", ReadOnly 0, Highlight Str[colorMValuesBC],
Visible (Flag_AnalysisType==0)}
- VTermPhase = {0, Name StrCat[mValuesBC, "1Bottom/1Potential on Term1, phase"], Units "rad", ReadOnly 0, Highlight Str[colorMValuesBC],
+ VTermPhase = {0, Name StrCat[mValuesBC, "1Bottom/1Potential on Term1, phase"],
+ Units "rad", ReadOnly 0, Highlight Str[colorMValuesBC],
Visible (Flag_AnalysisType==0)}
- ITermRMS = {1, Name StrCat[mValuesBC, "1Bottom/1Current through Term1 (enters the domain) (rms)"], Units "A" , ReadOnly 0, Highlight Str[colorMValuesBC],
+ ITermRMS = {1, Name StrCat[mValuesBC, "1Bottom/1Current through Term1 (enters the domain) (rms)"],
+ Units "A", ReadOnly 0, Highlight Str[colorMValuesBC],
Visible (Flag_AnalysisType==1) }
- ITermPhase = {0, Name StrCat[mValuesBC, "1Bottom/Current through Term1 (enters the domain), phase"], Units "rad", ReadOnly 0, Highlight Str[colorMValuesBC],
+ ITermPhase = {0, Name StrCat[mValuesBC, "1Bottom/Current through Term1 (enters the domain), phase"],
+ Units "rad", ReadOnly 0, Highlight Str[colorMValuesBC],
Visible (Flag_AnalysisType==1) }
-];
-
-
-
-DefineConstant[
//_use_transfinite = {0, Choices{0,1}, Name "{MeshSettings/0Transfinite mesh?", Highlight "Yellow", Closed !close_menu}
- s = {1, Name "{MeshSettings/1Mesh factor"}
+ s = {1, Name "{MeshSettings/1Mesh factor s"}
+ nbDelta = {1, Name "{MeshSettings/2Nb of elems per skin depth"}
];
-
-
-/*
-modelDim = 3; //
-Flag_Axi = 0; // 1 for AXI - it makes sense only for modelDim = 2
-
-If ((modelDim == 2)&&(Flag_Axi == 0))
- h2Ddepth = h;
-ElseIf ((modelDim == 2)&&(Flag_Axi == 1)) // 2D AXI
- h2Ddepth = 2*Pi;
-Else // 3D
- h2Ddepth = 1;
-EndIf
-*/
diff --git a/ECE_full_wave/LC.pro b/ECE_full_wave/LC.pro
index 71c7f1135c38b1391b70269fd35e86e63425174f..ef6fbe2ceb104c2e9b2750bcb133a48759b1158d 100644
--- a/ECE_full_wave/LC.pro
+++ b/ECE_full_wave/LC.pro
@@ -70,7 +70,7 @@ Group{
{ Name SetTerminalCurrent; Type Assign; // current excited terminals
Case {
If(Flag_AnalysisType==1)
- { Region Terminal; Value ITerminal1[]/h2Ddepth; } // here the depth is needed in 2D problems
+ { Region Terminal; Value ITerminal1[]; }
EndIf
}
}