ECE full wave formulation in E with two examples

ECE_full_wave/FullWave_E_ece_SISO.pro

+/* FullWave_E_ece_SISO.pro
+
+ Problem independent pro file
+ for passive, linear domains in Full Wave with ECE boundary conditions.
+
+ It uses a formulation with vec{E} strictly inside the domain and V strictly on the boundary.
+ The domain is simply connected, it has only electric terminals.
+ The terminals can be voltage or current excited.
+ The material inside is linear from all points of view: electric, magnetic, conduction
+ It is compulsory that there exists a ground terminal (its voltage is set to zero).
+ The post-operation assumes that there is one excited (active) terminal  (excitation equal to 1 - it can be in
+ voltage or current excited terminal).
+ Frequency domain simulations.
+*/
+
+DefineConstant[
+  h2Ddepth= 1 // Default value in 3D, to be adapted if axisymmetry or 2D planar
+];
+
+
+Jacobian {
+  { Name Vol;
+    Case {
+      { Region All; Jacobian Vol; }
+    }
+  }
+  { Name Sur;
+    Case {
+      { Region All; Jacobian Sur; }
+    }
+  }
+}
+
+Integration {
+  { Name Int;
+    Case {
+      { Type Gauss;
+        Case {
+          { GeoElement Point; NumberOfPoints  1; }
+          { GeoElement Line; NumberOfPoints  3; }
+          { GeoElement Triangle; NumberOfPoints  3; }
+          { GeoElement Quadrangle; NumberOfPoints  4; }
+          { GeoElement Tetrahedron; NumberOfPoints  4; }
+          { GeoElement Hexahedron; NumberOfPoints  6; }
+          { GeoElement Prism; NumberOfPoints  9; }
+          { GeoElement Pyramid; NumberOfPoints  8; }
+	}
+      }
+    }
+  }
+}
+
+FunctionSpace {
+
+  { Name Hcurl_E; Type Form1;
+    BasisFunction {
+      { Name se; NameOfCoef ee; Function BF_Edge;
+        Support Dom_FW ; Entity EdgesOf[All, Not Sur_FW]; }
+      { Name sn; NameOfCoef vn; Function BF_GradNode;
+        Support Dom_FW ; Entity NodesOf[Sur_FW, Not Sur_Terminals_FWece]; }
+      { Name sf; NameOfCoef vf; Function BF_GradGroupOfNodes;
+        Support Dom_FW ; Entity GroupsOfNodesOf[Sur_Terminals_FWece]; }
+    }
+    GlobalQuantity {
+      { Name TerminalPotential; Type AliasOf       ; NameOfCoef vf; }
+      { Name TerminalCurrent  ; Type AssociatedWith; NameOfCoef vf; }
+    }
+
+    SubSpace {
+      { Name dv ; NameOfBasisFunction {sn}; } // Subspace, it maybe use in equations or post-pro
+    }
+
+    Constraint {
+      { NameOfCoef TerminalPotential; EntityType GroupsOfNodesOf; NameOfConstraint SetTerminalPotential; }
+      { NameOfCoef TerminalCurrent;   EntityType GroupsOfNodesOf; NameOfConstraint SetTerminalCurrent; }
+    }
+  }
+
+}
+
+Formulation {
+
+  { Name FullWave_E_ece; Type FemEquation;
+    Quantity {
+      { Name e;  Type Local; NameOfSpace Hcurl_E; }
+      { Name dv; Type Local; NameOfSpace Hcurl_E[dv]; } // Just for post-processing issues
+      { Name V;  Type Global; NameOfSpace Hcurl_E[TerminalPotential]; }
+      { Name I;  Type Global; NameOfSpace Hcurl_E[TerminalCurrent]; }
+    }
+    Equation {
+      // \int_D curl(\vec{E}) \cdot  curl(\vec{e}) dv
+      Galerkin { [ nu[] * Dof{d e} , {d e} ]; In Vol_FW; Jacobian Vol; Integration Int; }
+
+      // \int_D j*\omega*(\sigma + j*\omega*\epsilon) \vec{E} \cdot \vec{e} dv
+      Galerkin { DtDof   [ sigma[]   * Dof{e} , {e} ]; In Vol_FW; Jacobian Vol; Integration Int; }
+      Galerkin { DtDtDof [ epsilon[] * Dof{e} , {e} ]; In Vol_FW; Jacobian Vol; Integration Int; }
+
+      //    alternative // j*\omega*sum (vk Ik);  for k - current excited terminals
+      //    GlobalTerm {DtDof [ -Dof{I} , {V} ]; In SurBCec; }
+
+      // j*\omega*sum (vk Ik);  for k - all terminals so that the currents through the terminals will be computed as well
+      GlobalTerm {DtDof [ -Dof{I} , {V} ]; In Sur_Terminals_FWece; }
+    }
+  }
+
+}
+
+
+Resolution {
+
+  { Name FullWave_E_ece;
+    System {
+      { Name Sys_Ele; NameOfFormulation FullWave_E_ece; Type Complex; Frequency Freq;}
+    }
+    Operation {
+      CreateDir[Dir];
+      SetTime[Freq];  // useful to save (and print) the current frequency
+      Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele];
+    }
+  }
+
+}
+
+PostProcessing {
+
+  { Name FW_E_ece; NameOfFormulation FullWave_E_ece;
+    Quantity {
+      { Name E;
+        Value {
+          Local { [ {e} ]; In Vol_FW; Jacobian Vol; }
+        }
+      }
+
+      { Name normE;
+        Value {
+          Local { [ Norm[{e}] ]; In Vol_FW; Jacobian Vol; }
+        }
+      }
+
+
+      { Name J;
+        Value {
+          Local { [ sigma[]*{e} ]; In Vol_FW; Jacobian Vol; }
+        }
+      }
+
+      { Name rmsJ;
+        Value {
+          Local { [ Norm[sigma[]*{e}] ]; In Vol_FW; Jacobian Vol; }
+        }
+      }
+
+      { Name gradV;
+        Value {
+          Local { [ {dv} ]; In Vol_FW; Jacobian Vol; }
+        }
+      }
+
+      { Name Vterminals;
+        Value {
+          // This is a global variable, you do not need to indicate the Jacobian
+          // It is only defined at the terminal and it is only one value for the whole line (surface in 3D)
+          // If the support you give here is wider, you will get 0 where it is not defined
+          Term { [ -1*{V} ]; In Sur_Terminals_FWece; } // => unknowns only in Sur_Terminals_FWece, outside it will be zero
+        }
+      }
+
+      { Name I;
+        Value {
+          Term { [ -1*{I}*h2Ddepth ]; In Sur_Terminals_FWece; }
+        }
+      }  // h2Ddepth is the depth for 2D problems,  and it has to be 1 for 3D problems
+
+      { Name VnotTerminals;
+        Value {
+          Local { [ -{dInv dv} ]; In Vol_FW; Jacobian Vol; }
+        }
+      }
+
+      { Name B;
+        Value {
+          Local { [ {d e}/(2*Pi*Freq*Complex[0,1]) ]; In Vol_FW; Jacobian Vol; }
+        }
+      }
+
+      { Name H;
+        Value {
+          Local { [ nu[]*{d e}/(2*Pi*Freq*Complex[0,1]) ]; In Vol_FW; Jacobian Vol; }
+        }
+      }
+
+
+
+      { Name Bz;
+        Value {
+          Local { [ CompZ[{d e}]/(2*Pi*Freq*Complex[0,1]) ]; In Vol_FW; Jacobian Vol; }
+        }
+      }
+
+      { Name Hz;
+        Value {
+          Local { [ nu[]*CompZ[{d e}]/(2*Pi*Freq*Complex[0,1]) ]; In Vol_FW; Jacobian Vol; }
+        }
+      }
+
+    }
+   }
+ }

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.
+*/
+
+
+Include "Ishape3d_data.pro";
+
+/* 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
+*/
+
+/*
+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
+*/
+
+Physical Volume ("MaterialX", 100) = {volDom} ;      /* MaterialX */
+
+Physical Surface ("Ground",   120) = {scLl} ;          /* Ground */
+Physical Surface ("Terminal", 121) = {scL} ;           /* Terminal */
+Physical Surface ("BoundaryNotTerminal", 131) = {scL12,scL23,scL34,scL41} ;  
+
+
+
+

ECE_full_wave/Ishape3d.pro

+/*  Ishape3d.pro
+Material description (for getdp) of Ishape3D test for ECE
+For details, see comments in the Ishape3d_data.pro file
+The Ishape3d geometry and meshing parameters are in the Ishape3d.geo file
+This file contains the setings that depend on the particular problem and
+it includes 1 general files
+*/
+
+Include "Ishape3d_data.pro"
+
+Group{
+  MaterialX = Region[100]; // MaterialX - homogenous domain
+  Ground = Region[120];    // Boundary - various parts
+  Terminal = Region[121];
+  BoundaryNotTerminal = Region[131];
+
+  Vol_FW = Region[{MaterialX}]; // domain without the boundary
+  Sur_Terminals_FWece = Region[{Ground, Terminal}]; // all terminals
+  Sur_FW = Region[{Sur_Terminals_FWece, BoundaryNotTerminal}]; //all boundary
+  Dom_FW = Region[ {Vol_FW, Sur_FW} ];
+ }
+
+
+ Function{
+   /*
+   Material properties are defined piecewise, for elementary groups
+   the values in the right hand side are defined in Ishape2d_data.pro
+   */
+
+   nu[#{MaterialX,Sur_FW}] = nu;             // magnetic reluctivity [m/H]
+   epsilon[#{MaterialX,Sur_FW}] = eps;       // permittivity [F/m]
+   sigma[#{MaterialX,Sur_FW}]  = sigma;      // conductivity [S/m]
+
+   // This is a problem with 2 terminals, out of which one is grounded.
+   // The non-grounded can be excited in voltage or in current.
+   // Depending on the excitation type, one of the following functions will be useful.
+   //
+   // When imposing voltage or current via a constraint, the region appears there
+   // If you have a complex, the value has to be a function => use square brackets []
+   If(Flag_AnalysisType==0)
+     VTerminal1[] = -VTermRMS*Complex[Cos[VTermPhase],Sin[VTermPhase]];
+   Else
+     ITerminal1[] = -ITermRMS*Complex[Cos[ITermPhase],Sin[ITermPhase]];
+   EndIf
+ }
+
+ Constraint{
+
+   // ece BC
+   { Name SetTerminalPotential; Type Assign; // voltage excited terminals
+     Case {
+       { Region Ground; Value 0.; }
+       If(Flag_AnalysisType==0)
+         { Region Terminal; Value VTerminal1[]; }
+       EndIf
+     }
+   }
+   { Name SetTerminalCurrent; Type Assign; // current excited terminals
+     Case {
+       If(Flag_AnalysisType==1)
+         { Region Terminal; Value ITerminal1[]/h2Ddepth; }  // here the depth is needed
+       EndIf
+     }
+   }
+
+}
+
+// For convenience, saving results in different directories
+// Dir can be empty: Dir ="";
+ If (Flag_AnalysisType==0)
+   Dir = "res/FWeceBC_voltExc/";
+ Else
+   Dir = "res/FWeceBC_crtExc/";
+ EndIf
+
+ // The formulation and its tools
+Include "FullWave_E_ece_SISO.pro"
+
+
+// Some Macros for changing Gmsh options
+// Strictly not necessary => for aestetics or easying result treatment
+Macro Change_post_options
+Echo[Str[ "nv = PostProcessing.NbViews;",
+    "For v In {0:nv-1}",
+    "View[v].NbIso = 25;",
+    "View[v].RangeType = 3;" ,// per timestep
+    "View[v].ShowTime = 3;",
+    "View[v].IntervalsType = 3;",
+    "EndFor"
+  ], 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" ]];
+
+       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
+     }
+   }
+
+}

ECE_full_wave/Ishape3d_data.pro

+/* Ishape3d_data.pro
+=========================== Ishape3D - test problem
+A cylindrical domain in space, one circular face in the plane z = 0, circle of radius a,
+axis of the cylinder is Oz, length = l
+
+The domain is filled with homogeneous and linear material, with epsilon, mu, sigma
+
+There is no no field source inside the domein (the "device" is passive).
+The field sources are only boundary conditions (BC) which are ECE.
+
+The problem has 2 terminals (on top and bottom), one has to be grounded, the other two can be either
+voltage excited or current excited.
+
+This is a test problem - used to validate the ECE implementation in onelab for a 3D problem
+
+===========================
+Names of input data - to be given by the user
+Geometry:
+a = radius of the cylinder [m]
+l = height of the cylinder [m]
+==
+Material
+epsr = relative permittivity
+mur = relative permeability
+sigma = conductivity [S/m]
+
+==
+Type of BC - always ece, top terminal is gnd; bottom terminal can be ev = electric voltage forced terminal; ec = electric current forced terminal
+typeBC = 0,1,
+0 = bottom ev
+1 = bottom ec
+==
+Excitation signals
+imposed by the user, on non-grounded terminals; default values - is 1 (V or A)
+=========================== */
+
+colorro = "LightGrey";
+
+/* new menu for the interface - holds geometrical data */
+cGUI= "{Ishape3D-ProblemData/";
+mGeo = StrCat[cGUI,"Geometrical parameters/0"];
+colorMGeo = "Ivory";
+close_menu = 1;
+DefineConstant[
+  a = {2.5e-6, Name StrCat[mGeo, "0a=0.5*Dimension along Ox"], Units "m", Highlight Str[colorMGeo], Closed close_menu},
+  l = {10e-6, Name StrCat[mGeo, "1l=Dimension along Oy"], Units "m", Highlight Str[colorMGeo]}
+];
+
+/* new menu for the interface - holds material data */
+mMat = StrCat[cGUI,"Material parameters/0"];
+colorMMat = "AliceBlue";
+DefineConstant[
+  epsr = {1, Name StrCat[mMat, "0Relative permittivity T part"], Units "-", Highlight Str[colorMMat], Closed close_menu},
+  mur = {1, Name StrCat[mMat, "1Relative permeability T part"], Units "-", Highlight Str[colorMMat]},
+  sigma = {6.6e7, Name StrCat[mMat, "2Conductivity T part"], Units "S/m", Highlight Str[colorMMat]}
+];
+
+eps0 = 8.854187818e-12; // F/m - absolute permitivity of vacuum
+mu0  = 4.e-7 * Pi;      // H/m - absolute permeability of vacuum
+nu0  = 1/mu0;           // m/H - absolute reluctivity of vacuum
+
+eps = eps0*epsr;        // F/m - absolute permitivity
+mu = mu0*mur;           // F/m - absolute permeability
+nu = 1/mu;              // m/H - absolute reluctivity
+
+
+/* new menu for the interface - type of BC */
+mTypeBC = StrCat[cGUI,"Type of BC/0"];
+colorMTypeBC = "Red";
+DefineConstant[
+  Flag_AnalysisType = { 0,
+    Choices{
+      0="bottomTerm ev",
+      1="bottomTerm ec"},
+    Name Str[mTypeBC], Highlight Str[colorMTypeBC], Closed !close_menu,
+    ServerAction Str["Reset", "GetDP/1ResolutionChoices"]}
+];
+
+
+
+/* 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
+
+DefineConstant[
+  Freq  = {fmin, Name StrCat[mValuesBC, "0Working Frequency"], Units "Hz", Highlight Str[colorMValuesBC], Closed !close_menu }
+];
+// You may also use UndefineConstant[]
+
+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],
+    Visible (Flag_AnalysisType==0) || (Flag_AnalysisType==1)}
+  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],
+    Visible (Flag_AnalysisType==0)}
+  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],
+    Visible (Flag_AnalysisType==1) }
+  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"}
+];
+
+
+/*
+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.geo

+Include "LC_data.pro";
+
+SetFactory("OpenCASCADE");
+
+// convert geometry read from STEP file to meters
+Geometry.OCCTargetUnit = "M";
+
+Geometry.NumSubEdges = 100; // Visualisation
+
+// We can activate the calculation of mesh element sizes based on curvature:
+Mesh.MeshSizeFromCurvature = 1; // Needed for spiral/helix
+
+// And we set the minimum number of elements per 2*Pi radians:
+Mesh.MinimumElementsPerTwoPi = 20/2;
+
+// We can constraint the min and max element sizes to stay within reasonnable
+// values (see `t10.geo' for more details):
+Mesh.MeshSizeMin = rw/4;
+Mesh.MeshSizeMax = rplateC;
+
+
+
+// read STEP file and store volume(s)
+v() = ShapeFromFile("LC_ostrowski.stp");
+
+vol_air() = v({10,11});
+vol_cond() = v({0,1,3:6,8,9});
+vol_diel() = v(7);
+vcore=v(2);
+
+// Doing intersections that are not in the initial step
+vol_aux()=BooleanFragments{Volume{vol_cond(),vol_diel()};Delete;}{}; //Plates and dielectric have to be intersected with connectors
+Coherence; // Air around
+
+// define physical groups
+Physical Volume("DomainAir", AIR) = vol_air();
+Physical Volume("DomainConductor",  CONDUCTOR) = vol_cond(); // 6 and 8 are the capacitor plates
+Physical Volume("DomainDielectric", DIELECTRIC) = vol_diel();
+Physical Volume("DomainCore", CORE) = vcore;
+
+// outside boundaries
+bnd_air() = Abs(CombinedBoundary{Volume{Volume{:}};});
+surf_elec_in = 130; //side L
+surf_elec_out = 144; // side C
+bnd_air() -= {surf_elec_in,surf_elec_out};
+
+Physical Surface("Terminal",  ELEC_IN)  = surf_elec_in;
+Physical Surface("Ground", ELEC_OUT) = surf_elec_out;
+Physical Surface("SurBCnotTerminalAir", BND_AIR) = bnd_air();
+
+
+// Adapting the mesh
+// get bounding box
+bb() = BoundingBox Volume {v()};
+
+// size along axes (useful for adapting mesh)
+sizeX = bb(3)-bb(0);
+sizeY = bb(4)-bb(1);
+sizeZ = bb(5)-bb(2);
+
+
+scale_factor=2;
+lc_ref = scale_factor*sizeX/200;
+ptos_cond() = PointsOf{Volume{vol_cond()};};
+ptos_air_in() = PointsOf{Volume{v(10)};};
+ptos_core() = PointsOf{Volume{vcore};};
+ptos_air_in() -= {ptos_cond(),ptos_core()};
+ptos_air_out() = PointsOf{Surface{152,153};}; // only the side
+
+MeshSize{ptos_air_out()} = 10*lc_ref;
+MeshSize{ptos_air_in()} = 3*lc_ref;
+MeshSize{ptos_core()}   = 1.5*lc_ref;
+MeshSize{ptos_cond()}   = lc_ref;
+
+// For aestetics
+Recursive Color SkyBlue { Volume{vol_air()};}
+Recursive Color Green  { Volume{vcore};}
+Recursive Color Red  { Volume{vol_cond()};}
+Recursive Color Yellow  { Volume{vol_diel()};}

ECE_full_wave/LC.pro

+/*  LC.pro
+Linked to step file from Joerg Ostrowski
+FW in E and V, ECE
+*/
+
+Include "LC_data.pro"
+
+Group{
+
+  Ground = Region[{ELEC_OUT}];
+  Terminal = Region[{ELEC_IN}];
+  SurBCnotTerminalAir = Region[{BND_AIR}];
+
+  DomainConductor = Region[{CONDUCTOR}];
+  DomainDielectric = Region[{DIELECTRIC}];
+  DomainCore = Region[{CORE}];
+  DomainAir = Region[{AIR}];
+
+  Vol_FW = Region[{DomainConductor,DomainDielectric,DomainCore,DomainAir}];
+  Sur_Terminals_FWece = Region[{Ground, Terminal}]; // all terminals
+  Sur_FW = Region[{Sur_Terminals_FWece, SurBCnotTerminalAir}]; //all boundary
+  Dom_FW = Region[ {Vol_FW, Sur_FW} ];
+ }
+
+
+ Function{
+   /*
+   Material properties are defined piecewise, for elementary groups
+   the values in the right hand side are defined in CL_data.pro
+   */
+
+   // conductor
+   nu[#{DomainConductor,Ground,Terminal}] = nuCond;             // magnetic reluctivity [m/H]
+   epsilon[#{DomainConductor,Ground,Terminal}] = epsCond;       // permittivity [F/m]
+   sigma[#{DomainConductor,Ground,Terminal}]  = sigmaCond;      // conductivity [S/m]
+
+   // dielectric
+   nu[#{DomainDielectric}] = nuDiel;
+   epsilon[#{DomainDielectric}] = epsDiel;
+   sigma[#{DomainDielectric}]  = sigmaDiel;
+
+   // core
+   nu[#{DomainCore}] = nuCore;
+   epsilon[#{DomainCore}] = epsCore;
+   sigma[#{DomainCore}]  = sigmaCore;
+
+   // air
+   nu[#{DomainAir,SurBCnotTerminalAir}] = nuAir;
+   epsilon[#{DomainAir,SurBCnotTerminalAir}] = epsAir;
+   sigma[#{DomainAir,SurBCnotTerminalAir}]  = sigmaAir;
+
+   If(Flag_AnalysisType==0)
+     VTerminal1[] = -VTermRMS*Complex[Cos[VTermPhase],Sin[VTermPhase]];
+   Else
+     ITerminal1[] = -ITermRMS*Complex[Cos[ITermPhase],Sin[ITermPhase]];
+   EndIf
+ }
+
+ Constraint{
+
+   // ece BC
+   { Name SetTerminalPotential; Type Assign;  // voltage excited terminals
+     Case {
+       { Region Ground; Value 0.; }
+       If(Flag_AnalysisType==0)
+         { Region Terminal; Value VTerminal1[]; }
+       EndIf
+     }
+   }
+   { 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
+       EndIf
+     }
+   }
+
+}
+
+
+If (Flag_AnalysisType==0)
+  Dir = "res/FWeceBC_voltExc/";
+Else
+  Dir = "res/FWeceBC_crtExc/";
+EndIf
+
+// The formulation and its tools
+Include "FullWave_E_ece_SISO.pro"
+
+// Output results
+Macro Change_post_options
+Echo[Str[ "nv = PostProcessing.NbViews;",
+    "For v In {0:nv-1}",
+    "View[v].NbIso = 25;",
+    "View[v].RangeType = 3;" ,// per timestep
+    "View[v].ShowTime = 3;",
+    "View[v].IntervalsType = 3;",
+    "EndFor"
+  ], File "post.opt"];
+Return
+
+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
+
+
+Macro Plugin_cuts
+  Echo[ Str[
+      "nview = PostProcessing.NbViews-1;",
+      "vname = Str[View[nview].Name];",
+      "Plugin(CutPlane).A        = 1;",
+      "Plugin(CutPlane).B        = 0;",
+      "Plugin(CutPlane).C        = 0;",
+      "Plugin(CutPlane).D        = -0.0289;", // half wire L side
+      "Plugin(CutPlane).View     = nview;",
+      "Plugin(CutPlane).Run;",
+      "View[nview+1].RangeType = 2;", // custom
+      "View[nview+1].CustomMax = View[nview+1].Max;", // max
+      "View[nview+1].CustomMin = View[nview+1].MinVisible;", // min
+      "View[nview+1].ScaleType = 2;", // logarithmic scale
+      "View[nview+1].Name = StrCat[vname, ' on YZ plane'];"
+    ], File StrCat[Dir,"cuts.geo"] ] ;
+Return
+
+PostOperation {
+
+  { Name Maps; NameOfPostProcessing FW_E_ece ;
+    Operation {
+      Print [ gradV, OnElementsOf ElementsOf[Vol_FW,OnOneSideOf Sur_FW], File StrCat[Dir,"map_gradV.pos" ]];
+
+      Print [ Vterminals, OnElementsOf ElementsOf[Sur_FW], File StrCat[Dir,"map_Vterminals.pos" ]];
+      Print [ VnotTerminals, OnElementsOf Vol_FW, File StrCat[Dir,"map_VnotTerminals.pos" ]];
+      Print [ H, OnElementsOf Vol_FW, File StrCat[Dir,"map_H.pos" ]];
+      Print [ E, OnElementsOf Vol_FW, File StrCat[Dir,"map_E.pos" ]];
+      Call Change_post_options;
+
+      Print [ normE, OnElementsOf Vol_FW, File StrCat[Dir,"map_normE.pos" ]];
+      Call Plugin_cuts;
+    }
+  }
+
+  { 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
+    }
+  }
+
+}
+
+// What follows only works if your solver is GetDP => adapt name if necessary
+
+DefineConstant[
+  // always solve this resolution (it's the only one provided)
+  // R_ = {"Analysis", Name "GetDP/1ResolutionChoices", Visible 1}
+  // set some command-line options for getdp
+  C_ = {"-solve -v 4 -pos -v2", Name "GetDP/9ComputeCommand", Visible 1}
+  // don't do the post-processing pass (some already included in Resolution)
+  P_ = {"Maps", Name "GetDP/2PostOperationChoices", Visible 1}
+];

ECE_full_wave/LC_data.pro

+DefineConstant[
+  cm = 1e-2,
+  mm = 1e-3
+  rw  = {5*mm/2, Name "Param. Geometry/03radius or wire", Highlight "Pink", ReadOnly}
+  rplateC = {6*cm, Name "Param. Geometry/10radius of capacitor plates", Highlight "LightCyan", ReadOnly}
+  dplateC = {3*mm, Name "Param. Geometry/11distance between plates", Highlight "LightCyan", ReadOnly}
+];
+
+sigma_cond = 1e7;
+mur_core = 1000;
+epsr_diel = 1e6;
+sigma_diel = 0;
+
+mu0 = 4.e-7 * Pi ;
+nu0  = 1/mu0;
+eps0 = 8.854187818e-12;
+
+// --------------------
+
+nuAir = nu0;
+epsAir = eps0;
+sigmaAir = 0;
+
+nuCore = 1/(mu0*mur_core);
+epsCore = eps0;
+sigmaCore = 0;
+
+nuDiel = nu0;
+epsDiel = eps0*epsr_diel;
+sigmaDiel = sigma_diel;
+
+nuCond = nu0;
+epsCond = eps0;
+sigmaCond = sigma_cond;
+
+
+//----------------------
+// Physical numbers...
+//----------------------
+CONDUCTOR     = 1000;  // volume, conductor touching the boundary at the excited termimnal, including the coil wire and one plate
+BND_CONDUCTOR = 1111;  // volume, conductor touching the boundary at the grounded terminal and including the other plate
+
+ELEC_IN  = 111; // surface, excited terminal
+ELEC_OUT = 112; // surface, grounded terminal
+
+DIELECTRIC   = 2000; // volume, the dielectric of the capacitor
+ELEC_DIEL_IN = 221;
+ELEC_DIEL_OUT = 222;
+
+CORE = 3000;  // volume, the core of the coil
+
+AIR = 4000;  // volume, the air
+BND_AIR = 4444; // surface, boundary touched by air
+
+// -----------------------
+cGUI= "{CL - Excitation/";
+close_menu = 1;
+colorro = "LightGrey";
+mTypeBC = StrCat[cGUI,"Type of ECE Excitation/0"];
+
+colorMTypeBC = "Red";
+DefineConstant[
+  Flag_AnalysisType = { 0,
+    Choices{
+      0="ev",
+      1="ec"},
+    Name Str[mTypeBC], Highlight Str[colorMTypeBC], Closed !close_menu,
+    ServerAction Str["Reset", "GetDP/1ResolutionChoices"]}
+];
+
+/* new menu for the interface - values of excitations - significance depends on the type of BC */
+mValuesBC = StrCat[cGUI,"Frequency values/0"];
+colorMValuesBC = "Ivory";
+
+fmin = 1000; // Hz
+fmax = 80000; // Hz
+
+nop = 10;
+freqs()=LogSpace[Log10[fmin],Log10[fmax],nop];
+//freqs()=LinSpace[fmin,fmax,nop];
+DefineConstant[
+  Freq  = {freqs(0), Choices{freqs()}, Loop 0, Name StrCat[mValuesBC, "0Working Frequency"], Units "Hz", Highlight  Str[colorMValuesBC], Closed !close_menu }
+];
+
+
+If (Flag_AnalysisType==1) // current exc
+  Printf("======================> Terminal is current => Z transfer function");
+  ITermRMS = 1;
+  ITermPhase = 0;
+Else  // voltage exc
+  Printf("======================> Terminal is voltage excited => Y transfer function");
+  VTermRMS = 1;
+  VTermPhase = 0;
+EndIf

ECE_full_wave/README.txt

+Electric circuit element  (ECE) boundary conditions in full-wave electromagnetic problems.
+
+Models developed by C. Ciuprina and R.V. Sabariego.
+
+This directory contains the pro and geo files necessary to run the test cases in article:
+
+G. Ciuprina, D. Ioan, and R. V. Sabariego. Electric circuit element boundary conditions
+in the finite element method for full-wave passive electromagnetic devices.
+Journal of Mathematics in Industry, 12:1-13, Art no. 7, 2022. DOI: 10.1186/s13362-022-00122-1.
+
+Common formulation files:
+-------------------------
+  "FullWave_E_ece_SISO.pro" - file with function spaces, formulation, resolution and postprocessing.
+
+Cylinder benchmark:
+-------------------
+  "Ishape3D_data.pro" - contains geometrical and physical parameters used by the geo and pro files.
+  "Ishape3D.geo" - contains the geometry description.
+  "Ishape3D.pro" - contains the problem definition.
+
+Inductor-Capacitor (LC) benchmark:
+-------------
+  "LC_data.pro" - contains geometrical and physical parameters used by the geo and pro files.
+  "LC_ostrowski.stp" - STEP geometry courtesy of J. Ostrowski.
+  "LC.geo" - adapts STEP file to use with Gmsh. Note that the computation of BooleanFragments is very slow.
+  "LC.pro" - contains the problem definition.
+
+Quick start
+-----------
+  Open "Ishape3D.pro" with Gmsh.
+  Open "LC.pro" with Gmsh.