first pass at tutorial 7

Magnetodynamics/Lib_MagStaDyn_av_2D_Cir.pro

+// This is a template .pro file containing a general formulation for 2D
+// magnetostatic and magnetodynamic problems in terms of the magnetic vector
+// potential a (potentially coupled with the electric scalar potential v), with
+// optional circuit coupling.
+
+// Below are definitions of the constants (inside "DefineConstant"), groups
+// (inside "DefineGroup") and functions (inside "DefineFunction") that can be
+// redefined from outside this template.
+
+DefineConstant[
+  Flag_FrequencyDomain = 1, // frequency-domain or time-domain simulation
+  Flag_CircuitCoupling = 0 // consider coupling with external electric circuit
+  CoefPower = 0.5, // coefficient for power calculations
+  Freq = 50, // frequency (for harmonic simulations)
+  TimeInit = 0, // intial time (for time-domain simulations)
+  TimeFinal = 1, // final time (for time-domain simulations)
+  DeltaTime = 0.01, // time step (for time-domain simulations)
+  InterpolationOrder = 1 // finite element order
+  Val_Rint = 0, // interior radius of annulus shell transformation region (VolInf_Mag)
+  Val_Rext = 0 // exterior radius of annulus shell  transformation region (VolInf_Mag)
+];
+
+Group {
+  DefineGroup[
+    VolCC_Mag, // the non-conducting part
+    VolC_Mag, // the conducting part
+    VolV_Mag, // a moving conducting part, with invariant mesh
+    VolM_Mag, // permanent magnets
+    VolS0_Mag, // current source domain with imposed current densities js0
+    VolS_Mag, // current source domain with imposed current, imposed voltage or
+              // circuit coupling
+    VolInf_Mag, // annulus where a infinite shell transformation is applied
+    SurFluxTube_Mag, // boundary with Neumann BC
+    SurPerfect_Mag, // boundary of perfect conductors (i.e. non-meshed)
+    SurImped_Mag // boundary of conductors approximated by a surface impedance
+                 // (i.e. non-meshed)
+  ];
+  If(Flag_CircuitCoupling)
+    DefineGroup[
+      SourceV_Cir, // voltage sources
+      SourceI_Cir, // current sources
+      Resistance_Cir, // resistances (linear)
+      Inductance_Cir, // inductances
+      Capacitance_Cir, // capacitances
+      Diode_Cir // diodes (nonlinear resistances)
+    ];
+  EndIf
+}
+
+Function {
+  DefineFunction[
+    nu, // reluctivity (in Vol_Mag)
+    sigma, // conductivity (in VolC_Mag and VolS_Mag)
+    br, // remanent magnetic flux density (in VolM_Mag)
+    js0, // source current density (in VolS0_Mag)
+    nxh, // n x magnetic field (on SurFluxTube_Mag)
+    Velocity, // velocity of moving part VolV_Moving
+    Ns, // number of turns (in VolS_Mag)
+    Sc, // cross-section of windings (in VolS_Mag)
+    CoefGeos, // geometrical coefficient for 2D or 2D axi model
+    Ysur // surface admittance (inverse of surface impedance Zsur) on SurImped_Mag
+  ];
+  If(Flag_CircuitCoupling)
+    DefineFunction[
+      Resistance, // resistance values
+      Inductance, // inductance values
+      Capacitance // capacitance values
+    ];
+  EndIf
+}
+
+// End of definitions.
+
+Group{
+  // all volumes
+  Vol_Mag = Region[ {VolCC_Mag, VolC_Mag, VolV_Mag, VolM_Mag, VolS0_Mag, VolS_Mag} ];
+  // all volumes + surfaces on which integrals will be computed
+  VolWithSur_Mag = Region[ {Vol_Mag, SurFluxTube_Mag, SurPerfect_Mag, SurImped_Mag} ];
+  If(Flag_CircuitCoupling)
+    // all circuit impedances
+    DomainZ_Cir = Region[ {Resistance_Cir, Inductance_Cir, Capacitance_Cir} ];
+    // all circuit sources
+    DomainSource_Cir = Region[ {SourceV_Cir, SourceI_Cir} ];
+    // all circuit elements
+    Domain_Cir = Region[ {DomainZ_Cir, DomainSource_Cir} ];
+  EndIf
+}
+
+Jacobian {
+  { Name Vol;
+    Case {
+      { Region VolInf_Mag ;
+        Jacobian VolSphShell {Val_Rint, Val_Rext} ; }
+      { Region All; Jacobian Vol; }
+    }
+  }
+  { Name Sur;
+    Case {
+      { Region All; Jacobian Sur; }
+    }
+  }
+}
+
+Integration {
+  { Name Gauss_v;
+    Case {
+      { Type Gauss;
+        Case {
+          { GeoElement Point; NumberOfPoints  1; }
+          { GeoElement Line; NumberOfPoints  5; }
+          { GeoElement Triangle; NumberOfPoints  7; }
+          { GeoElement Quadrangle; NumberOfPoints  4; }
+          { GeoElement Tetrahedron; NumberOfPoints 15; }
+          { GeoElement Hexahedron; NumberOfPoints 14; }
+          { GeoElement Prism; NumberOfPoints 21; }
+        }
+      }
+    }
+  }
+}
+
+// Same FunctionSpace for both static and dynamic formulations
+FunctionSpace {
+  { Name Hcurl_a_2D; Type Form1P; // 1-form (circulations) on edges
+                                  // perpendicular to the plane of study
+    BasisFunction {
+      // \vec{a}(x) = \sum_{n \in N(Domain)} a_n \vec{s}_n(x)
+      //   without nodes on perfect conductors (where a is constant)
+      { Name s_n; NameOfCoef a_n; Function BF_PerpendicularEdge;
+        Support VolWithSur_Mag; Entity NodesOf[All, Not SurPerfect_Mag]; }
+
+      // global basis function on boundary of perfect conductors
+      { Name s_skin; NameOfCoef a_skin; Function BF_GroupOfPerpendicularEdges;
+        Support VolWithSur_Mag; Entity GroupsOfNodesOf[SurPerfect_Mag]; }
+
+      // additional basis functions for 2nd order interpolation
+      If(InterpolationOrder == 2)
+        { Name s_e; NameOfCoef a_e; Function BF_PerpendicularEdge_2E;
+          Support Vol_Mag; Entity EdgesOf[All]; }
+      EndIf
+    }
+    GlobalQuantity {
+      { Name A; Type AliasOf; NameOfCoef a_skin; }
+      { Name I; Type AssociatedWith; NameOfCoef a_skin; }
+    }
+    Constraint {
+      { NameOfCoef a_n;
+        EntityType NodesOf; NameOfConstraint MagneticVectorPotential_2D; }
+
+      { NameOfCoef I;
+        EntityType GroupsOfNodesOf; NameOfConstraint Current_2D; }
+
+      If(InterpolationOrder == 2)
+        { NameOfCoef a_e;
+          EntityType EdgesOf; NameOfConstraint MagneticVectorPotential_2D_0; }
+      EndIf
+    }
+  }
+}
+
+FunctionSpace {
+  // Gradient of Electric scalar potential (2D)
+  { Name Hregion_u_2D; Type Form1P; // same as for \vec{a}
+    BasisFunction {
+      { Name sr; NameOfCoef ur; Function BF_RegionZ;
+        // constant vector (over the region) with nonzero z-component only
+        Support Region[{VolC_Mag, SurImped_Mag}];
+        Entity Region[{VolC_Mag, SurImped_Mag}]; }
+    }
+    GlobalQuantity {
+      { Name U; Type AliasOf; NameOfCoef ur; }
+      { Name I; Type AssociatedWith; NameOfCoef ur; }
+    }
+    Constraint {
+      { NameOfCoef U;
+        EntityType Region; NameOfConstraint Voltage_2D; }
+      { NameOfCoef I;
+        EntityType Region; NameOfConstraint Current_2D; }
+    }
+  }
+
+  // Current in stranded coil (2D)
+  { Name Hregion_i_2D; Type Vector;
+    BasisFunction {
+      { Name sr; NameOfCoef ir; Function BF_RegionZ;
+        Support VolS_Mag; Entity VolS_Mag; }
+    }
+    GlobalQuantity {
+      { Name Is; Type AliasOf; NameOfCoef ir; }
+      { Name Us; Type AssociatedWith; NameOfCoef ir; }
+    }
+    Constraint {
+      { NameOfCoef Us;
+        EntityType Region; NameOfConstraint Voltage_2D; }
+      { NameOfCoef Is;
+        EntityType Region; NameOfConstraint Current_2D; }
+    }
+  }
+}
+
+If(Flag_CircuitCoupling)
+  // UZ and IZ for impedances
+  FunctionSpace {
+    { Name Hregion_Z; Type Scalar;
+      BasisFunction {
+        { Name sr; NameOfCoef ir; Function BF_Region;
+          Support Domain_Cir; Entity Domain_Cir; }
+      }
+      GlobalQuantity {
+        { Name Iz; Type AliasOf; NameOfCoef ir; }
+        { Name Uz; Type AssociatedWith; NameOfCoef ir; }
+      }
+      Constraint {
+        { NameOfCoef Uz;
+          EntityType Region; NameOfConstraint Voltage_Cir; }
+        { NameOfCoef Iz;
+          EntityType Region; NameOfConstraint Current_Cir; }
+      }
+    }
+  }
+EndIf
+
+
+// Static Formulation
+Formulation {
+  { Name MagSta_a_2D; Type FemEquation;
+    Quantity {
+      { Name a; Type Local; NameOfSpace Hcurl_a_2D; }
+      { Name ir; Type Local; NameOfSpace Hregion_i_2D; }
+    }
+    Equation {
+      Galerkin { [ nu[] * Dof{d a} , {d a} ];
+        In Vol_Mag; Jacobian Vol; Integration Gauss_v; }
+
+      Galerkin { [ -nu[] * br[] , {d a} ];
+        In VolM_Mag; Jacobian Vol; Integration Gauss_v; }
+
+      Galerkin { [ -js0[] , {a} ];
+        In VolS0_Mag; Jacobian Vol; Integration Gauss_v; }
+
+      Galerkin { [ - (js0[]*Vector[0,0,1]) * Dof{ir} , {a} ];
+        In VolS_Mag; Jacobian Vol; Integration Gauss_v; }
+
+      Galerkin { [ nxh[] , {a} ];
+        In SurFluxTube_Mag; Jacobian Sur; Integration Gauss_v; }
+    }
+  }
+}
+
+// Dynamic Formulation (eddy currents)
+Formulation {
+  { Name MagDyn_a_2D; Type FemEquation;
+    Quantity {
+      { Name a; Type Local; NameOfSpace Hcurl_a_2D; }
+      { Name A_floating; Type Global; NameOfSpace Hcurl_a_2D [A]; }
+      { Name I_perfect; Type Global; NameOfSpace Hcurl_a_2D [I]; }
+
+      { Name ur; Type Local; NameOfSpace Hregion_u_2D; }
+      { Name I; Type Global; NameOfSpace Hregion_u_2D [I]; }
+      { Name U; Type Global; NameOfSpace Hregion_u_2D [U]; }
+
+      { Name ir; Type Local; NameOfSpace Hregion_i_2D; }
+      { Name Us; Type Global; NameOfSpace Hregion_i_2D [Us]; }
+      { Name Is; Type Global; NameOfSpace Hregion_i_2D [Is]; }
+
+      If(Flag_CircuitCoupling)
+        { Name Uz; Type Global; NameOfSpace Hregion_Z [Uz]; }
+        { Name Iz; Type Global; NameOfSpace Hregion_Z [Iz]; }
+      EndIf
+    }
+    Equation {
+      Galerkin { [ nu[] * Dof{d a} , {d a} ];
+        In Vol_Mag; Jacobian Vol; Integration Gauss_v; }
+      Galerkin { [ -nu[] * br[] , {d a} ];
+        In VolM_Mag; Jacobian Vol; Integration Gauss_v; }
+
+      // Electric field e = -Dt[{a}]-{ur},
+      // with {ur} = Grad v constant in each region of VolC
+      Galerkin { DtDof [ sigma[] * Dof{a} , {a} ];
+        In VolC_Mag; Jacobian Vol; Integration Gauss_v; }
+      Galerkin { [ sigma[] * Dof{ur} / CoefGeos[] , {a} ];
+        In VolC_Mag; Jacobian Vol; Integration Gauss_v; }
+
+      Galerkin { [ - sigma[] * (Velocity[] /\ Dof{d a}) , {a} ];
+        In VolV_Mag; Jacobian Vol; Integration Gauss_v; }
+
+      Galerkin { [ -js0[] , {a} ];
+        In VolS0_Mag; Jacobian Vol; Integration Gauss_v; }
+
+      Galerkin { [ nxh[] , {a} ];
+        In SurFluxTube_Mag; Jacobian Sur; Integration Gauss_v; }
+
+      Galerkin { DtDof [  Ysur[] * Dof{a} , {a} ];
+        In SurImped_Mag; Jacobian Sur; Integration Gauss_v; }
+      Galerkin { [ Ysur[] * Dof{ur} / CoefGeos[] , {a} ];
+        In SurImped_Mag; Jacobian Sur; Integration Gauss_v; }
+
+      // When {ur} act as a test function, one obtains the circuits relations,
+      // relating the voltage and the current of each region in VolC
+      Galerkin { DtDof [ sigma[] * Dof{a} , {ur} ];
+        In VolC_Mag; Jacobian Vol; Integration Gauss_v; }
+      Galerkin { [ sigma[] * Dof{ur} / CoefGeos[] , {ur} ];
+        In VolC_Mag; Jacobian Vol; Integration Gauss_v; }
+      GlobalTerm { [ Dof{I} *(CoefGeos[]/Fabs[CoefGeos[]]) , {U} ]; In VolC_Mag; }
+
+      Galerkin { DtDof [ Ysur[] * Dof{a} , {ur} ];
+        In SurImped_Mag; Jacobian Sur; Integration Gauss_v; }
+      Galerkin { [ Ysur[] * Dof{ur} / CoefGeos[] , {ur} ];
+        In SurImped_Mag; Jacobian Sur; Integration Gauss_v; }
+      GlobalTerm { [ Dof{I} *(CoefGeos[]/Fabs[CoefGeos[]]) , {U} ]; In SurImped_Mag; }
+
+      // js[0] should be of the form: Ns[]/Sc[] * Vector[0,0,1]
+      Galerkin { [ - (js0[]*Vector[0,0,1]) * Dof{ir} , {a} ];
+        In VolS_Mag; Jacobian Vol; Integration Gauss_v; }
+
+      Galerkin { DtDof [ Ns[]/Sc[] * Dof{a} , {ir} ];
+        In VolS_Mag; Jacobian Vol; Integration Gauss_v; }
+      Galerkin { [ Ns[]/Sc[] / sigma[] * (js0[]*Vector[0,0,1]) * Dof{ir} , {ir} ];
+        In VolS_Mag; Jacobian Vol; Integration Gauss_v; }
+      GlobalTerm { [ Dof{Us} / CoefGeos[] , {Is} ]; In VolS_Mag; }
+      // Attention: CoefGeo[.] = 2*Pi for Axi
+
+      GlobalTerm { [ -Dof{I_perfect} , {A_floating} ]; In SurPerfect_Mag; }
+
+      If(Flag_CircuitCoupling)
+	GlobalTerm { NeverDt[ Dof{Uz} , {Iz} ]; In Resistance_Cir; }
+        GlobalTerm { NeverDt[ Resistance[] * Dof{Iz} , {Iz} ]; In Resistance_Cir; }
+
+	GlobalTerm { [ Dof{Uz} , {Iz} ]; In Inductance_Cir; }
+	GlobalTerm { DtDof [ Inductance[] * Dof{Iz} , {Iz} ]; In Inductance_Cir; }
+
+	GlobalTerm { NeverDt[ Dof{Iz} , {Iz} ]; In Capacitance_Cir; }
+	GlobalTerm { DtDof [ Capacitance[] * Dof{Uz} , {Iz} ]; In Capacitance_Cir; }
+
+	GlobalTerm { NeverDt[ Dof{Uz} , {Iz} ]; In Diode_Cir; }
+	GlobalTerm { NeverDt[ Resistance[{Uz}] * Dof{Iz} , {Iz} ]; In Diode_Cir; }
+
+	GlobalTerm { [ 0. * Dof{Iz} , {Iz} ]; In DomainSource_Cir; }
+
+	GlobalEquation {
+	  Type Network; NameOfConstraint ElectricalCircuit;
+	  { Node {I};  Loop {U};  Equation {I};  In VolC_Mag; }
+	  { Node {Is}; Loop {Us}; Equation {Us}; In VolS_Mag; }
+	  { Node {Iz}; Loop {Uz}; Equation {Uz}; In Domain_Cir; }
+	}
+      EndIf
+
+    }
+  }
+}
+
+Resolution {
+  { Name MagDyn_a_2D;
+    System {
+      { Name Sys; NameOfFormulation MagDyn_a_2D;
+        If(Flag_FrequencyDomain)
+          Type ComplexValue; Frequency Freq;
+        EndIf
+      }
+    }
+    Operation {
+      If(Flag_FrequencyDomain)
+        Generate[Sys]; Solve[Sys]; SaveSolution[Sys];
+      Else
+        InitSolution[Sys]; // provide initial condition
+        TimeLoopTheta[Time0, TimeMax, DeltaTime, 1.]{
+          // Euler implicit (1) -- Crank-Nicolson (0.5)
+          Generate[Sys]; Solve[Sys]; SaveSolution[Sys];
+        }
+      EndIf
+    }
+  }
+  { Name MagSta_a_2D;
+    System {
+      { Name Sys; NameOfFormulation MagSta_a_2D; }
+    }
+    Operation {
+      Generate[Sys]; Solve[Sys]; SaveSolution[Sys];
+    }
+  }
+}
+
+// Same PostProcessing for both static and dynamic formulations (both refer to
+// the same FunctionSpace from which the solution is obtained)
+PostProcessing {
+  { Name MagDyn_a_2D; NameOfFormulation MagDyn_a_2D;
+    PostQuantity {
+      // In 2D, a is a vector with only a z-component: (0,0,az)
+      { Name a; Value { Term { [ {a} ]; In Vol_Mag; Jacobian Vol; } } }
+      // The equilines of az are field lines (giving the magnetic field direction)
+      { Name az; Value { Term { [ CompZ[{a}] ]; In Vol_Mag; Jacobian Vol; } } }
+      { Name b; Value { Term { [ {d a} ]; In Vol_Mag; Jacobian Vol; } } }
+      { Name h; Value {
+          Term { [ nu[] * {d a} ]; In Vol_Mag; Jacobian Vol; }
+          Term { [ -nu[] * br[] ]; In VolM_Mag; Jacobian Vol; }
+        }
+      }
+      { Name js; Value {
+          Term { [ js0[] ]; In VolS0_Mag; Jacobian Vol; }
+          Term { [  (js0[]*Vector[0,0,1])*{ir} ]; In VolS_Mag; Jacobian Vol; }
+          Term { [ Vector[0,0,0] ]; In Vol_Mag; Jacobian Vol; } // to force a vector result out of sources
+        }
+      }
+      { Name j;
+        // Only correct in MagDyn
+        Value {
+          Term { [ -sigma[] * (Dt[{a}]+{ur}/CoefGeos[]) ]; In VolC_Mag; Jacobian Vol; }
+          Term { [ js0[] ]; In VolS0_Mag; Jacobian Vol; }
+          Term { [  (js0[]*Vector[0,0,1])*{ir} ]; In VolS_Mag; Jacobian Vol; }
+          Term { [ Vector[0,0,0] ]; In Vol_Mag; Jacobian Vol; }
+          // Current density in A/m
+          Term { [ -Ysur[] * (Dt[{a}]+{ur}/CoefGeos[]) ]; In SurImped_Mag; Jacobian Sur; }
+        }
+      }
+
+      { Name JouleLosses;
+        Value {
+          Integral { [ CoefPower * sigma[]*SquNorm[Dt[{a}]+{ur}/CoefGeos[]] ];
+            In VolC_Mag; Jacobian Vol; Integration Gauss_v; }
+          Integral { [ CoefPower * 1./sigma[]*SquNorm[js0[]] ];
+            In VolS0_Mag; Jacobian Vol; Integration Gauss_v; }
+
+	  Integral { [ CoefPower * 1./sigma[]*SquNorm[(js0[]*Vector[0,0,1])*{ir}] ];
+            In VolS_Mag; Jacobian Vol; Integration Gauss_v; }
+
+          Integral { [ CoefPower * Ysur[]*SquNorm[Dt[{a}]+{ur}/CoefGeos[]] ];
+            In SurImped_Mag; Jacobian Sur; Integration Gauss_v; }
+	}
+      }
+
+      { Name U; Value {
+          Term { [ {U} ]; In VolC_Mag; }
+          Term { [ {Us} ]; In VolS_Mag; }
+          If(Flag_CircuitCoupling)
+            Term { [ {Uz} ]; In Domain_Cir; }
+          EndIf
+        }
+      }
+
+      { Name I; Value {
+          Term { [ {I} ]; In VolC_Mag; }
+          Term { [ {Is} ]; In VolS_Mag; }
+          If(Flag_CircuitCoupling)
+            Term { [ {Iz} ]; In Domain_Cir; }
+          EndIf
+        }
+      }
+
+    }
+  }
+
+  { Name MagSta_a_2D; NameOfFormulation MagSta_a_2D;
+    PostQuantity {
+      // In 2D, a is a vector with only a z-component: (0,0,az)
+      { Name a; Value { Term { [ {a} ]; In Vol_Mag; Jacobian Vol; } } }
+      // The equilines of az are field lines (giving the magnetic field direction)
+      { Name az; Value { Term { [ CompZ[{a}] ]; In Vol_Mag; Jacobian Vol; } } }
+      { Name b; Value { Term { [ {d a} ]; In Vol_Mag; Jacobian Vol; } } }
+      { Name h; Value { Term { [ nu[] * {d a} ]; In Vol_Mag; Jacobian Vol; } } }
+      { Name j;
+        Value {
+          Term { [ js0[] ]; In VolS0_Mag; Jacobian Vol; }
+          Term { [ Vector[0,0,0] ]; In Vol_Mag; Jacobian Vol; }
+        }
+      }
+    }
+  }
+}

Magnetodynamics/electromagnet.geo

+/* -------------------------------------------------------------------
+   File "electromagnet.geo"
+
+   This file is the geometrical description used by GMSH to produce
+   the file "electromagnet.msh".
+   ------------------------------------------------------------------- */
+
+dxCore =  50.e-3; dyCore = 100.e-3;
+xInd   =  75.e-3; dxInd  =  25.e-3; dyInd  =  50.e-3;
+
+Include "electromagnet_common.pro";
+
+s       = DefineNumber[1, Name "Model parameters/Global mesh size",
+		       Help "Reduce for finer mesh"]; 
+p0      = 12.e-3 *s;
+pCorex  =  4.e-3 *s; pCorey0 =  8.e-3 *s; pCorey  =  4.e-3 *s;
+pIndx   =  5.e-3 *s; pIndy   =  5.e-3 *s;
+pInt    = 12.5e-3*s; pExt    = 12.5e-3*s;
+
+Point(1) = {0,0,0,p0};
+Point(2) = {dxCore,0,0,pCorex};
+Point(3) = {dxCore,dyCore,0,pCorey};
+Point(4) = {0,dyCore,0,pCorey0};
+Point(5) = {xInd,0,0,pIndx};
+Point(6) = {xInd+dxInd,0,0,pIndx};
+Point(7) = {xInd+dxInd,dyInd,0,pIndy};
+Point(8) = {xInd,dyInd,0,pIndy};
+Point(9) = {rInt,0,0,pInt};
+Point(10) = {rExt,0,0,pExt};
+Point(11) = {0,rInt,0,pInt};
+Point(12) = {0,rExt,0,pExt};
+
+Line(1) = {1,2};  Line(2) = {2,5};   Line(3) = {5,6};
+Line(4) = {6,9};  Line(5) = {9,10};  Line(6) = {1,4};
+Line(7) = {4,11}; Line(8) = {11,12}; Line(9) = {2,3};
+Line(10) = {3,4}; Line(11) = {6,7};  Line(12) = {7,8};
+Line(13) = {8,5};
+
+Circle(14) = {9,1,11};  Circle(15) = {10,1,12};
+
+Line Loop(16) = {-6,1,9,10};                 Plane Surface(17) = {16};
+Line Loop(18) = {11,12,13,3};                Plane Surface(19) = {18};
+Line Loop(20) = {7,-14,-4,11,12,13,-2,9,10}; Plane Surface(21) = {-20}; // "-" for orientation
+Line Loop(22) = {8,-15,-5,14};               Plane Surface(23) = {-22};
+
+Physical Surface(101) = {21};  /* Air */
+Physical Surface(102) = {17};  /* Core */
+Physical Surface(103) = {19};  /* Ind */
+Physical Surface(111) = {23};  /* AirInf */
+
+Physical Line(1100) = {1,2,3,4,5}; /* Surface ht = 0 */
+Physical Line(1101) = {6,7,8};     /* Surface bn = 0 */
+Physical Line(1102) = {15};        /* Surface Inf */
+

Magnetodynamics/electromagnet.pro

+/* -------------------------------------------------------------------
+   Tutorial 7a : magnetostatic field of an electromagnet, bis
+
+   Features:
+   - Same as Tutorial 2, but using a generic template formulation library
+
+   To compute the solution in a terminal:
+       getdp electromagnet -solve MagSta_a_2D
+       getdp electromagnet -pos Map_a
+
+   To compute the solution interactively from the Gmsh GUI:
+       File > Open > electromagnet.pro
+       Run (button at the bottom of the left panel)
+   ------------------------------------------------------------------- */
+
+Include "electromagnet_common.pro";
+
+Group {
+  // Physical regions
+  Air    = Region[ 101 ];   Core   = Region[ 102 ];
+  Ind    = Region[ 103 ];   AirInf = Region[ 111 ];
+  Surface_ht0 = Region[ 1100 ];
+  Surface_bn0 = Region[ 1101 ];
+  Surface_Inf = Region[ 1102 ];
+
+  // Abstract regions used in the "Lib_MagStaDyn_av_2D_Cir.pro" template file
+  // that is included below:
+  VolCC_Mag = Region[{Air, Core, AirInf}]; // Non-conducting regions
+  VolS_Mag = Region[{Ind}]; // Stranded conductors, i.e., coils
+  VolInf_Mag = Region[{AirInf}]; // annulus for infinite shell transformation
+  Val_Rint = rInt; Val_Rext = rExt; // interior and exterior radii of annulus
+}
+
+Function {
+  DefineConstant[
+    murCore = {100, Name "Model parameters/Mur core"},
+    Current = {0.01, Name "Model parameters/Current"}
+  ];
+
+  mu0 = 4.e-7 * Pi;
+  nu[ Region[{Air, Ind, AirInf}] ]  = 1. / mu0;
+  nu[ Core ]  = 1. / (murCore * mu0);
+
+  Ns[ Ind ] = 1000 ; // number of turns in coil
+  Sc[ Ind ] = SurfaceArea[] ; // surface (cross section) of coil
+  // Current density in each coil portion for a unit current (will be multiplied
+  // by the actual total current in the coil)
+  js0[ Ind ] = Ns[]/Sc[] * Vector[0,0,-1];
+}
+
+Constraint {
+  { Name MagneticVectorPotential_2D;
+    Case {
+      { Region Surface_bn0; Value 0; }
+      { Region Surface_Inf; Value 0; }
+    }
+  }
+  { Name Current_2D;
+    Case {
+      { Region Ind; Value Current; }
+    }
+  }
+  { Name Voltage_2D;
+    Case {
+    }
+  }
+}
+
+Include "Lib_MagStaDyn_av_2D_Cir.pro";
+
+PostOperation {
+  { Name Map_a; NameOfPostProcessing MagSta_a_2D;
+    Operation {
+      Print[ a, OnElementsOf Vol_Mag, File "a.pos" ];
+    }
+  }
+}

Magnetodynamics/electromagnet_common.pro

+// Parameters shared by Gmsh and GetDP
+
+rInt   = 200.e-3; 
+rExt   = 250.e-3;

Magnetodynamics/transfo.geo

+
+Include "transfo_common.dat";
+
+Solver.AutoShowLastStep = 0; // don't automatically show the last time step
+
+// CORE
+
+p_Leg_1_L_0=newp; Point(newp) = {-width_Core/2., 0, 0, c_Core};
+p_Leg_1_R_0=newp; Point(newp) = {-width_Core/2.+width_Core_Leg, 0, 0, c_Core};
+
+p_Leg_1_L_1=newp; Point(newp) = {-width_Core/2., height_Core/2., 0, c_Core};
+p_Leg_1_R_1=newp; Point(newp) = {-width_Core/2.+width_Core_Leg, height_Core/2.-width_Core_Leg, 0, c_Core};
+
+
+p_Leg_2_L_0=newp; Point(newp) = {width_Core/2.-width_Core_Leg, 0, 0, c_Core};
+p_Leg_2_R_0=newp; Point(newp) = {width_Core/2., 0, 0, c_Core};
+
+p_Leg_2_L_1=newp; Point(newp) = {width_Core/2.-width_Core_Leg, height_Core/2.-width_Core_Leg, 0, c_Core};
+p_Leg_2_R_1=newp; Point(newp) = {width_Core/2., height_Core/2., 0, c_Core};
+
+
+l_Core_In[]={};
+l_Core_In[]+=newl; Line(newl) = {p_Leg_1_R_0, p_Leg_1_R_1};
+l_Core_In[]+=newl; Line(newl) = {p_Leg_1_R_1, p_Leg_2_L_1};
+l_Core_In[]+=newl; Line(newl) = {p_Leg_2_L_1, p_Leg_2_L_0};
+
+l_Core_Out[]={};
+l_Core_Out[]+=newl; Line(newl) = {p_Leg_2_R_0, p_Leg_2_R_1};
+l_Core_Out[]+=newl; Line(newl) = {p_Leg_2_R_1, p_Leg_1_L_1};
+l_Core_Out[]+=newl; Line(newl) = {p_Leg_1_L_1, p_Leg_1_L_0};
+
+l_Core_Leg_1_Y0[]={};
+l_Core_Leg_1_Y0[]+=newl; Line(newl) = {p_Leg_1_L_0, p_Leg_1_R_0};
+
+l_Core_Leg_2_Y0[]={};
+l_Core_Leg_2_Y0[]+=newl; Line(newl) = {p_Leg_2_L_0, p_Leg_2_R_0};
+
+
+ll_Core=newll; Line Loop(newll) = {l_Core_In[], l_Core_Leg_2_Y0[], l_Core_Out[], l_Core_Leg_1_Y0[]};
+s_Core=news; Plane Surface(news) = {ll_Core};
+
+Physical Surface("CORE", CORE) = {s_Core};
+
+// COIL_1_PLUS
+
+x_[]=Point{p_Leg_1_R_0};
+p_Coil_1_p_int_0=newp; Point(newp) = {x_[0]+gap_Core_Coil_1, 0, 0, c_Coil_1};
+p_Coil_1_p_ext_0=newp; Point(newp) = {x_[0]+gap_Core_Coil_1+width_Coil_1, 0, 0, c_Coil_1};
+p_Coil_1_p_int_1=newp; Point(newp) = {x_[0]+gap_Core_Coil_1, height_Coil_1/2, 0, c_Coil_1};
+p_Coil_1_p_ext_1=newp; Point(newp) = {x_[0]+gap_Core_Coil_1+width_Coil_1, height_Coil_1/2, 0, c_Coil_1};
+
+l_Coil_1_p[]={};
+l_Coil_1_p[]+=newl; Line(newl) = {p_Coil_1_p_ext_0, p_Coil_1_p_ext_1};
+l_Coil_1_p[]+=newl; Line(newl) = {p_Coil_1_p_ext_1, p_Coil_1_p_int_1};
+l_Coil_1_p[]+=newl; Line(newl) = {p_Coil_1_p_int_1, p_Coil_1_p_int_0};
+
+l_Coil_1_p_Y0[]={};
+l_Coil_1_p_Y0[]+=newl; Line(newl) = {p_Coil_1_p_int_0, p_Coil_1_p_ext_0};
+
+ll_Coil_1_p=newll; Line Loop(newll) = {l_Coil_1_p[], l_Coil_1_p_Y0[]};
+s_Coil_1_p=news; Plane Surface(news) = {ll_Coil_1_p};
+
+Physical Surface("COIL_1_PLUS", COIL_1_PLUS) = {s_Coil_1_p};
+
+
+// COIL_1_MINUS
+
+x_[]=Point{p_Leg_1_L_0};
+p_Coil_1_m_int_0=newp; Point(newp) = {x_[0]-gap_Core_Coil_1, 0, 0, c_Coil_1};
+p_Coil_1_m_ext_0=newp; Point(newp) = {x_[0]-(gap_Core_Coil_1+width_Coil_1), 0, 0, c_Coil_1};
+p_Coil_1_m_int_1=newp; Point(newp) = {x_[0]-gap_Core_Coil_1, height_Coil_1/2, 0, c_Coil_1};
+p_Coil_1_m_ext_1=newp; Point(newp) = {x_[0]-(gap_Core_Coil_1+width_Coil_1), height_Coil_1/2, 0, c_Coil_1};
+
+l_Coil_1_m[]={};
+l_Coil_1_m[]+=newl; Line(newl) = {p_Coil_1_m_ext_0, p_Coil_1_m_ext_1};
+l_Coil_1_m[]+=newl; Line(newl) = {p_Coil_1_m_ext_1, p_Coil_1_m_int_1};
+l_Coil_1_m[]+=newl; Line(newl) = {p_Coil_1_m_int_1, p_Coil_1_m_int_0};
+
+l_Coil_1_m_Y0[]={};
+l_Coil_1_m_Y0[]+=newl; Line(newl) = {p_Coil_1_m_int_0, p_Coil_1_m_ext_0};
+
+ll_Coil_1_m=newll; Line Loop(newll) = {l_Coil_1_m[], l_Coil_1_m_Y0[]};
+s_Coil_1_m=news; Plane Surface(news) = {ll_Coil_1_m};
+
+Physical Surface("COIL_1_MINUS", COIL_1_MINUS) = {s_Coil_1_m};
+
+
+// COIL_2_PLUS
+
+x_[]=Point{p_Leg_2_R_0};
+p_Coil_2_p_int_0=newp; Point(newp) = {x_[0]+gap_Core_Coil_2, 0, 0, c_Coil_2};
+p_Coil_2_p_ext_0=newp; Point(newp) = {x_[0]+gap_Core_Coil_2+width_Coil_2, 0, 0, c_Coil_2};
+p_Coil_2_p_int_1=newp; Point(newp) = {x_[0]+gap_Core_Coil_2, height_Coil_2/2, 0, c_Coil_2};
+p_Coil_2_p_ext_1=newp; Point(newp) = {x_[0]+gap_Core_Coil_2+width_Coil_2, height_Coil_2/2, 0, c_Coil_2};
+
+l_Coil_2_p[]={};
+l_Coil_2_p[]+=newl; Line(newl) = {p_Coil_2_p_ext_0, p_Coil_2_p_ext_1};
+l_Coil_2_p[]+=newl; Line(newl) = {p_Coil_2_p_ext_1, p_Coil_2_p_int_1};
+l_Coil_2_p[]+=newl; Line(newl) = {p_Coil_2_p_int_1, p_Coil_2_p_int_0};
+
+l_Coil_2_p_Y0[]={};
+l_Coil_2_p_Y0[]+=newl; Line(newl) = {p_Coil_2_p_int_0, p_Coil_2_p_ext_0};
+
+ll_Coil_2_p=newll; Line Loop(newll) = {l_Coil_2_p[], l_Coil_2_p_Y0[]};
+s_Coil_2_p=news; Plane Surface(news) = {ll_Coil_2_p};
+
+Physical Surface("COIL_2_PLUS", COIL_2_PLUS) = {s_Coil_2_p};
+
+
+// COIL_2_MINUS
+
+x_[]=Point{p_Leg_2_L_0};
+p_Coil_2_m_int_0=newp; Point(newp) = {x_[0]-gap_Core_Coil_2, 0, 0, c_Coil_2};
+p_Coil_2_m_ext_0=newp; Point(newp) = {x_[0]-(gap_Core_Coil_2+width_Coil_2), 0, 0, c_Coil_2};
+p_Coil_2_m_int_1=newp; Point(newp) = {x_[0]-gap_Core_Coil_2, height_Coil_2/2, 0, c_Coil_2};
+p_Coil_2_m_ext_1=newp; Point(newp) = {x_[0]-(gap_Core_Coil_2+width_Coil_2), height_Coil_2/2, 0, c_Coil_2};
+
+l_Coil_2_m[]={};
+l_Coil_2_m[]+=newl; Line(newl) = {p_Coil_2_m_ext_0, p_Coil_2_m_ext_1};
+l_Coil_2_m[]+=newl; Line(newl) = {p_Coil_2_m_ext_1, p_Coil_2_m_int_1};
+l_Coil_2_m[]+=newl; Line(newl) = {p_Coil_2_m_int_1, p_Coil_2_m_int_0};
+
+l_Coil_2_m_Y0[]={};
+l_Coil_2_m_Y0[]+=newl; Line(newl) = {p_Coil_2_m_int_0, p_Coil_2_m_ext_0};
+
+ll_Coil_2_m=newll; Line Loop(newll) = {l_Coil_2_m[], l_Coil_2_m_Y0[]};
+s_Coil_2_m=news; Plane Surface(news) = {ll_Coil_2_m};
+
+Physical Surface("COIL_2_MINUS", COIL_2_MINUS) = {s_Coil_2_m};
+
+
+// AIR_WINDOW
+
+l_Air_Window_Y0[]={};
+l_Air_Window_Y0[]+=newl; Line(newl) = {p_Leg_1_R_0, p_Coil_1_p_int_0};
+l_Air_Window_Y0[]+=newl; Line(newl) = {p_Coil_1_p_ext_0, p_Coil_2_m_ext_0};
+l_Air_Window_Y0[]+=newl; Line(newl) = {p_Coil_2_m_int_0, p_Leg_2_L_0};
+
+ll_Air_Window=newll; Line Loop(newll) = {-l_Core_In[], -l_Coil_1_p[], l_Coil_2_m[], l_Air_Window_Y0[]};
+s_Air_Window=news; Plane Surface(news) = {ll_Air_Window};
+
+Physical Surface("AIR_WINDOW", AIR_WINDOW) = {s_Air_Window};
+
+
+// AIR_EXT
+
+x_[]=Point{p_Leg_2_R_1};
+p_Air_Ext_1_R_0=newp; Point(newp) = {x_[0]+gap_Core_Box_X, 0, 0, c_Box};
+p_Air_Ext_1_R_1=newp; Point(newp) = {x_[0]+gap_Core_Box_X, x_[1]+gap_Core_Box_Y, 0, c_Box};
+
+x_[]=Point{p_Leg_1_L_1};
+p_Air_Ext_1_L_0=newp; Point(newp) = {x_[0]-gap_Core_Box_X, 0, 0, c_Box};
+p_Air_Ext_1_L_1=newp; Point(newp) = {x_[0]-gap_Core_Box_X, x_[1]+gap_Core_Box_Y, 0, c_Box};
+
+l_Air_Ext[]={};
+l_Air_Ext[]+=newl; Line(newl) = {p_Air_Ext_1_R_0, p_Air_Ext_1_R_1};
+l_Air_Ext[]+=newl; Line(newl) = {p_Air_Ext_1_R_1, p_Air_Ext_1_L_1};
+l_Air_Ext[]+=newl; Line(newl) = {p_Air_Ext_1_L_1, p_Air_Ext_1_L_0};
+
+
+l_Air_Ext_Y0[]={};
+l_Air_Ext_Y0[]+=newl; Line(newl) = {p_Leg_2_R_0, p_Coil_2_p_int_0};
+l_Air_Ext_Y0[]+=newl; Line(newl) = {p_Coil_2_p_ext_0, p_Air_Ext_1_R_0};
+
+l_Air_Ext_Y0[]+=newl; Line(newl) = {p_Air_Ext_1_L_0, p_Coil_1_m_ext_0};
+l_Air_Ext_Y0[]+=newl; Line(newl) = {p_Coil_1_m_int_0, p_Leg_1_L_0};
+
+
+ll_Air_Ext=newll; Line Loop(newll) = {-l_Core_Out[], -l_Coil_2_p[], l_Coil_1_m[], l_Air_Ext[], l_Air_Ext_Y0[]};
+s_Air_Ext=news; Plane Surface(news) = {ll_Air_Ext};
+
+Physical Surface("AIR_EXT", AIR_EXT) = {s_Air_Ext};
+Physical Line("SUR_AIR_EXT", SUR_AIR_EXT) = {l_Air_Ext[]};

Magnetodynamics/transfo.pro

+/* -------------------------------------------------------------------
+   Tutorial 7b : magnetodyamic model of a single-phase transformer
+
+   Features:
+   - Time-domain and frequency-domain dynamical problems
+   - Use of a generic template formulation library
+   - Circuit coupling used as a black-box (see Tutorial 8 for more)
+
+   To compute the solution in a terminal:
+       getdp transfo -solve MagDyn_a_2D
+       getdp electromagnet -pos Map
+
+   To compute the solution interactively from the Gmsh GUI:
+       File > Open > transfo.pro
+       Run (button at the bottom of the left panel)
+   ------------------------------------------------------------------- */
+
+Include "transfo_common.pro";
+
+DefineConstant[
+  type_Conds = {2, Choices{1 = "Massive", 2 = "Coil"}, Highlight "Blue",
+    Name "Parameters/01Conductor type"}
+  type_Source = {2, Choices{1 = "Current", 2 = "Voltage"}, Highlight "Blue",
+    Name "Parameters/02Source type"}
+  type_Analysis = {1, Choices{1 = "Frequency-domain", 2 = "Time-domain"}, Highlight "Blue",
+    Name "Parameters/03Analysis type"}
+  Freq = {50, Min 0, Max 1e3, Step 1,
+    Name "Parameters/Frequency"}
+];
+
+Group {
+  /* Abstract regions that will be used in the "Lib_MagStaDyn_av_2D_Cir.pro"
+  template file included below; the regions are first intialized as empty,
+  before being filled with physical groups */
+
+  VolCC_Mag = Region[{}]; // Non-conducting regions
+  VolC_Mag = Region[{}]; // Massive conductors
+  VolS_Mag = Region[{}]; // Stranded conductors, i.e., coils
+
+  // air physical groups
+  Air = Region[{AIR_WINDOW, AIR_EXT}];
+  VolCC_Mag += Region[Air];
+
+  // exterior boundary
+  Sur_Air_Ext = Region[{SUR_AIR_EXT}];
+
+  // magnetic core of the transformer, assumed to be non-conducting
+  Core = Region[CORE];
+  VolCC_Mag += Region[Core];
+
+  Coil_1_P = Region[COIL_1_PLUS];
+  Coil_1_M = Region[COIL_1_MINUS];
+  Coil_1 = Region[{Coil_1_P, Coil_1_M}];
+
+  Coil_2_P = Region[COIL_2_PLUS];
+  Coil_2_M = Region[COIL_2_MINUS];
+  Coil_2 = Region[{Coil_2_P, Coil_2_M}];
+
+  Coils = Region[{Coil_1, Coil_2}];
+
+  If (type_Conds == 1)
+    VolC_Mag += Region[{Coils}];
+  ElseIf (type_Conds == 2)
+    VolS_Mag += Region[{Coils}];
+    VolCC_Mag += Region[{Coils}];
+  EndIf
+}
+
+
+Function {
+  mu0 = 4e-7*Pi;
+
+  mu[Air] = 1 * mu0;
+
+  mur_Core = 1;
+  mu[Core] = mur_Core * mu0;
+
+  mu[Coils] = 1 * mu0;
+  sigma[Coils] = 1e7;
+
+  // For a correct definition of the voltage
+  CoefGeo = thickness_Core;
+
+  // To be defined separately for each coil portion
+  Sc[Coil_1_P] = SurfaceArea[];
+  SignBranch[Coil_1_P] = 1; // To fix the convention of positive current (1:
+                            // along Oz, -1: along -Oz)
+
+  Sc[Coil_1_M] = SurfaceArea[];
+  SignBranch[Coil_1_M] = -1;
+
+  Sc[Coil_2_P] = SurfaceArea[];
+  SignBranch[Coil_2_P] = 1;
+
+  Sc[Coil_2_M] = SurfaceArea[];
+  SignBranch[Coil_2_M] = -1;
+
+  // Number of turns (same for PLUS and MINUS portions) (half values because
+  // half coils are defined)
+  Ns[Coil_1] = 1;
+  Ns[Coil_2] = 1;
+
+  // Global definitions (nothing to change):
+
+  // Current density in each coil portion for a unit current (will be multiplied
+  // by the actual total current in the coil)
+  js0[Coils] = Ns[]/Sc[] * Vector[0,0,SignBranch[]];
+  CoefGeos[Coils] = SignBranch[] * CoefGeo;
+
+  // The reluctivity will be used
+  nu[] = 1/mu[];
+}
+
+If(type_Analysis == 1)
+  Flag_FrequencyDomain = 1;
+Else
+  Flag_FrequencyDomain = 0;
+EndIf
+
+If (type_Source == 1) // current
+
+  Flag_CircuitCoupling = 0;
+
+ElseIf (type_Source == 2) // voltage
+
+  // PLUS and MINUS coil portions to be connected in series, with applied
+  // voltage on the resulting branch
+  Flag_CircuitCoupling = 1;
+
+  // Here is the definition of the circuits on primary and secondary sides:
+  Group {
+    // Empty Groups to be filled
+    Resistance_Cir  = Region[{}];
+    Inductance_Cir  = Region[{}] ;
+    Capacitance_Cir = Region[{}] ;
+    SourceV_Cir = Region[{}]; // Voltage sources
+    SourceI_Cir = Region[{}]; // Current sources
+
+    // Primary side
+    E_in = Region[10001]; // arbitrary region number (not linked to the mesh)
+    SourceV_Cir += Region[{E_in}];
+
+    // Secondary side
+    R_out = Region[10101]; // arbitrary region number (not linked to the mesh)
+    Resistance_Cir += Region[{R_out}];
+  }
+
+  Function {
+    deg = Pi/180;
+    // Input RMS voltage (half of the voltage because of symmetry; half coils
+    // are defined)
+    val_E_in = 1.;
+    phase_E_in = 90 *deg; // Phase in radian (from phase in degree)
+    // High value for an open-circuit test; Low value for a short-circuit test;
+    // any value in-between for any charge
+    Resistance[R_out] = 1e6;
+  }
+
+  Constraint {
+
+    { Name Current_Cir ;
+      Case {
+      }
+    }
+
+    { Name Voltage_Cir ;
+      Case {
+        { Region E_in; Value val_E_in; TimeFunction F_Cos_wt_p[]{2*Pi*Freq, phase_E_in};}
+      }
+    }
+
+    { Name ElectricalCircuit ; Type Network ;
+      Case Circuit_1 {
+        // PLUS and MINUS coil portions to be connected in series, together with
+        // E_in (an additional resistance should be defined to represent the
+        // Coil_1 end-winding (not considered in the 2D model))
+        { Region E_in; Branch {1,2}; }
+
+        { Region Coil_1_P; Branch {2,3} ; }
+        { Region Coil_1_M; Branch {3,1} ; }
+      }
+      Case Circuit_2 {
+        // PLUS and MINUS coil portions to be connected in series, together with
+        // R_out (an additional resistance should be defined to represent the
+        // Coil_2 end-winding (not considered in the 2D model))
+        { Region R_out; Branch {1,2}; }
+
+        { Region Coil_2_P; Branch {2,3} ; }
+        { Region Coil_2_M; Branch {3,1} ; }
+      }
+    }
+
+  }
+
+EndIf
+
+Constraint {
+  { Name MagneticVectorPotential_2D;
+    Case {
+      { Region Sur_Air_Ext; Value 0; }
+    }
+  }
+  { Name Current_2D;
+    Case {
+     If (type_Source == 1)
+      // Current in each coil (same for PLUS and MINUS portions)
+      { Region Coil_1; Value 1; }
+      { Region Coil_2; Value 0; }
+     EndIf
+    }
+  }
+  { Name Voltage_2D;
+    Case {
+    }
+  }
+}
+
+Include "Lib_MagStaDyn_av_2D_Cir.pro";
+
+PostOperation {
+  { Name Map_a; NameOfPostProcessing MagDyn_a_2D;
+    Operation {
+      Print[ j, OnElementsOf Region[{VolC_Mag, VolS_Mag}], Format Gmsh, File "j.pos" ];
+      Print[ b, OnElementsOf Vol_Mag, Format Gmsh, File "b.pos" ];
+      Print[ az, OnElementsOf Vol_Mag, Format Gmsh, File "az.pos" ];
+
+      If (type_Source == 1) // current
+        // In text file UI.txt: voltage and current for each coil portion (note
+        // that the voltage is not equally distributed in PLUS and MINUS
+        // portions, which is the reason why we must apply the total voltage
+        // through a circuit -> type_Source == 2)
+        Echo[ "Coil_1_P", Format Table, File "UI.txt" ];
+        Print[ U, OnRegion Coil_1_P, Format FrequencyTable, File > "UI.txt" ];
+        Print[ I, OnRegion Coil_1_P, Format FrequencyTable, File > "UI.txt"];
+        Echo[ "Coil_1_M", Format Table, File > "UI.txt" ];
+        Print[ U, OnRegion Coil_1_M, Format FrequencyTable, File > "UI.txt" ];
+        Print[ I, OnRegion Coil_1_M, Format FrequencyTable, File > "UI.txt"];
+
+        Echo[ "Coil_2_P", Format Table, File > "UI.txt" ];
+        Print[ U, OnRegion Coil_2_P, Format FrequencyTable, File > "UI.txt" ];
+        Print[ I, OnRegion Coil_2_P, Format FrequencyTable, File > "UI.txt"];
+        Echo[ "Coil_2_M", Format Table, File > "UI.txt" ];
+        Print[ U, OnRegion Coil_2_M, Format FrequencyTable, File > "UI.txt" ];
+        Print[ I, OnRegion Coil_2_M, Format FrequencyTable, File > "UI.txt"];
+
+      ElseIf (type_Source == 2)
+        // In text file UI.txt: voltage and current of the primary coil (from E_in)
+        // (real and imaginary parts!)
+        Echo[ "E_in", Format Table, File "UI.txt" ];
+        Print[ U, OnRegion E_in, Format FrequencyTable, File > "UI.txt" ];
+        Print[ I, OnRegion E_in, Format FrequencyTable, File > "UI.txt"];
+
+        // In text file UI.txt: voltage and current of the secondary coil (from R_out)
+        Echo[ "R_out", Format Table, File > "UI.txt" ];
+        Print[ U, OnRegion R_out, Format FrequencyTable, File > "UI.txt" ];
+        Print[ I, OnRegion R_out, Format FrequencyTable, File > "UI.txt"];
+      EndIf
+    }
+  }
+}

Magnetodynamics/transfo_common.pro

+
+// Dimensions
+
+width_Core = 1.;
+height_Core = 1.;
+
+// Thickness along Oz (to be considered for a correct definition of voltage)
+thickness_Core = 1.;
+
+width_Window = 0.5;
+height_Window = 0.5;
+
+width_Core_Leg = (width_Core-width_Window)/2.;
+
+width_Coil_1 = 0.10;
+height_Coil_1 = 0.25;
+gap_Core_Coil_1 = 0.05;
+
+width_Coil_2 = 0.10;
+height_Coil_2 = 0.25;
+gap_Core_Coil_2 = 0.05;
+
+gap_Core_Box_X = 1.;
+gap_Core_Box_Y = 1.;
+
+// Characteristic lenghts (for mesh sizes)
+
+s = 1;
+
+c_Core = width_Core_Leg/10. *s;
+
+c_Coil_1 = height_Coil_1/2/5 *s;
+c_Coil_2 = height_Coil_2/2/5 *s;
+
+c_Box = gap_Core_Box_X/6. *s;
+
+// Physical regions
+
+AIR_EXT = 1001;
+SUR_AIR_EXT = 1002;
+AIR_WINDOW = 1011;
+
+CORE = 1050;
+
+COIL_1_PLUS = 1101;
+COIL_1_MINUS = 1102;
+
+COIL_2_PLUS = 1201;
+COIL_2_MINUS = 1202;

Magnetostatics/electromagnet.pro

@@ -8,7 +8,7 @@
 
    To compute the solution in a terminal:
        getdp electromagnet -solve MagSta_a
-       getdp electromagnet -pos Map
+       getdp electromagnet -pos Map_a
 
    To compute the solution interactively from the Gmsh GUI:
        File > Open > electromagnet.pro
@@ -66,13 +66,13 @@ Group {
      The abstract regions in this model have the following interpretation:
      - Vol_Nu_Mag  = region where the term [ nu[] * Dof{d a} , {d a} ] is assembled
      - Vol_Js_Mag  = region where the term [ - Dof{js} , {a} ] is assembled
-     - Vol_Inf     = region where the infinite ring geometric transformation is applied
+     - Vol_Inf_Mag = region where the infinite ring geometric transformation is applied
      - Sur_Dir_Mag = Homogeneous Dirichlet part of the model's boundary;
      - Sur_Neu_Mag = Homogeneous Neumann part of the model's boundary;
   */
   Vol_Nu_Mag  = Region[ {Air, AirInf, Core, Ind} ];
   Vol_Js_Mag  = Region[ Ind ];
-  Vol_Inf     = Region[ {AirInf} ];
+  Vol_Inf_Mag = Region[ {AirInf} ];
   Sur_Dir_Mag = Region[ {Surface_bn0, Surface_Inf} ];
   Sur_Neu_Mag = Region[ {Surface_ht0} ];
 }
@@ -168,7 +168,7 @@ Include "electromagnet_common.pro";
 Val_Rint = rInt; Val_Rext = rExt;
 Jacobian {
   { Name Vol ;
-    Case { { Region Vol_Inf ;
+    Case { { Region Vol_Inf_Mag ;
              Jacobian VolSphShell {Val_Rint, Val_Rext} ; }
            { Region All ; Jacobian Vol ; }
     }