OCC geo for Ishape3D + removing not used 2D parameters

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 {

ECE_full_wave/Ishape3d.geo

 /*  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;

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

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
-*/

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