GetDP ¶

// Gmsh script describing the geometry of the microstrip line and the associated
// meshing constraints.
//
// The microstrip line is a three-dimensional structure, but here we exploit the
// translation invariance along the line and model a 2D cross-section (a
// two-dimensional cut in a plane of constant "z"). All results are thus
// understood as values "per unit length along z".
//
// In this first tutorial we use the built-in Gmsh CAD kernel. Later tutorials
// (starting with tutorial 2) will use the more powerful OpenCASCADE kernel to
// build geometries from boolean operations on solid primitives.

// Dimensions (there are no units in Gmsh -- the following values should be
// interpreted as dimensions in meters due to the definition of the physical
// parameters in MKS units in "microstrip.pro"):
h = 1.e-3; // thickness of dielectric substrate
w = 4.72e-3; // width of microstrip line
t = 0.035e-3; // thickness of microstrip line
xBox = w/2. * 6.; // width of air box
yBox = h * 12.; // height of air box

// Global mesh size factor (that can be modified interactively in the Gmsh
// graphical interface):
s = DefineNumber[1., Name "Parameters/Global mesh size factor"];

// Target mesh sizes:
p0 = h / 10 * s;
pLine0 = w / 20 * s;
pLine1 = w / 100 * s;
pxBox = xBox / 10 * s;
pyBox = yBox / 8 * s;

// We create this first geometry in a bottom-up manner with the built-in Gmsh
// CAD kernel, successively defining model points, lines, loops and surfaces:
Point(1) = {0, 0, 0, p0};
Point(2) = {xBox, 0, 0, pxBox};
Point(3) = {xBox, h, 0, pxBox};
Point(4) = {0, h, 0, pLine0};
Point(5) = {w/2, h, 0, pLine1};
Point(6) = {0, h+t, 0, pLine0};
Point(7) = {w/2, h+t, 0, pLine1};
Point(8) = {0, yBox, 0, pyBox};
Point(9) = {xBox, yBox, 0, pyBox};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 9};
Line(4) = {9, 8};
Line(5) = {8, 6};
Line(7) = {4, 1};
Line(8) = {5, 3};
Line(9) = {4, 5};
Line(10) = {6, 7};
Line(11) = {5, 7};
Curve Loop(12) = {1, 2, -8, -9, 7};
Plane Surface(13) = {12};
Curve Loop(14) = {10, -11, 8, 3, 4, 5};
Plane Surface(15) = {14};

// The last step is to define physical groups, assigning unique tags (and names)
// to groups of model entities. Physical groups serve two purposes: when they
// are defined they tell Gmsh which elements to save in the mesh file (by
// default only elements belonging to at least one physical group are saved),
// and they provide the tags that GetDP uses in "Region[]" to identify parts of
// the domain.
Physical Surface("Air", 1) = {15};
Physical Surface("Dielectric", 2) = {13};
Physical Curve("Ground", 10) = {1};
Physical Curve("Microstrip boundary", 11) = {9, 10, 11};
Physical Curve("Inf", 12) = {2, 3, 4};

// We consider the calculation of the electric field given a static distribution
// of electric potential. This corresponds to an electrostatic physical model,
// obtained by combining the time-invariant Faraday equation ("curl e = 0", with
// "e" the electric field) with Gauss' law ("div d = rho", with "d" the
// displacement field and "rho" the charge density) and the dielectric
// constitutive law ("d = epsilon e", with "epsilon" the dielectric
// permittivity).
//
// Since "curl e = 0", "e" can be derived from a scalar electric potential "v",
// such that "e = -grad v". Plugging this potential in Gauss' law and using the
// constitutive law leads to a scalar (generalized) Poisson equation in terms of
// the electric potential: "-div(epsilon grad v) = rho".
//
// We consider here the special case where "rho = 0" to model a conducting
// microstrip line on top of a dielectric substrate. A Dirichlet boundary
// condition sets the potential to 1 mV on the boundary of the microstrip line
// and to 0 V on the ground. A homogeneous Neumann boundary condition (zero
// normal component of the displacement field, i.e. "n . d = 0") is imposed on
// the left boundary of the domain to account for the symmetry of the problem. A
// homogeneous Neumann condition is also called a "natural" boundary condition:
// it appears in the weak formulation as a vanishing boundary integral, and
// therefore requires no explicit treatment in the model -- we simply do
// nothing. The domain is truncated on the top and right with a homogeneous
// Dirichlet boundary condition ("v = 0"), assumed to be imposed sufficiently
// far away from the microstrip.

Group {
  // Create region groups associated with the physical groups defined in the
  // "microstrip.msh" mesh file:
  Air = Region[ 1 ];
  Dielectric = Region[ 2 ];
  Ground = Region[ 10 ];
  Microstrip = Region[ 11 ];
  Inf = Region[ 12 ];

  // Define abstract regions to be used below in the definition of the scalar
  // electric potential formulation:
  //  - "Vol_Ele": overall volume domain where "-div(epsilon grad v) = 0" is
  //    solved; it contains only "volume" elements of the mesh (triangles here)
  //  - "Sur_Neu_Ele": surface where non homogeneous Neumann boundary conditions
  //    (on "n . d = -n . (epsilon grad v)") are imposed; it contain only
  //    "surface" elements of the mesh (lines here).
  //
  // The purpose of abstract regions is to allow a generic definition of the
  // FunctionSpace, Formulation and PostProcessing objects with no reference to
  // model-specific physical groups. We will show in tutorial 9 how abstract
  // formulations can then be isolated in geometry-independent template files,
  // thanks to an appropriate declaration mechanism (using "DefineConstant[]",
  // "DefineGroup[]" and "DefineFunction[]").
  //
  // Since there are no non-homogeneous Neumann conditions in this particular
  // example, "Sur_Neu_Ele" is defined as empty.
  //
  // Note that volume elements are those that correspond to the higher dimension
  // of the model at hand (2D elements here), surface elements correspond to the
  // higher dimension of the model minus one (1D elements here).
  Vol_Ele = Region[ {Air, Dielectric} ];
  Sur_Neu_Ele = Region[ {} ];
}

Function {
  // Material laws (here the dielectric permittivity) are defined as piecewise
  // functions (note the square brackets), in terms of the above defined groups:
  eps0 = 8.854187818e-12;
  epsilon[ Air ] = 1. * eps0;
  epsilon[ Dielectric ] = 9.8 * eps0;
}

Constraint {
  // The Dirichlet boundary condition is also defined piecewise, through the
  // following "v_Ele" constraint, invoked in the FunctionSpace below. The
  // constraint type "Assign" means that the coefficients in the finite element
  // expansion will be assigned the prescribed values:
  { Name v_Ele; Type Assign;
    Case {
      { Region Ground; Value 0.; }
      { Region Microstrip; Value 1.e-3; }
      { Region Inf; Value 0; }
    }
  }
}

Group{
  // The domain of definition of a FunctionSpace lists all regions on which a
  // field is defined. Domains of definitions may contain both volume and
  // surface regions, so we use a "Dom_" prefix to avoid confusion:
  Dom_H1_v_Ele =  Region[ {Vol_Ele, Sur_Neu_Ele} ];
}

FunctionSpace {
  // The function space in which we pick the electric scalar potential "v"
  // solution is defined by
  //  - a domain of definition (the "Support": "Dom_H1_v_Ele")
  //  - a type ("Form0", meaning "scalar field")
  //  - a set of basis functions ("BF_Node" for scalar nodal basis functions,
  //    i.e. isoparametric Lagrange elements)
  //  - a set of entities to which the basis functions are associated ("Entity":
  //    here all the nodes of the domain of definition "NodesOf[All]")
  //  - a constraint (here the Dirichlet boundary conditions)
  //
  // The finite element expansion of the unknown field "v" reads
  //
  //   v(x, y) = Sum_k vn_k sn_k(x, y)
  //
  // where the "vn_k" coefficients are the nodal values and "sn_k(x, y)" are the
  // nodal basis functions. Not all coefficients vn_k are unknowns of the finite
  // element problem, due to the Dirichlet constraint, which assigns particular
  // values to the nodes of the "Ground" and "Microstrip" regions.
  //
  // The default mesh "microstrip.msh" is made of first order (3-node)
  // triangles: the nodal basis functions sn_k(x, y) are then piece-wise linear
  // on each triangle. If second order (6-node) triangles are used instead
  // ("gmsh microstrip.geo -order 2 -2"), the basis functions will be piece-wise
  // quadratic on each triangle. In all cases, with Lagrange elements, we have
  // "sn_k(x_l, y_l) = \delta_kl" (the Kronecker delta, which is 1 if "k = l"
  // and 0 otherwise) if "(x_l, y_l)" denotes the coordinates of node "l".
  { Name H1_v_Ele; Type Form0;
    BasisFunction {
      { Name sn; NameOfCoef vn; Function BF_Node;
        Support Dom_H1_v_Ele; Entity NodesOf[ All ]; }
      // Using "NodesOf[All]" instead of "NodesOf[Dom_H1_v_Ele]" is an
      // optimization, which avoids explicitly building the list of all the
      // nodes. It is always safe here: GetDP restricts the actual basis
      // functions to those that have support in "Dom_H1_v_Ele". In cases where
      // different basis functions must be associated with disjoint subsets of
      // nodes (see e.g. tutorial 6), the explicit form "NodesOf[ <subset> ]"
      // becomes necessary.
    }
    Constraint {
      { NameOfCoef vn; EntityType NodesOf; NameOfConstraint v_Ele; }
    }
  }
}

Jacobian {
  // Jacobians specify the mapping between elements in the mesh and the
  // reference elements over which integration is performed.
  //
  // "Vol" stands for the 1-to-1 mapping between identical spatial dimensions,
  // i.e. in this case a reference triangle (2D) onto triangles in the "z = 0"
  // plane (2D):
  { Name Vol;
    Case {
      { Region All; Jacobian Vol; }
    }
  }
  // "Sur" is used to map the reference line segment (1D) onto lines in the
  // plane (2D). It is not strictly needed in this particular example (because
  // "Sur_Neu_Ele" is empty), but is defined here for completeness, so that
  // any surface term that might be added later to the Formulation (for
  // example the optional non-homogeneous Neumann term shown below) can
  // readily refer to it:
  { Name Sur;
    Case {
      { Region All; Jacobian Sur; }
    }
  }
}

Integration {
  // A Gauss numerical quadrature is used for the integrations:
  { Name Int;
    Case {
      { Type Gauss;
        Case {
          // One quadrature point is sufficient to integrate exactly the
          // stiffness matrix with first order elements (constant value per
          // element):
          { GeoElement Triangle; NumberOfPoints 1; }
          // But we need 3 quadrature points for second order elements:
          { GeoElement Triangle2; NumberOfPoints 3; }
        }
      }
    }
  }
}

Formulation {
  // The Formulation encodes the weak formulation of the original partial
  // differential equation, i.e. of "-div(epsilon grad v) = 0". This weak
  // formulation involves finding "v" such that
  //
  //   (-div(epsilon grad v) , v')_Vol_Ele = 0
  //
  // holds for all test functions "v'", where "(. , .)_D" denotes the inner
  // product over a domain "D". Integrating by parts (i.e. using Green's first
  // identity: "(div(f), u)_D = -(f, grad u)_D + (n . f, u)_Bnd_D"), the weak
  // formulation becomes: find "v" such that
  //
  //   (epsilon grad v, grad v')_Vol_Ele
  //     - (n . (epsilon grad v), v')_Bnd_Vol_Ele = 0
  //
  // holds for all "v'", where "Bnd_Vol_Ele" is the boundary of "Vol_Ele". In
  // our microstrip example this surface term vanishes, since
  //  - on the "Ground", "Microstrip" and "Inf" parts of the boundary we impose a
  //    Dirichlet boundary condition on "v" (through the "Assign" constraint)
  //    and the test functions "v'" vanish (there are actually none, as the
  //    corresponding unknowns have been fixed);
  //  - on the remaining part of the boundary (on the left side, at y = 0) by
  //    symmetry "n . d = - n. (epsilon grad v) = 0", i.e. we have a homogeneous
  //    Neumann ("natural") boundary condition.
  //
  // We are thus eventually looking for functions "v" in the function space
  // "H1_v_Ele" such that
  //
  //   (epsilon grad v, grad v')_Vol_Ele = 0
  //
  // holds for all "v'". Finally, our choice here is to use a Galerkin method,
  // where the test functions "v'" are the same basis functions ("sn_k") as the
  // ones used to interpolate the unknown potential "v".
  //
  // The "Integral" statement in the Formulation is a symbolic representation of
  // this weak formulation. It has 4 semicolon separated arguments:
  //  - the density "[. , .]" to be integrated (note the square brackets instead
  //    of the parentheses), with the test functions (always) after the comma;
  //  - the domain of integration;
  //  - the Jacobian of the transformation between the reference element and the
  //    element in the mesh;
  //  - the integration method.
  //
  // In the density, braces around a quantity (such as "{v}") indicate that this
  // quantity belongs to a FunctionSpace. Differential operators can be applied
  // within braces (such as "{Grad v}"); in particular the symbol "d" represents
  // the exterior derivative, and it is a synonym of "Grad" when applied to a
  // scalar function, declared with type "Form0" in the FunctionSpace.
  //
  // As the Galerkin method uses as test functions the same basis functions
  // "sn_k" as for the unknown field "v", the second term in the density should
  // be something like [ ... , basis_functions_of {d v} ]. However, since the
  // second term is always devoted to test functions, the operator
  // "basis_functions_of" would always be there.  It can therefore be made
  // implicit and, in the GetDP syntax, it is omitted. So, one simply writes "[
  // ... , {d v} ]".
  //
  // On the other hand, the first term can contain a much wider variety of
  // expressions. In our case it should be expressed in terms of the finite
  // element expansion of "v" at the present system solution, i.e. when the
  // coefficients "vn_k" in the expansion "v = Sum_k vn_k sn_k" are
  // unknown. This is indicated by prefixing the braces with "Dof" ("degree of
  // freedom"), which leads to the following density:
  //
  //   [ epsilon[] * Dof{d v} , {d v} ],
  //
  // a bilinear term that will contribute to the stiffness matrix of the
  // electrostatic problem at hand.
  //
  // Another option, which does not work here (because it would make the
  // problem trivial: there would be no unknowns left), is to evaluate the
  // first argument with the last available computed solution, i.e. simply
  // perform the interpolation with known coefficients "vn_k". The "Dof"
  // prefix is then omitted and we would have:
  //
  //   [ epsilon[] * {d v} , {d v} ],
  //
  // a linear term that would contribute to the right-hand side of the linear
  // system. This mechanism is not useful for a linear, one-shot problem like
  // the present one, but it is essential whenever a term depends on a
  // previously computed field: see tutorial 3 (where the reluctivity "nu" is
  // evaluated at the current iterate "{d a}" inside a Newton-Raphson loop)
  // and tutorial 7 (where the electrical conductivity is evaluated at the
  // temperature field "{T}" computed by a separate formulation).

  { Name Electrostatics_v; Type FemEquation;
    Quantity {
      { Name v; Type Local; NameOfSpace H1_v_Ele; }
    }
    Equation {
      Integral { [ epsilon[] * Dof{d v} , {d v} ];
        In Vol_Ele; Jacobian Vol; Integration Int; }

      // Additional "Integral" terms can be added here. For example, the
      // following term may account for non-homogeneous Neumann boundary
      // conditions, provided that the function "nd[]" is defined:
      //
      // Integral { [ nd[] , {v} ];
      //   In Sur_Neu_Ele; Jacobian Sur; Integration Int; }

      // All the terms in the "Equation" are added, and an implicit "= 0" is
      // considered at the end.
    }
  }
}

// In "Resolution" we specify what to do with a weak formulation: here we simply
// generate a linear system, solve it and save the solution to disk (in a ".res"
// file).
Resolution {
  { Name Ele;
    System {
      { Name Sys_Ele; NameOfFormulation Electrostatics_v; }
    }
    Operation {
      Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele];
    }
  }
}

// Post-processing is done in two parts.
//
// The first part defines, in terms of the Formulation, which itself refers to
// the FunctionSpace, a number of quantities that can be evaluated at the
// post-processing stage. The three quantities defined here are:
//  - the scalar electric potential;
//  - the electric field;
//  - the electric displacement.
PostProcessing {
  { Name Ele; NameOfFormulation Electrostatics_v;
    Quantity {
      { Name v; Value {
          Term { [ {v} ]; In Vol_Ele; Jacobian Vol; }
        }
      }
      { Name e; Value {
          Term { [ -{d v} ]; In Vol_Ele; Jacobian Vol; }
        }
      }
      { Name d; Value {
          Term { [ -epsilon[] * {d v} ]; In Vol_Ele; Jacobian Vol; }
        }
      }
    }
  }
}

// The second part consists in defining post-processing operations, which can be
// invoked separately. (The first PostOperation is invoked by default when Gmsh
// is run interactively. The generated post-processing files are automatically
// displayed by Gmsh if the "Merge result automatically" option is enabled in
// the Gmsh "gear" menu.)

tol = 1.e-7; // small offset to ensure the cut is inside the simulation domain
yCut = 2.e-3; // vertical position of the cut

PostOperation {
  { Name Map; NameOfPostProcessing Ele;
     Operation {
       Print [ v, OnElementsOf Vol_Ele, File "v.pos" ];
       Print [ e, OnElementsOf Vol_Ele, File "e.pos" ];
       Print [ d, OnElementsOf Vol_Ele, File "d.pos" ];
       Print [ e, OnLine {{tol, yCut, 0}{14.e-3, yCut, 0}}{60}, File "e_cut.pos" ];
     }
  }
  { Name Cut; NameOfPostProcessing Ele;
    // Same cut as above, with more points and exported in raw text format.
    // "Format Table" requests an ASCII file (one row of numbers per evaluation
    // point) instead of the default Gmsh ".pos" format, which carries mesh
    // information and is automatically rendered by Gmsh. A ".txt" file is
    // primarily intended for post-processing by external tools (spreadsheets,
    // plotting scripts, etc.).
    //
    // Note that PostOperations other than the default one are invoked from the
    // Gmsh GUI through the "Modules > GetDP > Post-processing" menu (or from
    // the command line, e.g. with "getdp microstrip.pro -solve Ele -pos Cut"):
    Operation {
      Print [ e, OnLine {{tol, tol, 0}{14.e-3, tol, 0}} {500}, Format Table,
        File "e_cut.txt" ];
    }
  }
}

bw = 8e-3; // base plate width (for one fin)
DefineConstant[
  N = {1, Min 1, Max 100, Step 1, Name "Parameters/Number of fins"}
];

// Gmsh script describing the geometry of the radiator and the associated
// meshing constraints.

// Instead of the built-in Gmsh geometry kernel we use the OpenCASCADE kernel,
// in order to benefit from boolean operations:
SetFactory("OpenCASCADE");

// Definition of a global mesh size factor that can be modified interactively in
// the Gmsh graphical interface (in the same way as with DefineNumber[] - see
// tutorial 1), but that can also be overridden beforehand, e.g. on the command
// line with "-setnumber s value":
DefineConstant[
  s = {1., Name "Parameters/Global mesh size factor"}
];

fh = 40e-3; // fin height
fw = 2e-3; // fin width
bh = 10e-3; // base plate height
lc = fw / 3 * s; // mesh size

// We include a file that defines the number of fins and the width of one cell:
// the same file will be included in radiator.pro to guarantee consistency
// between Gmsh and GetDP when defining the periodic boundary condition
Include "radiator_common.pro";

DefineConstant[
  position = {0.5, Min 0.1, Max 0.9, Step 0.1,
    Name "Parameters/Fin position in [0,1]"}
  dim = {2, Choices{2="2D", 3="3D"},
    Name "Parameters/0Model dimension"}
  zh = {fh / 20, Min 1e-3, Max 2 * fh, Step 1e-3, Visible dim == 3,
    Name "Parameters/Fin thickness [m]"}
];

For i In {0 : N - 1}
  x0 = i * bw;
  Rectangle(2 * i + 1) = {x0 + position * (bw - fw), bh, 0, fw, fh}; // fin
  Rectangle(2 * i + 2) = {x0, 0, 0, bw, bh}; // base plate
EndFor

surf() = BooleanUnion{ Surface{1}; Delete; }{ Surface{2 : 2 * N}; Delete; };

MeshSize{:} = lc;

// Note that a number prefix in a parameter name is invisible in the graphical
// user interface ("0Model dimension" will appear as "Model dimension"). It
// enables sorting the menu entries, i.e. in this case making sure that "Model
// dimension" is the first entry in the "Parameters" menu.

e = 1e-6; // small tolerance for bounding box searches below

If(dim == 2)
  interior = surf(0);
  bottom() = Curve In BoundingBox{-e, -e, -e, N * bw + e, e, e};
  left() = Curve In BoundingBox{-e, -e, -e, e, bh + e, e};
  right() = Curve In BoundingBox{N * bw - e, -e, -e, N * bw + e, bh + e, e};
  top() = Boundary{ Surface{interior}; };
  top() -= {bottom(), left(), right()};

  // Set a periodic mesh constraint on the right curve: this forces Gmsh to
  // generate identical meshes on the "right" and "left" boundaries (the nodes
  // and elements of one are obtained from those of the other by the prescribed
  // translation). This matching of the two meshes is what enables the "Link"
  // constraint defined in "radiator.pro" to connect the degrees of freedom of
  // the two boundaries node-by-node:
  Periodic Curve { right() } = { left() } Translate {N * bw, 0, 0};
Else
  ex() = Extrude {0, 0, zh}{ Surface{surf()}; };
  interior = ex(1);
  bottom() = Surface In BoundingBox{-e, -e, -e, N * bw + e, e, zh + e};
  left() = Surface In BoundingBox{-e, -e, -e, e, bh + e, zh + e};
  right() = Surface In BoundingBox{N * bw - e, -e, -e, N * bw + e, bh + e, zh + e};
  top() = Boundary{ Volume{interior}; };
  top() -= {bottom(), left(), right()}; // try adding e.g. surf() and ex(0)

  // Set a periodic mesh constraint on the right surface (same purpose as in
  // the 2D case above):
  Periodic Surface { right() } = { left() } Translate {N * bw, 0, 0};
EndIf

// Define physical groups (GeoEntity{dim} can be used as a synonym for Curve,
// Surface and Volume respectively for dim == 1, 2 and 3):
Physical GeoEntity{dim} ("Radiator", 1) = interior;
Physical GeoEntity{dim - 1} ("Bottom base plate", 10) = bottom();
Physical GeoEntity{dim - 1} ("Left base plate", 11) = left();
Physical GeoEntity{dim - 1} ("Right base plate", 12) = right();
Physical GeoEntity{dim - 1} ("Top fins", 13) = top();

// This model computes the temperature profile in a radiator with fins.
// Starting at an initial temperature "T = T0", a heat flux "qn = qn0" is
// imposed through the bottom surface of the radiator base plate. Convection on
// the fins evacuates the heat based on Newton's law of cooling "qn = h (T -
// T0)", with "h" the convection coefficient.
//
// The governing partial differential equation links the temperature "T" and the
// heat flux density vector "q" through
//
//   rho * cp * \partial_t T + div(q) = 0,
//
// where "rho" is the mass density and "cp" the specific heat. Introducing
// Fourier's law "q = -k grad T", with "k" the thermal conductivity, we get the
// heat diffusion equation in terms of the temperature:
//
//   rho * cp * \partial_t T - div(k grad T) = 0.
//
// We will see below that the heat flux will be imposed naturally through a
// non-homogeneous Neumann boundary condition, while the convection condition
// will be imposed through a Robin boundary condition.

Group {
  // Physical regions:
  Radiator = Region[ 1 ];
  Bottom = Region[ 10 ];
  Left = Region[ 11 ];
  Right = Region[ 12 ];
  Top = Region[ 13 ];

  // The abstract regions in this model have the following interpretation:
  //  - "Vol_The": overall domain
  //  - "Sur_Neu_The": part of the boundary with non-homogeneous Neumann
  //    condition
  //  - "Sur_Rob_The": part of the boundary with Robin condition
  Vol_The = Region[ {Radiator} ];
  Sur_Neu_The = Region[ {Bottom} ];
  Sur_Rob_The = Region[ {Top} ];
}

Function {
  // "DefineConstant[]" declares parameters that can be modified interactively
  // in the Gmsh GUI, or overridden on the command line with "-setnumber name
  // value":
  DefineConstant[
    AnalysisType = {0, Choices{0="Steady-state", 1="Transient"},
      Name "Parameters/Analysis"}
    qn0 = {3e4, Min 1, Max 1e6, Step 1e3,
      Name "Parameters/Heat flux [W m^-2]"}
    h = {60, Min 0, Max 120, Step 1,
      Name "Parameters/Convection coefficient [W m^-2 K^-1]"}
    tmax = {360, Min 1, Max 3600, Step 10,
      Name "Parameters/Simulation time [s]", Visible AnalysisType}
    dt = {5, Min 0.1, Max 100, Step 0.1,
      Name "Parameters/Time step [s]", Visible AnalysisType}
  ];

  rho_cp[] = 2700 * 900;
  k[] = 170;
  qn0[] = -qn0; // negative sign to have flux coming into "Vol_The"
  h[] = h;
  T0[] = 20;
}

// The weak formulation reads: find "T" such that
//
//   (rho * cp * \partial_t T, T')_Vol_The - (div(k grad T), T')_Vol_The = 0
//
// holds for all test functions "T'". After integration by parts it reads: find
// "T" such that
//
//   (rho * cp * \partial_t T, T')_Vol_The + (k grad T, grad T')_Vol_The
//     - (k grad T . n, T')_Bnd_Vol_The = 0
//
// holds for all test functions "T'". The boundary term is split to handle the
// imposed flux ("-k grad T . n = qn0") as a non-homogeneous Neumann boundary
// condition on "Sur_Neu_The", and the convection condition ("-k grad T . n = h
// (T - T0)") as a Robin condition on "Sur_Rob_The". The final weak formulation
// then reads: find "T" such that
//
//   (rho * cp * \partial_t T, T')_Vol_The + (k grad T, grad T')_Vol_The
//     + (qn0, T')_Sur_Neu_The + (h (T - T0), T')_Sur_Rob_The = 0
//
// holds for all test functions "T'".
//
// In the steady-state, the first term vanishes. Note that this problem is still
// well-posed, even without Dirichlet boundary conditions, thanks to the Robin
// condition (with "h > 0"). Indeed, the bilinear form is coercive: "a(T, T) =
// \int_Vol_The k |grad T|^2 + \int_Sur_Rob_The h T^2 >= C ||T||^2_H1".

Jacobian {
  { Name Vol;
    Case {
      { Region All; Jacobian Vol; }
    }
  }
  { Name Sur;
    Case {
      { Region All; Jacobian Sur; }
    }
  }
}

Integration {
  { Name Int;
    Case {
      { Type Gauss;
        Case {
          // We need to integrate on lines and triangles for the 2D model, and
          // on triangles and tetrahedra for the 3D model. We need sufficiently
          // many quadrature points to integrate quadratic functions (inner
          // product "(T , T')" of piecewise linear functions):
          { GeoElement Line; NumberOfPoints 2; }
          { GeoElement Triangle; NumberOfPoints 3; }
          { GeoElement Tetrahedron; NumberOfPoints 4; }
        }
      }
    }
  }
}

// Include definitions needed for the Link constraint below:
Include "radiator_common.pro";

Constraint {
  { Name T_The;
    Case {
      // Initial condition for transient solution (different constraint types
      // can be specified for each region within the "Case" -- compare with the
      // global "Type Assign" used in tutorial 1):
      If(AnalysisType == 1)
        { Region Vol_The; Type Init; Value T0[]; }
      EndIf
      // For a symmetric geometry (with position = 0.5 in "radiator.geo") a
      // homogeneous Neumann boundary condition on the Left and Right boundaries
      // is sufficient to ensure periodicity. Otherwise an explicit periodicity
      // condition is required. This can be achieved in GetDP with a "Link"
      // constraint, which links degrees of freedom from two regions ("Region"
      // and "RegionRef", geometrically mapped onto each other with "Function"),
      // with a coefficient. This requires the mesh to match on both regions.
      { Region Right; Type Link; RegionRef Left; Coefficient 1;
        // The function maps Region onto RegionRef, using the built-in "X[]",
        // "Y[]", "Z[]" coordinate functions to define the translation vector:
        Function Vector[X[] - N * bw, Y[], Z[]]; }
    }
  }
}

Group{
  Dom_H1_T_The =  Region[ {Vol_The, Sur_Neu_The, Sur_Rob_The} ];
}

FunctionSpace {
  { Name H1_T_The; Type Form0;
    BasisFunction {
      { Name sn; NameOfCoef Tn; Function BF_Node; Support Dom_H1_T_The;
        Entity NodesOf[All]; }
    }
    Constraint {
      { NameOfCoef Tn; EntityType NodesOf; NameOfConstraint T_The; }
    }
  }
}

Formulation {
  { Name Thermal_T; Type FemEquation;
    Quantity {
      { Name T; Type Local; NameOfSpace H1_T_The; }
    }
    Equation {
      // A "DtDof" term enables GetDP to handle time discretization
      // automatically (depending on the operation) in a "Resolution":
      Integral { DtDof [ rho_cp[] * Dof{T} , {T} ];
        In Vol_The; Jacobian Vol; Integration Int; }

      Integral { [ k[] * Dof{d T} , {d T} ];
        In Vol_The; Jacobian Vol; Integration Int; }

      Integral { [ qn0[] , {T} ];
        In Sur_Neu_The; Jacobian Sur; Integration Int; }

      // In GetDP equation terms with "Dof" must be linear with respect to the
      // value of the degree of freedom, so a "Dof" term with an affine part
      // such as "[ h[] * (Dof{T} - T0[]), {T} ]" is invalid and must be split
      // into two separate terms:
      Integral { [ h[] * Dof{T} , {T} ];
        In Sur_Rob_The; Jacobian Sur; Integration Int; }
      Integral { [ -h[] * T0[] , {T} ];
        In Sur_Rob_The; Jacobian Sur; Integration Int; }
    }
  }
}

Resolution {
  { Name The;
    System {
      { Name Sys_The; NameOfFormulation Thermal_T; }
    }
    Operation {
      If(AnalysisType == 1) // Transient
        // Initialize the solution:
        InitSolution[Sys_The];
        // Perform a time loop with an implicit Euler scheme, from time == 0 to
        // time == "tmax", with fixed time step "dt":
        TimeLoopTheta[0, tmax, dt, 1] {
          Generate[Sys_The]; Solve[Sys_The]; SaveSolution[Sys_The];
        }
        // More advanced time discretizations are available as built-in
        // resolution operations (see e.g. "TimeLoopAdaptive"). They could also be
        // implemented manually by replacing the "DtDof" term in the formulation
        // with the chosen approximation, e.g.
        //
        //   Integral { [ rho_cp[] * Dof{T} / dt, {T} ];
        //     In Vol_The; Integration Int; Jacobian Vol;  }
        //   Integral { [ - rho_cp[] * {T}[1] / dt, {T} ];
        //     In Vol_The; Integration Int; Jacobian Vol;  }
        //
        // for implicit Euler ("{T}[1]" denotes the value of "{T}" at the
        // previous step). A custom time integration loop can then be
        // constructed in the resolution with the "While[]" operation.
      Else // Steady-state
        // The "DtDof" term is automatically disregarded by GetDP if there is no
        // time loop, so we can generate and solve the steady-state system
        // directly:
        Generate[Sys_The]; Solve[Sys_The]; SaveSolution[Sys_The];
      EndIf
    }
  }
}

PostProcessing {
  { Name The; NameOfFormulation Thermal_T;
    Quantity {
      { Name T; Value{ Local{ [ {T} ]; In Vol_The; Jacobian Vol;} } }
      { Name q; Value{ Local{ [ -k[] * {d T} ]; In Vol_The; Jacobian Vol; } } }
    }
  }
}

PostOperation {
  { Name Map; NameOfPostProcessing The;
    Operation {
      Print[ T, OnElementsOf Vol_The, File "T.pos"];
      Print[ q, OnElementsOf Vol_The, File "q.pos"];
    }
  }
}

// Parameters shared by Gmsh and GetDP.

mm = 1e-3;

dxCore = 50 * mm;
dyCore = 100 * mm;
xCoil = 75 * mm;
dxCoil = 25 * mm;
dyCoil = 100 * mm;
rInt = 200 * mm;
rExt = 250 * mm;

DefineConstant[
  SymmetryType = {0, Choices{0="None", 1="X-axis", 2="Y-axis", 3="All"},
    Name "Parameters/0Symmetry"}
];

// Gmsh script describing the geometry of an electromagnet with an iron core and
// two coil sides, surrounded by air. The 3D setup is invariant along the
// out-of-plane direction "z" (the axis of the coils), so we model a 2D
// cross-section in the "(x, y)" plane. The two additional symmetries of this
// cross-section -- with respect to the X-axis (separating the two coil sides)
// and to the Y-axis (passing through the core) -- make it possible to model
// only a half (or even a quarter) of the cut. The "SymmetryType" parameter
// defined in "electromagnet_common.pro" selects which part is meshed.
//
// Note: tutorial 1 used the built-in Gmsh CAD kernel. Here we switch to
// OpenCASCADE, which offers boolean operations on solid primitives.

// We include a file with the dimensions as well as the symmetry type; the same
// file will be included in electromagnet.pro to guarantee consistency between
// Gmsh and GetDP:
Include "electromagnet_common.pro";

SetFactory("OpenCASCADE");
Rectangle(1) = {-dxCore, -dyCore, 0, 2 * dxCore, 2 * dyCore};
Rectangle(2) = {xCoil, -dyCoil / 2, 0, dxCoil, dyCoil};
Rectangle(3) = {-xCoil - dxCoil, -dyCoil / 2, 0, dxCoil, dyCoil};
Disk(4) = {0, 0, 0, rInt};
Disk(5) = {0, 0, 0, rExt};

If(SymmetryType == 0) // no symmetry
  // With the OpenCASCADE kernel the "Coherence" command is a shortcut for
  // "BooleanFragments" applied to all entities of the highest dimension. It
  // computes the intersections of all overlapping entities and replaces them
  // with a conformal (non-overlapping) set of entities. This is essential when
  // building a geometry from overlapping primitives with the OpenCASCADE
  // kernel: without it, entities would overlap and the resulting mesh would not
  // be conformal at region interfaces. After "Coherence" each original primitive
  // has been fragmented into the pieces that actually belong to each physical
  // region (core, coils, air, infinite shell), which is what the "Closest"
  // queries below rely on to retrieve them by position.
  Coherence;
Else
  // A bounding rectangle is built and intersected with the previously defined
  // primitives. The three non-trivial values of "SymmetryType" are handled
  // uniformly by the following formulas:
  //   - SymmetryType == 1: keep the upper half (y >= 0), X-axis is a symmetry
  //     plane, so "bot()" later carries the "n x h = 0" boundary condition;
  //   - SymmetryType == 2: keep the right half (x >= 0), Y-axis is a symmetry
  //     plane, so "left()" later carries the "b . n = 0" boundary condition;
  //   - SymmetryType == 3: keep the first quadrant (x >= 0 and y >= 0), both
  //     axes are symmetry planes, so both "bot()" and "left()" are defined
  //     below.
  d = 1.1 * rExt;
  x = (SymmetryType == 1) ? -d : 0; // symmetry w.r.t. X-axis ?
  dx = (SymmetryType == 1) ? 2 * d : d;
  y = (SymmetryType == 2) ? -d : 0; // symmetry w.r.t. Y-axis ?
  dy = (SymmetryType == 2) ? 2 * d : d;
  Rectangle(6) = {x, y, 0, dx, dy};
  BooleanIntersection{ Surface{6}; Delete; }{ Surface{1:5}; Delete; }
EndIf

// After boolean operations the original surface tags may have changed. The
// "Closest" command retrieves the entities closest to a given point, sorted by
// distance, which is a convenient way to identify the desired surfaces by their
// geometric location -- it is enough that the probe point lies inside the target
// region to retrieve it unambiguously:
core() = Closest{0, 0, 0}{ Surface{:}; };
air() = Closest{dxCore + mm, 0, 0}{ Surface{:}; };
airinf() = Closest{rInt + mm, 0, 0}{ Surface{:}; };
indr() = Closest{xCoil + mm, 0, 0}{ Surface{:}; };
indl() = Closest{-xCoil - mm, 0, 0}{ Surface{:}; };
Physical Surface("Core", 1) = core(0);
Physical Surface("Air", 2) = air(0);
Physical Surface("AirInf", 3) = airinf(0);
Physical Surface("CoilRight", 4) = indr(0);
If(SymmetryType == 0 || SymmetryType == 1)
  Physical Surface("CoilLeft", 5) = indl(0);
EndIf

bot() = {};
left() = {};
If(SymmetryType == 1 || SymmetryType == 3)
  bot() = Curve In BoundingBox{-d, -mm, -mm, d, mm, mm};
  Physical Curve("Bottom", 10) = bot();
EndIf
If(SymmetryType == 2 || SymmetryType == 3)
  left() = Curve In BoundingBox{-mm, -d, -mm, mm, d, mm};
  Physical Curve("Left", 11) = left();
EndIf

inf() = CombinedBoundary{ Surface{:}; };
inf() -= {bot(), left()};
Physical Curve("Inf", 12) = inf();

// Mesh size constraints:
DefineConstant[
  s = {1, Name "Parameters/}Global mesh size factor"}
];

MeshSize{:} = 12.5 * mm * s;
MeshSize{ PointsOf{ Surface{indr(0)}; } } = 5 * mm * s;
If(SymmetryType == 0 || SymmetryType == 1)
  MeshSize{ PointsOf{ Surface{indl(0)}; } } = 5 * mm * s;
EndIf
MeshSize{ PointsOf{ Surface{core(0)}; } } = 4 * mm * s;

// This model computes the static magnetic field produced by a DC current in an
// electromagnet. This corresponds to a magnetostatic physical model, obtained
// by combining the time-invariant Maxwell-Ampere equation ("curl h = js", with
// "h" the magnetic field and "js" the source current density) with the magnetic
// Gauss law ("div b = 0", with "b" the magnetic flux density) and the magnetic
// constitutive law ("b = mu h", with "mu" the magnetic permeability).
//
// Since "div b = 0", "b" can be derived from a vector magnetic potential "a",
// such that "b = curl a". Plugging this vector potential in Maxwell-Ampere's
// law and using the constitutive law leads to a generalized vector Poisson
// equation in terms of the magnetic vector potential: "curl(nu curl a) = js",
// where "nu = 1/mu" is the magnetic reluctivity.
//
// In the general (3D) case the vector potential is not unique, since for all
// scalar functions "f", "b = curl a = curl (a + grad f)". We will introduce
// explicit gauge conditions in tutorial 8 and 10 to ensure uniqueness. In the
// 2D setting however, with a source current density along the z-axis, the
// magnetic vector potential "a" is a vector field with a single non-zero
// component along the z-axis, i.e.:
//
//   a = (0, 0, az(x, y))
//
// which automatically removes gradient fields from the admissible space (the
// kernel reduces to constants, which a Dirichlet boundary condition is
// sufficient to fix).
//
// Magnetostatic fields expand to infinity. The corresponding boundary condition
// can be imposed rigorously by means of a geometrical transformation that maps
// a ring (a "shell") of finite elements to the complementary of its interior.
// As this is a mere geometric transformation, it is enough in the model
// description to attribute the special Jacobian "VolSphShell" to the ring
// region "AirInf". This Jacobian takes as arguments the inner ("rInt") and
// outer ("rOut") radii of the transformed ring region: to ensure consistency
// they are defined in a separate file, which is included in both the .geo and
// .pro file.  For a derivation and discussion of this transformation technique,
// see e.g. F. Henrotte et al., "Finite element modelling with transformation
// techniques", IEEE Transactions on Magnetics, 35(3), 1434-1437, 1999.

Include "electromagnet_common.pro";

// With this information, GetDP is able to correctly transform all quantities
// involved in the model. The model exhibits two symmetries: with respect to the
// X-axis and with respect to the Y-axis. The "SymmetryType" constant defined in
// "electromagnet_common.pro" allows to solve the problem with or without taking
// advantage of them.

Group {
  // Physical regions:
  Core = Region[ 1 ];
  Air = Region[ 2 ];
  AirInf = Region[ 3 ];
  CoilRight = Region[ 4 ];
  // Left coil only if no symmetry, or if single symmetry with respect to X-axis:
  CoilLeft = Region[ {} ];
  If(SymmetryType == 0 || SymmetryType == 1)
    CoilLeft += Region[ 5 ];
  EndIf
  Coils = Region[ {CoilLeft, CoilRight} ];
  // Bottom boundary only if symmetry with respect to X-axis:
  Bottom = Region[ {} ];
  If(SymmetryType == 1 || SymmetryType == 3)
    Bottom += Region[ 10 ];
  EndIf
  // Left boundary only if symmetry with respect to Y-axis:
  Left = Region[ {} ];
  If(SymmetryType == 2 || SymmetryType == 3)
    Left += Region[ 11 ];
  EndIf
  Inf = Region[ 12 ];

  // The abstract regions in this model have the following interpretation:
  //  - "Vol_Mag": overall domain
  //  - "Vol_S_Mag": region with imposed current source js
  //  - "Vol_NL_Mag": region with nonlinear magnetic constitutive law
  //  - "Vol_L_Mag": region with linear magnetic constitutive law
  //  - "Vol_Inf_Mag": region with infinite shell geometric transformation
  //  - "Sur_Neu_Mag": part of the boundary with non-homogeneous Neumann
  //    conditions
  Vol_Mag = Region[ {Air, AirInf, Core, Coils} ];
  Vol_S_Mag = Region[ Coils ];
  Vol_Inf_Mag = Region[ AirInf ];
  Sur_Neu_Mag = Region[ {} ]; // empty

  DefineConstant[
    NonlinearCore = {0, Choices {0, 1},
      Help "Is the magnetic law in the core nonlinear?",
      Name "Parameters/Nonlinear core"}
  ];

  If(NonlinearCore)
    Vol_NL_Mag = Region[ Core ];
  Else
    Vol_NL_Mag = Region[ {} ];
  EndIf
  // "-" to subtract "Vol_NL_Mag" from "Vol_Mag". "Vol_L_Mag" is therefore the
  // linear part of the magnetic domain: in the linear case it coincides with
  // the full domain, and in the nonlinear case it is everything except the
  // (nonlinear) core:
  Vol_L_Mag = Region[ {Vol_Mag, -Vol_NL_Mag} ];
}

// The weak formulation is derived in a similar way as for the electrostatic or
// thermal problems from previous tutorials. The main differences are that the
// unknown field is vector-valued, and that we want to handle the nonlinear case
// when nu is not constant. The weak formulation reads: find "a" such that
//
//   (curl(nu curl a), a')_Vol_Mag = (js, a')_Vol_S_Mag
//
// holds for all test functions "a'". After integration by parts it reads: find
// "a" such that
//
//   (nu curl a, curl a')_Vol_Mag + (n x (nu curl a), a')_Bnd_Vol_Mag =
//      (js, a')_Vol_S_Mag
//
// holds for all test functions "a'". In the electromagnet model the second
// (boundary) term vanishes:
//
//  - On the "Left" and "Inf" parts of the boundary, the normal component of "b"
//    must vanish by symmetry and decay, respectively. Imposing "a = 0" there
//    (homogeneous Dirichlet) enforces "b . n = 0" (since the tangential "a" is
//    zero, the flux of "b" through any surface patch on the boundary is zero by
//    Stokes' theorem). The boundary term vanishes because the test functions
//    "a'" are zero there -- GetDP automatically removes the corresponding test
//    functions when an "Assign" constraint is prescribed. For the "Left" part
//    of the 2D planar model, the "a = 0" condition can also be understood
//    directly: by symmetry with respect to the Y-axis the magnetic vector
//    potential has opposite signs on either side (the current density that
//    sources it is antisymmetric); continuity then forces "az = 0" on the axis
//    itself.
//
//  - On the "Bottom" boundary, the tangential component of "h" is zero by
//    symmetry, i.e. "n x h = n x (nu curl a) = 0". This is a homogeneous
//    Neumann condition: it is satisfied naturally (the boundary integrand is
//    zero).
//
// Note that depending on the value of "SymmetryType", the "Left" and/or
// "Bottom" region might be empty.
//
// We are thus eventually looking for functions a such that
//
//   (nu curl a, curl a')_Vol_Mag = (js, a')_Vol_S_Mag
//
// holds for all test functions "a'".
//
// If the magnetic constitutive law is nonlinear, i.e. if "nu" is not a constant
// but a function of "b = curl a", we need to linearize the formulation before
// we can solve it in GetDP. Several linearization methods are possible:
//
// 1) Newton-Raphson method - at iteration "k" we approximate
//
//   h(b_k) \approx h(b_k-1) + (dh/db)_k-1 (b_k - b_k-1)
//
// i.e.
//
//   (nu(curl a_k) curl a_k, curl a')
//     \approx (nu(curl a_k-1) curl a_k-1, curl a')
//           + (dh/db(curl a_k-1) curl a_k, curl a')
//           - (dh/db(curl a_k-1) curl a_k-1, curl a')
//
// 2) Picard method - at iteration "k" we approximate
//
//   (nu(curl a_k) curl a_k, curl a')
//     \approx (nu(curl a_k-1) curl a_k, curl a')

Function {
  DefineConstant[
    NewtonRaphson = {1, Choices {0="Picard", 1="Newton-Raphson"},
      Visible NonlinearCore,
      Name "Parameters/Linearization method"}
    Current = {1, Min 0.01, Max 100, Step 0.1,
      Name "Parameters/Current [A]"}
  ];

  mu0 = 4.e-7 * Pi;
  nu0 = 1 / mu0;
  nu [ Air ] = nu0;
  nu [ AirInf ] = nu0;
  nu [ Coils ] = nu0;

  If(!NonlinearCore)
    // Linear magnetic law in the core
    mur = 1000;
    nu [ Core ]  = 1 / (mur * mu0);
  Else
    // Nonlinear magnetic law in the core, defined by interpolating b-h curve
    // samples, provided in the following "data_h()" and "data_b()" lists:
    data_h() = {
      0.0000e+00, 5.5023e+00, 1.1018e+01, 1.6562e+01, 2.2149e+01, 2.7798e+01,
      3.3528e+01, 3.9363e+01, 4.5335e+01, 5.1479e+01, 5.7842e+01, 6.4481e+01,
      7.1470e+01, 7.8906e+01, 8.6910e+01, 9.5644e+01, 1.0532e+02, 1.1620e+02,
      1.2868e+02, 1.4322e+02, 1.6050e+02, 1.8139e+02, 2.0711e+02, 2.3932e+02,
      2.8028e+02, 3.3314e+02, 4.0231e+02, 4.9395e+02, 6.1678e+02, 7.8320e+02,
      1.0110e+03, 1.3257e+03, 1.7645e+03, 2.3819e+03, 3.2578e+03, 4.5110e+03,
      6.3187e+03, 8.9478e+03, 1.2802e+04, 1.8500e+04, 2.6989e+04, 3.9739e+04,
      5.9047e+04, 8.8520e+04, 1.3388e+05, 2.0425e+05, 3.1434e+05, 4.8796e+05,
      7.6403e+05
    };
    data_b() = {
      0.0000e+00, 5.0000e-02, 1.0000e-01, 1.5000e-01, 2.0000e-01, 2.5000e-01,
      3.0000e-01, 3.5000e-01, 4.0000e-01, 4.5000e-01, 5.0000e-01, 5.5000e-01,
      6.0000e-01, 6.5000e-01, 7.0000e-01, 7.5000e-01, 8.0000e-01, 8.5000e-01,
      9.0000e-01, 9.5000e-01, 1.0000e+00, 1.0500e+00, 1.1000e+00, 1.1500e+00,
      1.2000e+00, 1.2500e+00, 1.3000e+00, 1.3500e+00, 1.4000e+00, 1.4500e+00,
      1.5000e+00, 1.5500e+00, 1.6000e+00, 1.6500e+00, 1.7000e+00, 1.7500e+00,
      1.8000e+00, 1.8500e+00, 1.9000e+00, 1.9500e+00, 2.0000e+00, 2.0500e+00,
      2.1000e+00, 2.1500e+00, 2.2000e+00, 2.2500e+00, 2.3000e+00, 2.3500e+00,
      2.4000e+00
    };

    // We first compute a list of reluctivity values for each b-h sample (and
    // fix the first value, when "h = b = 0"):
    data_nu = data_h() / data_b();
    data_nu(0) = data_nu(1);

    // We then create a function "nu(|b|^2)" by spline interpolation. Note
    // that we interpolate against "|b|^2" (the squared norm) rather than "|b|"
    // itself: this avoids computing a square root at each evaluation, and
    // makes the derivative "d(nu)/d(|b|^2)" (needed for Newton-Raphson below)
    // directly available:
    data_b2_nu = ListAlt[data_b()^2, data_nu()];
    nu[ Core ] = InterpolationAkima[ SquNorm[$1] ]{ data_b2_nu() };

    // "nu[]" is a piecewise-defined function: its left-hand side above is
    // "nu[ Core ]", meaning that for elements in "Core" the formula on the
    // right is used. When called inside the Formulation as "nu[{d a}]", the
    // argument "{d a}" (= curl a = b) is substituted for "$1" in the formula.
    // With no argument -- "nu[]" as it appears in the linear branch and in
    // the Formulation in the linear regions Vol_L_Mag -- the formula is
    // simply evaluated as written (here, a constant).

    // The function "nu[]" is expected to take b (a vector value) as its first
    // (and only) argment, "$1". As we will also see below in the Resolution,
    // the $ sign identifies runtime quantities (function arguments, global
    // variables, user runtime variables), whose values are not yet known when
    // GetDP parses the ".pro" file. Here "nu[]" computes the square of the norm
    // of its first argument ("SquNorm[$1]") and passes it as argument to the
    // built-in "InterpolationAkima[]" function. In addition to its single
    // argument, "InterpolationAkima[]" takes a list of values as parameters -
    // provided after the arguments between braces - here the list of values to
    // interpolate.

    // The derivative of "nu" with respect to "|b|^2" is directly available
    // through the built-in "dInterpolationAkima[]" function:
    dnudb2[ Core ] = dInterpolationAkima[ SquNorm[$1] ]{ data_b2_nu() };

    // The components of the Jacobian "dh/db", for "h = nu b", are
    //
    //   (dh/db)_ij = dh_i/db_j = nu delta_ij + dnu/db_j b_i
    //
    // which, in matrix form, reads (with "I" the identity tensor -- i.e. the
    // 3x3 identity matrix in Cartesian components):
    //
    //   dh/db = I nu + b (grad nu)^T
    //
    // With "nu = nu(|b|^2)", we have "grad nu = 2 nu'(|b|^2) b", and thus
    //
    //   dh/db = I nu + 2 nu'(|b|^2) b b^T
    //
    dhdb[ Core ] = TensorDiag[1, 1, 1] * nu[$1#1] + 2 * dnudb2[#1] *
      SquDyadicProduct[#1];

    // A small optimization has been applied in the above function definition to
    // avoid multiple evaluations of argument "$1" at runtime: "$1#1" stores the
    // value of the argument in a register ("#1"), which is then reused twice.

    // We also define absolute and relative tolerances on the residual, as well
    // as the maximum number of iterations, to control the iterative loop:
    NLTolAbs = 1e-10;
    NLTolRel = 1e-6;
    NLIterMax = 20;
  EndIf

  // Number of turns in the coil:
  NbTurns = 1000;

  // Current density in the inductor, along the z-axis:
  js[ CoilRight ] = Vector[0, 0, -NbTurns * Current / (dxCoil * dyCoil)];
  js[ CoilLeft ] = Vector[0, 0, NbTurns * Current / (dxCoil * dyCoil)];
}

Constraint {
  // Homogenous Dirichlet condition on the "Left" and "Inf" parts of the
  // boundary. When constraint type is given, "Type Assign" is assumed by
  // default:
  { Name a_Mag_2D;
    Case {
      { Region Left; Value 0.; }
      { Region Inf; Value 0.; }
    }
  }
}

// The magnetic vector potential has a single non-zero component along the
// z-axis: "a = (0, 0, az(x, y))". This is reflected in the "Type",
// "BasisFunction" and "Entity" specified in the "Hcurl_a_Mag_2D" FunctionSpace
// below: the vector potential is a "perpendicular 1-form" interpolated with
// "BF_PerpendicularEdge" basis functions associated to nodes of the mesh. With
// this information GetDP is able to correctly apply geometrical transformations
// and differentiate the vector potential in the Formulation and PostProcessing
// terms.

Group {
  Dom_Hcurl_a_Mag_2D = Region[ {Vol_Mag, Sur_Neu_Mag} ];
}

FunctionSpace {
  { Name Hcurl_a_Mag_2D; Type Form1P;
    BasisFunction {
      { Name se; NameOfCoef ae; Function BF_PerpendicularEdge;
        Support Dom_Hcurl_a_Mag_2D; Entity NodesOf[ All ]; }
    }
    Constraint {
      { NameOfCoef ae; EntityType NodesOf; NameOfConstraint a_Mag_2D; }
    }
  }
}

// The special Jacobian "VolSphShell" for the infinite elements takes 2
// parameters in this case, "rInt" and "rExt", which represent the inner and
// outer radii of the transformed ring region. The value of the parameters must
// match those used in the geometrical description of the model, so to ensure
// consistency there are defined in a separate file, which is included in both
// the .geo and .pro file:
Include "electromagnet_common.pro";

Jacobian {
  { Name Vol;
    Case {
      // Use the special infinite ring Jacobian in "Vol_Inf_Mag"...
      { Region Vol_Inf_Mag; Jacobian VolSphShell {rInt, rExt}; }
      // ... and the standard "Vol" Jacobian everywhere else:
      { Region All; Jacobian Vol; }
    }
  }
}

Integration {
  { Name Int;
    Case {
      { Type Gauss;
        Case {
          { GeoElement Triangle; NumberOfPoints 1; }
        }
      }
    }
  }
}

Formulation {
  { Name Magnetostatics_a_2D; Type FemEquation;
    Quantity {
      { Name a; Type Local; NameOfSpace Hcurl_a_Mag_2D; }
    }
    Equation {
      // Since a is a "Form1P", the differential operator "d" means "curl":
      Integral { [ nu[] * Dof{d a} , {d a} ];
        In Vol_L_Mag; Jacobian Vol; Integration Int; }

      If(NonlinearCore && NewtonRaphson)
        // Newton-Raphson linearization:
        Integral { [ nu[{d a}] * {d a} , {d a} ];
          In Vol_NL_Mag; Jacobian Vol; Integration Int; }
        Integral { [ dhdb[{d a}] * Dof{d a} , {d a} ];
          In Vol_NL_Mag; Jacobian Vol; Integration Int; }
        Integral { [ - dhdb[{d a}] * {d a} , {d a} ];
          In Vol_NL_Mag; Jacobian Vol; Integration Int; }

        // [The following block can be skipped on a first reading.]
        //
        // Note that this implementation of the Newton-Raphson method can be
        // simplified (and accelerated) by directly defining the incremental
        // contribution to the system matrix through a "JacNL[]" term.
        //
        // With "nu = nu(|b|^2)", we have
        //
        //   h(b_k) \approx h(b_k-1) + (dh/db)_k-1 (b_k - b_k-1)
        //      = nu_k-1 b_k-1 + (I nu + 2 nu'(|b|^2) b b^T)_k-1 (b_k - b_k-1)
        //      = nu_k-1 b_k + 2 (nu'(|b|^2) b b^T)_k-1 (b_k - b_k-1)
        //
        // Defining the function
        //
        //   dhdb_NL[] = 2 * dnudb2[$1#1] * SquDyadicProduct[#1]
        //
        // we can then rewrite the Newton-Raphson linearization as follows:
        //
        //   Integral { [ nu[{d a}] * Dof{d a} , {d a} ]; // note the Dof{}!
        //     In Vol_NL_Mag; Jacobian Vol; Integration Int; }
        //   Integral { JacNL [ dhdb_NL[{d a}] * Dof{d a} , {d a} ];
        //     In Vol_NL_Mag; Jacobian Vol; Integration Int; }
        //
        // The system should then be assembled and solved with "GenerateJac[]"
        // and "SolveJac[]" (instead of "Generate[]" and "Solve[]"), which will
        // automatically handle the increment "b_k - b_k-1".

      ElseIf(NonlinearCore)
        // Picard linearization:
        Integral { [ nu[{d a}] * Dof{d a}, {d a} ];
          In Vol_NL_Mag; Jacobian Vol; Integration Int; }
      EndIf

      Integral { [ - js[] , {a} ];
        In Vol_S_Mag; Jacobian Vol; Integration Int; }
    }
  }
}

Resolution {
  { Name Mag;
    System {
      { Name Sys_Mag; NameOfFormulation Magnetostatics_a_2D; }
    }
    Operation {
      // Generate matrix "Mat_0" and right-hand-side "rhs_0" of the linear
      // system "Sys_Mag":
      Generate[Sys_Mag];

      // Solve "Mat_0 x_0 = rhs_0", i.e. compute "x_0 = Mat_0^-1 rhs_0":
      Solve[Sys_Mag];

      If(NonlinearCore)
        // Re-generate the system: "Generate[]" re-assembles a matrix "Mat_1"
        // and right-hand side "rhs_1", using the latest available solution to
        // evaluate the nonlinear coefficients (here "nu(|b_0|^2)" with "b_0 =
        // curl a_0"). The source "js" is unchanged since it is prescribed. The
        // resulting "Mat_1" thus reflects the new operating point of the
        // nonlinear law:
        Generate[Sys_Mag];

        // Compute residual "rhs_1 - Mat_1 x_0" and store its 2-norm in the user
        // runtime variable "$res0":
        GetResidual[Sys_Mag, $res0];

        // Initialize runtime variables to track the residual and the iteration
        // count, then print out the absolute and relative residual. User
        // runtime variables (prefixed by "$") must be declared and updated
        // through an "Evaluate[]" operation, rather than assigned directly like
        // ordinary ".pro" symbols (which are constants evaluated once at parse
        // time). "Evaluate[]" is the mechanism that lets the resolution create
        // or update such variables during the operation sequence:
        Evaluate[ $res = $res0, $iter = 0 ];
        Print[{$iter, $res, $res / $res0},
          Format "Residual %03g: abs %14.12e rel %14.12e"];

        // Iterate until convergence:
        While[$res > NLTolAbs && $res / $res0 > NLTolRel &&
              $res / $res0 <= 1 && $iter < NLIterMax]{
          // Solve "Mat_k x_k = rhs_k", generate "Mat_k+1" and "rhs_k+1", and
          // compute the residual "rhs_k+1 - Mat_k+1 x_k":
          Solve[Sys_Mag]; Generate[Sys_Mag]; GetResidual[Sys_Mag, $res];

          Evaluate[ $iter = $iter + 1 ];
          Print[{$iter, $res, $res / $res0},
            Format "Residual %03g: abs %14.12e rel %14.12e"];
        }
        // Note that this explicit implementation of the iterative loop can be
        // replaced by the built-in "IterativeLoop[]" or "IterativeLoopN[]"
        // resolution operations, which offer several refinements. The simplest
        // implementation would be (this requires using "JacNL[]" for
        // Newton-Raphson, as explained above):
        //
        //   IterativeLoop[NLIterMax, NLTolRel, 1] {
        //     GenerateJac[Sys_Mag]; SolveJac[Sys_Mag];
        //   }
      EndIf

      SaveSolution[Sys_Mag];
    }
  }
}

PostProcessing {
  { Name Mag; NameOfFormulation Magnetostatics_a_2D;
    Quantity {
      { Name a;
        Value {
          Term { [ {a} ]; In Vol_Mag; Jacobian Vol; }
        }
      }
      { Name az; // z-component of the vector potential
        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 js;
        Value {
          Term { [ js[] ]; In Vol_S_Mag; Jacobian Vol; }
        }
      }
    }
  }
}

PostOperation {
  { Name Map; NameOfPostProcessing Mag;
    Operation {
      Print[ a, OnElementsOf Vol_Mag, File "a.pos" ];
      Print[ js, OnElementsOf Vol_S_Mag, File "js.pos" ];
      Print[ az, OnElementsOf Vol_Mag, File "az.pos" ];
      Print[ b, OnElementsOf Vol_Mag, File "b.pos" ];
      Print[ b, OnLine{{mm, mm, 0}{rInt, mm, 0}}{50}, File "cutb.pos" ];
    }
  }
}

// Parameters shared by Gmsh and GetDP.

mm = 1e-3;

dxCore = 50 * mm;
dyCore = 100 * mm;
xCoil = 75 * mm;
dxCoil = 25 * mm;
dyCoil = 100 * mm;
rInt = 200 * mm;
rExt = 250 * mm;

// Gmsh script describing the geometry of the electromagnet. This is a pared
// down version of the geometry from tutorial 3, where all symmetries are
// taken advantage of.

Include "electromagnet_common.pro";

SetFactory("OpenCASCADE");
Rectangle(1) = {-dxCore, -dyCore, 0, 2 * dxCore, 2 * dyCore};
Rectangle(2) = {xCoil, -dyCoil / 2, 0, dxCoil, dyCoil};
Rectangle(3) = {-xCoil - dxCoil, -dyCoil / 2, 0, dxCoil, dyCoil};
Disk(4) = {0, 0, 0, rInt};
Disk(5) = {0, 0, 0, rExt};
d = 1.1 * rExt;
Rectangle(6) = {0, 0, 0, d, d};
BooleanIntersection{ Surface{6}; Delete; }{ Surface{1:5}; Delete; }

core() = Closest{0, 0, 0}{ Surface{:}; };
air() = Closest{dxCore + mm, 0, 0}{ Surface{:}; };
airinf() = Closest{rInt + mm, 0, 0}{ Surface{:}; };
indr() = Closest{xCoil + mm, 0, 0}{ Surface{:}; };
Physical Surface("Core", 1) = core(0);
Physical Surface("Air", 2) = air(0);
Physical Surface("AirInf", 3) = airinf(0);
Physical Surface("CoilRight", 4) = indr(0);

bot() = Curve In BoundingBox{-mm, -mm, -mm, d, mm, mm};
left() = Curve In BoundingBox{-mm, -mm, -mm, mm, d, mm};
inf() = CombinedBoundary{ Surface{:}; };
inf() -= {bot(), left()};
Physical Curve("Bottom", 10) = bot();
Physical Curve("Left", 11) = left();
Physical Curve("Inf", 12) = inf();

// Mesh size constraints:
DefineConstant[
  s = {1, Name "Parameters/}Global mesh size factor"}
];

MeshSize{:} = 12.5 * mm * s;
MeshSize{ PointsOf{ Surface{indr(0)}; } } = 5 * mm * s;
MeshSize{ PointsOf{ Surface{core(0)}; } } = 4 * mm * s;

// Refine to capture the skin depth:
skin() = Curve In BoundingBox{dxCore - mm, -mm, -mm, dxCore + mm, d, mm};
MeshSize{ PointsOf{ Curve{skin(0)}; } } = 0.5 * mm * s;

// This model computes the magnetic field produced by an AC current in the same
// electromagnet as the one considered in tutorial 3 with a DC current.
//
// This corresponds to a magneto-quasistatic (or "magnetodynamic") physical
// model, obtained by combining the Maxwell-Ampere equation where the
// displacement current is neglected ("curl h = j", with "h" the magnetic field
// and "j" the current density), the Faraday equation ("curl e = -\partial_t b",
// with "e" the electric field and "b" the magnetic flux density), the magnetic
// Gauss law ("div b = 0"), the magnetic constitutive law ("b = mu h", with "mu"
// the magnetic permeability) and Ohm's law ("j = sigma e", with "sigma" the
// electric conductivity). In the coil of the electromagnet we assume that eddy
// currents can be neglected and we simply impose a source current density "j =
// js". This is a standard modeling choice for "stranded" coils made of many
// thin insulated wires: each wire is thin compared to the skin depth, which
// prevents eddy loops from forming inside the coil. The assumption would no
// longer be valid for a massive ("solid") conductor, or for a coil connected to
// a voltage source (see tutorial 7), potentially through a circuit (see
// tutorial 8).
//
// We can introduce a vector magnetic potential "a" such that "b = curl a"
// and "e = -\partial_t a", that satisfies both Gauss' law and Faraday's
// equation. We call this the "modified" vector potential because the usual
// definition "e = -\partial_t a - grad v" has been simplified by implicitly
// fixing the electric scalar potential "v" to zero -- equivalently, we have
// chosen a gauge that absorbs "grad v" into "a". This is a valid choice
// whenever the current is imposed (as here) rather than driven by a voltage.
// In tutorial 7 we will reintroduce "v" to allow prescribing a voltage drop
// or computing it from an imposed current, using the so-called "a-v"
// formulation.
//
// Plugging this vector potential into Ampere's equation and using the
// constitutive laws leads to
//
//   curl(nu curl a) + sigma \partial_t a = js
//
// where "nu = 1/mu" is the magnetic reluctivity. The right-hand side is the
// imposed (source) current density; on the left-hand side, "sigma \partial_t
// a" is the induced (eddy) current density, due to the time-varying magnetic
// flux. The total current density in a conducting region is thus
// "j = -sigma \partial_t a" (purely induced) and in a source region it is
// "j = js" (purely imposed). See the PostProcessing section below, which
// defines the full "j" piecewise on the two kinds of regions.
//
// We consider both time-domain (transient) and frequency-domain (steady-state)
// resolutions, with a sinusoidal source current density at fixed frequency f:
//
//   js(t) = (0, 0, Js sin(2 pi f t))
//
// As in tutorial 3, with a source current density along the z-axis, the
// magnetic vector potential a is a vector field with a single non-zero
// component along the z-axis, i.e.,
//
//   a(x, y, t) = (0, 0, az(x, y, t))
//
// and we can use the same "Form1P" function space with the same
// "BF_PerpendicularEdge" basis functions.
//
// In the time-domain an implicit Euler time stepping scheme is used, whereas in
// the frequency domain, we solve for the phasor "A(x, y)" such that
//
//   a(x, y, t) = Re( A(x, y) exp(i 2 pi f t) )
//
// Note that with this definition of phasors in terms of peak values, a factor
// 1/2 appears in time-averaged quadratic quantities such as electromagnetic
// power (due to the time-averaged cosine squared over a period). This factor
// dissapears for phasors defined in terms of RMS values, or for the calculation
// of instantaneous power in the time-domain.

Include "electromagnet_common.pro";

Group {
  // Physical regions:
  Core = Region[ 1 ];
  Air = Region[ 2 ];
  AirInf = Region[ 3 ];
  CoilRight = Region[ 4 ];
  Bottom = Region[ 10 ];
  Left = Region[ 11 ];
  Inf = Region[ 12 ];

  // Abstract regions:
  //  - "Vol_Mag": full volume domain
  //  - "Vol_S_Mag": region with imposed current source js
  //  - "Vol_C_Mag": conducting regions with eddy currents
  //  - "Vol_Inf_Mag": region with infinite ring geometric transformation
  //  - "Sur_Neu_Mag": part of the boundary with non-homogeneous Neumann
  //    conditions
  Vol_Mag     = Region[ {Air, AirInf, Core, CoilRight} ];
  Vol_S_Mag   = Region[ CoilRight ];
  Vol_C_Mag   = Region[ Core ];
  Vol_Inf_Mag = Region[ AirInf ];
  Sur_Neu_Mag = Region[ {} ]; // empty
}

Function {
  DefineConstant[
    AnalysisType = {1, Choices{0="Time-domain", 1="Frequency-domain"},
      Name "Parameters/Analysis type"}
    // In addition to a 2D model (invariant along the z-axis) we also consider
    // an axisymmetric model (invariant by rotation around the y-axis). This is
    // appropriate when the 3D device is itself a solid of revolution -- as with
    // a cylindrical coil. It would not be appropriate for, e.g., a
    // rotating-machine cross-section, which is only piecewise-invariant by
    // rotation:
    Axisymmetric = {0, Choices{0, 1},
      Name "Parameters/Axisymmetric model?"}
    Current = {1, Min 0.01, Max 100, Step 0.1,
      Name "Parameters/Current [A]"}
    f = {5, Min 1, Max 1000, Step 1,
      Name "Parameters/Frequency [Hz]"}
    mur = {1000, Min 1, Max 5000, Step 1,
      Name "Parameters/Core relative permeability"}
    sigma = {1e6, Min 1, Max 1e6, Step 100,
      Name "Parameters/Core conductivity [S m^-1]"}
  ];

  T = 1 / f;
  dt = T / 20;
  tmax = 2 * T;

  mu0 = 4.e-7 * Pi;
  nu0 = 1 / mu0;
  nu [ Region[{Air, AirInf, CoilRight}] ] = nu0;
  nu [ Core ]  = 1 / (mur * mu0);

  sigma[ Core ] = sigma;

  SkinDepth = DefineNumber[ 1e3 * Sqrt[2 / (2 * Pi * f * mur * mu0 * sigma)],
    Name "Parameters/Skin depth [mm]", ReadOnly];

  NbTurns = 1000;

  // The built-in function "F_Sin_wt_p[]{2 * Pi * f, 0}" will be interpreted as
  // "Sin[2 * Pi * f * $Time]" in the time domain, where "$Time" denotes the
  // current value of time, and as "Complex[0, -1]" in the frequency domain.
  //
  // The syntax "Complex[Re, Im]" builds a complex number from its real and
  // imaginary parts: "Complex[0, 1]" is the imaginary unit "i", and "Complex[0,
  // -1] = -i" is the phasor of "sin(omega t)" under the convention "x(t) = Re[X
  // exp(i omega t)]" (indeed, "sin(omega t) = Re[-i exp(i omega t)]").
  js[ CoilRight ] = Vector[0, 0, -NbTurns * Current / (dxCoil * dyCoil)] *
    F_Sin_wt_p[]{2 * Pi * f, 0};

  // Coefficient for computing electromagnetic power: 1 for time-domain
  // (instantaneous) power or for RMS phasors, 1/2 for peak-valued phasors:
  CoefPower = (AnalysisType == 0) ? 1 : 0.5;
}

Constraint {
  { Name a_Mag_2D;
    Case {
      { Region Left; Value 0.; }
      { Region Inf; Value 0.; }
    }
  }
}

Group {
  Dom_Hcurl_a_Mag_2D = Region[ {Vol_Mag, Sur_Neu_Mag} ];
}

FunctionSpace {
  { Name Hcurl_a_Mag_2D; Type Form1P;
    BasisFunction {
      { Name se; NameOfCoef ae; Function BF_PerpendicularEdge;
        Support Dom_Hcurl_a_Mag_2D; Entity NodesOf[ All ]; }
    }
    Constraint {
      { NameOfCoef ae; EntityType NodesOf; NameOfConstraint a_Mag_2D; }
    }
  }
}

Jacobian {
  { Name Vol;
    Case {
      // Solving an axisymmetric problem in GetDP is achieved by specifying the
      // axisymmetric variants of the Jacobian transformations, i.e. replacing
      // "Vol" with "VolAxi" and "VolSphShell" with "VolAxiSphShell":
      If(Axisymmetric)
        { Region Vol_Inf_Mag; Jacobian VolAxiSphShell {rInt, rExt}; }
        { Region All; Jacobian VolAxi; }
        // The axisymmetric Jacobians assume that the domain is in the z == 0
        // plane, and that the rotation axis is the y-axis.
      Else
        { Region Vol_Inf_Mag; Jacobian VolSphShell {rInt, rExt}; }
        { Region All; Jacobian Vol; }
      EndIf
    }
  }
}

Integration {
  { Name Int;
    Case {
      { Type Gauss;
        Case {
          // With pieceweise linear basis functions we must be able to integrate
          // quadratic polynomials in "Vol_C_Mag", hence we use 4 integration
          // points instead of 1 as in tutorial 3:
          { GeoElement Triangle; NumberOfPoints 3; }
        }
      }
    }
  }
}

// The weak formulation is derived in a similar way as for the magnetostatic
// problem from tutorial 3: find "a" such that
//
//   (curl(nu curl a), a')_Vol_Mag + (sigma \partial_t a, a')_Vol_C_Mag
//     = (js, a')_Vol_S_Mag
//
// holds for all test functions "a'". After integration by parts, and with the
// same boundary conditions as in tutorial 3, it reads: find a such that
//
//   (nu curl a, curl a')_Vol_Mag + (sigma \partial_t a, a')_Vol_C_Mag
//     = (js, a')_Vol_S_Mag
//
// holds for all test functions "a'".
Formulation {
  { Name Magnetoquasistatics_a_2D; Type FemEquation;
    Quantity {
      { Name a; Type Local; NameOfSpace Hcurl_a_Mag_2D; }
    }
    Equation {
      Integral { [ nu[] * Dof{d a} , {d a} ];
        In Vol_Mag; Jacobian Vol; Integration Int; }

      Integral { DtDof [ sigma[] * Dof{a} , {a} ];
        In Vol_C_Mag; Jacobian Vol; Integration Int; }

      Integral { [ - js[] , {a} ];
        In Vol_S_Mag; Jacobian Vol; Integration Int; }
    }
  }
}

Resolution {
  { Name Mag;
    System {
      { Name Sys_Mag; NameOfFormulation Magnetoquasistatics_a_2D;
        If(AnalysisType == 1)
          // With "Type Complex" and "Frequency f", GetDP transforms the "DtDof"
          // operator as a multiplication by "Complex[0, 2 * Pi * f]", and
          // handles all fields as complex-valued phasors:
          Type Complex; Frequency f;
        EndIf
      }
    }
    Operation {
      If(AnalysisType == 1)
        Generate[Sys_Mag]; Solve[Sys_Mag]; SaveSolution[Sys_Mag];
      Else
        // If no "Init" constraint is given, "InitSolution[]" initializes the
        // solution to 0:
        InitSolution[Sys_Mag];
        TimeLoopTheta[0, tmax, dt, 1] {
          Generate[Sys_Mag]; Solve[Sys_Mag]; SaveSolution[Sys_Mag];
        }
      EndIf
    }
  }
}

PostProcessing {
  { Name Mag; NameOfFormulation Magnetoquasistatics_a_2D;
    Quantity {
      { Name a;
        Value {
          Term { [ {a} ]; In Vol_Mag; Jacobian Vol; }
        }
      }
      { Name az; // z-component of the vector potential
        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; }
        }
      }
      // The current density can be defined piecewise: as the induced (eddy)
      // current density in the conducting regions, and as the source current
      // density in the inductor regions:
      { Name j;
        Value {
          Term { [ - sigma[] * Dt[{a}] ]; In Vol_C_Mag; Jacobian Vol; }
          Term { [ js[] ]; In Vol_S_Mag; Jacobian Vol; }
        }
      }
      // We can compute Joule losses by integrating "1 / sigma |j|^2" over the
      // domain (beware of time-averaging with phasors, cf. the definition of
      // "CoefPower" above), which can be done by defining an "Integral"
      // post-processing quantity, again piecewise. In the conducting regions
      // the current is induced, "j = -sigma \partial_t a", so "|j|^2 / sigma =
      // sigma |\partial_t a|^2". In the source regions the current is imposed,
      // and we keep "|js|^2 / sigma" directly:
      { Name JouleLosses;
        Value {
          Integral { [ CoefPower * sigma[] * SquNorm[Dt[{a}]] ];
            In Vol_C_Mag; Jacobian Vol; Integration Int; }
          Integral { [ CoefPower / sigma[] * SquNorm[js[]] ];
            In Vol_S_Mag; Jacobian Vol; Integration Int; }
        }
      }
    }
  }
}

PostOperation {
  { Name Map; NameOfPostProcessing Mag;
    Operation {
      // In the time domain, the output files will contain the values at each
      // time step. In the frequency domain, they will contain two values: the
      // real and imaginary parts.
      Print[ a, OnElementsOf Vol_Mag, File "a.pos" ];
      Print[ az, OnElementsOf Vol_Mag, File "az.pos" ];
      Print[ b, OnElementsOf Vol_Mag, File "b.pos" ];
      Print[ b, OnLine{{mm, mm, 0}{rInt, mm, 0}}{50}, File "cutb.pos" ];
      Print[ j, OnElementsOf Vol_S_Mag, File "js.pos" ];
      Print[ j, OnElementsOf Vol_C_Mag, File "j.pos" ];
      // Integral post-processing quantities do not make sense on elements on
      // the mesh, as they are global values. Here we evaluate the Joule losses
      // over the conducting regions:
      Print[ JouleLosses[Vol_C_Mag], OnGlobal, Format Table, File "losses.txt" ];
    }
  }
}

DefineConstant[
  dim = {2, Choices{2="2D", 3="3D"},
    Name "Parameters/0Model dimension"}
  Wx = {4, Min 1, Max 100, Step 1,
    Name "Parameters/Waveguide width [cm]"}
  Wz = {3, Min 1, Max 100, Step 1, Visible (dim == 3),
    Name "Parameters/Waveguide height [cm]"}
];

// back to meters:
Wx = Wx / 100;
Wz = Wz / 100;

SetFactory("OpenCASCADE");

Include "full_wave_common.pro";

DefineConstant[
  L = {9, Min 1, Max 100, Step 1,
    Name "Parameters/Waveguide length [cm]"}
  res = {0.4, Min 1e-2, Max 50e-1, Step 1e-2,
    Name "Parameters/Mesh size [cm]"}
  airRadius = {6, Min 0, Max 100, Step 1,
    Name "Parameters/Air radius [cm]"}
];

// Convert dimensions back to meters:
L = L / 100;
res = res / 100;
airRadius = airRadius / 100;

// Inner part of waveguide:
xMinIn = -Wx / 2;
yMinIn = -L;
zMinIn = -Wz / 2;

// Wall thickness:
t = 2e-3;

// Outer part of waveguide:
xMinOut = xMinIn - t;
yMinOut = yMinIn - t;
zMinOut = zMinIn - t;
WxOut = Wx + 2 * t;
WzOut = Wz + 2 * t;

If(dim == 2)
  // 2D geometry in z == 0 plane:
  zMinIn = 0;
  zMinOut = 0;
  If(airRadius > 0)
    Disk(1) = {0, 0, 0, airRadius};
    Rectangle(2) = {xMinOut, yMinOut, zMinOut, WxOut, L + t};
    vAir() = BooleanDifference{ Surface{1}; Delete; }{ Surface{2}; Delete; };
  Else
    vAir() = {};
  EndIf
  vGuideIn = news; // next available surface tag
  Rectangle(vGuideIn) = {xMinIn, yMinIn, zMinIn, Wx, L};
Else
  // 3D geometry:
  If(airRadius > 0)
    Sphere(1) = {0, 0, 0, airRadius};
    Box(2) = {xMinOut, yMinOut, zMinOut, WxOut, L + t, WzOut};
    vAir() = BooleanDifference{ Volume{1}; Delete; }{ Volume{2}; Delete; };
  Else
    vAir() = {};
  EndIf
  vGuideIn = newv; // next available volume tag
  Box(vGuideIn) = {xMinIn, yMinIn, zMinIn, Wx, L, Wz};
EndIf

If(airRadius > 0)
  // Make all entities conformal (see tutorial 3):
  Coherence;
EndIf

// Retrieve surfaces for boundary conditions using a bounding box search, with a
// small tolerance:
tol = 1e-6;

// Excitation port at yMinIn:
sPort() = GeoEntity{dim - 1} In BoundingBox {
  xMinIn - tol, yMinIn - tol, zMinIn - tol,
  xMinIn + Wx + tol, yMinIn + tol, zMinIn + Wz + tol};

// Output port at yMinIn + L:
sOut() = GeoEntity{dim - 1} In BoundingBox {
  xMinIn - tol, yMinIn + L - tol, zMinIn - tol,
  xMinIn + Wx + tol, yMinIn + L + tol, zMinIn + Wz + tol};

// Waveguide walls excluding ports for PEC conditions:
sPEC() = GeoEntity{dim - 1} In BoundingBox {
  xMinOut - tol, yMinOut - tol, zMinOut - tol,
  xMinOut + WxOut + tol, yMinOut + L + t + tol, zMinOut + WzOut + tol};
sPEC() -= {sPort(), sOut()};

// Infinite boundary for ABC. "CombinedBoundary" returns the exterior boundary
// curves (resp. surfaces) of a set of surfaces (resp. volumes), combining
// (i.e. removing) internal boundaries shared by adjacent entities.  Here it
// returns the outermost boundary of the entire model:
sInf() = CombinedBoundary{ GeoEntity{dim}{:}; };
sInf() -= {sPort(), sOut(), sPEC()};

// Uniform mesh size constraint:
MeshSize{:} = res;

// Physical groups:
Physical GeoEntity{dim}("Waveguide", 1000) = {vGuideIn};
Physical GeoEntity{dim}("Air", 1001) = {vAir()};
Physical GeoEntity{dim - 1}("PEC", 2000) = {sPEC()};
Physical GeoEntity{dim - 1}("Port", 2001) = {sPort()};
Physical GeoEntity{dim - 1}("Inf", 2002) = {sInf()};

// In this tutorial we introduce the full-wave Maxwell problem in the frequency
// domain, with lowest order Whitney edge elements to discretize the vector
// electric field in both 2D and 3D. The vector-valued basis functions are
// associated with the edges of the mesh, i.e. the electric field "e(x, y, z)"
// is expanded as
//
//   e(x, y, z) = Sum_{k in E} ee_k se_k(x, y, z)
//
// where the sum runs over all the edges of the mesh and the vector-valued edge
// basis functions "se_k(x, y, z)" are tangentially continuous across the mesh
// elements. The coefficients "ee_k" (i.e. the values of the degrees of
// freedom), correspond to the circulations of the electric field on the edges
// of the mesh. This identification between DoFs and edge circulations is a
// defining property of the lowest-order Whitney edge element: each basis
// function "se_k" is designed so that its line integral along edge "k" is
// equal to 1, and its line integral along any other edge is equal to 0. By
// linearity, "\int_{edge k} e . dl = ee_k", so each coefficient "ee_k" is
// literally the circulation of "e" along edge "k".
//
// Edge elements are one member of a larger family known as Whitney
// elements. Nodal "BF_Node" functions (associated with nodes, one pointwise
// scalar DoF per node) discretize scalar field and belong to "H1"; edge
// "BF_Edge" functions (associated with edges, one circulation DoF per edge)
// discretize vector fields whose curl must remain square-integrable and belong
// to "H(curl)"; facet "BF_Facet" functions (associated with faces, one flux DoF
// per face) discretize vector fields whose divergence must remain
// square-integrable and belong to "H(div)"; and volume "BF_Volume" functions
// (associated with elements, one DoF per element) discretize scalar densities
// in "L2". These spaces are linked by the gradient, curl and divergence
// operators: "grad" maps "H1" into "H(curl)", "curl" maps "H(curl)" into
// "H(div)", and "div" maps "H(div)" into "L2", reflecting the identities "curl
// grad = 0" and "div curl = 0" at the discrete level.  Each physical quantity
// in encountered in the tutorials has a natural home in this hierarchy: "v" and
// "T" live in "H1", "e" and "h" live in "H(curl)", "b" and "d" live in
// "H(div)", and charge or energy densities live in "L2".
//
// We model a rectangular waveguide, excited by a given electric field mode on
// an input port. Since the degrees of freedom correspond to the circulation of
// the electric field, imposing a Dirichlet boundary condition is not as trivial
// as with Lagrange finite elements (where the coefficients correspond to the
// value of the field at the nodes). With edge elements, the coefficients will
// change depending on the length and the orientation of the mesh edges: in this
// tutorial we will use a pre-resolution to compute the coefficients on the
// waveguide port.
//
// We assume that the waveguide is a good conductor, so we impose a "Perfect
// Electric Conductor" (PEC) boundary condition on the rest of the waveguide
// boundary, which amounts to imposing that the tangential component of the
// electric field, "e_t", is zero. Since the degrees of freedom are the
// circulation of the electric field, this amounts to setting the coefficients
// "ee_k" to 0 on the boundary edges of the waveguide.
//
// To truncate the computational domain, we cannot use the same infinite
// elements as in tutorial 3: the oscillatory nature of the wave solution
// prevents simply compressing an outside ring around the computational
// domain. Instead, we use a Silver-Muller Absorbing Boundary Condition (ABC),
// which expresses that the field behaves locally like an outgoing plane wave:
//
//   n x h = - sqrt(epsilon / mu) e_t,
//
// with "n" the exterior normal on the boundary, "h" the magnetic field,
// "epsilon" the dielectric permittivity, "mu" the magnetic permeability and
// "e_t = e - (e . n) n = (n x e) x n" the tangential electric field.
//
// We use the same phasor convention as in tutorial 4, i.e. "e(x, y, z, t) = Re(
// E(x, y, z) exp(i omega t) )", with the angular frequency "omega = 2 pi f".

Group {
  // Physical regions:
  Waveguide = Region[ 1000 ];
  Air = Region[ 1001 ];
  BndPEC = Region[ 2000 ]; // Boundary with PEC condition
  BndPort = Region[ 2001 ]; // Input port
  BndInf = Region[ 2002 ]; // Boundary with ABC

  // Abstract regions:
  //  - "Vol_Wav": overall domain
  //  - "Sur_Port_Wav": port with imposed electric field mode
  //  - "Sur_Inf_Wav": truncation surface with absorbing boundary condition
  Vol_Wav = Region[ {Waveguide, Air} ];
  Sur_Port_Wav = Region[ BndPort ];
  Sur_Inf_Wav = Region[ BndInf ];
}

Include "full_wave_common.pro";

Function {
  eps0 = 8.8541878176e-12;
  mu0 = 4 * Pi * 1e-7;

  DefineConstant[
    m = { 1, Min 1, Max 10, Step 1,
      Name "Parameters/Excitation mode (x)"},
    n = { 1, Min 1, Max 10, Step 1, Visible (dim == 3),
      Name "Parameters/Excitation mode (z)"},
    epsr = { 1, Min 1, Max 10, Step 0.1,
      Name "Parameters/Relative permittivity in waveguide"},
    f = { 8e9, Min 1e8, Max 1e10, Step 1e8,
      Name "Parameters/Frequency [Hz]"}
  ];

  epsilon[ Region[ {Air, BndInf} ] ] = eps0;
  epsilon[ Region[ {Waveguide, BndPort} ] ] = epsr * eps0;
  mu[] = mu0;
  nu[] = 1 / mu0;

  i[] = Complex[0, 1];
  omega = 2 * Pi * f;

  // The expressions that follow give the transverse electric-field profile of
  // the mode excited at the port, "ePort". Both in 2D and in 3D we excite a TM
  // mode (transverse magnetic with respect to the propagation direction y,
  // i.e. "H_y = 0"). In a TM mode, all six field components can be recovered
  // from the single longitudinal scalar "E_y" by algebra and one transverse
  // derivative (see the explicit formulas below).
  //
  // In our FEM setup we do not impose "E_y" at the port: only the tangential
  // components -- "E_x" in 2D, and "(E_x, E_z)" in 3D -- are prescribed there,
  // and "E_y" emerges in the interior as part of the full vector solve. The "0"
  // written for the y-component of "ePort" below is normal to the port and is
  // silently dropped by the tangential L2 projection on "Sur_Port_Wav".
  //
  // The generating scalar "E_y" satisfies a Helmholtz equation with Dirichlet
  // BC at the PEC walls ("E_y" being tangential to them):
  //   - In 3D,
  //
  //         E_y = E_0 sin(kx (x - Wx / 2)) sin(kz (z - Wz / 2)),
  //
  //     giving the TM_mn mode of the rectangular waveguide. Both indices must
  //     satisfy "m >= 1" and "n >= 1", otherwise "E_y" vanishes identically and
  //     the mode degenerates (the family TM_m0 or TM_0n does not exist).
  //   - In 2D the geometry is invariant in z (no PEC walls at z = const), so
  //     only "kx = m Pi / Wx" is quantized:
  //
  //        E_y = E_0 sin(kx (x - Wx / 2)),
  //
  //     giving the TM_m mode of the parallel-plate waveguide (m >= 1; the m = 0
  //     case would be the TEM mode).
  //
  // The transverse components follow from
  //
  //   E_x = -gamma / kc^2 * \partial_x E_y
  //
  // and (in 3D)
  //
  //   E_z = -gamma / kc^2 * \partial_z E_y,
  //
  // with "gamma" the longitudinal propagation constant and "kc^2 = kx^2 + kz^2"
  // in 3D, reducing to "kc^2 = kx^2" in 2D. From the dispersion relation
  //
  //   k^2 = kx^2 + kz^2 + (-i gamma)^2
  //
  // with "k = omega sqrt(mu eps)", it follows that
  //
  //   gamma = sqrt(k^2 - (kx^2 + kz^2)) * i
  //
  // i.e. "gamma" is purely imaginary above cutoff (propagating mode) and purely
  // real below cutoff (evanescent mode).

  k[] = omega * Sqrt[epsilon[] * mu[]];
  kx = m * Pi / Wx; // transverse wavevector (x)
  kz = n * Pi / Wz; // transverse wavevector (z)
  If(dim == 2)
    gamma[] = Sqrt[k[]^2 - kx^2] * i[];
    ePort[] = Vector[-gamma[] / kx * Cos[kx * (X[] - Wx / 2)], 0, 0];
  Else
    kc = Sqrt[kx^2 + kz^2];
    gamma[] = Sqrt[k[]^2 - kc^2] * i[];
    ePort[] = Vector[
      -gamma[] * kx / kc^2 * Cos[kx * (X[] - Wx / 2)] * Sin[kz * (Z[] - Wz / 2)],
      0.,
      -gamma[] * kz / kc^2 * Sin[kx * (X[] - Wx / 2)] * Cos[kz * (Z[] - Wz / 2)]
    ];
  EndIf

  c = 1 / Sqrt[epsr * eps0 * mu0];
  DefineConstant[
    lambda = {c / f * 100, ReadOnly 1, Highlight "LightGrey",
      Name "Parameters/}Wavelength [cm]"},
    f_cutoff = { (dim == 2) ? (0.5 * c * m / Wx) :
      0.5 * c * Sqrt[(m * m) / (Wx * Wx) + (n * n) / (Wz * Wz)],
      ReadOnly 1, Highlight "LightGrey",
      Name "Parameters/}Cutoff frequency [Hz]"}
  ];
}

Jacobian {
  // Instead of explicitly defining a Jacobian method for volumes and surfaces,
  // a single Jacobian method can be defined piecewise:
  { Name Jac;
    Case {
      { Region Sur_Port_Wav; Jacobian Sur; }
      { Region Sur_Inf_Wav; Jacobian Sur; }
      { Region Vol_Wav; Jacobian Vol; }
    }
  }
}

Integration {
  { Name Int;
    Case {
      { Type Gauss;
        Case {
          { GeoElement Line; NumberOfPoints 4; }
          { GeoElement Triangle; NumberOfPoints 4; }
          { GeoElement Tetrahedron; NumberOfPoints 4; }
        }
      }
    }
  }
}

Constraint {
  { Name e_Wav;
    Case {
      // On the PEC boundary, the circulation of the electric field is imposed
      // to 0 on the edges:
      { Region BndPEC; Type Assign; Value 0.; }
      // On the input port, the circulation of the electric field is imposed
      // through the resolution of a linear system, performing the L2 projection
      // of a known mode:
      { Region BndPort; Type AssignFromResolution; NameOfResolution Proj_e; }
    }
  }
}

Group{
  Dom_Hcurl_e_Wav =  Region[ {Vol_Wav, Sur_Port_Wav, Sur_Inf_Wav} ];
}

FunctionSpace {
  // The electric "e" field belongs to H(curl), meaning that both "e" and "curl
  // e" are square integrable. Basis functions "BF_Edge" are associated with the
  // edges of the mesh, and ensure the tangential continuity of the electric
  // field across mesh elements:
  { Name Hcurl_e_Wav; Type Form1;
    BasisFunction {
      { Name se; NameOfCoef ee; Function BF_Edge;
        Support Dom_Hcurl_e_Wav; Entity EdgesOf[All]; }
    }
    Constraint {
      { NameOfCoef ee; EntityType EdgesOf; NameOfConstraint e_Wav; }
    }
  }
}

Formulation {
  // To derive the finite element formulation in terms of the electric field e
  // in the frequency domain, we start from Faraday's law:
  //
  //   curl e = - i omega b,
  //
  // with "b = mu h" the magnetic flux density. Taking its curl and introducing
  // the reluctivity "nu = 1 / mu", we get
  //
  //   curl (nu curl e) = -i omega curl h.
  //
  // Assuming no conduction currents (we treat the waveguide walls as perfect
  // conductors), Maxwell-Ampere's equation reduces to
  //
  //   curl h = i omega d,
  //
  // with "d = epsilon e" the displacement field. Combining the last two
  // equations leads to the strong form
  //
  //   curl (nu curl e) - omega^2 epsilon e = 0.
  //
  // The weak formulation is obtained by multiplying by test functions "e'" in
  // H(curl) and integrating over the domain "Vol_Wav". This leads to finding
  // "e" such that
  //
  //   (curl (nu curl e), e')_Vol_Wav - (omega^2 epsilon e, e')_Vol_Wav = 0.
  //
  // holds for all test functions "e'". After integration by parts, the weak form
  // becomes: find e such that, for all "e'":
  //
  //   (nu curl e, curl e')_Vol_Wav + (n x nu curl e, e')_Bnd_Vol_Wav
  //      - (omega^2 epsilon e, e')_Vol_Wav = 0.
  //
  // The boundary term vanishes where Dirichlet boundary conditions are applied
  // (PEC and input port). On the outside boundary, the Silver-Muller ABC leads
  // to a Robin boundary condition:
  //
  //   (n x nu curl e, e')_Sur_Inf
  //     = (n x nu (- i omega b), e')_Sur_Inf
  //     = - (i omega n x h, e')_Sur_Inf
  //     = (i omega sqrt(epsilon / mu) e_t, e')_Sur_Inf
  //
  { Name FullWave_e; Type FemEquation;
    Quantity {
      { Name e; Type Local; NameOfSpace Hcurl_e_Wav; }
    }
    Equation {
      Integral { [ nu[] * Dof{d e} , {d e} ];
        In Vol_Wav; Integration Int; Jacobian Jac; }

      Integral { [ -omega^2 * epsilon[] * Dof{e} , {e} ];
        In Vol_Wav; Integration Int; Jacobian Jac; }

      Integral { [ i[] * omega * Sqrt[epsilon[] / mu[]] * Dof{e} , {e} ];
        In Sur_Inf_Wav; Integration Int; Jacobian Jac; }
    }
  }

  // In order to compute the coefficients of the imposed electric field mode on
  // the input port, we define an auxiliary weak formulation to perform the
  // L2-projection of the "ePort[]" function onto the finite element basis: find
  // "e" such that
  //
  //   (e, e')_Sur_Port_Wav - (ePort, e')_Sur_Port_Wav = 0
  //
  // holds for all test functions "e'".
  //
  // A direct Dirichlet assignment (as we did for the PEC boundary above, where
  // the imposed circulation was simply zero) is impractical here: each
  // coefficient "ee_k" is an integral of "ePort" along an individual mesh edge,
  // which would have to be computed by hand for every edge on the port --
  // depending on its length, orientation, and the local value of "ePort".  The
  // L2 projection does this automatically and consistently: it solves a small
  // auxiliary system whose solution is, by construction, the best
  // tangentially-continuous edge-element approximation of "ePort" on the port.
  { Name Projection_e; Type FemEquation;
    Quantity {
      { Name e; Type Local; NameOfSpace Hcurl_e_Wav; }
    }
    Equation {
      Integral { [ Dof{e} , {e} ];
        In Sur_Port_Wav; Integration Int; Jacobian Jac; }
      Integral { [ -ePort[] , {e} ];
        In Sur_Port_Wav; Integration Int; Jacobian Jac; }
    }
  }
}

Resolution {
  { Name Wav;
    System {
      { Name Sys_Wav; NameOfFormulation FullWave_e; Type Complex; }
    }
    Operation {
      Generate[Sys_Wav]; Solve[Sys_Wav]; SaveSolution[Sys_Wav];
    }
  }

  // One constraint associated to the function space used in the "FullWave_e"
  // formulation is of type "AssignFromResolution Proj_e". GetDP will
  // automatically perform the "Proj_e" resolution during the pre-processing
  // stage of the main "Wav" resolution:
  { Name Proj_e;
    System {
      { Name Sys_Proj; NameOfFormulation Projection_e; Type Complex;
        DestinationSystem Sys_Wav; }
    }
    Operation {
      // Note the "TransferSolution[]" operation, which will copy the ee_k
      // coefficients computed on "Sur_Port_Wav" into the "DestinationSystem"
      // "Sys_Wav", defined in the "Wav" resolution:
      Generate[Sys_Proj]; Solve[Sys_Proj]; TransferSolution[Sys_Proj];
    }
  }
}

PostProcessing {
  { Name Wav; NameOfFormulation FullWave_e;
    Quantity {
      { Name e;
        Value {
          Term{ [ {e} ]; In Vol_Wav; Jacobian Jac; }
        }
      }
      { Name eNorm;
        Value {
          Term{ [ Norm[{e}] ]; In Vol_Wav; Jacobian Jac; }
        }
      }
      { Name h;
        Value{
          Term{ [ i[] * nu[] / omega * {d e} ]; In Vol_Wav; Jacobian Jac; }
        }
      }
      // The (complex) Poynting vector "s = e x h*" gives the time-averaged
      // power flux density (with factor 1/2 for peak-value phasors). The symbol
      // "/\" is GetDP's cross-product operator; "Conj[...]" takes the complex
      // conjugate:
      { Name s;
        Value{
          Term{ [ 0.5 * {e} /\ Conj[i[] * nu[] / omega * {d e}] ]; In Vol_Wav;
            Jacobian Jac; }
        }
      }
    }
  }
}

PostOperation {
  { Name Map; NameOfPostProcessing Wav;
    Operation {
      Print [ e, OnElementsOf Vol_Wav, File "e.pos"];
      Print [ eNorm, OnElementsOf Vol_Wav, File "eNorm.pos"];
      Print [ h, OnElementsOf Vol_Wav, File "h.pos"];
      Print [ s, OnElementsOf Vol_Wav, File "s.pos"];
    }
  }
}

// Gmsh script describing the geometry of the microstrip line and the associated
// meshing constraints.

// Dimensions (there are no units in Gmsh -- the following values should be
// interpreted as dimensions in meters due to the definition of the physical
// parameters in MKS units in "microstrip.pro"):
h = 1.e-3; // thickness of dielectric substrate
w = 4.72e-3; // width of microstrip line
t = 0.035e-3; // thickness of microstrip line
xBox = w/2. * 6.; // width of air box
yBox = h * 12.; // height of air box

// Global mesh size factor (that can be modified interactively in the Gmsh
// graphical interface):
s = DefineNumber[1., Name "Parameters/Global mesh size factor"];

// Target mesh sizes:
p0 = h / 10 * s;
pLine0 = w / 20 * s;
pLine1 = w / 100 * s;
pxBox = xBox / 10 * s;
pyBox = yBox / 8 * s;

// We create this first geometry in a bottom-up manner with the built-in Gmsh
// CAD kernel, successively defining model points, lines, loops and surfaces:
Point(1) = {0, 0, 0, p0};
Point(2) = {xBox, 0, 0, pxBox};
Point(3) = {xBox, h, 0, pxBox};
Point(4) = {0, h, 0, pLine0};
Point(5) = {w/2, h, 0, pLine1};
Point(6) = {0, h+t, 0, pLine0};
Point(7) = {w/2, h+t, 0, pLine1};
Point(8) = {0, yBox, 0, pyBox};
Point(9) = {xBox, yBox, 0, pyBox};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 9};
Line(4) = {9, 8};
Line(5) = {8, 6};
Line(7) = {4, 1};
Line(8) = {5, 3};
Line(9) = {4, 5};
Line(10) = {6, 7};
Line(11) = {5, 7};
Curve Loop(12) = {1, 2, -8, -9, 7};
Plane Surface(13) = {12};
Curve Loop(14) = {10, -11, 8, 3, 4, 5};
Plane Surface(15) = {14};

// The last step is to define physical groups, assigning unique tags (and names)
// to groups of model entities. Physical groups serve two purposes: when they
// are defined they tell Gmsh which elements to save in the mesh file (by
// default only elements belonging to at least one physical group are saved),
// and they provide the tags that GetDP uses in "Region[]" to identify parts of
// the domain.
Physical Surface("Air", 1) = {15};
Physical Surface("Dielectric", 2) = {13};
Physical Curve("Ground", 10) = {1};
Physical Curve("Microstrip boundary", 11) = {9, 10, 11};
Physical Curve("Inf", 12) = {2, 3, 4};

// In a two-dimensional electrostatic model (see tutorial 1), the finite element
// expansion of the electric scalar potential "v(x, y)" associates basis
// functions with individual nodes in the mesh:
//
//   v(x, y) = Sum_{k in N} vn_k sn_k(x, y)
//
// with "N" the set of all the nodes in the mesh and "sn_k(x, y)" the Lagrange
// basis function associated with node "k". Now consider the situation where
// "v(x, y)" is constant over a region of the problem, e.g. "v(x, y) =
// vf_electrode" on an electrode. By factorizing the identical nodal values
// "vf_electrode", a global (non-local) basis function "sf_electrode(x, y)" is
// obtained as factor:
//
//   v(x, y) = Sum_{k in N} vn_k sn_k(x, y)
//     = Sum_{k in N \ E} vn_k sn_k(x, y) + Sum_{k in E} vn_k sn_k(x, y)
//     = Sum_{k in N \ E} vn_k sn_k(x, y) + vf_electrode Sum_{k in E} sn_k(x, y)
//     = Sum_{k in N \ E} vn_k sn_k(x, y) + vf_electrode sf_electrode(x, y)
//
// with "E" the set of nodes on the electrode. The global basis function
//
//   sf_electrode(x, y) = Sum_{k in E} sn_k(x, y)
//
// is the sum of the Lagrange (nodal) basis functions associated with the nodes
// of the electrode. This function
//  - is a continuous fonction, scalar in this case;
//  - is equal to 1 at the nodes of the electrode, and to 0 at all other nodes;
//  - decreases from 1 to 0 over the one-element-thick layer of elements sharing
//    at least one node with the electrode region.
//
// One such global basis function can be associated with each electrode in the
// system, so that the finite element expansion of the electric scalar potential
// reads:
//
//   v(x, y) = Sum_k vn_k sn_k(x, y)
//     + Sum_electrode vf_electrode sf_electrode(x, y)
//
// where "Sum_electrode" runs over all the electrodes in the model (one term per
// electrode, so two in the present example: the boundary of the microstrip line
// and the ground), and "Sum_k" runs over all nodes except those of the
// electrode regions.
//
// In this tutorial we modify the electrostatic problem from tutorial 1 to show
// how GetDP takes advantage of global basis functions
//  - to efficiently compute the electrode charges "Q_electrode", which are the
//    energy duals of the global "vf_electrode" potentials;
//  - to deal with floating potentials, which are the computed electrode
//    potentials when the electrode charge is imposed;
//  - to provide output global quantities (charges, voltages, capacitances) that
//    can be used in an external circuit.

Group {
  // Physical regions:
  Air = Region[ 1 ];
  Dielectric = Region[ 2 ];
  Ground = Region[ 10 ];
  Microstrip = Region[ 11 ];
  Inf = Region[ 12 ];

  // Abstract regions:
  //  - "Vol_Ele": overall domain
  //  - "Sur_Neu_Ele": surface with non homogeneous Neumann boundary conditions
  //  - "Sur_Electrodes_Ele": electrode regions
  Vol_Ele = Region[ {Air, Dielectric} ];
  Sur_Neu_Ele = Region[ {} ];
  Sur_Electrodes_Ele = Region [ {Ground, Microstrip} ];
}

Function {
  DefineConstant[
    TypeBC = {0, Choices{0="Fixed voltage", 1="Fixed charge"},
      Name "Parameters/0Type"}
    ValueBC = {1e-3,
      Name "Parameters/1Value"}
    epsr = {9.8,
      Name "Parameters/2Relative permittivity"}
    DisplayGlobalBasisFunctions = {0, Choices {0, 1},
      Name "Parameters/3Display global basis functions"}
    OverwriteOutput = {1, Choices {0, 1},
      Name "Parameters/4Overwrite output.txt file"}
  ];

  eps0 = 8.854187818e-12;
  epsilon[ Air ] = eps0;
  epsilon[ Dielectric ] = epsr * eps0;
}

Constraint {
  // The Dirichlet boundary condition on the local electric potential is only
  // used for the homogeneous condition on the top and right boundaries. The
  // boundary of the microstrip line and the ground are now treated as
  // electrodes, on which either the potential or the charge is imposed through
  // the "Voltage_Ele" or "Charge_Ele" constraints below.
  { Name v_Ele;
    Case {
      { Region Inf; Value 0; }
    }
  }

  { Name Voltage_Ele;
    Case {
      // If the "Fixed voltage" option is enabled (i.e. when "TypeBC == 0"),
      // impose the global potential on the boundary of the microstrip line (and
      // thus the voltage, i.e. the potential difference between the microstrip
      // and the ground) to "ValueBC":
      If(TypeBC == 0)
        { Region Microstrip; Value ValueBC; }
      EndIf
      // Impose the global (i.e. region-wise) potential to 0 on the ground
      // electrode:
      { Region Ground; Value 0; }
    }
  }
  { Name Charge_Ele;
    Case {
      // If the "Fixed charge" option is enabled (i.e. when "TypeBC == 1"),
      // impose the global charge on the boudary of the microstrip line:
      If(TypeBC == 1)
        { Region Microstrip; Value ValueBC; }
      EndIf
      // Since "Ground" is always constrained to zero potential by "Voltage_Ele"
      // above, its charge cannot be prescribed; it will be computed as a global
      // output, cf. "Q" in the PostProcessing below.
    }
  }
}

Group{
  Dom_H1_v_Ele =  Region[ {Vol_Ele, Sur_Neu_Ele, Sur_Electrodes_Ele} ];
}

FunctionSpace {
  // The natural treatment of global quantities in GetDP stems from the fact
  // that nearly all the work is done at the level of the FunctionSpace
  // definition. As seen above, the finite element expansion of the potential
  // "v" is
  //
  //   v(x, y) = Sum_k vn_k sn_k(x, y)
  //     + Sum_electrode vf_electrode sf_electrode(x, y)
  //
  // with "Sum_k" running over all nodes except those of the electrode
  // regions. The global basis function "sf_electrode(x, y)" is the sum of the
  // nodal basis functions associated with the nodes of the electrode. This is
  // precisely what one finds in the FunctionSpace definition below, where:
  //  - "BF_Node" is the standard Lagrange basis function, that we associate
  //    with all the nodes of the domain except those on the electrode regions
  //    ("Entity NodesOf[ All, Not Sur_Electrodes_Ele ]"). One degree of freedom
  //    is created per node of the region.
  //  - "BF_GroupOfNodes" is the sum of Lagrange basis functions associated with
  //    the group of nodes of the electrodes ("Entity GroupsOfNodesOf[
  //    Sur_Electrodes_Ele ]"). A single degree of freedom is shared by all the
  //    nodes of each listed region -- whose value is the "vf_electrode"
  //    coefficient attached to the global basis function "sf_electrode".
  //
  // By activating the "Display global basis functions" option, you can
  // visualize the two "sf_electrode" basis functions in the model: one
  // associated with the boundary of the microstrip line, the other with the
  // ground.
  //
  // The global quantities are attributed an explicit name in the
  // "GlobalQuantity" section; these names are used in the corresponding
  // "GlobalTerm" in the Formulation. Such global terms are the equivalent of a
  // "Integral" term, but where no integration needs to be performed. The
  // "AssociatedWith" statement refers to the fact that the global potential of
  // an electrode is the (electrostatic) energy dual of the electric charge
  // carried by that electrode. Indeed, let us consider the electrostatic weak
  // formulation derived in tutorial 1: find "v" in H1_v_Ele such that
  //
  //   (epsilon grad v, grad v')_Vol_Ele
  //     - (n . (epsilon grad v), v')_Bnd_Vol_Ele = 0
  //
  // holds for all test functions "v'". When the test-function "v'" is
  // "sf_electrode", the boundary term reduces to
  //
  //   - (n . (epsilon grad v), sf_electrode)_Sur_Electrodes_Ele.
  //
  // Since "sf_electrode == 1" on the electrode, and since "e = -grad v" and "d
  // = epsilon e", the boundary term is equal to
  //
  //   - \int_Sur_Electrodes_Ele n . (epsilon grad v)
  //     = \int_Sur_Electrodes_Ele n . (epsilon e)
  //     = \int_Sur_Electrodes_Ele n . d
  //     = - \int_Electrodes (div d)
  //     = - \int_Electrodes rho
  //     = - Q_electrode
  //
  // i.e. the charge carried by the electrode.
  //
  // Constraints can be set on either component of the FunctionSpace. Besides
  // the usual Dirichlet boundary condition on the local field, one may fix
  // either the "Voltage" or the "Charge" of each indidual electrode (never
  // both, of course). When the "Charge" is fixed, the computed "Voltage" for
  // that electrode is the so-called floating potential.
  { Name H1_v_Ele; Type Form0;
    BasisFunction {
      { Name sn; NameOfCoef vn; Function BF_Node;
        Support Dom_H1_v_Ele; Entity NodesOf[ All, Not Sur_Electrodes_Ele ]; }
      { Name sf; NameOfCoef vf; Function BF_GroupOfNodes;
        Support Dom_H1_v_Ele; Entity GroupsOfNodesOf[ Sur_Electrodes_Ele ]; }
    }
    GlobalQuantity {
      { Name Voltage; Type AliasOf; NameOfCoef vf; }
      { Name Charge; Type AssociatedWith; NameOfCoef vf; }
    }
    Constraint {
      { NameOfCoef vn; EntityType NodesOf;
        NameOfConstraint v_Ele; }
      { NameOfCoef Voltage; EntityType GroupsOfNodesOf;
        NameOfConstraint Voltage_Ele; }
      { NameOfCoef Charge; EntityType GroupsOfNodesOf;
	NameOfConstraint Charge_Ele; }
    }
    // Subspace definition only needed to display the global basis functions in
    // post-processing:
    SubSpace {
      { Name vf; NameOfBasisFunction sf; }
    }
  }
}

Jacobian {
  { Name Vol;
    Case {
      { Region All; Jacobian Vol; }
    }
  }
}

Integration {
  { Name Int;
    Case {
      { Type Gauss;
        Case {
          { GeoElement Triangle; NumberOfPoints 1; }
        }
      }
    }
  }
}

Formulation {
  // The formulation contains two changes compared to the formulation from the
  // first tutorial: the global quantities are declared as "Global" in the
  // "Quantity" section, and a "GlobalTerm" is added that triggers the assembly
  // of the additional equation per electrode (the "pre-integrated" boundary
  // term) in the system to compute the charge "Q_electrode". "Q_electrode"
  // satisfies (just consider the equation corresponding to the test function
  // "sf_electrode"):
  //
  //   Q_electrode = (epsilon grad v, grad sf_electrode)_Vol_Ele
  //
  // Note that, although "Vol_Ele" appears as the integration domain, "grad
  // sf_electrode" is nonzero only in the one-element-thick layer adjacent to
  // the electrode.
  { Name Electrostatics_v; Type FemEquation;
    Quantity {
      { Name v; Type Local; NameOfSpace H1_v_Ele; }
      { Name V; Type Global; NameOfSpace H1_v_Ele [Voltage]; }
      { Name Q; Type Global; NameOfSpace H1_v_Ele [Charge]; }
      // The following line is only needed to display the global basis functions
      // in post-processing:
      { Name vf; Type Local; NameOfSpace H1_v_Ele [vf]; }
    }
    Equation {
      Integral { [ epsilon[] * Dof{d v} , {d v} ];
        In Vol_Ele; Jacobian Vol; Integration Int; }
      // As with "Integral" terms, the second argument of "GlobalTerm [. , .]"
      // is the test function (the "AliasOf" quantity, here "{V}"), and the
      // first argument is the term multiplying the basis function (involving
      // the "AssociatedWith" Dof, here "{Q}"). The "In" keyword specifies the
      // region on which the GlobalQuantity is defined. No integration is
      // performed (hence "Global", not "Integral"):
      GlobalTerm { [ -Dof{Q} , {V} ]; In Sur_Electrodes_Ele; }
    }
  }
}

Resolution {
  { Name Ele;
    System {
      { Name Sys_Ele; NameOfFormulation Electrostatics_v; }
    }
    Operation {
      Generate[Sys_Ele]; Solve[Sys_Ele]; SaveSolution[Sys_Ele];
    }
  }
}

PostProcessing {
 { Name Ele; NameOfFormulation Electrostatics_v;
    Quantity {
      // Local post-processing quantities are defined as usual:
      { Name v; Value {
          Term { [ {v} ]; In Vol_Ele; Jacobian Vol; }
        }
      }
      { Name e; Value {
          Term { [ -{d v} ]; In Vol_Ele; Jacobian Vol; }
        }
      }
      { Name d; Value {
          Term { [ -epsilon[] * {d v} ]; In Vol_Ele; Jacobian Vol; }
        }
      }
      // Global quantities defined in a similar way:
      { Name Q; Value {
          Term { [ {Q} ]; In Sur_Electrodes_Ele; }
        }
      }
      { Name V; Value {
          Term { [ {V} ]; In Sur_Electrodes_Ele; }
        }
      }
      // Since we have access to the voltage and the charge we can directly
      // compute the capacitance:
      { Name C; Value {
          Term { [ {Q} / {V} ]; In Sur_Electrodes_Ele; }
        }
      }
      // We could also compute the capacitance by way of the energy ("C = 2 *
      // energy / V^2"), obtained by integrating ("epsilon / 2 * |grad v|^2")
      // over the domain, which can be done by defining an "Integral"
      // post-processing quantity:
      { Name energy; Value {
          Integral { [ epsilon[] / 2. * SquNorm[{d v}] ]; In Vol_Ele;
            Jacobian Vol; Integration Int; }
	}
      }
      // Finally we can use the "BF[]" operator to retrieve the global basis
      // functions on the two electrodes, in order to visualize them:
      { Name sf_Ground; Value {
          Term { [ BF{vf} ]; In Vol_Ele; SubRegion Ground; Jacobian Vol; }
        }
      }
      { Name sf_Microstrip; Value {
          Term { [ BF{vf} ]; In Vol_Ele; SubRegion Microstrip; Jacobian Vol; }
        }
      }
    }
  }
}

PostOperation {
  { Name Map; NameOfPostProcessing Ele;
    Operation {
      Print[ v, OnElementsOf Vol_Ele, File "v.pos" ];

      If(DisplayGlobalBasisFunctions)
        // If "DisplayGlobalBasisFunctions" is set, save the global basis
        // function associated with the ground in a Gmsh post-processing file:
        Print[ sf_Ground, OnElementsOf Vol_Ele, File "sf_Ground.pos" ];
        // To better visualize the global function, print some Gmsh options in
        // an additional option file, that will also be automatically merged by
        // Gmsh and that will modify the view's display options:
        Echo[ StrCat[
            "l = PostProcessing.NbViews - 1;", // get last view index
            "View[l].IntervalsType = 1;", // display iso-curves
            "View[l].NbIso = 40;", // display 40 iso-curves
            "View[l].ShowElement = 1;" // display the element outlines
          ], File "sf.opt"];

        // Do the same with the global function associated with the microstrip:
        Print[ sf_Microstrip, OnElementsOf Vol_Ele,
          File "sf_Microstrip.pos" ];
        // Instead of creating the same option file, just append (using ">") an
        // empty string to the existing file, which will also trigger Gmsh to
        // merge it:
        Echo[ "", File > "sf.opt"];
      EndIf

      If(OverwriteOutput)
        // Remove the "output.txt" file if the "OverwriteOutput" option is set:
        DeleteFile[ "output.txt" ];
      EndIf

      // Output global quantities in the "output.txt" file, and also send them
      // to the graphical user interface using the "SendToServer" option:
      Echo[ "Microstrip charge [C]:", Format Table, File > "output.txt"];
      Print[ Q, OnRegion Microstrip, File > "output.txt", Color "AliceBlue",
        Format Table, SendToServer "}Output/Microstrip/Charge [C]" ];

      Echo[ "Microstrip potential [V]:", Format Table, File > "output.txt"];
      Print[ V, OnRegion Microstrip, File > "output.txt", Color "AliceBlue",
        Format Table, SendToServer "}Output/Microstrip/Potential [V]" ];

      Echo[ "Ground charge [C]:", Format Table, File > "output.txt"];
      Print[ Q, OnRegion Ground, File > "output.txt", Color "AliceBlue",
        Format Table, SendToServer "}Output/Ground/Charge [C]" ];

      Echo[ "Microstrip capacitance [F]:", Format Table, File > "output.txt"];
      Print[ C, OnRegion Microstrip, File > "output.txt", Color "AliceBlue",
        Format Table, SendToServer "}Output/Global/Capacitance [F]" ];

      Echo[ "Electrostatic energy [J]:", Format Table, File > "output.txt"];
      Print[ energy[Vol_Ele], OnGlobal, File > "output.txt", Color "AliceBlue",
        Format Table, SendToServer "}Output/Global/Energy [J]" ];
    }
  }
}

mm = 1e-3;

DefineConstant[
  d = {10, Min 1, Max 100, Step 1,
    Name "Parameters/Conductor spacing [mm]"}
  e = {4, Min 1, Max 10, Step 1,
    Name "Parameters/Conductor thickness [mm]"}
  h = {40, Min 5, Max 100, Step 1,
    Name "Parameters/Conductor height [mm]"}
  g = {15, Min 1, Max 10, Step 1,
    Name "Parameters/Dielectric thickness [mm]"}
  t = {1.5, Min 1, Max 10, Step 1,
    Name "Parameters/Shield thickness [mm]"}
];

// Go back to meters:
d = d * mm;
e = e * mm;
h = h * mm;
g = g * mm;
t = t * mm;

// Internal and external radii of "infinite" region
rInt = 4 * (d + e + g + t);
rExt = 5.5 * (d + e + g + t);

SetFactory("OpenCASCADE");

// Include the file with shared parameters, so that the same values are used
// here and in "busbar.pro":
Include "busbar_common.pro";

// Three rectangular conductors carrying balanced three-phase currents. The
// conductors are placed side by side, separated by a distance d, each with
// width e and height h:
For i In {1:3}
  Rectangle(i) = {-d - 1.5 * e + (i - 1) * (e + d), -h / 2, 0, e, h};
EndFor

// Insulation surrounding the three conductors, with a gap g on each side:
xx = -d - 1.5 * e - g;
yy = -h / 2 - g;
ww = 2 * d + 3 * e + 2 * g;
Rectangle(4) = {xx, yy, 0, ww, h + 2 * g};

// Metallic shield (thickness t) around the insulation:
Rectangle(5) = {xx - t, yy - t, 0, ww + 2 * t, h + 2 * g + 2 * t};

// Air region (disk) and infinite ring for the VolSphShell transformation (see
// tutorial 3); rInt and rExt are defined in "busbar_common.pro":
Disk(6) = {0, 0, 0, rInt};
Disk(7) = {0, 0, 0, rExt};

// Make all overlapping surfaces conformal (see tutorial 3):
Coherence;

// Mesh size constraints:
MeshSize { PointsOf{ Surface{6, 7}; } } = h / 5;
MeshSize { PointsOf{ Surface{4, 5}; } } = e / 3;
MeshSize { PointsOf{ Surface{1:3}; } } = e / 5;

// Physical groups:
Physical Surface("Conductor 1", 100) = {1};
Physical Surface("Conductor 2", 101) = {2};
Physical Surface("Conductor 3", 102) = {3};
Physical Surface("Insulator", 200) = {4};
Physical Surface("Shield", 201) = {5};
Physical Surface("Air", 202) = {6};
Physical Surface("Infinite shell", 203) = {7};

// Outermost boundary for the homogeneous Dirichlet boundary condition on the
// magnetic vector potential:
Physical Curve("Infinity", 300) = CombinedBoundary{ Surface{:}; };

// Boundary of all solid regions (conductors, insulator and shield) for the
// convection (Robin) boundary condition in the thermal problem:
Physical Curve("Convection", 301) = CombinedBoundary{ Surface{1:5}; };

// This tutorial models the coupled magneto-thermal behavior of a three-phase
// busbar system, consisting of three conductors carrying balanced AC currents
// inside an insulating enclosure surrounded by a grounded metallic shield.
//
// The magnetic part solves a magneto-quasistatic problem (as in tutorial 4) in
// the frequency domain, using an "a-v" formulation that combines the magnetic
// vector potential "a" with an electric scalar potential "v". The electric
// scalar potential is introduced to handle conductors with imposed currents or
// voltages -- this extends the approach from tutorial 4 (which only had an
// imposed source current density) and from tutorial 6 (which introduced global
// quantities in an electrostatic context).
//
// More precisely, since "div b = 0", we start by introducing the vector
// potential "a" such that "b = curl a". Faraday's equation leads to "curl e =
// -\partial_t b = -\partial_t curl a", i.e. "curl (e + \partial_t a) = 0". We
// can thus derive "e + \partial_t a" from a scalar electric potential "v" such
// that "e + \partial_t a = -grad v".
//
// Whereas in tutorial 4 we fixed the gradient to zero (leading to the modified
// vector potential "a-formulation"), here we keep it in the conductors so that,
// in each conducting region, the electric field is expressed as
//
//   e = -\partial_t a - grad v,
//
// where "grad v" is assumed to be a constant (region-wise) vector along the
// z-axis in this 2D model. The current density in each conductor is then
//
//   j = sigma e = sigma (-\partial_t a - grad v).
//
// Integrating "j" over the cross-section of a conductor yields the total
// current, which can be imposed as a constraint. Similarly, the voltage drop
// per unit length along "z" can be imposed or computed. This is the "a-v"
// formulation, where "v" is discretized with a single global degree of freedom
// per conducting region.
//
// The thermal part solves the steady-state heat equation (as in tutorial 2),
// with time-averaged Joule losses from the magnetic solution as a volumetric
// heat source. Convection boundary conditions evacuate the heat on the outer
// surfaces of the solid regions.
//
// The coupling arises because the electrical conductivity "sigma" depends on
// temperature, and the Joule losses depend on the current distribution (which
// depends on "sigma"). A fixed-point (Picard) iteration (see tutorial 3)
// alternates between the magnetic and thermal solves until convergence.

Group {
  // Physical regions:
  Cond_1 = Region[ 100 ];
  Cond_2 = Region[ 101 ];
  Cond_3 = Region[ 102 ];
  Conds = Region[ {Cond_1, Cond_2, Cond_3} ];
  Insulator = Region[ 200 ];
  Shield = Region[ 201 ];
  Air = Region[ 202 ];
  AirInf = Region[ 203 ];
  Surface_Inf = Region[ 300 ];
  Surface_Conv = Region[ 301 ];

  // Abstract magnetic regions:
  //  - "Vol_Mag": full domain for the magnetic problem
  //  - "Vol_C_Mag": conducting regions where eddy currents and/or imposed
  //    currents flow (both the three-phase conductors and the grounded shield)
  //  - "Vol_Inf_Mag": region with infinite shell geometric transformation (see
  //    tutorial 3)
  //  - "Sur_Neu_Mag": surface with non-homogeneous Neumann conditions (empty
  //    here)
  Vol_Mag = Region[ {Conds, Insulator, Shield, Air, AirInf} ];
  Vol_C_Mag = Region[ {Conds, Shield} ];
  Vol_Inf_Mag = Region[ AirInf ];
  Sur_Neu_Mag = Region[ {} ];

  // Abstract thermal regions:
  //  - "Vol_The": domain for the thermal problem (solid regions only)
  //  - "Sur_Neu_The": surface with non-homogeneous Neumann conditions (empty)
  //  - "Sur_Rob_The": surface with Robin (convection) boundary condition
  Vol_The = Region[ {Conds, Insulator, Shield} ];
  Sur_Neu_The = Region[ {} ];
  Sur_Rob_The = Region[ {Surface_Conv} ];
}

Function {
  DefineConstant[
    Current = {1400, Min 1, Max 10000, Step 1,
      Name "Parameters/0Current [A]"}
    f = {50, Min 1, Max 1000, Step 1,
      Name "Parameters/1Frequency"}
    NonLinear = {1, Choices{0, 1},
      Name "Parameters/3Temperature-dependent conductivity" }
  ];

  // Magnetic parameters:
  mu0 = 4 * Pi * 1e-7;
  nu[] = 1. / mu0;

  If(NonLinear)
    // Temperature-dependent electrical conductivity of the copper conductors:
    // the argument "$1" will receive the local temperature value at runtime,
    // which is passed through "sigma[{T}]" in the Formulation (see below). This
    // is the mechanism through which the thermal solution influences the
    // magnetic problem:
    sigma[ Conds ] = 5.8e7 / (1 + 0.00393 * ($1 - 20));
  Else
    // Temperature independent electrical conductivity, for comparison purposes:
    sigma[ Conds ] = 5.8e7;
  EndIf

  // The electrical conductivity of the shield is assumed to be temperature
  // independent:
  sigma[ Shield ] = 3.5e7;

  // "CoefGeo[]" encodes the out-of-plane extent implicit in a 2D model:
  //   - for a planar 2D model, "CoefGeo == L", the out-of-plane length along
  //     the z-axis;
  //   - for an axisymmetric model, "CoefGeo == 2 * Pi" (the full revolution).
  // Its role in the a-v formulation, and why its sign matters, will be
  // explained below. In the present 2D model we set "L = +1 m", so
  // per-unit-length quantities are recovered directly:
  CoefGeo[] = 1;

  // Thermal parameters -- "k" is the thermal conductivity, "h" the convection
  // heat-transfer coefficient and "T0" the ambient temperature:
  k[ Conds ] = 385;
  k[ Shield ] = 205;
  k[ Insulator ] = 1;
  h[] = 60;
  T0[] = 20;

  // Tolerances and maximum iterations for the magneto-thermal fixed-point loop,
  // in the nonlinear case:
  NLTolAbs = 1e-12;
  NLTolRel = 1e-6;
  NLIterMax = 20;
}

Constraint {
  // Homogeneous Dirichlet condition on the magnetic vector potential at infinity
  // (as in tutorials 3 and 4):
  { Name a_Mag_2D;
    Case {
      { Region Surface_Inf; Value 0; }
    }
  }

  // Imposed total current in each conductor. The three-phase balanced currents
  // are specified with 120 degree phase shifts using the built-in function
  // "F_Cos_wt_p[]{omega, phi}", which is interpreted as "Complex[Cos[phi],
  // Sin[phi]]" in the frequency domain (and as "Cos[omega * $Time + phi]" in
  // the time domain -- see tutorial 4 for a similar mechanism with
  // "F_Sin_wt_p[]"):
  { Name Current_Mag_2D;
    Case {
      { Region Cond_1; Value Current;
        TimeFunction F_Cos_wt_p[]{2 * Pi * f, 0}; }
      { Region Cond_2; Value Current;
        TimeFunction F_Cos_wt_p[]{2 * Pi * f, -2 * Pi / 3}; }
      { Region Cond_3; Value Current;
        TimeFunction F_Cos_wt_p[]{2 * Pi * f, 2 * Pi / 3}; }
    }
  }

  // The shield is grounded, i.e. zero voltage drop along the z-axis. This
  // allows induced (eddy) currents to flow freely in the shield, driven by the
  // time-varying magnetic flux:
  { Name Voltage_Mag_2D;
    Case {
      { Region Shield; Value 0; }
    }
  }

  // Initial temperature for the thermal problem:
  { Name T_The; Type Init;
    Case {
      { Region Vol_The; Value T0[]; }
    }
  }
}

// Import the values of "rInt" and "rExt" for the infinite shell Jacobian:
Include "busbar_common.pro";

Jacobian {
  { Name Vol;
    Case {
      // Use the infinite shell Jacobian on the annular region:
      { Region Vol_Inf_Mag; Jacobian VolSphShell{rInt, rExt}; }
      { Region All; Jacobian Vol; }
    }
  }
  { Name Sur;
    Case {
      { Region All; Jacobian Sur; }
    }
  }
}

Integration {
  { Name Int;
    Case {
      { Type Gauss;
        Case {
          { GeoElement Line; NumberOfPoints 4; }
          { GeoElement Triangle; NumberOfPoints 4; }
        }
      }
    }
  }
}

Group {
  Dom_Hcurl_a_Mag_2D = Region[ {Vol_Mag, Sur_Neu_Mag} ];
  Dom_H1_T_The = Region[ {Vol_The, Sur_Neu_The, Sur_Rob_The} ];
}

FunctionSpace {
  // Magnetic vector potential, discretized as in tutorials 3 and 4 with
  // perpendicular edge basis functions "BF_PerpendicularEdge" associated with
  // the nodes of the mesh:
  { Name Hcurl_a_Mag_2D; Type Form1P;
    BasisFunction {
      { Name se; NameOfCoef ae; Function BF_PerpendicularEdge;
        Support Dom_Hcurl_a_Mag_2D; Entity NodesOf[All]; }
    }
    Constraint {
      { NameOfCoef ae; EntityType NodesOf; NameOfConstraint a_Mag_2D; }
    }
  }

  // The gradient of the electric scalar potential, "grad v", is discretized
  // with one global degree of freedom per conducting region using the
  // "BF_RegionZ" basis function. "BF_RegionZ" is the unit vector along the
  // z-axis -- constant and equal to "(0, 0, 1)" -- on each region in
  // "Vol_C_Mag", and zero everywhere else. It is the vector counterpart of the
  // scalar "BF_GroupOfNodes" basis function introduced in tutorial 6 for global
  // quantities.
  //
  // The finite element expansion of "grad v" thus reads
  //
  //   grad v = Sum_{r in Vol_C_Mag} ur_r sr_r(x, y)
  //
  // where "ur_r" is the (constant, complex-valued) z-component of the gradient
  // of the electric potential in region r, and "sr_r(x, y)" is the unit vector
  // along z, equal to 1 in region r and 0 elsewhere.
  //
  // As with the electrostatic global quantities in tutorial 6, the voltage
  // ("AliasOf") and the current ("AssociatedWith") are energy duals. The rule
  // of thumb when choosing which gets which type: the quantity that plays the
  // role of the primary degree of freedom -- the one the variational
  // formulation uses as test function -- is the "AliasOf"; its energy-dual (the
  // quantity that appears in the bilinear form paired with the test function)
  // is the "AssociatedWith". Accordingly, the GlobalTerm below pairs "Dof{I}"
  // (AssociatedWith) with "{U}" (AliasOf), exactly the way "GlobalTerm
  // [-Dof{Q}, {V}]" pairs charge with voltage in tutorial 6.
  { Name Hcurl_u_Mag_2D; Type Form1P;
    BasisFunction {
      // constant vector over the region, with only nonzero z-component
      { Name sr; NameOfCoef ur; Function BF_RegionZ;
        Support Vol_C_Mag; Entity Vol_C_Mag; }
    }
    GlobalQuantity {
      { Name Voltage; Type AliasOf; NameOfCoef ur; }
      { Name Current; Type AssociatedWith; NameOfCoef ur; }
    }
    Constraint {
      { NameOfCoef Voltage; EntityType Region; NameOfConstraint Voltage_Mag_2D; }
      { NameOfCoef Current; EntityType Region; NameOfConstraint Current_Mag_2D; }
    }
  }

  // Temperature, discretized with standard nodal (Lagrange) basis functions (as
  // in tutorial 2):
  { Name H1_T_The; Type Form0;
    BasisFunction {
      { Name sn; NameOfCoef Tn; Function BF_Node; Support Dom_H1_T_The;
        Entity NodesOf[All]; }
    }
    Constraint {
      { NameOfCoef Tn; EntityType NodesOf; NameOfConstraint T_The; }
    }
  }

}

Formulation {
  // Magneto-quasistatic a-v formulation.
  //
  // Using the constitutive law "h = nu b" together with "b = curl a" in
  // Ampere's equation, and plugging Ohm's law "j = sigma e" with "e =
  // -\partial_t a - grad v", we obtain
  //
  //   curl(nu curl a) = j = sigma e = -sigma \partial_t a - sigma grad v.
  //
  // In each conducting region, "grad v" is assumed to be constant, and zero
  // everywhere else. The weak a-v formulation then reads: find "a" and "v" such
  // that
  //
  //   (nu curl a, curl a')_Vol_Mag + (sigma \partial_t a, a')_Vol_C_Mag
  //     + (sigma grad v, a')_Vol_C_Mag = 0
  //
  //   - (sigma \partial_t a, ur')_Vol_C_Mag - (sigma grad v, ur')_Vol_C_Mag
  //     = I(ur')
  //
  // hold for all test functions "a'" and "ur'", with "I(ur') = I_r", the
  // current in conductor "r", for each "r" in "Vol_C_Mag". For each test
  // function "sr_r" (the global basis function for conductor "r"), the second
  // equation yields the circuit relation linking the voltage drop and the
  // current in region "r". In the first equation, as in tutorial 4, the
  // boundary term resulting from integration by parts vanishes since "a" is set
  // to 0 on the boundary of "Vol_Mag".
  //
  // In order to account for the out-of-plane length "L" of the model, we define
  // "grad v = ur / CoefGeo[]". Here "{ur}" is the discrete DoF attached to
  // "BF_RegionZ" (a constant z-vector per conducting region); dividing by
  // "CoefGeo[]" converts the region-wise voltage drop to the physical gradient
  // that appears in the PDE. "CoefGeo[]" plays two roles:
  //  - Everywhere "grad v" enters an integrand, it is written as "{ur} /
  //    CoefGeo[]", with "CoefGeo == L" in planar 2D and "CoefGeo == 2 * Pi" in
  //    axisymmetric. This is a modelling convention, not a Jacobian correction:
  //    the "Vol" (resp. "VolAxi") Jacobian already rescales each surface
  //    integrand to its 3D counterpart. So "CoefGeo" divides "{ur}" in
  //    "Integral" terms only where "grad v" appears explicitly.
  //  - The sign of "CoefGeo[]" lets the user flip the orientation of the
  //    current in a region (e.g. for a return conductor). This sign enters the
  //    circuit relation via the "GlobalTerm" with "Sign[CoefGeo[]]" below. In
  //    the present 2D model, "CoefGeo[] = 1", so the sign has no numerical
  //    effect -- it is kept so that the same formulation can be reused in
  //    models with oppositely-oriented conductors.
  { Name Magnetoquasistatics_av_2D; Type FemEquation;
    Quantity {
      { Name a; Type Local; NameOfSpace Hcurl_a_Mag_2D; }
      { Name ur; Type Local; NameOfSpace Hcurl_u_Mag_2D; }
      { Name I; Type Global; NameOfSpace Hcurl_u_Mag_2D [Current]; }
      { Name U; Type Global; NameOfSpace Hcurl_u_Mag_2D [Voltage]; }

      // Declaring "T" in this formulation gives access to the temperature field
      // from the thermal FunctionSpace, which is needed to evaluate the
      // temperature-dependent conductivity sigma:
      { Name T; Type Local; NameOfSpace H1_T_The; }
    }
    Equation {
      Integral { [ nu[] * Dof{d a} , {d a} ];
        In Vol_Mag; Jacobian Vol; Integration Int; }
      // The coefficients in the magneto-quasistatic formulation will be
      // complex-valued, while those in the thermal formulation will be
      // real-valued. The "<T>[ ... ]" syntax instructs GetDP to evaluate the
      // expression inside in real arithmetic (in the FunctionSpace of "T"),
      // even though the formulation is complex-valued:
      Integral { DtDof [ sigma[<T>[{T}]] * Dof{a} , {a} ];
        In Vol_C_Mag; Jacobian Vol; Integration Int; }
      Integral { [ sigma[<T>[{T}]] * Dof{ur} / CoefGeo[] , {a} ];
        In Vol_C_Mag; Jacobian Vol; Integration Int; }

      Integral { DtDof [ sigma[<T>[{T}]] * Dof{a} , {ur} ];
        In Vol_C_Mag; Jacobian Vol; Integration Int; }
      Integral { [ sigma[<T>[{T}]] * Dof{ur} / CoefGeo[] , {ur} ];
        In Vol_C_Mag; Jacobian Vol; Integration Int; }
      // "Sign[CoefGeo[]]" accounts for a possible sign reversal of the current
      // direction (when "CoefGeo[] < 0"):
      GlobalTerm { [ Dof{I} * Sign[CoefGeo[]] , {U} ]; In Vol_C_Mag; }
    }
  }

  // Thermal formulation.
  //
  // We proceed as in tutorial 2, but in steady state and with a volume heat
  // source "Q". The steady-state thermal conduction problem is governed by
  //
  //    - div(k grad T) = Q,
  //
  // which, with the same convection boundary condition as in tutorial 2, leads
  // to the following weak form: find "T" such that
  //
  //   (k grad T, grad T')_Vol_The + (h (T - T0), T')_Sur_Rob_The
  //     = (Q, T')_Vol_C_Mag
  //
  // holds for all test functions "T'". In our magneto-thermal setting, "Q" is
  // the time-averaged Joule loss density over a period "[0, Tp]":
  //
  //   Q = 1/Tp \int_0^Tp sigma |e(t)|^2 dt = sigma/2 |E|^2
  //
  // with "E" the complex phasor of the electric field (the factor 1/2 comes
  // from time-averaging the squared sinusoid).

  { Name Thermal_T; Type FemEquation;
    Quantity {
      { Name T; Type Local; NameOfSpace H1_T_The; }

      // Declaring "{a}" and "{ur}" here gives access to the magnetic solution
      // (computed in the frequency domain) in order to evaluate Joule losses:
      { Name a; Type Local; NameOfSpace Hcurl_a_Mag_2D; }
      { Name ur; Type Local; NameOfSpace Hcurl_u_Mag_2D; }
    }
    Equation {
      Integral { [ k[] * Dof{d T} , {d T} ];
        In Vol_The; Jacobian Vol; Integration Int; }

      Integral { [ h[] * Dof{T} , {T} ];
        In Sur_Rob_The; Jacobian Sur; Integration Int; }
      Integral { [ -h[] * T0[] , {T} ];
        In Sur_Rob_The; Jacobian Sur; Integration Int; }

      // The "<a>[ ... ]" syntax instructs GetDP to evaluate the expression
      // inside in complex arithmetic, even though the thermal formulation is
      // real-valued. Without it, only the real part of "{a}" and "{ur}" would
      // be used, which would give incorrect results. Note that, at the point
      // the thermal system is assembled, "{a}" and "{ur}" refer to the magnetic
      // solution that has already been computed (there is no "Dof" prefix),
      // i.e. they refer to known phasor fields. The expression "SquNorm[...]"
      // computes the squared modulus "|...|^2 = Re[...]^2 + Im[...]^2":
      Integral { [ -0.5 * sigma[{T}] * <a>[SquNorm[Dt[{a}] + {ur} / CoefGeo[]]],
          {T} ]; In Vol_C_Mag; Jacobian Vol; Integration Int; }
    }
  }
}

Resolution {
  // The coupled magneto-thermal resolution alternates between the magnetic and
  // thermal solves. The coupling is through:
  //  1) "sigma[{T}]": the conductivity in the magnetic problem depends on "T";
  //  2) Joule losses: the heat source "Q" in the thermal problem depends on the
  //     magnetic solution ("a", "ur").
  //
  // This is a fixed-point (Picard) iteration scheme: at each iteration, the
  // magnetic system is solved with the current temperature distribution, then
  // the thermal system is solved with the resulting Joule losses. The process
  // repeats until the magnetic residual converges (indicating self-consistency
  // of the two coupled problems).
  //
  // A monolithic coupling could also be implemented, where a single coupled
  // system is solved.
  //
  // If the NonLinear flag is not set, we disregard the temperature dependency
  // of copper conductivity and the coupled problem can be solved without
  // iterating.
  { Name MagThe;
    System {
      // The magnetic system is complex-valued and solved at a single frequency
      // (as in the frequency-domain case of tutorial 4):
      { Name Sys_Mag; NameOfFormulation Magnetoquasistatics_av_2D;
        Type Complex; Frequency f; }
      // The thermal system is real-valued:
      { Name Sys_The; NameOfFormulation Thermal_T; }
    }
    Operation {
      // Initialize the temperature to the initial condition "T0[]":
      InitSolution[Sys_The];

      // First solve: magnetic with the initial temperature, then thermal:
      Generate[Sys_Mag]; Solve[Sys_Mag];
      Generate[Sys_The]; Solve[Sys_The];

      If(NonLinear)
        // Re-generate the magnetic system with the updated temperature (which
        // changes sigma), and compute the initial residual:
        Generate[Sys_Mag];
        GetResidual[Sys_Mag, $res0];

        // Initialize runtime variables to track the residual and the iteration
        // count, then print out the absolute and relative residual:
        Evaluate[ $res = $res0, $iter = 0 ];
        Print[{$iter, $res, $res / $res0},
          Format "Residual %03g: abs %14.12e rel %14.12e"];

        // Iterate until convergence (same loop structure as in tutorial 3):
        While[$res > NLTolAbs && $res / $res0 > NLTolRel &&
          $res / $res0 <= 1 && $iter < NLIterMax]{
          Solve[Sys_Mag];
          Generate[Sys_The]; Solve[Sys_The];
          Generate[Sys_Mag]; GetResidual[Sys_Mag, $res];
          Evaluate[ $iter = $iter + 1 ];
          Print[{$iter, $res, $res / $res0},
            Format "Residual %03g: abs %14.12e rel %14.12e"];
        }
      EndIf

      SaveSolution[Sys_Mag];
      SaveSolution[Sys_The];
    }
  }
}

PostProcessing {
  // Post-processing for the magnetic formulation. Since the magnetic solution
  // is complex-valued, all quantities defined here are complex phasors:
  { Name Mag; NameOfFormulation Magnetoquasistatics_av_2D;
    Quantity {
      { Name a;
        Value {
          Term { [ {a} ]; In Vol_Mag; Jacobian Vol; }
        }
      }
      { Name az;
        Value {
          Term { [ CompZ[{a}] ]; In Vol_Mag; Jacobian Vol; }
        }
      }
      { Name b;
        Value {
          Term { [ {d a} ]; In Vol_Mag; Jacobian Vol; }
        }
      }
      { Name j;
        Value {
          Term { [ -sigma[<T>[{T}]] * (Dt[{a}] + {ur} / CoefGeo[]) ];
            In Vol_C_Mag; Jacobian Vol; }
        }
      }
      { Name JouleLosses;
        Value {
          Integral { [ 0.5 * sigma[<T>[{T}]] * SquNorm[Dt[{a}] + {ur} /
                CoefGeo[]] ]; In Vol_C_Mag; Jacobian Vol; Integration Int; }
        }
      }
      { Name U;
        Value {
          Term { [ {U} ]; In Vol_C_Mag; }
        }
      }
      { Name I;
        Value {
          Term { [ {I} ]; In Vol_C_Mag; }
        }
      }
    }
  }

  // Post-processing for the thermal formulation:
  { Name The; NameOfFormulation Thermal_T;
    Quantity {
      { Name T;
        Value {
          Term { [ {T} ]; In Vol_The; Jacobian Vol; }
        }
      }
    }
  }
}

PostOperation {
  { Name MapMag; NameOfPostProcessing Mag;
    Operation {
      Print[ az, OnElementsOf Vol_Mag, File "az.pos" ];
      Print[ b, OnElementsOf Vol_Mag, File "b.pos" ];
      Print[ j, OnElementsOf Vol_C_Mag, File "j.pos" ];
      // "RegionTable" outputs one row per region (instead of one row per
      // element), which is appropriate for global quantities:
      Print[ U, OnRegion Vol_C_Mag, File "U.txt", Format RegionTable ];
      Print[ I, OnRegion Vol_C_Mag, File "I.txt", Format RegionTable ];
    }
  }
  { Name MapThe; NameOfPostProcessing The;
    Operation {
      Print[ T, OnElementsOf Vol_The, File "T.pos" ];
    }
  }
}

// Setting a variable with the reserved name "GetDP/2PostOperationChoices"
// allows to control which PostOperations are performed by default by Gmsh
// when the model is run interactively. Here we request both "MapMag" and
// "MapThe", instead of just "MapMag" (the default):
DefineConstant[
  PostOp = {"MapMag, MapThe", Name "GetDP/2PostOperationChoices"}
];

DefineConstant[
  dim = {2, Choices{2="2D", 3="3D"},
    Name "Parameters/0Model dimension"}
  thick = {5, Min 0.1, Max 10, Step 0.1,
    Name "Parameters/Thickness [mm]"}
];

// back to meters:
mm = 1e-3;
thick = thick * mm;

// Gmsh script describing the geometry of two rectangular conductors. In 2D the
// model consists of two rectangles (cross-sections of the conductors). In 3D
// the rectangles are extruded along z to form volumes, and the bottom and top
// surfaces of each conductor are identified as separate physical groups: these
// will serve as electrodes in the 3D circuit coupling.

SetFactory("OpenCASCADE");

Include "circuit_common.pro";

DefineConstant[
  w1 = {1, Min 0.1, Max 10, Step 0.1,
    Name "Parameters/Cond1 width [mm]"}
  w2 = {2, Min 0.1, Max 10, Step 0.1,
    Name "Parameters/Cond2 width [mm]"}
  s = {1, Min 0.1, Max 10, Step 0.1,
    Name "Parameters/}Global mesh size factor"}
];

w1 = w1 * mm;
w2 = w2 * mm;
h = 5 * mm;

Rectangle(1) = {0, 0, 0, w1, h};
Rectangle(2) = {w1 + 1 * mm, 0, 0, w2, h};

If(dim == 3)
  // Extrude both rectangles along z. The "Extrude" command returns a list: the
  // first element (index 0) is the top surface, and the second (index 1) is the
  // newly created volume:
  left() = Extrude{0, 0, thick}{ Surface{1}; };
  cond1 = left(1);
  right() = Extrude{0, 0, thick}{ Surface{2}; };
  cond2 = right(1);
Else
  cond1 = 1;
  cond2 = 2;
EndIf

MeshSize{:} = 0.2 * mm * s;

bnd() = CombinedBoundary{ GeoEntity{dim}{:}; };

Physical GeoEntity{dim}("Conductor 1", 100) = cond1;
Physical GeoEntity{dim}("Conductor 2", 101) = cond2;
Physical GeoEntity{dim - 1}("Boundary", 200) = bnd();

If(dim == 3)
  // The bottom and top surfaces of each conductor are needed as electrode
  // surfaces for the 3D a-v formulation with circuit coupling. Surface tag "1"
  // is the original rectangle (bottom of cond1), and "left(0)" is the top
  // surface created by the extrusion. Similarly for cond2:
  Physical Surface("Bottom left", 300) = 1;
  Physical Surface("Top left", 301) = left(0);
  Physical Surface("Bottom right", 302) = 2;
  Physical Surface("Top right", 303) = right(0);
EndIf

// This tutorial shows how to couple a finite element model with a lumped
// element circuit. We consider two conductors (modelled with finite elements)
// connected in parallel and driven by a voltage source through a series
// resistance and inductance (modelled as lumped circuit elements). The circuit
// topology is described as follows (with circled numbers representing the
// circuit nodes):
//
//    (2)     +-----+   (3)   +-----+       (4)
//     +------+ R_1 +----+----+ L_1 +--------+------------+
//     |      +-----+         +-----+        |            |
//     |                                +----+----+  +----+----+
//  +--+--+                             |         |  |         |
//  | V_1 |                             | Cond_1  |  | Cond_2  |
//  +--+--+                             |         |  |         |
//     |                                +----+----+  +----+----+
//     |                                     |            |
//     +-------------------------------------+------------+
//    (1)
//
// The two conductors "Cond_1" and "Cond_2" are discretized with finite
// elements, while "V_1" (voltage source), "R_1" (resistance) and "L_1"
// (inductance) are lumped circuit elements.
//
// Both the 2D and the 3D cases are treated using an "a-v" magneto-quasistatic
// formulation in the conductors, which takes into account skin effect.
// Proximity effects are not modelled, as the magnetic problem is only solved in
// the conductors: there is no surrounding air region, so the two conductors are
// magnetically uncoupled.
//
// In 2D, the magnetic vector potential "a" is a perpendicular 1-form
// ("Form1P"), as in tutorials 4 and 7, and the electric scalar potential "v" is
// discretized with a single global degree of freedom per conducting region, as
// in tutorial 7. In 3D, "a" is a 1-form ("Form1") discretized with Whitney edge
// elements, as in tutorial 5. The electric scalar potential "v" is now a full
// scalar field inside each conductor, whose gradient is discretized with
// "BF_GradNode" (the gradient of the standard Lagrange basis functions) and
// "BF_GradGroupOfNodes" (the gradient of the global basis functions introduced
// in tutorial 6). An explicit gauge condition (tree-cotree gauging) is needed
// to remove the gradient fields from "a", which would otherwise be redundant
// with "grad v" (see below).
//
// Both frequency-domain and time-domain analyses are supported (as in tutorial
// 4).

Include "circuit_common.pro";

Group {
  // Physical regions (in the mesh):
  Cond_1 = Region[ 100 ];
  Cond_2 = Region[ 101 ];
  Bnd = Region[ 200 ];
  If(dim == 3)
    BotCond_1 = Region[ 300 ];
    TopCond_1 = Region[ 301 ];
    BotCond_2 = Region[ 302 ];
    TopCond_2 = Region[ 303 ];
  EndIf

  // Circuit regions. These are fictitious regions used to represent lumped
  // circuit elements. Their tags must not conflict with any physical group tag
  // in the mesh, so that GetDP can distinguish them from mesh-based regions:
  V_1 = Region[ 1000 ];
  R_1 = Region[ 1001 ];
  L_1 = Region[ 1002 ];

  // Abstract magnetic regions:
  //  - "Vol_Mag": full domain for the magnetic problem
  //  - "Vol_C_Mag": conducting regions
  //  - "Sur_Electrodes_Mag": electrode surfaces with voltages and currents in
  //    3D
  Vol_Mag = Region[ {Cond_1, Cond_2} ];
  Vol_C_Mag = Region[ {Cond_1, Cond_2} ];
  If(dim == 3)
    Sur_Electrodes_Mag = Region[ {BotCond_1, TopCond_1, BotCond_2, TopCond_2} ];
  EndIf

  // Abstract circuit regions. These groups allow a generic treatment of the
  // lumped element constitutive relations in the Formulation (see the
  // GlobalTerms below):
  //  - "Resistance_Cir": resistors (U = R I)
  //  - "Inductance_Cir": inductors (U = L dI/dt)
  //  - "Capacitance_Cir": capacitors (I = C dU/dt)
  //  - "Dom_Cir": all lumped circuit elements
  Resistance_Cir = Region[ {R_1} ];
  Inductance_Cir = Region[ {L_1} ];
  Capacitance_Cir = Region[ {} ]; // empty, but kept for completeness
  Dom_Cir = Region[ {V_1, R_1, L_1} ];
}

Function {
  DefineConstant[
    AnalysisType = {1, Choices{0="Time-domain", 1="Frequency-domain"},
      Name "Parameters/1Analysis type"}
    Voltage = {0.1, Min 1e-3, Max 10, Step 1,
      Name "Parameters/2Voltage source [V]"}
    // In the frequency domain we consider an AC voltage source with frequency
    // "f":
    f = {1000, Min 1, Max 10000, Step 10, Visible AnalysisType == 1,
      Name "Parameters/3Frequency [Hz]"}
    // In the time domain we consider a step voltage with given rise time "tr",
    // and run the simulation for "10 * tr":
    tr = {1e-3, Min 1e-5, Max 1, Step 1e-3, Visible AnalysisType == 0,
      Name "Parameters/3Rise time"}
    Rval = {10e-3, Min 0, Max 100, Step 0.1,
      Name "Parameters/4Series resistance [Ohm]"}
    Lval = {5e-6, Min 0, Max 100e-6, Step 1e-6,
      Name "Parameters/4Series inductance [H]"}
  ];

  If(AnalysisType == 1)
    tfct[] = F_Cos_wt_p[]{2 * Pi * f, 0};
  Else
    tmax = 10 * tr;
    dt = tr / 10;
    // A simple ramp function for the step voltage:
    tfct[] = ($Time < tr) ? $Time / tr : 1;
  EndIf

  // Magnetic parameters:
  mu0 = 4 * Pi * 1e-7;
  nu[] = 1. / mu0;
  sigma[ Cond_1 ] = 5.8e7;
  sigma[ Cond_2 ] = 5.8e7;

  If(dim == 2)
    // In 2D, "CoefGeo[]" represents the out-of-plane thickness (as in tutorial
    // 7):
    CoefGeo[] = thick;
  Else
    // In 3D the thickness is part of the geometrical model and the mesh, so
    // "CoefGeo[] = 1":
    CoefGeo[] = 1;
  EndIf

  // Circuit parameters: the "Resistance[]", "Inductance[]" and "Capacitance[]"
  // functions are piecewise constants defined on the circuit regions, used in
  // the lumped element GlobalTerms of the Formulation:
  Resistance[ R_1 ] = Rval;
  Inductance[ L_1 ] = Lval;
  Capacitance[ ] = 0; // unused
}

Constraint {
  // Homogeneous Dirichlet condition on the magnetic vector potential on the
  // outer boundaries:
  { Name a_Mag;
    Case {
      { Region Bnd; Value 0.; }
    }
  }

  // Tree-cotree gauge condition for the 3D vector potential (see the
  // FunctionSpace definition below for details):
  { Name a_Gauge_Mag;
    Case {
      { Region Vol_C_Mag ; SubRegion Bnd; Value 0.; }
    }
  }

  // No imposed current constraint on the FE conductors (the current will be
  // determined by the circuit):
  { Name Current_Mag;
    Case {
    }
  }

  { Name Voltage_Mag;
    Case {
      If(dim == 3)
        // In 3D, set the absolute voltage on the bottom electrode of "Cond_1"
        // to 0 as a reference:
        { Region BotCond_1; Value 0; }
      EndIf
    }
  }

  // Initial condition on the circuit inductor current (for the time-domain
  // analysis):
  { Name Current_Cir ;
    Case {
      { Region L_1; Type Init; Value 0.; }
    }
  }

  // Imposed voltage on the voltage source:
  { Name Voltage_Cir ;
    Case {
      { Region V_1; Value Voltage; TimeFunction tfct[]; }
    }
  }

  // The circuit topology is described as a constraint of type "Network". Each
  // entry specifies a circuit element (a "Region") and the two circuit nodes
  // ("Branch") it connects, defining a directed branch from the first node to
  // the second. The circuit node numbers (1, 2, 3, 4, ...) are arbitrary labels
  // that refer to the ASCII art diagram at the top of this file.
  //
  // The sign convention is the same for all elements (sources and loads alike):
  // for a branch "Branch {a, b}", the current "I" is positive from node "a" to
  // node "b", and the voltage "U" is the potential rise from "a" to "b", i.e.
  // "U = V(b) - V(a)". This uniform "source" (or "generator") convention
  // differs from the "passive sign convention" commonly used in circuit
  // textbooks, where the current through a load is defined as flowing into the
  // positive voltage terminal.
  //
  // The constitutive law of each element determines the relation between "U"
  // and "I". With the source convention:
  //   - Resistor: U + R I = 0 (voltage drops in the direction of current)
  //   - Inductor: U + L dI/dt = 0
  //   - Capacitor: I + C dU/dt = 0
  //   - Voltage source: U = Voltage * tfct[] (via the "Voltage_Cir" constraint)
  //
  // The power delivered by any element is "P = U I": sources have "P > 0" and
  // loads have "P < 0" (e.g. for a resistor, "P = U I = -R I^2 < 0").
  { Name ElectricalCircuit ; Type Network ;
    Case Circuit_1 {
      // Lumped elements:
      { Region V_1; Branch {1, 2}; }
      { Region R_1; Branch {2, 3}; }
      { Region L_1; Branch {3, 4}; }

      If(dim == 2)
        // In 2D, each conducting region has a single voltage drop and current
        // flowing through it (in the direction perpendicular to the plane). It
        // behaves as a two-terminal circuit element and can be directly
        // inserted between two circuit nodes. Here both conductors are
        // connected in parallel between nodes 4 and 1:
        { Region Cond_1; Branch {4, 1} ; }
        { Region Cond_2; Branch {4, 1} ; }
      Else
        // In 3D, each electrode surface carries a voltage (absolute potential)
        // and a current. To impose a voltage drop across a conductor (say
        // "Cond_1") in the circuit -- rather than an absolute voltage on each
        // electrode surface independently -- we use an "anti-serial"
        // connection.  The key idea is to connect the top and bottom electrodes
        // with opposite node orderings, so that their absolute voltages are
        // subtracted:
        //
        //   V_drop(Cond_1) = V(TopCond_1) - V(BotCond_1).
        //
        // This is achieved by connecting "TopCond_1" from node 4 to an
        // intermediate node 10, and "BotCond_1" from node 1 to the same node 10
        // (note the reversed direction). The voltage between nodes 4 and 1 then
        // corresponds to the voltage drop over "Cond_1", regardless of the
        // absolute voltage values. The current convention is also consistent:
        // the current flowing into "TopCond_1" equals the current flowing out
        // of "BotCond_1".
        //
        // Note: in this particular circuit, since "BotCond_1" is fixed at 0 V
        // through the "Voltage_Mag" constraint, the anti-serial connection for
        // the bottom face is not strictly necessary. However, it is included
        // here for generality. Beware that if the bottom face is removed from
        // the Network, the corresponding region "BotCond_1" has to be removed
        // from "Sur_Electrodes_Mag" as well!
        { Region TopCond_1; Branch {4, 10} ; }
        { Region BotCond_1; Branch {1, 10} ; }

        { Region TopCond_2; Branch {4, 11} ; }
        { Region BotCond_2; Branch {1, 11} ; }
      EndIf
    }
  }
}

Jacobian {
  { Name Vol;
    Case {
      { Region All; Jacobian Vol; }
    }
  }
  { Name Sur;
    Case {
      { Region All; Jacobian Sur; }
    }
  }
}

Integration {
  { Name Int;
    Case {
      { Type Gauss;
        Case {
          { GeoElement Line; NumberOfPoints 4; }
          { GeoElement Triangle; NumberOfPoints 4; }
          { GeoElement Tetrahedron; NumberOfPoints 4; }
        }
      }
    }
  }
}

FunctionSpace {
  If(dim == 2)
    // In 2D, the magnetic vector potential and the gradient of the electric
    // scalar potential are discretized as in tutorial 7:
    { Name Hcurl_a_Mag_2D; Type Form1P;
      BasisFunction {
        { Name se; NameOfCoef ae; Function BF_PerpendicularEdge;
          Support Vol_Mag; Entity NodesOf[All]; }
      }
      Constraint {
        { NameOfCoef ae; EntityType NodesOf; NameOfConstraint a_Mag; }
      }
    }
    { Name Hcurl_u_Mag_2D; Type Form1P;
      BasisFunction {
        { Name sr; NameOfCoef ur; Function BF_RegionZ;
          Support Vol_C_Mag; Entity Vol_C_Mag; }
      }
      GlobalQuantity {
        { Name Voltage; Type AliasOf; NameOfCoef ur; }
        { Name Current; Type AssociatedWith; NameOfCoef ur; }
      }
      Constraint {
        { NameOfCoef Voltage; EntityType Region; NameOfConstraint Voltage_Mag; }
        { NameOfCoef Current; EntityType Region; NameOfConstraint Current_Mag; }
      }
    }
  Else
    // In 3D, the magnetic vector potential is a true 1-form (Form1), as in
    // tutorial 5, discretized with Whitney edge basis functions "BF_Edge"
    // associated with the edges of the mesh:
    { Name Hcurl_a_Mag_3D; Type Form1;
      BasisFunction {
        { Name se; NameOfCoef ae; Function BF_Edge;
          Support Vol_Mag; Entity EdgesOf[All]; }
      }
      Constraint {
        { NameOfCoef ae; EntityType EdgesOf; NameOfConstraint a_Mag; }
        // Tree-cotree gauge condition for 3D. In the conducting regions, the
        // electric field "e = -\partial_t a - grad v" involves two unknowns,
        // "a" and "grad v", but only their sum matters physically. Since the
        // edge element space of "a" can represent arbitrary gradient fields,
        // any gradient could be moved freely between "a" and "grad v" without
        // changing "e": the system is thus singular without an additional
        // constraint.
        //
        // Tree-cotree gauging removes gradient fields from "a" by setting to
        // zero the degrees of freedom of "a" on the edges of a spanning tree of
        // the mesh graph in "Vol_C_Mag". The "EdgesOfTreeIn" entity type
        // selects such a tree, rooted on the boundary "Bnd" (specified by
        // "EntitySubType StartingOn") where "a" is already constrained to
        // zero. This action eliminates the irrotational (gradient) components
        // of "a", leaving only the solenoidal part; "ur" alone will carry the
        // gradient part of the electric field.
        //
        // See tutorial 10 for more details about gauging.
        { NameOfCoef ae; EntityType EdgesOfTreeIn; EntitySubType StartingOn;
          NameOfConstraint a_Gauge_Mag; }
      }
    }

    // In 3D, the electric scalar potential "v" is a full field (not just a
    // single DOF per region as in 2D). Its gradient "grad v" is discretized in
    // H(curl) with two types of basis functions:
    //  - "BF_GradNode": the gradient of the standard Lagrange (nodal) basis
    //    functions, associated with the interior nodes of each conductor
    //    (excluding the electrode surfaces). These capture the spatial
    //    variation of "v" inside the conductor.
    //  - "BF_GradGroupOfNodes": the gradient of the global basis functions
    //    associated with the electrode surfaces (as in tutorial 6 with
    //    "BF_GroupOfNodes", but here through the gradient). Each electrode
    //    carries a single global voltage and current.
    //
    // The resulting expansion of "grad v" in conductor "r" reads:
    //
    //   grad v(x,y,z) = Sum_k un_k grad(sn_k(x, y, z))
    //     + Sum_electrode ur_electrode grad(sr_electrode(x,y,z))
    //
    // where the first sum runs over interior nodes and the second over
    // electrode surfaces. Note that we could alternatively have defined the
    // scalar potential "v" directly (as a Form0 with "BF_Node" and
    // "BF_GroupOfNodes"), but then the Formulation would need to be written
    // differently in 2D and 3D. By defining "grad v" directly as a Form1, the
    // same Formulation can be shared between 2D and 3D (see below).
    { Name Hcurl_u_Mag_3D; Type Form1;
      BasisFunction {
        { Name sn; NameOfCoef un; Function BF_GradNode;
          Support Vol_C_Mag; Entity NodesOf[Vol_C_Mag, Not Sur_Electrodes_Mag]; }
        { Name sr; NameOfCoef ur; Function BF_GradGroupOfNodes;
          Support Vol_C_Mag; Entity GroupsOfNodesOf[Sur_Electrodes_Mag]; }
      }
      GlobalQuantity {
        { Name Voltage; Type AliasOf; NameOfCoef ur; }
        { Name Current; Type AssociatedWith; NameOfCoef ur; }
      }
      Constraint {
        { NameOfCoef Voltage; EntityType GroupsOfNodesOf;
          NameOfConstraint Voltage_Mag; }
        { NameOfCoef Current; EntityType GroupsOfNodesOf;
          NameOfConstraint Current_Mag; }
      }
    }
  EndIf

  // The circuit currents and voltages are discretized with a scalar ("Type
  // Scalar") function space using "BF_Region" basis functions. "BF_Region" is a
  // region-wise constant scalar function (equal to 1 in its associated region
  // and 0 elsewhere), which is the scalar counterpart of the vector
  // "BF_RegionZ" used in 2D. Each lumped circuit element carries one current
  // DoF and one voltage DoF (its energy dual):
  { Name Hregion_i_Cir; Type Scalar;
    BasisFunction {
      { Name sr; NameOfCoef ir; Function BF_Region;
        Support Dom_Cir; Entity Dom_Cir; }
    }
    GlobalQuantity {
      { Name Current; Type AliasOf; NameOfCoef ir; }
      { Name Voltage; Type AssociatedWith; NameOfCoef ir; }
    }
    Constraint {
      { NameOfCoef Voltage; EntityType Region; NameOfConstraint Voltage_Cir; }
      { NameOfCoef Current; EntityType Region; NameOfConstraint Current_Cir; }
    }
  }

}

Formulation {
  // The formulation combines the FE magneto-quasistatic a-v weak formulation (as
  // in tutorial 7) with the lumped element constitutive relations and the
  // circuit topology (Kirchhoff's laws). All these contributions are assembled
  // into a single linear system.
  { Name Magnetoquasistatics_av; Type FemEquation;
    Quantity {
      // FE degrees of freedom (2D or 3D):
      If(dim == 2)
        { Name a; Type Local; NameOfSpace Hcurl_a_Mag_2D; }
        { Name ur; Type Local; NameOfSpace Hcurl_u_Mag_2D; }
        { Name I; Type Global; NameOfSpace Hcurl_u_Mag_2D [Current]; }
        { Name U; Type Global; NameOfSpace Hcurl_u_Mag_2D [Voltage]; }
      Else
        { Name a; Type Local; NameOfSpace Hcurl_a_Mag_3D; }
        { Name ur; Type Local; NameOfSpace Hcurl_u_Mag_3D; }
        { Name I; Type Global; NameOfSpace Hcurl_u_Mag_3D [Current]; }
        { Name U; Type Global; NameOfSpace Hcurl_u_Mag_3D [Voltage]; }
      EndIf

      // Circuit degrees of freedom:
      { Name Iz; Type Global; NameOfSpace Hregion_i_Cir [Current]; }
      { Name Uz; Type Global; NameOfSpace Hregion_i_Cir [Voltage]; }
    }
    Equation {
      // FE terms for the a-v formulation (as in tutorial 7):
      Integral { [ nu[] * Dof{d a} , {d a} ];
        In Vol_Mag; Jacobian Vol; Integration Int; }
      Integral { DtDof [ sigma[] * Dof{a} , {a} ];
        In Vol_C_Mag; Jacobian Vol; Integration Int; }
      Integral { [ sigma[] * Dof{ur} / CoefGeo[] , {a} ];
        In Vol_C_Mag; Jacobian Vol; Integration Int; }

      Integral { DtDof [ sigma[] * Dof{a} , {ur} ];
        In Vol_C_Mag; Jacobian Vol; Integration Int; }
      Integral { [ sigma[] * Dof{ur} / CoefGeo[] , {ur} ];
        In Vol_C_Mag; Jacobian Vol; Integration Int; }

      // In 2D the global current-voltage relation is per conducting region; in
      // 3D it is per electrode surface:
      If(dim == 2)
        GlobalTerm { [ Dof{I} * Sign[CoefGeo[]] , {U} ]; In Vol_C_Mag; }
      Else
        GlobalTerm { [ Dof{I} * Sign[CoefGeo[]] , {U} ]; In Sur_Electrodes_Mag; }
      EndIf

      // Lumped element constitutive relations, expressed as GlobalTerms. For
      // each lumped element, we have two unknowns: the current "Iz" and its
      // energy dual, the voltage "Uz". Two equations are needed per element:
      // one is the constitutive relation (below), the other comes from
      // Kirchhoff's laws (enforced by the GlobalEquation).
      //
      // All constitutive relations use the source sign convention defined in the
      // "ElectricalCircuit" network constraint above:
      //
      // Resistor: Uz + R * Iz = 0
      GlobalTerm { [ Dof{Uz} , {Iz} ]; In Resistance_Cir; }
      GlobalTerm { [ Resistance[] * Dof{Iz} , {Iz} ]; In Resistance_Cir; }

      // Inductor: Uz + L * dIz/dt = 0
      GlobalTerm { [ Dof{Uz} , {Iz} ]; In Inductance_Cir; }
      GlobalTerm { DtDof [ Inductance[] * Dof{Iz} , {Iz} ]; In Inductance_Cir; }

      // Capacitor: Iz + C * dUz/dt = 0 (unused here, but included for
      // completeness)
      GlobalTerm { [ Dof{Iz} , {Iz} ]; In Capacitance_Cir; }
      GlobalTerm { DtDof [ Capacitance[] * Dof{Uz} , {Iz} ]; In Capacitance_Cir; }

      // The "GlobalEquation" block enforces Kirchhoff's current and voltage
      // laws ("KCL" and "KVL", respectively), as specified by the
      // "ElectricalCircuit" network constraint. It couples the FE global
      // quantities (I, U) with the circuit global quantities (Iz, Uz) through
      // Kirchhoff's laws at each circuit node.
      //
      // Each entry in the GlobalEquation has:
      //  - "Node {x}": the current variable contributing to KCL at each node
      //  - "Loop {x}": the voltage variable contributing to KVL around each
      //    loop
      //  - "Equation {x}": the variable whose constitutive equation has been
      //    written (in the GlobalTerms above)
      //  - "In ...": the region(s) where this entry applies
      GlobalEquation {
        Type Network; NameOfConstraint ElectricalCircuit;
        If(dim == 2)
          { Node {I};  Loop {U};  Equation {I};  In Vol_C_Mag; }
        Else
          { Node {I};  Loop {U};  Equation {I};  In Sur_Electrodes_Mag; }
        EndIf
        { Node {Iz}; Loop {Uz}; Equation {Uz}; In Dom_Cir; }
      }
    }
  }
}

Resolution {
  { Name Mag;
    System {
      { Name Sys_Mag; NameOfFormulation Magnetoquasistatics_av;
        If(AnalysisType == 1)
          Type Complex; Frequency f;
        EndIf
      }
    }
    Operation {
      InitSolution[Sys_Mag];
      If(AnalysisType == 1)
        Generate[Sys_Mag]; Solve[Sys_Mag]; SaveSolution[Sys_Mag];
      Else
        TimeLoopTheta[0, tmax, dt, 1] {
          Generate[Sys_Mag]; Solve[Sys_Mag]; SaveSolution[Sys_Mag];
        }
      EndIf
    }
  }
}

PostProcessing {
  { Name Mag; NameOfFormulation Magnetoquasistatics_av;
    Quantity {
      { Name a;
        Value {
          Term { [ {a} ]; In Vol_Mag; Jacobian Vol; }
        }
      }
      { Name b;
        Value {
          Term { [ {d a} ]; In Vol_Mag; Jacobian Vol; }
        }
      }
      { Name j;
        Value {
          Term { [ -sigma[] * (Dt[{a}] + {ur} / CoefGeo[]) ];
            In Vol_C_Mag; Jacobian Vol; }
        }
      }
      // Global quantities can be evaluated both on FE regions and on circuit
      // regions, using piecewise definitions:
      { Name U;
        Value {
          If(dim == 2)
            Term { [ {U} ]; In Vol_C_Mag; }
          Else
            Term { [ {U} ]; In Sur_Electrodes_Mag; }
          EndIf
          Term { [ {Uz} ]; In Dom_Cir; }
        }
      }
      { Name I;
        Value {
          If(dim == 2)
            Term { [ {I} ]; In Vol_C_Mag; }
          Else
            Term { [ {I} ]; In Sur_Electrodes_Mag; }
          EndIf
          Term { [ {Iz} ]; In Dom_Cir; }
        }
      }
    }
  }
}

PostOperation {
  { Name MapMag; NameOfPostProcessing Mag;
    Operation {
      Print[ a, OnElementsOf Vol_Mag, File "a.pos" ];
      Print[ b, OnElementsOf Vol_Mag, File "b.pos" ];
      Print[ j, OnElementsOf Vol_C_Mag, File "j.pos" ];
      If(dim == 2)
        Print[ U, OnRegion Vol_C_Mag, File "U.txt", Format RegionTable ];
        Print[ I, OnRegion Vol_C_Mag, File "I.txt", Format RegionTable ];
      Else
        Print[ U, OnRegion Sur_Electrodes_Mag, File "U.txt", Format RegionTable ];
        Print[ I, OnRegion Sur_Electrodes_Mag, File "I.txt", Format RegionTable ];
      EndIf
      Print[ U, OnRegion Dom_Cir, File "Uz.txt", Format RegionTable ];
      Print[ I, OnRegion Dom_Cir, File "Iz.txt", Format RegionTable ];

      // "SendToServer" sends the value to the Gmsh graphical user interface for
      // display in the left panel. The "}" prefix in the server name places
      // these entries at the bottom of the parameter tree:
      Print[ U, OnRegion V_1, File "", Format Table,
        SendToServer "}Output/V (V_1) [V]" ];
      Print[ I, OnRegion V_1, File "", Format Table,
        SendToServer "}Output/I (V_1) [A]" ];
      Print[ U, OnRegion L_1, File "", Format Table,
        SendToServer "}Output/V (L_1) [V]" ];
      Print[ I, OnRegion L_1, File "", Format Table,
        SendToServer "}Output/I (L_1) [A]" ];
    }
  }
}

// Parameters shared by Gmsh and GetDP: core and coil dimensions.

mm = 1e-3;
wCore = 100 * mm;
hCore = 100 * mm;
wCoil1 = 10 * mm;
hCoil1 = 40 * mm;
wCoil2 = 20 * mm;
hCoil2 = 40 * mm;

DefineConstant[
  thickCore = {30, Min 5, Max 60, Step 0.1,
    Name "Parameters/Core thickness [mm]", Visible 0}
  gapCoil = {2, Min 0.5, Max 7.5, Step 0.1,
    Name "Parameters/Coil gap [mm]"}
  wCoreLeg = {20, Min 10, Max 27, Step 0.1,
    Name "Parameters/Core width [mm]"}
];

thickCore = thickCore * mm;
gapCoil = gapCoil * mm;
wCoreLeg = wCoreLeg * mm;

// Inner and outer radii for the infinite shell Jacobian (see tutorial 3),
// chosen to be large enough to enclose the full geometry:
r = Max[(wCore + 2 * gapCoil + wCoil1 + wCoil2) / 2, hCore / 2];
JacRadiusInt = r * Sqrt[2] + 30 * mm;
JacRadiusExt = JacRadiusInt + 50 * mm;

// Gmsh script describing the geometry of a single-phase transformer: a
// ferromagnetic core with a primary winding (two coil sides) and a secondary
// winding (two coil sides), surrounded by an air region with an infinite shell
// for the unbounded domain (see tutorial 3).

SetFactory("OpenCASCADE");

Include "transfo_common.pro";

// The core is is created by subtracting the inner rectangle (the window) from
// the outer rectangle:
Rectangle(1) = {-wCore / 2, -hCore / 2, 0,
  wCore, hCore};
Rectangle(2) = {-wCore / 2 + wCoreLeg, -hCore / 2 + wCoreLeg, 0,
  wCore - 2 * wCoreLeg, hCore - 2 * wCoreLeg};

// Primary winding: two coil sides (minus and plus) on either side of the left
// core leg. "Coil1M" carries the current in the -z direction and "Coil1P" in
// the +z direction (the go and return paths of the same winding):
Rectangle(3) = {-wCore / 2 - gapCoil - wCoil1, -hCoil1 / 2, 0,
  wCoil1, hCoil1};
Rectangle(4) = {-wCore / 2 + wCoreLeg + gapCoil, -hCoil1 / 2, 0,
  wCoil1, hCoil1};

// Secondary winding: two coil sides on either side of the right core leg:
Rectangle(5) = {wCore / 2 - wCoreLeg - gapCoil - wCoil2, -hCoil2 / 2, 0,
  wCoil2, hCoil1};
Rectangle(6) = {wCore / 2 + gapCoil, -hCoil2 / 2, 0,
  wCoil2, hCoil2};

// Air region (disk) and infinite ring:
Disk(7) = {0, 0, 0, JacRadiusInt};
Disk(8) = {0, 0, 0, JacRadiusExt};

// Make all overlapping surfaces conformal (see tutorial 3):
Coherence;

// Use the "Closest" command to retrieve the regions that were split by the
// boolean operation in "Coherence":
airinf() = Closest {JacRadiusInt + 1e-6, 0, 0} { Surface{:}; };
core() = Closest {wCore / 2 - 1e-6, 0, 0} { Surface{:}; };

// The air region is everything except the infinite shell, the core and the
// coils:
air() = Surface{:};
air() -= {airinf(0), core(0), 3, 4, 5, 6};

bnd() = CombinedBoundary{ Surface{:}; };

DefineConstant[
  s = {1, Min 0.1, Max 10, Step 0.1,
    Name "Parameters/}Global mesh size factor"}
];

MeshSize{:} = 2 * mm * s;
MeshSize{ PointsOf{ Surface{airinf(0)}; } } = 8 * mm * s;

Physical Surface("Air", 1) = air();
Physical Surface("AirInf", 2) = airinf(0);
Physical Surface("Coil1M", 3) = 3;
Physical Surface("Coil1P", 4) = 4;
Physical Surface("Coil2M", 5) = 5;
Physical Surface("Coil2P", 6) = 6;
Physical Surface("Core", 7) = core(0);
Physical Curve("Bnd", 10) = bnd();

// This tutorial models a single-phase transformer with two windings (primary
// and secondary) wound around a laminated ferromagnetic core, using a
// frequency-domain magneto-quasistatic a-v formulation coupled with two
// external circuits.
//
// Whereas previous tutorials built the complete formulation, function spaces,
// resolution and post-processing from scratch, this tutorial illustrates the
// use of a generic template library ("Lib_Magnetoquasistatics_av_2D_Cir.pro").
// The library provides a ready-made a-v formulation that handles massive
// conductors, stranded coils, nonlinear materials, permanent magnets, impedance
// boundary conditions, and circuit coupling. The user only needs to define the
// physical and abstract groups, the material properties, the constraints, and a
// few configuration variables -- the rest is taken care of by the library.
//
// The main new modeling concepts introduced in the tutorial are the modeling of
// stranded coils and of the ferromagnetic laminations in the core.
//
// Unlike the massive (solid) conductors of tutorials 7 and 8, where eddy
// currents are resolved inside the conductor cross-section, a stranded coil is
// modeled as a homogenized region carrying a uniform current density
//
//   j = Ns / Sc * I * z_hat,
//
// where "Ns" is the number of turns, "Sc" the cross-sectional area of the
// winding, "I" the total current (the same in each turn), and "z_hat" the unit
// vector along the conductor direction. This approximation is valid when the
// coil consists of a large number of turns in series with each other (so that
// the each turn carries the same current), where each individual turn is thin
// compared to the skin depth (so that the current distribution within each turn
// is essentially uniform).

// Each winding consists of two coil sides (the "go" and "return" halves of the
// same winding), modeled as separate regions with opposite current
// directions. The two windings are connected to two independent circuits:
//
//   Circuit 1 (primary):   voltage source E_in -> Coil_1_M -> Coil_1_P
//   Circuit 2 (secondary): Coil_2_M -> Coil_2_P -> load impedance Z_out
//
// The magnetic coupling between the two circuits occurs through the shared
// magnetic flux in the core -- there is no explicit electrical connection
// between them, i.e., they are galvanically isolated. The load connected at the
// terminals of the secondary can be a resistance, an inductance, or a
// capacitance, selected at runtime through a "DefineConstant" parameter.
//
// We define an effective complex reluctivity "nu = nu_r + i * nu_i = nu_r + i *
// (nu_H + nu_EC)" to model the hysteresis and eddy current losses in the
// laminated ferromagnetic core. The imaginary part of the reluctivity
// introduces a dissipative contribution "(1/2) * omega * (nu_H + nu_EC) *
// |B_peak|^2" per unit volume.
//
//  - Hysteresis losses: Each magnetisation cycle dissipates approximately
//
//      p_H = kh * |B_peak|^2,
//
//    with an empirical coefficient kh. Matching with "(1/2) * omega * nu_H *
//    |B_peak|^2p_H" gives
//
//      nu_H = kh / Pi.
//
//    This equivalent reluctivity is frequency-independent and hysteresis losses
//    p_H are proportional to the frequency "f" (through "omega").
//
//  - Eddy current losses: in a lamination of thickness "d" and bulk
//    conductivity "sigma_iron", eddy-current losses are approximately given by
//    (see e.g. J. Gyselinck et al., "Calculation of eddy currents and
//    associated losses in electrical steel laminations", IEEE Transactions on
//    Magnetics, 35(3), 1191-1194, 1999):
//
//      p_EC = (1/2) * sigma_iron * (d * omega * |B_peak|)^2 / 12.
//
//    Matching with "(1/2) * omega * nu_EC * |B_peak|^2" gives
//
//      nu_EC = sigma_iron * omega * d^2 / 12.
//
//    This equivalent reluctivity is proportional to the frequency, and
//    eddy-current losses "p_EC" thus scale as "f^2".

Include "transfo_common.pro";

Group {
  // Physical regions:
  Air = Region[ 1 ];
  AirInf = Region[ 2 ];
  Coil_1_M = Region[ 3 ];
  Coil_1_P = Region[ 4 ];
  Coil_2_M = Region[ 5 ];
  Coil_2_P = Region[ 6 ];
  Coil_1 = Region[ {Coil_1_M, Coil_1_P} ];
  Coil_2 = Region[ {Coil_2_M, Coil_2_P} ];
  Coils = Region[ {Coil_1, Coil_2} ];
  Core = Region[ 7 ];
  Bnd = Region[ 10 ]; // exterior boundary

  // Abstract regions expected by the template library. The role of each group
  // is documented in "Lib_Magnetoquasistatics_av_2D_Cir.pro":
  //  - "Vol_Mag": full magnetic domain
  //  - "Vol_S_Mag": stranded conductors (windings with Ns turns)
  //  - "Vol_Inf_Mag": region with infinite shell transformation
  Vol_Mag = Region[ {Air, AirInf, Core, Coils} ];
  Vol_S_Mag = Region[ Coils ];
  Vol_Inf_Mag = Region[ AirInf ];

  // Circuit regions (fictitious, as in tutorial 8): the voltage source on the
  // primary side and the load impedance on the secondary side:
  E_in = Region[ 100 ];
  Z_out = Region[ 102 ];

  // Abstract circuit regions expected by the library:
  //  - "Source_Cir": all voltage and current sources
  //  - "Resistance_Cir": all resistors
  //  - "Inductance_Cir": all inductors
  //  - "Capacitance_Cir": all capacitors
  Source_Cir = Region[ {E_in} ];
  DefineConstant[
    LoadType = {0, Choices {0="Resistance", 1="Inductance", 2="Capacitance"},
      Name "Parameters/2Load type"}
  ];
  // "LoadType" determines which abstract circuit group "Z_out" belongs to.
  // Each group has a different constitutive relation in the library formulation
  // (see tutorial 8 for the lumped element GlobalTerms):
  If(LoadType == 0)
    Resistance_Cir = Region[ {Z_out} ];
  ElseIf(LoadType == 1)
    Inductance_Cir = Region[ {Z_out} ];
  Else
    Capacitance_Cir = Region[ {Z_out} ];
  EndIf
}

Function {
  DefineConstant[
    Voltage = {1, Min 1e-3, Max 10,
      Name "Parameters/0Primary voltage [V]"}
    f = {50, Min 0, Max 1e3, Step 1,
      Name "Parameters/1Frequency [Hz]"}
    Rval = {50, Min 1e-3, Max 50, Step 0.25, Visible (LoadType == 0),
      Name "Parameters/3Load resistance [Ohm]"}
    Lval = {1e-6, Min 1e-9, Max 0.1, Step 1e-3, Visible (LoadType == 1),
      Name "Parameters/3Load inductance [H]"}
    Cval = {1e-6, Min 1e-9, Max 0.1, Step 1e-4, Visible (LoadType == 2),
      Name "Parameters/3Load capacitance [F]"}
    // Typical values for grain-oriented Si-steel (M3, 0.35 mm laminations):
    murCore = {2000, Min 1, Max 10000, Step 1,
      Name "Parameters/4Core relative permeability"}
    dlam = {0.35e-3, Min 0.1e-3, Max 1e-3, Step 0.05e-3,
      Name "Parameters/5Lamination thickness [m]"}
    kh = {100, Min 10, Max 500, Step 10,
      Name "Parameters/6Hysteresis coeff [W m^-3 T^-2 s]"}
  ];

  mu0 = 4e-7 * Pi;
  nu[ Air ] = 1 / mu0;
  nu[ AirInf ] = 1 / mu0;

  // Homogenized effective reluctivity to take into account losses in the
  // laminated ferromagnetic core:
  nu_r = 1 / (murCore * mu0);
  sigma_iron = 2e6;
  nu_i = kh / Pi + sigma_iron * (2 * Pi * f) * dlam^2 / 12;
  nu[ Core ] = Complex[nu_r, nu_i];

  // Stranded coil parameters expected by the library:
  //  - "Ns[]": number of turns in each winding
  //  - "Sc[]": cross-sectional area of each coil side. "SurfaceArea[]" is a
  //    built-in GetDP function that returns the area of the region on which it
  //    is evaluated
  Ns[ Coil_1 ] = 100;
  Ns[ Coil_2 ] = 200;
  Sc[ Coil_1_M ] = SurfaceArea[];
  Sc[ Coil_1_P ] = SurfaceArea[];
  Sc[ Coil_2_M ] = SurfaceArea[];
  Sc[ Coil_2_P ] = SurfaceArea[];
  nu[ Coils ] = 1 / mu0;
  // The conductivity of the stranded coils enters the Joule losses calculation
  // in the library (via the DC resistance term "Ns^2 / (sigma * Sc^2)"):
  sigma[ Coils ] = 5.8e7;

  // Geometrical coefficient in all the regions with global voltages and
  // currents, to specify the out-of-plane thickness and control the current
  // direction (as in tutorials 7 and 8):
  CoefGeo[ Coil_1_M ] = -1 * thickCore;
  CoefGeo[ Coil_1_P ] = 1 * thickCore;
  CoefGeo[ Coil_2_M ] = -1 * thickCore;
  CoefGeo[ Coil_2_P ] = 1 * thickCore;
  CoefGeo[ Core ] = thickCore;

  // Lumped circuit element parameters (as in tutorial 8):
  Resistance[ Z_out ] = Rval;
  Inductance[ Z_out ] = Lval;
  Capacitance[ Z_out ] = Cval;
}

Constraint {
  // Homogeneous Dirichlet condition on the vector potential at infinity:
  { Name a_Mag_2D;
    Case {
      { Region Bnd; Value 0; }
    }
  }

  // No imposed current or voltage on the FE coils (determined by the circuit):
  { Name Current_Mag_2D;
    Case {
    }
  }
  { Name Voltage_Mag_2D;
    Case {
    }
  }

  // No imposed current on circuit elements (the inductor initial condition is
  // only needed in the time domain, which is not used here):
  { Name Current_Cir ;
    Case {
    }
  }

  // Imposed voltage on the primary voltage source (as in tutorial 8):
  { Name Voltage_Cir ;
    Case {
      { Region E_in; Value Voltage; TimeFunction F_Cos_wt_p[]{2 * Pi * f, 0}; }
    }
  }

  // The circuit topology uses two independent circuits. "Circuit_1" and
  // "Circuit_2" define separate connected networks, each with its own set of
  // circuit nodes (the node numbers are local to each case). This allows
  // modeling magnetically coupled but electrically isolated windings -- the
  // coupling between the two circuits happens through the shared magnetic flux
  // in the finite element model, not through any electrical connection.
  { Name ElectricalCircuit ; Type Network ;
    Case Circuit_1 {
      // Primary circuit: voltage source driving the primary winding. The two
      // coil sides are connected in series (nodes 1 -> 2 -> 3 -> back to 1):
      { Region E_in; Branch {1, 2}; }
      { Region Coil_1_M; Branch {2, 3} ; }
      { Region Coil_1_P; Branch {3, 1} ; }
    }
    Case Circuit_2 {
      // Secondary circuit: the secondary winding drives the load. Note that the
      // secondary has no external source -- the voltage is induced by the
      // time-varying flux from the primary:
      { Region Coil_2_M; Branch {1, 2} ; }
      { Region Coil_2_P; Branch {2, 3} ; }
      { Region Z_out; Branch {3, 1}; }
    }
  }
}

// Configuration of the template library. The following variables are read by
// "Lib_Magnetoquasistatics_av_2D_Cir.pro" to select the analysis type and
// enable the appropriate features (the full list of configuration variables is
// documented at the beginning of the library file):
CircuitCoupling = 1;
AnalysisType = 2; // frequency-domain
CoefPower = 0.5; // peak-valued phasors
Freq = f;
Include "Lib_Magnetoquasistatics_av_2D_Cir.pro";

// The library is made generic through three GetDP declaration commands:
//  - "DefineConstant[]": defines a constant only if it has not already been
//    defined. The library uses it to set default values for all its
//    configuration variables (AnalysisType, Axisymmetric, Freq, etc.). If the
//    user sets a variable before the "Include", it takes precedence; otherwise
//    the library's default applies.
//  - "DefineGroup[]": similarly, declares a group only if it does not already
//    exist. The library declares all the abstract groups it knows about
//    (Vol_Mag, Vol_C_Mag, Vol_S_Mag, Vol_NL_Mag, Vol_M_Mag, etc.) with empty
//    defaults. The user populates only the groups relevant to the problem at
//    hand; unused groups remain empty and the corresponding terms in the
//    formulation are silently skipped (since an integral over an empty region
//    contributes nothing).
//  - "DefineFunction[]": declares a function only if it does not already exist.
//    The library declares all the material functions it may use (nu[], sigma[],
//    br[], js0[], etc.); the user defines only those needed for the problem.
//
// This "define-if-absent" mechanism allows a single library file to handle
// magnetostatics, eddy currents, permanent magnets, stranded coils, circuit
// coupling, nonlinear materials and impedance boundary conditions -- all
// controlled by which groups and functions the user chooses to populate before
// the "Include".

// The PostOperation is specific to this model: the library provides the
// PostProcessing (with quantities like "az", "b", "j", "U", "I", "NormU",
// "NormI", etc.), but the user decides what to print and where.
PostOperation {
  { Name Map; NameOfPostProcessing Mag;
    Operation {
      Print[ j, OnElementsOf Vol_Mag, Format Gmsh, File "j.pos" ];
      Print[ b, OnElementsOf Vol_Mag, Format Gmsh, File "b.pos" ];
      Print[ NormB, OnElementsOf Vol_Mag, Format Gmsh, File "normb.pos" ];
      Print[ az, OnElementsOf Vol_Mag, Format Gmsh, File "az.pos" ];
      Print[ h, OnElementsOf Core, Format Gmsh, File "h.pos" ];

      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"];
      Echo[ "Z_out", Format Table, File > "UI.txt" ];
      Print[ U, OnRegion Z_out, Format FrequencyTable, File > "UI.txt"];
      Print[ I, OnRegion Z_out, Format FrequencyTable, File > "UI.txt"];

      // The "{0}" suffix after the parameter name selects only the first
      // component to be sent to Gmsh (instead of both the real and imaginary
      // part, as all values are complex):
      Print[ NormU, OnRegion E_in, Format Table, File "",
        SendToServer "}Output/0|V_in| [V]"{0} ];
      Print[ NormI, OnRegion E_in, Format Table, File "",
        SendToServer "}Output/1|I_in| [A]"{0} ];
      Print[ NormU, OnRegion Z_out, Format Table, File "",
        SendToServer "}Output/2|V_out| [V]"{0} ];
      Print[ NormI, OnRegion Z_out, Format Table, File "",
        SendToServer "}Output/3|I_out| [A]"{0} ];

      // To interactively plot the exterior characteristic |V_out|(|I_out|) with
      // Gmsh, click the "Loop" arrow button next to the load value in the
      // "Parameters" menu (this sweeps the load over its allowed range), then
      // press "Run". Click the "Draw" graph button next to |I_out| and select
      // e.g. "Top left -> X", then click the "Draw" graph button next to
      // |V_out| and select "Top left -> Y" -- you can do this while the sweep
      // is still running to watch the curve build up live. The two quantities
      // are then plotted against each other in a graph anchored to the top-left
      // corner of the Gmsh window, with |I_out| on the X axis and |V_out| on
      // the Y axis.

      Print[ JouleLosses[Vol_Mag], OnGlobal, Format Table,
        File "losses.txt", SendToServer "}Output/p_Joule [W m^-3]"{0} ];
      Print[ HarmonicMagneticLosses[Vol_Mag], OnGlobal, Format Table,
        File > "losses.txt", SendToServer "}Output/p_H+EC [W m^-3]"{0} ];
    }
  }
}

// Lib_Magnetoquasistatics_av_2D_Cir.pro
//
// Template library for 2D magneto-quasistatic problems in terms of the magnetic
// vector potential "a" and the electric scalar potential "v", with optional
// circuit coupling.

// Default definitions of constants, groups and functions that can (should) be
// redefined from outside the template:

DefineConstant[
  AnalysisType = 1, // static (0), time-domain (1) or frequency-domain (2)
  Axisymmetric = 0, // axisymmetric model?
  CircuitCoupling = 0, // consider coupling with external electric circuit
  NewtonRaphson = 1, // Newton-Raphson or Picard method for nonlinear iterations
  CoefPower = 1, // power calculation coefficient (0.5 for peak-valued phasors)
  Freq = 50, // frequency (for harmonic simulations)
  TimeInit = 0, // initial time (for time-domain simulations)
  TimeFinal = 1/50, // final time (for time-domain simulations)
  DeltaTime = 1/500, // time step (for time-domain simulations)
  FEOrder = 1, // finite element order
  JacRadiusInt = 0, // interior radius for spherical shell (Vol_Inf_Mag)
  JacRadiusExt = 0, // exterior radius of spherical shell (Vol_Inf_Mag)
  JacCenterX = 0, // x-coordinate of center of Vol_Inf_Mag
  JacCenterY = 0, // y-coordinate of center of Vol_Inf_Mag
  JacCenterZ = 0, // z-coordinate of center of Vol_Inf_Mag
  NLTolAbs = 1e-10, // absolute tolerance on residual for nonlinear iterations
  NLTolRel = 1e-6, // relative tolerance on residual for nonlinear iterations
  NLIterMax = 20 // maximum number of nonlinear iterations
];

Group {
  DefineGroup[
    // The full magnetic domain:
    Vol_Mag,

    // Subsets of Vol_Mag:
    Vol_C_Mag, // massive conductors
    Vol_S_Mag, // stranded conductors
    Vol_S0_Mag, // stranded conductors with imposed current density js0
    Vol_NL_Mag, // nonlinear magnetic materials
    Vol_V_Mag, // moving massive conducting parts (with invariant mesh)
    Vol_M_Mag, // permanent magnets
    Vol_Inf_Mag, // annulus where an infinite shell transformation is applied

    // Boundaries:
    Sur_Neu_Mag, // non-homogeneous Neumann BC (flux tube with n x h = nxh[])
    Sur_Imped_Mag // conductors approximated by an impedance (non-meshed)
  ];
  If(CircuitCoupling)
    DefineGroup[
      Source_Cir, // voltage and current sources
      Resistance_Cir, // resistors
      Inductance_Cir, // inductors
      Capacitance_Cir // capacitors
    ];
  EndIf
}

Function {
  DefineFunction[
    nu, // reluctivity (in Vol_Mag)
    sigma, // conductivity (in Vol_C_Mag and Vol_S_Mag)
    br, // remanent magnetic flux density (in Vol_M_Mag)
    js0, // source current density (in Vol_S0_Mag)
    dhdb, // Jacobian for Newton-Raphson method (in Vol_NL_Mag)
    nxh, // n x magnetic field (on Sur_Neu_Mag)
    Velocity, // velocity of moving part (in Vol_V_Mag)
    Ns, // number of turns (in Vol_S_Mag)
    Sc, // cross-section of windings (in Vol_S_Mag)
    CoefGeo, // thickness (2D) or 2 * Pi (axisymmetric) geometrical coefficient
             // (in Vol_Mag); negative to set positive current along -z
    Ysur // surface admittance (on Sur_Imped_Mag)
  ];
  If(CircuitCoupling)
    DefineFunction[
      Resistance, // lumped resistance values
      Inductance, // lumped inductance values
      Capacitance // lumped capacitance values
    ];
  EndIf
}

// End of default definitions.

Jacobian {
  { Name JacVol_Mag;
    Case {
      If(Axisymmetric)
        { Region Vol_Inf_Mag;
          Jacobian VolAxiSquSphShell{JacRadiusInt, JacRadiusExt,
                                     JacCenterX, JacCenterY, JacCenterZ}; }
        { Region All; Jacobian VolAxiSqu; }
      Else
        { Region Vol_Inf_Mag;
          Jacobian VolSphShell{JacRadiusInt, JacRadiusExt,
                               JacCenterX, JacCenterY, JacCenterZ}; }
        { Region All; Jacobian Vol; }
      EndIf
    }
  }
  { Name JacSur_Mag;
    Case {
      If(Axisymmetric)
        { Region All; Jacobian SurAxi; }
      Else
        { Region All; Jacobian Sur; }
      EndIf
    }
  }
}

Integration {
  { Name Int_Mag;
    Case {
      { Type Gauss;
        Case {
          { GeoElement Point; NumberOfPoints  1; }
          { GeoElement Line; NumberOfPoints (FEOrder == 1) ? 2 : 4; }
          { GeoElement Triangle; NumberOfPoints (FEOrder == 1) ? 3 : 6; }
          { GeoElement Quadrangle; NumberOfPoints (FEOrder == 1) ? 4 : 16; }
        }
      }
    }
  }
}

Group{
  Dom_Hcurl_a_Mag_2D = Region[ {Vol_Mag, Sur_Neu_Mag, Sur_Imped_Mag} ];
  Dom_Hcurl_u_Mag_2D = Region[ {Vol_C_Mag, Sur_Imped_Mag} ];
  If(CircuitCoupling)
    Dom_Cir = Region[ {Source_Cir, Resistance_Cir, Inductance_Cir, Capacitance_Cir} ];
  EndIf
}

FunctionSpace {
  // Magnetic vector potential
  { Name Hcurl_a_Mag_2D; Type Form1P;
    BasisFunction {
      { Name se; NameOfCoef ae; Function BF_PerpendicularEdge;
        Support Dom_Hcurl_a_Mag_2D; Entity NodesOf[All]; }
      If(FEOrder == 2)
        { Name se2; NameOfCoef ae2; Function BF_PerpendicularEdge_2E;
          Support Dom_Hcurl_a_Mag_2D; Entity EdgesOf[All]; }
      EndIf
    }
    Constraint {
      { NameOfCoef ae; EntityType NodesOf; NameOfConstraint a_Mag_2D; }
      If(FEOrder == 2)
        { NameOfCoef ae2; EntityType EdgesOf; NameOfConstraint a0_Mag_2D; }
      EndIf
    }
  }

  // Gradient of scalar electric potential
  { Name Hcurl_u_Mag_2D; Type Form1P;
    BasisFunction {
      { Name sr; NameOfCoef ur; Function BF_RegionZ;
        Support Dom_Hcurl_u_Mag_2D; Entity Dom_Hcurl_u_Mag_2D; }
    }
    GlobalQuantity {
      { Name Voltage; Type AliasOf; NameOfCoef ur; }
      { Name Current; Type AssociatedWith; NameOfCoef ur; }
    }
    Constraint {
      { NameOfCoef Voltage; EntityType Region; NameOfConstraint Voltage_Mag_2D; }
      { NameOfCoef Current; EntityType Region; NameOfConstraint Current_Mag_2D; }
    }
  }

  // Current in stranded coils
  { Name Hregion_i_Mag_2D; Type Vector;
    BasisFunction {
      { Name sr; NameOfCoef ir; Function BF_RegionZ;
        Support Vol_S_Mag; Entity Vol_S_Mag; }
    }
    GlobalQuantity {
      { Name Current; Type AliasOf; NameOfCoef ir; }
      { Name Voltage; Type AssociatedWith; NameOfCoef ir; }
    }
    Constraint {
      { NameOfCoef Voltage; EntityType Region; NameOfConstraint Voltage_Mag_2D; }
      { NameOfCoef Current; EntityType Region; NameOfConstraint Current_Mag_2D; }
    }
  }

  If(CircuitCoupling)
    // Current in lumped circuit elements
    { Name Hregion_i_Cir; Type Scalar;
      BasisFunction {
        { Name sr; NameOfCoef ir; Function BF_Region;
          Support Dom_Cir; Entity Dom_Cir; }
      }
      GlobalQuantity {
        { Name Current; Type AliasOf; NameOfCoef ir; }
        { Name Voltage; Type AssociatedWith; NameOfCoef ir; }
      }
      Constraint {
        { NameOfCoef Voltage; EntityType Region; NameOfConstraint Voltage_Cir; }
        { NameOfCoef Current; EntityType Region; NameOfConstraint Current_Cir; }
      }
    }
  EndIf
}

Group{
  // all linear materials
  Vol_L_Mag = Region[ {Vol_Mag, -Vol_NL_Mag} ];
}

Formulation {
  { Name Magnetoquasistatics_av_2D; Type FemEquation;
    Quantity {
      { Name a; Type Local; NameOfSpace Hcurl_a_Mag_2D; }
      { Name ir; Type Local; NameOfSpace Hregion_i_Mag_2D; }
      { Name Us; Type Global; NameOfSpace Hregion_i_Mag_2D [Voltage]; }
      { Name Is; Type Global; NameOfSpace Hregion_i_Mag_2D [Current]; }
      If(AnalysisType > 0)
        { Name ur; Type Local; NameOfSpace Hcurl_u_Mag_2D; }
        { Name U; Type Global; NameOfSpace Hcurl_u_Mag_2D [Voltage]; }
        { Name I; Type Global; NameOfSpace Hcurl_u_Mag_2D [Current]; }
        If(CircuitCoupling)
          { Name Uz; Type Global; NameOfSpace Hregion_i_Cir [Voltage]; }
          { Name Iz; Type Global; NameOfSpace Hregion_i_Cir [Current]; }
        EndIf
      EndIf
    }
    Equation {
      Integral { [ nu[] * Dof{d a} , {d a} ];
        In Vol_L_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }

      If(NewtonRaphson)
        Integral { [ nu[{d a}] * {d a} , {d a} ];
          In Vol_NL_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
        Integral { [ dhdb[{d a}] * Dof{d a} , {d a} ];
          In Vol_NL_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
        Integral { [ - dhdb[{d a}] * {d a} , {d a} ];
          In Vol_NL_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
      Else
        Integral { [ nu[{d a}] * Dof{d a}, {d a} ];
          In Vol_NL_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
      EndIf

      Integral { [ - nu[] * br[] , {d a} ];
        In Vol_M_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }

      Integral { [ - js0[] , {a} ];
        In Vol_S0_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
      Integral { [ - Ns[] / Sc[] * Sign[CoefGeo[]] * Dof{ir} , {a} ];
        In Vol_S_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }

      Integral { [ Ns[]^2 / Sc[]^2 / sigma[] * Sign[CoefGeo[]] * Dof{ir} , {ir} ];
        In Vol_S_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
      GlobalTerm { [ Dof{Us} / CoefGeo[] , {Is} ]; In Vol_S_Mag; }

      Integral { [ nxh[] , {a} ];
        In Sur_Neu_Mag; Jacobian JacSur_Mag; Integration Int_Mag; }

      If(AnalysisType > 0)
        Integral { DtDof [ sigma[] * Dof{a} , {a} ];
          In Vol_C_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
        Integral { [ sigma[] * Dof{ur} / CoefGeo[] , {a} ];
          In Vol_C_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }

        Integral { [ - sigma[] * (Velocity[] /\ Dof{d a}) , {a} ];
          In Vol_V_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }

        Integral { DtDof [  Ysur[] * Dof{a} , {a} ];
          In Sur_Imped_Mag; Jacobian JacSur_Mag; Integration Int_Mag; }
        Integral { [ Ysur[] * Dof{ur} / CoefGeo[] , {a} ];
          In Sur_Imped_Mag; Jacobian JacSur_Mag; Integration Int_Mag; }

        Integral { DtDof [ sigma[] * Dof{a} , {ur} ];
          In Vol_C_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
        Integral { [ sigma[] * Dof{ur} / CoefGeo[] , {ur} ];
          In Vol_C_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
        GlobalTerm { [ Dof{I} *(CoefGeo[]/Fabs[CoefGeo[]]) , {U} ]; In Vol_C_Mag; }

        Integral { DtDof [ Ysur[] * Dof{a} , {ur} ];
          In Sur_Imped_Mag; Jacobian JacSur_Mag; Integration Int_Mag; }
        Integral { [ Ysur[] * Dof{ur} / CoefGeo[] , {ur} ];
          In Sur_Imped_Mag; Jacobian JacSur_Mag; Integration Int_Mag; }
        GlobalTerm { [ Dof{I} * (CoefGeo[] / Fabs[CoefGeo[]]) , {U} ]; In Sur_Imped_Mag; }

        Integral { DtDof [ Ns[] / Sc[] * Dof{a} , {ir} ];
          In Vol_S_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }

        If(CircuitCoupling)
          GlobalTerm { NeverDt[ Dof{Uz} , {Iz} ]; In Resistance_Cir; }
          GlobalTerm { NeverDt[ Resistance[{Uz}] * 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 { [ 0. * Dof{Iz} , {Iz} ]; In Source_Cir; }

          GlobalEquation {
            Type Network; NameOfConstraint ElectricalCircuit;
            { Node {I};  Loop {U};  Equation {I};  In Vol_C_Mag; }
            { Node {Is}; Loop {Us}; Equation {Us}; In Vol_S_Mag; }
            { Node {Iz}; Loop {Uz}; Equation {Uz}; In Dom_Cir; }
          }
        EndIf
      EndIf
    }
  }
}

Resolution {
  { Name Mag;
    System {
      { Name A; NameOfFormulation Magnetoquasistatics_av_2D;
        If(AnalysisType == 2)
          Type ComplexValue; Frequency Freq;
        EndIf
      }
    }
    Operation {
      If(AnalysisType == 0 || AnalysisType == 2)
        Generate[A]; Solve[A]; SaveSolution[A];
      Else
        InitSolution[A];
        TimeLoopTheta[TimeInit, TimeFinal, DeltaTime, 1.]{
          Generate[A]; Solve[A];
          If(NbrRegions[Vol_NL_Mag])
            Generate[A]; GetResidual[A, $res0];
            Evaluate[ $res = $res0, $iter = 0 ];
            Print[{$iter, $res, $res / $res0},
              Format "Residual %03g: abs %14.12e rel %14.12e"];
            While[$res > NLTolAbs && $res / $res0 > NLTolRel &&
                  $res / $res0 <= 1 && $iter < NLIterMax]{
              Solve[A]; Generate[A]; GetResidual[A, $res];
              Evaluate[ $iter = $iter + 1 ];
              Print[{$iter, $res, $res / $res0},
                Format "Residual %03g: abs %14.12e rel %14.12e"];
            }
          EndIf
          SaveSolution[A];
        }
      EndIf
    }
  }
}

PostProcessing {
  { Name Mag; NameOfFormulation Magnetoquasistatics_av_2D;
    Quantity {
      { Name a; Value {
          Term { [ {a} ]; In Vol_Mag; Jacobian JacVol_Mag; }
        }
      }
      { Name az; Value {
          Term { [ CompZ[{a}] ]; In Vol_Mag; Jacobian JacVol_Mag; }
        }
      }
      { Name xaz; Value {
          Term { [ X[] * CompZ[{a}] ]; In Vol_Mag; Jacobian JacVol_Mag; }
        }
      }
      { Name b; Value {
          Term { [ {d a} ]; In Vol_Mag; Jacobian JacVol_Mag; }
        }
      }
      { Name NormB; Value {
          Term { [ Norm[{d a}] ]; In Vol_Mag; Jacobian JacVol_Mag; }
        }
      }
      { Name h; Value {
          Term { [ nu[] * {d a} ]; In Vol_L_Mag; Jacobian JacVol_Mag; }
          Term { [ nu[{d a}] * {d a} ]; In Vol_NL_Mag; Jacobian JacVol_Mag; }
          Term { [ -nu[] * br[] ]; In Vol_M_Mag; Jacobian JacVol_Mag; }
        }
      }
      { Name j; Value {
          If(AnalysisType > 0)
            Term { [ -sigma[] * (Dt[{a}] + {ur} / CoefGeo[]) ];
              In Vol_C_Mag; Jacobian JacVol_Mag; }
            // Current density in A/m
            Term { [ -Ysur[] * (Dt[{a}] + {ur} / CoefGeo[]) ];
              In Sur_Imped_Mag; Jacobian JacSur_Mag; }
          EndIf
          Term { [ js0[] ];
            In Vol_S0_Mag; Jacobian JacVol_Mag; }
          Term { [ Ns[] / Sc[] * Sign[CoefGeo[]] * {ir} ];
            In Vol_S_Mag; Jacobian JacVol_Mag; }
          Term { [ Vector[0, 0, 0] ];
            In Vol_Mag; Jacobian JacVol_Mag; }
        }
      }
      { Name U; Value {
          If(AnalysisType > 0)
            Term { [ {U} ]; In Vol_C_Mag; }
            If(CircuitCoupling)
              Term { [ {Uz} ]; In Dom_Cir; }
            EndIf
          EndIf
          Term { [ {Us} ]; In Vol_S_Mag; }
        }
      }
      { Name I; Value {
          If(AnalysisType > 0)
            Term { [ {I} ]; In Vol_C_Mag; }
            If(CircuitCoupling)
              Term { [ {Iz} ]; In Dom_Cir; }
            EndIf
          EndIf
          Term { [ {Is} ]; In Vol_S_Mag; }
        }
      }
      { Name NormU; Value {
          If(AnalysisType > 0)
            Term { [ Norm[{U}] ]; In Vol_C_Mag; }
            If(CircuitCoupling)
              Term { [ Norm[{Uz}] ]; In Dom_Cir; }
            EndIf
          EndIf
          Term { [ Norm[{Us}] ]; In Vol_S_Mag; }
        }
      }
      { Name NormI; Value {
          If(AnalysisType > 0)
            Term { [ Norm[{I}] ]; In Vol_C_Mag; }
            If(CircuitCoupling)
              Term { [ Norm[{Iz}] ]; In Dom_Cir; }
            EndIf
          EndIf
          Term { [ Norm[{Is}] ]; In Vol_S_Mag; }
        }
      }
      { Name JouleLosses; Value {
          If(AnalysisType > 0)
            Integral { [ CoefPower * sigma[] * SquNorm[Dt[{a}] + {ur} / CoefGeo[]] ];
              In Vol_C_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
            Integral { [ CoefPower * Ysur[] * SquNorm[Dt[{a}] + {ur} / CoefGeo[]] ];
              In Sur_Imped_Mag; Jacobian JacSur_Mag; Integration Int_Mag; }
          EndIf
          Integral { [ CoefPower * 1. / sigma[] * SquNorm[js0[]] ];
            In Vol_S0_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
          Integral { [ CoefPower * 1. / sigma[] *
              SquNorm[Ns[] / Sc[] * Sign[CoefGeo[]] * {ir}] ];
            In Vol_S_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
        }
      }
      { Name HarmonicMagneticLosses; Value {
          If(AnalysisType == 2)
            Integral { [ CoefPower * 2 * Pi * Freq * Im[nu[]] * SquNorm[{d a}] ];
              In Vol_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
          EndIf
	}
      }
      { Name LaplaceForce; Value {
          If(AnalysisType == 2)
            // DC component (phasor at 2 * Freq in HarmonicLaplaceForce2F)
            Integral { [
                Re[-sigma[] * (Dt[{a}] + {ur} / CoefGeo[])] /\ Re[{d a}] / 2. +
                Im[-sigma[] * (Dt[{a}] + {ur} / CoefGeo[])] /\ Im[{d a}] / 2. ];
              In Vol_C_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
            Integral { [
                Re[js0[]] /\ Re[{d a}] / 2. +
                Im[js0[]] /\ Im[{d a}] / 2. ];
              In Vol_S0_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
            Integral { [
                Re[(js0[] *Vector[0, 0, 1]) * {ir}] /\ Re[{d a}] / 2. +
                Im[(js0[] *Vector[0, 0, 1]) * {ir}] /\ Im[{d a}] / 2. ];
              In Vol_S_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
          Else
            If(AnalysisType == 1)
              Integral { [ -sigma[] * (Dt[{a}] + {ur} / CoefGeo[]) /\ {d a} ];
                In Vol_C_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
            EndIf
            Integral { [ js0[] /\ {d a} ];
              In Vol_S0_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
            Integral { [ Ns[] / Sc[] * Sign[CoefGeo[]] * {ir} /\ {d a} ];
              In Vol_S_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
          EndIf
        }
      }
      { Name HarmonicLaplaceForce2F; Value {
          If(AnalysisType == 2)
            // 2 * Freq component
            Integral { [ Complex[
                  Re[-sigma[] * (Dt[{a}] + {ur} / CoefGeo[])] /\ Re[{d a}] / 2. -
                  Im[-sigma[] * (Dt[{a}] + {ur} / CoefGeo[])] /\ Im[{d a}] / 2. ,
                  Im[-sigma[] * (Dt[{a}] + {ur} / CoefGeo[])] /\ Re[{d a}] / 2. +
                  Re[-sigma[] * (Dt[{a}] + {ur} / CoefGeo[])] /\ Im[{d a}] / 2. ] ];
              In Vol_C_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
            Integral { [ Complex [
                  Re[js0[]] /\ Re[{d a}] / 2. -
                  Im[js0[]] /\ Im[{d a}] / 2. ,
                  Im[js0[]] /\ Re[{d a}] / 2. +
                  Re[js0[]] /\ Im[{d a}] / 2. ] ];
              In Vol_S0_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
            Integral { [ Complex [
                  Re[Ns[] / Sc[] * Sign[CoefGeo[]] * {ir}] /\ Re[{d a}] / 2. -
                  Im[Ns[] / Sc[] * Sign[CoefGeo[]] * {ir}] /\ Im[{d a}] / 2. ,
                  Im[Ns[] / Sc[] * Sign[CoefGeo[]] * {ir}] /\ Re[{d a}] / 2. +
                  Re[Ns[] / Sc[] * Sign[CoefGeo[]] * {ir}] /\ Im[{d a}] / 2. ] ];
              In Vol_S_Mag; Jacobian JacVol_Mag; Integration Int_Mag; }
          EndIf
        }
      }
    }
  }
}

// Parameters shared by Gmsh and GetDP.

mm = 1e-3;

dxCore = 50 * mm;
dyCore = 100 * mm;
xCoil = 75 * mm;
dxCoil = 25 * mm;
dyCoil = 100 * mm;
rInt = 200 * mm;
rExt = 250 * mm;

// Gmsh script describing the geometry of the electromagnet. This is a pared
// down version of the geometry from tutorial 3, where all symmetries are
// taken advantage of.

Include "electromagnet_common.pro";

SetFactory("OpenCASCADE");
Rectangle(1) = {-dxCore, -dyCore, 0, 2 * dxCore, 2 * dyCore};
Rectangle(2) = {xCoil, -dyCoil / 2, 0, dxCoil, dyCoil};
Rectangle(3) = {-xCoil - dxCoil, -dyCoil / 2, 0, dxCoil, dyCoil};
Disk(4) = {0, 0, 0, rInt};
Disk(5) = {0, 0, 0, rExt};
d = 1.1 * rExt;
Rectangle(6) = {0, 0, 0, d, d};
BooleanIntersection{ Surface{6}; Delete; }{ Surface{1:5}; Delete; }

core() = Closest{0, 0, 0}{ Surface{:}; };
air() = Closest{dxCore + mm, 0, 0}{ Surface{:}; };
airinf() = Closest{rInt + mm, 0, 0}{ Surface{:}; };
indr() = Closest{xCoil + mm, 0, 0}{ Surface{:}; };
Physical Surface("Core", 1) = core(0);
Physical Surface("Air", 2) = air(0);
Physical Surface("AirInf", 3) = airinf(0);
Physical Surface("CoilRight", 4) = indr(0);

bot() = Curve In BoundingBox{-mm, -mm, -mm, d, mm, mm};
left() = Curve In BoundingBox{-mm, -mm, -mm, mm, d, mm};
inf() = CombinedBoundary{ Surface{:}; };
inf() -= {bot(), left()};
Physical Curve("Bottom", 10) = bot();
Physical Curve("Left", 11) = left();
Physical Curve("Inf", 12) = inf();

// Mesh size constraints:
DefineConstant[
  s = {1, Name "Parameters/}Global mesh size factor"}
];

MeshSize{:} = 12.5 * mm * s;
MeshSize{ PointsOf{ Surface{indr(0)}; } } = 5 * mm * s;
MeshSize{ PointsOf{ Surface{core(0)}; } } = 4 * mm * s;

// Refine to capture the skin depth:
skin() = Curve In BoundingBox{dxCore - mm, -mm, -mm, dxCore + mm, d, mm};
MeshSize{ PointsOf{ Curve{skin(0)}; } } = 0.5 * mm * s;

// This file reimplements the electromagnet models from tutorials 3 and 4 using
// the "Lib_Magnetoquasistatics_av_2D_Cir.pro" template library. The library
// acts as a black box that hides the complexity of the formulation, function
// spaces, resolution and post-processing: the user only needs to define groups,
// material properties and constraints. The resulting model file is compact, yet
// supports static, time-domain and frequency-domain analyses, with or without
// axisymmetry.

Include "electromagnet_common.pro";

Group {
  // Physical regions:
  Core = Region[ 1 ];
  Air = Region[ 2 ];
  AirInf = Region[ 3 ];
  CoilRight = Region[ 4 ];
  Bottom = Region[ 10 ];
  Left = Region[ 11 ];
  Inf = Region[ 12 ];

  // Abstract regions expected by the library. The inductor is a stranded coil
  // ("Vol_S_Mag"), the core is a massive conductor with eddy currents
  // ("Vol_C_Mag"):
  Vol_Mag     = Region[ {Air, AirInf, Core, CoilRight} ];
  Vol_S_Mag   = Region[ CoilRight ];
  Vol_C_Mag   = Region[ Core ];
  Vol_Inf_Mag = Region[ AirInf ];
}

Function {
  DefineConstant[
    AnalysisType = {0, Choices{0="Static", 1="Time-domain", 2="Frequency-domain"},
      Name "Parameters/Analysis type"}
    Axisymmetric = {0, Choices{0, 1},
      Name "Parameters/Axisymmetric model?"}
    Current = {1, Min 0.01, Max 100, Step 0.1,
      Name "Parameters/Current [A]"}
    f = {5, Min 1, Max 1000, Step 1, Visible AnalysisType > 0,
      Name "Parameters/Frequency [Hz]"}
    mur = {1000, Min 1, Max 5000, Step 1,
      Name "Parameters/Core relative permeability"}
    sigma = {1e6, Min 1, Max 1e6, Step 100,
      Name "Parameters/Core conductivity [S m^-1]"}
  ];

  T = 1 / f;
  dt = T / 20;
  tmax = 2 * T;

  mu0 = 4.e-7 * Pi;
  nu0 = 1 / mu0;
  nu [ Region[{Air, AirInf, CoilRight}] ] = nu0;
  nu [ Core ]  = 1 / (mur * mu0);

  sigma[ Core ] = sigma;
  sigma[ CoilRight ] = sigma;

  If(AnalysisType > 0)
    time_fct[] = F_Sin_wt_p[]{2 * Pi * f, 0};
    SkinDepth = DefineNumber[ 1e3 * Sqrt[2 / (2 * Pi * f * mur * mu0 * sigma)],
      Name "Parameters/Skin depth [mm]", ReadOnly];
  Else
    time_fct[] = 1;
  EndIf

  // Stranded coil parameters: these four definitions are all the library needs
  // to handle the 1000-turn inductor with imposed current:
  Ns[ CoilRight ] = 1000;
  Sc[ CoilRight ] = dxCoil * dyCoil;
  CoefGeo[ CoilRight ] = (Axisymmetric ? 2 * Pi : 1) * (-1);
  CoefGeo[ Core ] = (Axisymmetric ? 2 * Pi : 1);
}

Constraint {
  { Name a_Mag_2D;
    Case {
      { Region Left; Value 0.; }
      { Region Inf; Value 0.; }
    }
  }
  { Name Voltage_Mag_2D;
    Case {
      // Zero voltage on the core
      { Region Core; Value 0; }
    }
  }
  { Name Current_Mag_2D;
    Case {
      // Imposed current in the inductor:
      { Region CoilRight; Value Current; TimeFunction time_fct[]; }
    }
  }
}

// Configure and include the library:
JacRadiusInt = rInt;
JacRadiusExt = rExt;
Freq = f;
TimeFinal = tmax;
DeltaTime = dt;
CoefPower = (AnalysisType == 2) ? 0.5 : 1;
Include "../Lib_Magnetoquasistatics_av_2D_Cir.pro";

PostOperation {
  { Name Map; NameOfPostProcessing Mag;
    Operation {
      Print[ a, OnElementsOf Vol_Mag, File "a.pos" ];
      Print[ az, OnElementsOf Vol_Mag, File "az.pos" ];
      Print[ b, OnElementsOf Vol_Mag, File "b.pos" ];
      Print[ b, OnLine{{mm, mm, 0}{rInt, mm, 0}}{50}, File "cutb.pos" ];
      Print[ j, OnElementsOf Vol_S_Mag, File "js.pos" ];
      Print[ j, OnElementsOf Vol_C_Mag, File "j.pos" ];
      Print[ JouleLosses[Vol_C_Mag], OnGlobal, Format Table, File "losses.txt" ];
    }
  }
}

// Parameters shared by Gmsh and GetDP: core and coil dimensions.

mm = 1e-3;
wCore = 100 * mm;
hCore = 100 * mm;
wCoil1 = 10 * mm;
hCoil1 = 40 * mm;
wCoil2 = 20 * mm;
hCoil2 = 40 * mm;

DefineConstant[
  thickCore = {30, Min 5, Max 60, Step 0.1,
    Name "Parameters/Core thickness [mm]"}
  gapCoil = {2, Min 0.5, Max 7.5, Step 0.1,
    Name "Parameters/Coil gap [mm]"}
  wCoreLeg = {20, Min 10, Max 27, Step 0.1,
    Name "Parameters/Core width [mm]"}
];

thickCore = thickCore * mm;
gapCoil = gapCoil * mm;
wCoreLeg = wCoreLeg * mm;

// Inner and outer radii for the infinite shell Jacobian (see tutorial 3),
// chosen to be large enough to enclose the full geometry:
r = Max[(wCore + 2 * gapCoil + wCoil1 + wCoil2) / 2, hCore / 2];
r = Max[r, (thickCore + 2 * gapCoil + Max[wCoil1, wCoil2]) / 2];
JacRadiusInt = r * Sqrt[2] + 30 * mm;
JacRadiusExt = JacRadiusInt + 50 * mm;

// Gmsh script describing the geometry of the same single-phase transformer as
// in tutorial 9, but in 3D.

SetFactory("OpenCASCADE");

Include "transfo_common.pro";

// The core is created by subtracting the inner box (the window) from the
// outer box:
Box(1) = {-wCore / 2, -hCore / 2, -thickCore / 2,
  wCore, hCore, thickCore};
Box(2) = {-wCore / 2 + wCoreLeg, -hCore / 2 + wCoreLeg, -thickCore / 2,
  wCore - 2 * wCoreLeg, hCore - 2 * wCoreLeg, thickCore};
BooleanDifference(3) = { Volume{1}; Delete; }{ Volume{2}; Delete; };

// We define a macro to create a coil, which expects three input parameters to
// be defined:
//  - "xCoil": leftmost position of the coil
//  - "wCoil": width of the coil
//  - "hCoil": height of the coil
Macro Coil
  r1 = news; // "news" returns the next available surface tag
  Rectangle(r1) = {xCoil, -hCoil / 2,
    thickCore / 2 + gapCoil + wCoil,
    2 * wCoil + 2 * gapCoil + wCoreLeg,
    thickCore + 2 * gapCoil + 2 * wCoil,
    wCoil / 2}; // rounded corner with radius "wCoil / 2"
  r2 = news;
  Rectangle(r2) = {xCoil + wCoil, -hCoil / 2 + wCoil,
    thickCore / 2 + gapCoil + wCoil,
    wCoreLeg + 2 * gapCoil,
    thickCore + 2 * gapCoil};
  // Subtract the inner rectangle from the outer rounded rectangle:
  d() = BooleanDifference{ Surface{r1}; Delete; }{ Surface{r2}; Delete; };
  // Rotate the resulting surface to create the base of the coil:
  Rotate {{1, 0, 0}, {0, -hCoil / 2, thickCore / 2 + gapCoil + wCoil}, -Pi / 2} {
    Surface{d(0)};
  }
  // Extrude the base along the y-axis:
  e() = Extrude{0, hCoil, 0}{ Surface{d(0)}; };
  // Create a "cut" in the coil to define an electrode. The cut surface
  // is needed in 3D so that the source current density inside the coil
  // can be derived from a single-valued scalar potential, with a unit
  // jump across the cut (see the .pro file).
  r3 = news;
  Rectangle(r3) = {xCoil, -hCoil / 2, 0, wCoil, hCoil};
  BooleanFragments { Surface{r3}; Delete; }{ Volume{e(1)}; Delete; }
Return

// Call the "Coil" macro twice to create the two coils:
xCoil = -wCore / 2 - gapCoil - wCoil1;
wCoil = wCoil1;
hCoil = hCoil1;
Call Coil;
xCoil = wCore / 2 - wCoreLeg - gapCoil - wCoil2;
wCoil = wCoil2;
hCoil = hCoil2;
Call Coil;

// Create the spherical shell region:
Sphere(10) = {0, 0, 0, JacRadiusInt};
Sphere(11) = {0, 0, 0, JacRadiusExt};

// Fragment all entities to get a conformal geometry:
Coherence;

// Use the "Closest" command to retrieve the regions of interest:
coil1() = Closest {-wCore / 2 - gapCoil - wCoil1 / 2, 0, 0} { Volume{:}; };
elec1() = Closest {-wCore / 2 - gapCoil - wCoil1 / 2, 0, 0} { Surface{:}; };
coil2() = Closest {wCore / 2 + gapCoil + wCoil2 / 2, 0, 0} { Volume{:}; };
elec2() = Closest {wCore / 2 - wCoreLeg - gapCoil - wCoil2 / 2, 0, 0} { Surface{:}; };
core() = Closest {0, hCore / 2 - mm, 0} { Volume{:}; };
airinf() = Closest {JacRadiusInt + mm, 0, 0} { Volume{:}; };

air() = Volume{:};
air() -= {airinf(0), core(0), coil1(0), coil2(0)};
bnd() = CombinedBoundary{ Volume{:}; };

// Assign mesh sizes:
DefineConstant[
  s = {3, Min 0.1, Max 10, Step 0.1,
    Name "Parameters/}Global mesh size factor"}
];
MeshSize{:} = 2 * mm * s;
MeshSize{ PointsOf{ Volume{airinf(0)}; } } = 8 * mm * s;

// Create physical groups:
Physical Volume("Air", 1) = air();
Physical Volume("AirInf", 2) = airinf(0);
Physical Volume("Coil1", 3) = coil1(0);
Physical Volume("Coil2", 4) = coil2(0);
Physical Volume("Core", 5) = core(0);
Physical Surface("Bnd", 10) = bnd();
Physical Surface("Electrode1", 11) = elec1(0);
Physical Surface("Electrode2", 12) = elec2(0);

// GetDP formulation for same transformer model as in tutorial 9, but in 3D. The
// 3D vector potential "a" requires gauging for uniqueness in 3D; this tutorial
// presents two strategies, selectable via the "GaugeType" parameter:
// tree-cotree gauging and Coulomb gauging (through a Lagrange multiplier).

// Resolution sequence: we first pre-compute a normalized source current density
// in each coil, then link it to global variables for field-circuit coupling in
// the 3D a-v formulation. The actual source current density is enforced in each
// coil by multiplying the pre-computed distribution by the global current
// quantity computed from the circuit equations.
Include "transfo_common.pro";

Group {
  // Physical regions:
  Air = Region[ 1 ];
  AirInf = Region[ 2 ];
  Coil_1 = Region[ 3 ];
  Coil_2 = Region[ 4 ];
  Coils = Region[ {Coil_1, Coil_2} ];
  Core = Region[ 5 ];
  Bnd = Region[ 10 ];
  Electrode_1 = Region[ 11 ];
  Electrode_2 = Region[ 12 ];

  // Abstract regions:
  //  - "Vol_Mag": full magnetic domain
  //  - "Vol_S_Mag": stranded conductors (windings with Ns turns)
  //  - "Vol_Inf_Mag": region with infinite shell transformation
  //  - "Vol_Ele~{i}": i-th coil for source current density pre-calculation
  //  - "Sur_Electrode_Ele~{i}": electrode in i-th coil
  Vol_Mag = Region[ {Air, AirInf, Core, Coil_1, Coil_2} ];
  Vol_S_Mag = Region[ {Coil_1, Coil_2} ];
  Vol_Inf_Mag = Region[ AirInf ];
  NumCoils = 2;
  For i In {1 : NumCoils}
    // "~{i}" in the group names below will be transformed into "_1" when "i ==
    // 1" and "_2" when "i == 2" in the "For" loop. This is a general mechanism
    // in GetDP that can be applied not only to groups but to most object names
    // (constants, functions, constraints, function spaces, ...) in order to
    // create their names programmatically:
    Vol_Ele~{i} = Region[ Coil~{i} ];
    Sur_Electrode_Ele~{i} = Region[ Electrode~{i} ];
  EndFor

  // Circuit regions:
  E_in = Region[ 100 ];
  Z_out = Region[ 102 ];

  // Abstract circuit regions:
  //  - "Source_Cir": all voltage and current sources
  //  - "Resistance_Cir": all resistors
  //  - "Inductance_Cir": all inductors
  //  - "Capacitance_Cir": all capacitors
  Source_Cir = Region[ {E_in} ];
  DefineConstant[
    LoadType = {0, Choices {0="Resistance", 1="Inductance", 2="Capacitance"},
      Name "Parameters/2Load type"}
  ];
  If(LoadType == 0)
    Resistance_Cir = Region[ {Z_out} ];
  ElseIf(LoadType == 1)
    Inductance_Cir = Region[ {Z_out} ];
  Else
    Capacitance_Cir = Region[ {Z_out} ];
  EndIf
}

Function {
  DefineConstant[
    Voltage = {1, Min 1e-3, Max 10,
      Name "Parameters/0Primary voltage [V]"}
    f = {50, Min 0, Max 1e3, Step 1,
      Name "Parameters/1Frequency [Hz]"}
    Rval = {50, Min 1e-3, Max 50, Step 0.25, Visible (LoadType == 0),
      Name "Parameters/3Load resistance [Ohm]"}
    Lval = {1e-6, Min 1e-9, Max 0.1, Step 1e-3, Visible (LoadType == 1),
      Name "Parameters/3Load inductance [H]"}
    Cval = {1e-6, Min 1e-9, Max 0.1, Step 1e-4, Visible (LoadType == 2),
      Name "Parameters/3Load capacitance [F]"}
    murCore = {2000, Min 1, Max 10000, Step 1,
      Name "Parameters/4Core relative permeability"}
    dlam = {0.35e-3, Min 0.1e-3, Max 1e-3, Step 0.05e-3,
      Name "Parameters/5Lamination thickness [m]"}
    kh = {100, Min 10, Max 500, Step 10,
      Name "Parameters/6Hysteresis coeff [W m^-3 T^-2 s]"}
    GaugeType = {0, Choices {0="Tree-Cotree", 1="Coulomb"},
      Name "Parameters/7Gauge type"}
  ];

  mu0 = 4e-7 * Pi;
  nu[ Air ] = 1 / mu0;
  nu[ AirInf ] = 1 / mu0;

  // Same modelling of ferromagnetic laminations as in tutorial 9:
  nu_r = 1 / (murCore * mu0);
  sigma_iron = 2e6;
  nu_i = kh / Pi + sigma_iron * (2 * Pi * f) * dlam^2 / 12;
  nu[ Core ] = Complex[nu_r, nu_i];

  // Number of turns "Ns" of each coil:
  Ns[ Coil_1 ] = 100;
  Ns[ Coil_2 ] = 200;

  nu[ Coils ] = 1 / mu0;
  sigma[ Coils ] = 5.8e7;

  Resistance[ Z_out ] = Rval;
  Inductance[ Z_out ] = Lval;
  Capacitance[ Z_out ] = Cval;
}

Constraint {
  { Name a_Mag;
    Case {
      { Region Bnd; Value 0; }
    }
  }

  // Tree-cotree gauging constraint on the vector potential "a" (only used when
  // "GaugeType == 0"); "a" is set to zero on the edges of a spanning tree of
  // the mesh; the "SubRegion Bnd" anchors the tree on the outer boundary, where
  // the homogeneous Dirichlet boundary condition "a_Mag" is imposed:
  { Name a_Gauge_Mag;
    Case {
      { Region Vol_Mag ; SubRegion Bnd; Value 0.; }
    }
  }

  // Boundary condition for the Coulomb gauge multiplier "xi" (only used when
  // "GaugeType == 1"):
  { Name xi_Mag;
    Case {
      // N.B.: we are only interested in the gradient of "xi", so we can set "xi"
      // to any value on the boundary; here we choose zero as it allows for an
      // easier mixed formulation (see explanations in the "Formulation" below):
      { Region Bnd; Value 0; }
    }
  }

  // No imposed current or voltage in the magnetic formulation -- the voltage
  // source is imposed in the circuit as in tutorial 9:
  { Name Current_Mag;
    Case {
    }
  }
  { Name Voltage_Mag;
    Case {
    }
  }

  { Name Current_Cir ;
    Case {
    }
  }

  { Name Voltage_Cir ;
    Case {
      { Region E_in; Value Voltage; TimeFunction F_Cos_wt_p[]{2 * Pi * f, 0}; }
    }
  }

  // Contrary to the massive conductor case (see tutorial 8), the global basis
  // function for circuit coupling is not associated with the electrodes
  // anymore, but directly with the regions in "Vol_S_Mag", i.e., with "Coil_1"
  // and "Coil_2", through the use of the pre-computed source current density
  // with "BF_Global" in the "Hregion_vs_Mag" function space (see below):
  { Name ElectricalCircuit ; Type Network ;
    Case Circuit_1 {
      { Region E_in; Branch {1, 2}; }
      { Region Coil_1; Branch {2, 1} ; }
    }
    Case Circuit_2 {
      { Region Coil_2; Branch {1, 2} ; }
      { Region Z_out; Branch {2, 1}; }
    }
  }
}

Jacobian {
  { Name Vol;
    Case {
      { Region All; Jacobian Vol; }
    }
  }
  { Name Sur;
    Case {
      { Region All; Jacobian Sur; }
    }
  }
}

Integration {
  { Name Int;
    Case {
      { Type Gauss;
        Case {
          { GeoElement Line; NumberOfPoints 4; }
          { GeoElement Triangle; NumberOfPoints 4; }
          { GeoElement Tetrahedron; NumberOfPoints 4; }
        }
      }
    }
  }
}

Group {
  DefineGroup[ Source_Cir, Resistance_Cir, Inductance_Cir, Capacitance_Cir ];
  Dom_Cir = Region[ {Source_Cir, Resistance_Cir, Inductance_Cir,
      Capacitance_Cir} ];
}

// Pre-computation of the source current density distribution in each coil. The
// scalar potential "vs" solves a homogeneous Laplace problem on the coil
// volume, with a unit jump across the electrode (the cut introduced in the
// .geo). The resulting field "js = -grad(vs)" is therefore divergence-free and,
// by construction, integrates to one across any cross-section. Below, each "js"
// field will be multiplied by the net current in Coil~{i} to model the actual
// current density distribution in the coil.
//
// Note that "js" is not exactly uniform across the cross-section: like a DC
// current in a homogeneous solid conductor of the same shape, it concentrates
// on the inner side of bends, where the current path is shorter. For a stranded
// winding (whose ideal current density would be uniform along the wire axis,
// "js = tau / Sc" with "tau" the unit tangent), this is therefore an
// approximation. Alternatives include prescribing "js" analytically (only
// practical for simple shapes like straight or circular coils), computing it
// from a magnetostatic problem with a known driving current (e.g. via
// Biot-Savart for line segments), or homogenizing the winding. The Laplace
// approach used here is fully general and divergence-free by construction.

For i In {1 : NumCoils}

  Constraint {
    // We enforce unit current across each coil's electrode
    { Name UnitCurrent~{i};
      Case {
        { Region Sur_Electrode_Ele~{i}; Type Assign; Value 1; }
      }
    }
  }

  FunctionSpace {
    { Name H1_vs~{i}; Type Form0;
      BasisFunction {
        { Name sn; NameOfCoef vn; Function BF_Node;
          Support Vol_Ele~{i};
          Entity NodesOf[ All, Not Sur_Electrode_Ele~{i} ]; }
        // We restrict the support of the global basis function to the positive
        // side of the electrode, i.e. the side where the electrode normal
        // points, using "ElementsOf[ ..., OnPositiveSideOf ...]". Since the
        // electrode surface "Sur_Electrode_Ele~{i}" is connected on both sides
        // to the same coil volume "Coil~{i}", this is necessary to properly set
        // up the discontinuity of "vs": the negative side of the electrode in
        // "Coil~{i}" will not be part of the support, hence the global function
        // will evaluate to zero there.
        //
        // As an alternative, we could have cut each "Coil~{i}" volume into two
        // parts, and restricted the support of the global function to one of
        // the parts.
        { Name sf; NameOfCoef vf; Function BF_GroupOfNodes;
          Support ElementsOf[Vol_Ele~{i}, OnPositiveSideOf Sur_Electrode_Ele~{i}];
          Entity GroupsOfNodesOf[ Sur_Electrode_Ele~{i} ]; }
      }
      GlobalQuantity {
        { Name Voltage; Type AliasOf; NameOfCoef vf; }
        { Name Current; Type AssociatedWith; NameOfCoef vf; }
      }
      Constraint {
        { NameOfCoef Current; EntityType GroupsOfNodesOf;
          NameOfConstraint UnitCurrent~{i}; }
      }
    }
  }

  Formulation {
    { Name Electrokinetics_vs~{i}; Type FemEquation;
      Quantity {
        { Name vs; Type Local; NameOfSpace H1_vs~{i}; }
        { Name I; Type Global; NameOfSpace H1_vs~{i} [Current]; }
        { Name V; Type Global; NameOfSpace H1_vs~{i} [Voltage]; }
      }
      Equation {
        Integral { [ Dof{d vs} , {d vs} ];
          In Vol_Ele~{i}; Jacobian Vol; Integration Int; }
        GlobalTerm { [ Dof{I} , {V} ]; In Sur_Electrode_Ele~{i}; }
      }
    }
  }

  Resolution {
    // The "Ele~{i}" resolutions will be triggered automatically in a
    // pre-resolution phase, due to the "BF_Global" definition below; we hide
    // these resolutions from the list of available resolutions with the "Hidden
    // 1" option:
    { Name Ele~{i}; Hidden 1;
      System {
        { Name Sys_Ele; NameOfFormulation Electrokinetics_vs~{i}; }
      }
      Operation {
        Generate[Sys_Ele]; Solve[Sys_Ele];
      }
    }
  }

EndFor

FunctionSpace {

  // Function space for the magnetic vector potential "a". In 3D, the curl-curl
  // operator has a non-trivial kernel: since "curl(grad(phi)) = 0", any "a +
  // grad(phi)" gives the same flux density "b = curl(a)" and the same curl-curl
  // energy, so "a" is indeterminate up to gradients of arbitrary scalars "phi
  // in H1_0". A gauge condition is needed to obtain a unique solution:
  //  - Tree-cotree gauging ("GaugeType == 0"): pin "ae" to zero on the edges of
  //    a spanning tree of the mesh, rooted on the outer boundary
  //    ("EntitySubType StartingOn"). The gradient null space is spanned exactly
  //    by the tree edges (modulo boundary), so this kills it without affecting
  //    the physical solution. This constraint is directly implemented in the
  //    "Hcurl_a_Mag" space.
  //  - Coulomb gauging (GaugeType == 1): the gauge is imposed weakly in the
  //    formulation (see explanations in the "Formulation" below), thanks to a
  //    Lagrange multiplier "xi" in the "H1_xi_Mag" space.

  { Name Hcurl_a_Mag; Type Form1;
    BasisFunction {
      { Name se; NameOfCoef ae; Function BF_Edge;
        Support Vol_Mag; Entity EdgesOf[All]; }
    }
    Constraint {
      { NameOfCoef ae; EntityType EdgesOf; NameOfConstraint a_Mag; }
      If (GaugeType == 0)
        { NameOfCoef ae; EntityType EdgesOfTreeIn; EntitySubType StartingOn;
          NameOfConstraint a_Gauge_Mag; }
      EndIf
    }
  }

  If (GaugeType == 1)
    // We introduce a scalar Lagrange multiplier "xi" in H1_0 to enforce the
    // Coulomb gauge constraint; see explanations in the "Formulation" below.
    { Name H1_xi_Mag ; Type Form0 ;
      BasisFunction {
        { Name sn ; NameOfCoef xin ; Function BF_Node ;
          Support Vol_Mag; Entity NodesOf[All]; }
      }
      Constraint {
        { NameOfCoef xin; EntityType NodesOf; NameOfConstraint xi_Mag; }
      }
    }
  EndIf

  // Function space for the stranded-conductor global DoFs: one current (and its
  // dual voltage) per coil. Each basis function is the entire spatial
  // distribution of "vs" pre-computed by "Electrokinetics_vs~{i}" on the
  // corresponding "Coil~{i}", and the associated coefficient "ir" is the actual
  // current flowing through that coil.
  { Name Hregion_vs_Mag; Type Form0;
    BasisFunction {
      { Name sr; NameOfCoef ir;
        // "BF_Global" builds a basis function from a quantity computed by
        // another formulation -- here the unit-current source distribution
        // "-grad(vs)" pre-computed for each coil. The number of global basis
        // functions ("NumCoils == 2") should coincide with the number of
        // regions in "Group Vol_S_Mag" ("Coil_1" and "Coil_2"). The special
        // syntax "Electrokinetics_vs {NumCoils}" for the "Formulation" expands
        // to "Electrokinetics_vs_1" and "Electrokinetics_vs_2", matched
        // one-to-one with the regions of "Group Vol_S_Mag". "Ele {NumCoils}"
        // expands in the same way to "Ele_1" and "Ele_2" for the "Resolution":
        Function BF_Global {
          Quantity vs; Formulation Electrokinetics_vs {NumCoils};
          Group Vol_S_Mag; Resolution Ele {NumCoils};
        };
        // Each global basis function is directly associated with a "Global"
        // entity, i.e. each region in "Vol_S_Mag":
        Support Vol_S_Mag; Entity Global [Vol_S_Mag]; }
    }
    // "Current" is an alias of the coefficient "ir" (so "Current" is directly
    // the coil's current). "Voltage" is the dual global quantity associated
    // with "ir", representing the coil's terminal voltage; the link between
    // them is established by the "GlobalTerm" in the formulation and closed by
    // the "Network" constraint.
    GlobalQuantity {
      { Name Current; Type AliasOf; NameOfCoef ir; }
      { Name Voltage; Type AssociatedWith; NameOfCoef ir; }
    }
    // Both globals are left free here (cf. the empty "Current_Mag" and
    // "Voltage_Mag" constraints above): their runtime values are set by the
    // circuit equations.
    Constraint {
      { NameOfCoef Current; EntityType Global; NameOfConstraint Current_Mag; }
      { NameOfCoef Voltage; EntityType Global; NameOfConstraint Voltage_Mag; }
    }
  }

  { Name Hregion_i_Cir; Type Scalar;
    BasisFunction {
      { Name sr; NameOfCoef ir; Function BF_Region;
        Support Dom_Cir; Entity Dom_Cir; }
    }
    GlobalQuantity {
      { Name Current; Type AliasOf; NameOfCoef ir; }
      { Name Voltage; Type AssociatedWith; NameOfCoef ir; }
    }
    Constraint {
      { NameOfCoef Voltage; EntityType Region; NameOfConstraint Voltage_Cir; }
      { NameOfCoef Current; EntityType Region; NameOfConstraint Current_Cir; }
    }
  }

}

Formulation {
  { Name Magnetoquasistatics_av; Type FemEquation;
    Quantity {
      { Name a; Type Local; NameOfSpace Hcurl_a_Mag; }
      If(GaugeType == 1)
        { Name xi ; Type Local ; NameOfSpace H1_xi_Mag; }
      EndIf

      { Name vs; Type Local ; NameOfSpace Hregion_vs_Mag; }
      { Name Is; Type Global; NameOfSpace Hregion_vs_Mag [Current]; }
      { Name Us; Type Global; NameOfSpace Hregion_vs_Mag [Voltage]; }

      { Name Iz; Type Global; NameOfSpace Hregion_i_Cir [Current]; }
      { Name Uz; Type Global; NameOfSpace Hregion_i_Cir [Voltage]; }
    }
    Equation {
      Integral { [ nu[] * Dof{d a} , {d a} ];
        In Vol_Mag; Jacobian Vol; Integration Int; }

      // The Coulomb gauge weakly imposes div(a) = 0 via a Lagrange multiplier
      // saddle-point formulation. It picks, within each equivalence class "a +
      // grad(H1_0)", the unique representative L^2-orthogonal to the kernel:
      //
      //   (a, grad(xi')) = 0,   for all xi' in H1_0                 (*)
      //
      // By integration by parts ("xi' = 0" on "Bnd" cancels the boundary term),
      // (*) is equivalent to the variational form of "div(a) = 0", hence the
      // name.
      //
      // The constraint (*) is enforced by introducing a Lagrange multiplier "xi
      // in H1_0" and forming the Lagrangian
      //
      //   L(a, xi) = (1/2) * (nu * curl(a), curl(a)) - (j, a) + (a, grad(xi)).
      //
      // The two stationarity conditions give the saddle-point system. Variation
      // in "a" (replace "a" by "a + eps * a'", differentiate "L" with respect
      // to "eps" and set to zero at "eps = 0") gives, for any test direction
      // "a'", the augmented curl-curl equation
      //
      //   (nu * curl(a), curl(a')) + (grad(xi), a') = (j, a')       (1)
      //
      // The same procedure for "xi" affects only the multiplier-coupling term
      // and gives the gauge constraint
      //
      //   (a, grad(xi')) = 0                                        (2)
      //
      // (The term "(grad(xi), a')" in (1) is the transpose of "(a, grad(xi'))"
      // in (2), making the off-diagonal blocks of the block system transposes
      // of each other and the system symmetric.)
      //
      // Note: an "augmented Lagrangian" variant adds a penalty term "(gamma /
      // 2) * (div(a), div(a))" to "L", contributing an extra "gamma * (div(a),
      // div(a'))" to equation (1). With "gamma > 0" this improves the
      // conditioning of the otherwise indefinite system at the cost of tuning
      // "gamma"; the implementation below uses the bare Lagrangian ("gamma ==
      // 0"):
      If(GaugeType == 1)
        Integral { [ Dof{a}, {d xi} ] ;
          In Vol_Mag; Jacobian Vol; Integration Int; }
        Integral { [ Dof{d xi}, {a} ] ;
          In Vol_Mag; Jacobian Vol; Integration Int; }
      EndIf

      // The pre-computed source distribution "-grad(vs)" is normalized to a
      // unit current across the cross-section; the "Ns" factor scales it up to
      // the actual current density of the Ns-turn winding. The minus sign on
      // the dual term reflects "j_unit = -grad(vs)": the flux linkage "Lambda =
      // Ns * (a, j_unit) = -Ns * (a, grad(vs))" and its time derivative (the
      // back-EMF in the voltage equation) thus carry the opposite sign of the
      // source term:
      Integral { [ Ns[] * Dof{d vs} , {a} ];
        In Vol_S_Mag; Jacobian Vol; Integration Int; }

      Integral { DtDof[ -Ns[] * Dof{a} , {d vs} ];
        In Vol_S_Mag; Jacobian Vol; Integration Int; }

      // DC resistance of the winding. The "Ns^2" factor comes from Ns turns in
      // series, each one with cross-section Sc/Ns and length proportional to
      // Ns, i.e. R_total = Ns^2 * L / (sigma * Sc):
      Integral { [ Ns[]^2 / sigma[] * Dof{d vs} , {d vs} ];
        In Vol_S_Mag; Jacobian Vol; Integration Int; }

      GlobalTerm { [ Dof{Us} , {Is} ]; In Vol_S_Mag; }

      // The circuit coupling is handled in the same way as in tutorials 8 and
      // 9:
      GlobalTerm { [ Dof{Uz} , {Iz} ]; In Resistance_Cir; }
      GlobalTerm { [ Resistance[] * Dof{Iz} , {Iz} ]; In Resistance_Cir; }

      GlobalTerm { [ Dof{Uz} , {Iz} ]; In Inductance_Cir; }
      GlobalTerm { DtDof [ Inductance[] * Dof{Iz} , {Iz} ]; In Inductance_Cir; }

      GlobalTerm { [ Dof{Iz} , {Iz} ]; In Capacitance_Cir; }
      GlobalTerm { DtDof [ Capacitance[] * Dof{Uz} , {Iz} ]; In Capacitance_Cir; }

      GlobalEquation {
        Type Network; NameOfConstraint ElectricalCircuit;
        { Node {Is}; Loop {Us}; Equation {Us}; In Vol_S_Mag; }
        { Node {Iz}; Loop {Uz}; Equation {Uz}; In Dom_Cir; }
      }
    }
  }
}

Resolution {
  { Name Mag;
    System {
      { Name Sys_Mag; NameOfFormulation Magnetoquasistatics_av;
        Type Complex; Frequency f;
      }
    }
    Operation {
      Generate[Sys_Mag]; Solve[Sys_Mag]; SaveSolution[Sys_Mag];
    }
  }
}

PostProcessing {
  { Name Mag; NameOfFormulation Magnetoquasistatics_av;
    Quantity {
      { Name a;
        Value {
          Term { [ {a} ]; In Vol_Mag; Jacobian Vol; }
        }
      }
      { Name vs;
        Value {
          Term { [ {vs} ]; In Vol_S_Mag; Jacobian Vol; }
        }
      }
      { Name js;
        Value {
          Term { [ - {d vs} ]; In Vol_S_Mag; Jacobian Vol; }
        }
      }
      { Name b;
        Value {
          Term { [ {d a} ]; In Vol_Mag; Jacobian Vol; }
        }
      }
      { Name U; Value {
          Term { [ {Uz} ]; In Dom_Cir; }
          Term { [ {Us} ]; In Vol_S_Mag; }
        }
      }
      { Name I; Value {
          Term { [ {Iz} ]; In Dom_Cir; }
          Term { [ {Is} ]; In Vol_S_Mag; }
        }
      }
      { Name NormU; Value {
          Term { [ Norm[{Uz}] ]; In Dom_Cir; }
          Term { [ Norm[{Us}] ]; In Vol_S_Mag; }
        }
      }
      { Name NormI; Value {
          Term { [ Norm[{Iz}] ]; In Dom_Cir; }
          Term { [ Norm[{Is}] ]; In Vol_S_Mag; }
        }
      }
    }
  }
}


PostOperation {
  { Name Map; NameOfPostProcessing Mag;
    Operation {
      Print[ vs, OnElementsOf Vol_S_Mag, Format Gmsh, File "vs.pos" ];
      Print[ js, OnElementsOf Vol_S_Mag, Format Gmsh, File "js.pos" ];
      Print[ a, OnElementsOf Vol_Mag, Format Gmsh, File "a.pos" ];
      Print[ b, OnElementsOf Vol_Mag, Format Gmsh, File "b.pos" ];
      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"];
      Echo[ "Z_out", Format Table, File > "UI.txt" ];
      Print[ U, OnRegion Z_out, Format FrequencyTable, File > "UI.txt"];
      Print[ I, OnRegion Z_out, Format FrequencyTable, File > "UI.txt"];
      Print[ NormU, OnRegion E_in, Format Table, File "",
        SendToServer "}Output/0|V_in| [V]"{0} ];
      Print[ NormI, OnRegion E_in, Format Table, File "",
        SendToServer "}Output/1|I_in| [A]"{0} ];
      Print[ NormU, OnRegion Z_out, Format Table, File "",
        SendToServer "}Output/2|V_out| [V]"{0} ];
      Print[ NormI, OnRegion Z_out, Format Table, File "",
        SendToServer "}Output/3|I_out| [A]"{0} ];
    }
  }
}