Magnetostatic Problem on a Topologically Closed Conductor

Summary

Solves for the magnetic field around a topologically closed toroidal conductor carrying a net current, using a magnetic vector potential discretized using $var element = document.getElementById("moose-equation-b4dbc5ec-590a-42d8-94f9-346dde239126");katex.render("H(\\mathrm{curl})", element, {displayMode:false,throwOnError:false});$ conforming Nédélec elements, following the method described in P. Dular, International Compumag Society Newsletter, 7, no. 2 (2000): 4-7. for enforcing global constraints.

Description

This problem solves the magnetostatic problem with strong form:

$var element = document.getElementById("moose-equation-12e360b3-08dc-4276-af4c-9a594315d490");katex.render("\\begin{split} \\vec \\nabla \\cdot \\vec B &= 0 \\\\ \\vec \\nabla \\times \\vec H &= \\vec J \\\\ \\vec \\nabla \\times \\vec E &= - \\partial_t \\vec B = \\vec 0 \\end{split}", element, {displayMode:true,throwOnError:false});$

with the constitutive relations

$var element = document.getElementById("moose-equation-ce812f65-bc74-4c87-a876-3a6aab206977");katex.render("\\begin{split} \\vec B &= \\mu \\vec H \\\\ \\vec E &= \\rho \\vec J + \\vec E_{ext}\\\\ \\vec J &= \\sigma (\\vec E - \\vec E_{ext}) \\end{split}", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-d8c092c7-330f-4a92-af16-01b8e8ea936f");katex.render("\\vec B", element, {displayMode:false,throwOnError:false});$ is the magnetic flux density, $var element = document.getElementById("moose-equation-b976e049-f366-4d02-96e6-b7bfd0a49ed0");katex.render("\\vec H", element, {displayMode:false,throwOnError:false});$ is the magnetic field, $var element = document.getElementById("moose-equation-0627c407-428a-49d0-9961-d1d47c09985e");katex.render("\\vec J", element, {displayMode:false,throwOnError:false});$ is the current density, $var element = document.getElementById("moose-equation-09b19dba-64de-4fdf-ab89-e1c0edc9d888");katex.render("\\vec E", element, {displayMode:false,throwOnError:false});$ is the electric field, $var element = document.getElementById("moose-equation-f103ca36-9f10-44e1-ba84-faae68fa377f");katex.render("\\vec E_{ext}", element, {displayMode:false,throwOnError:false});$ is the externally applied electric field, $var element = document.getElementById("moose-equation-eb943ed4-67fb-4b1c-8c9a-8f72e98c5376");katex.render("\\mu", element, {displayMode:false,throwOnError:false});$ is the material permeability, $var element = document.getElementById("moose-equation-4f4c6ca1-6878-44c0-b7b8-036782dc428e");katex.render("\\rho", element, {displayMode:false,throwOnError:false});$ is the material resistivity, and $var element = document.getElementById("moose-equation-1eee8b8c-d474-4256-8a94-d6bbcead9a6d");katex.render("\\sigma", element, {displayMode:false,throwOnError:false});$ is the material conductivity. To enforce the divergence free condition on $var element = document.getElementById("moose-equation-09f73964-6f1a-4e0e-ab2d-d2c01fe66b07");katex.render("\\vec B", element, {displayMode:false,throwOnError:false});$ , we shall also introduce the magnetic vector potential $var element = document.getElementById("moose-equation-6abd6fb5-d594-4cac-9714-4c74a7f43e2a");katex.render("\\vec A", element, {displayMode:false,throwOnError:false});$ defined such that $var element = document.getElementById("moose-equation-5a820baf-1b07-4567-aca2-bbb2caf497d2");katex.render("\\vec \\nabla \\times \\vec A = \\vec B", element, {displayMode:false,throwOnError:false});$ .

In addition, we wish to impose a global voltage constraint on a toroidal conductor domain $var element = document.getElementById("moose-equation-c7adafe9-d4e0-4558-a68d-f76d842f3c89");katex.render("\\Omega_c", element, {displayMode:false,throwOnError:false});$ of the form

$var element = document.getElementById("moose-equation-26b405b7-01a2-403b-bc60-4ca93a478969");katex.render(" \\oint_{\\mathcal C} \\vec E \\cdot d\\vec l = -\\mathcal V", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-c9b8f849-2db6-4a3c-8031-abd42649da42");katex.render("\\mathcal C", element, {displayMode:false,throwOnError:false});$ is any path traversing around the conductor lying entirely within $var element = document.getElementById("moose-equation-299274a5-4cce-4b66-ae91-74ca10bb82a3");katex.render("\\Omega_c", element, {displayMode:false,throwOnError:false});$ with winding number equal to one, and $var element = document.getElementById("moose-equation-2294c335-509f-4ba3-a9a5-07e64bc98fe5");katex.render("\\mathcal V", element, {displayMode:false,throwOnError:false});$ is the loop voltage imposed on the conductor.

In a simply connected domain, satisfying $var element = document.getElementById("moose-equation-fd1f958d-9328-4de5-b941-87af6057c6d5");katex.render("\\vec \\nabla \\times \\vec E = \\vec 0", element, {displayMode:false,throwOnError:false});$ can be achieved by introducing a continuous single-valued scalar electric potential $var element = document.getElementById("moose-equation-13a80bbc-5d79-4054-88c6-fd25f860b1a0");katex.render("\\phi", element, {displayMode:false,throwOnError:false});$ such that $var element = document.getElementById("moose-equation-c3444883-0365-4757-b9f5-820c6cdef5a7");katex.render("\\vec E = -\\vec\\nabla \\phi", element, {displayMode:false,throwOnError:false});$ . However, in order for the electric field $var element = document.getElementById("moose-equation-e55e9668-72f8-46d4-a690-9a074719a0c4");katex.render("\\vec E", element, {displayMode:false,throwOnError:false});$ to be nonzero circulating the closed conductor, this is no longer sufficient, since $var element = document.getElementById("moose-equation-a0b0df52-5ab3-40cd-8325-907c5da334a6");katex.render("\\oint_{\\mathcal C} \\vec\\nabla \\phi \\cdot d\\vec l = 0", element, {displayMode:false,throwOnError:false});$ . We shall therefore represent the electric field by

$var element = document.getElementById("moose-equation-43490969-85d3-488d-a562-ee58ed17b7c9");katex.render(" \\vec E = - \\vec \\nabla \\phi + \\vec E_{ext}.", element, {displayMode:true,throwOnError:false});$

With this, we can arrive at the weak form for $var element = document.getElementById("moose-equation-5beb3342-0daf-4aa7-b6cf-feafdacb7bdd");katex.render("\\vec \\nabla \\times \\vec H = \\vec J", element, {displayMode:false,throwOnError:false});$ we intend to solve, of the form

$var element = document.getElementById("moose-equation-c98d66d6-5a83-478f-976f-c64ec05a57f9");katex.render("\\begin{split} (\\mu^{-1} \\vec \\nabla \\times \\vec A, \\vec \\nabla \\times \\vec A')_\\Omega - (\\vec H \\times \\vec n, \\vec A')_{\\partial \\Omega} +(\\sigma \\vec \\nabla \\phi, \\vec A')_{\\Omega_c} + (\\sigma \\vec E_{ext}, \\vec A')_\\Omega &= 0 \\,\\,\\, \\forall \\vec A' \\in V_A \\\\ (\\sigma \\vec \\nabla \\phi, \\vec \\nabla \\phi')_{\\Omega_c} +(\\sigma \\vec E_{ext}, \\vec \\nabla \\phi')_{\\Omega_c} &= 0 \\,\\,\\, \\forall \\phi' \\in V_\\phi \\\\ (\\sigma \\vec \\nabla \\phi, \\vec E_{ext}')_{\\Omega_c} +(\\sigma \\vec E_{ext}, \\vec E_{ext}')_{\\Omega_c} &= IV(\\vec E_{ext}') \\,\\,\\, \\forall \\vec E_{ext}' \\in V_E \\end{split}", element, {displayMode:true,throwOnError:false});$

where

$var element = document.getElementById("moose-equation-4a45eff5-2ba4-4de7-93f7-d2339dab270d");katex.render("\\begin{split} \\vec A, \\vec A' \\in H(\\mathrm{curl})(\\Omega) &: \\vec A \\times \\hat n= \\vec 0 \\,\\,\\, \\mathrm{on}\\,\\, \\partial\\Omega \\\\ \\phi, \\phi' \\in H^1(\\Omega_c) &: \\vec \\nabla \\phi \\cdot \\hat n= 0 \\,\\,\\, \\mathrm{on}\\,\\, \\partial\\Omega_c \\\\ \\vec E_{ext}, \\vec E_{ext}' \\in H(\\mathrm{curl})(\\Omega_c) &: \\vec \\nabla \\times \\vec E_{ext} = \\vec 0 \\,\\,\\, \\mathrm{in}\\,\\, \\Omega_c, \\, \\oint_{\\mathcal C} \\vec E_{ext} \\cdot d\\vec l = \\mathcal V \\end{split}", element, {displayMode:true,throwOnError:false});$

Discretizing the External Electric Field

We must now find a suitable choice of basis $var element = document.getElementById("moose-equation-17d0df5c-de9f-466f-800d-293ced5c988f");katex.render("{\\vec c_i \\in H(\\mathrm{curl})(\\Omega_c)}", element, {displayMode:false,throwOnError:false});$ for the external electric field, such that

$var element = document.getElementById("moose-equation-4b400f70-05f9-4649-9adb-0db603cae100");katex.render(" \\vec E_{ext} = - \\sum_{i \\in \\Omega_c} V_i \\vec c_i", element, {displayMode:true,throwOnError:false});$

for some unique set of degrees of freedom $var element = document.getElementById("moose-equation-e34b1f33-e25b-41d3-9fd2-62c0a93b7fb4");katex.render("\\{V_i\\}", element, {displayMode:false,throwOnError:false});$ . We shall also wish to restrict the support of $var element = document.getElementById("moose-equation-883235a4-68ad-4756-94d9-1c34dfec9b60");katex.render("\\{\\vec c_i\\}", element, {displayMode:false,throwOnError:false});$ to the smallest number of elements possible in $var element = document.getElementById("moose-equation-8672e7fc-7ce5-4d83-8fb8-c173c318e7c0");katex.render("\\Omega_c", element, {displayMode:false,throwOnError:false});$ to minimise the number of degrees of freedom required to represent $var element = document.getElementById("moose-equation-a73cfcc4-120a-4b22-8e10-bfa81b0893d2");katex.render("\\vec E_{ext}", element, {displayMode:false,throwOnError:false});$ .

To do this, we first identify a 'cut' surface $var element = document.getElementById("moose-equation-4a9ab667-a27a-4118-b79d-ec2cc829ac2b");katex.render("\\Gamma_c", element, {displayMode:false,throwOnError:false});$ of the toroidal conductor that is required to make it simply connected. Next, we define a one element wide 'transition' region $var element = document.getElementById("moose-equation-6c115d7b-6b0f-41dc-8304-1cd9f889f29e");katex.render("\\Omega_t \\subset \\Omega_c", element, {displayMode:false,throwOnError:false});$ on one side of the cut surface, using an MFEMCutTransitionSubMesh object, comprised of all elements with at least one vertex lying on the cut surface that lie on one side of the cut.

Next, to enforce $var element = document.getElementById("moose-equation-0b421530-f1b1-4873-b960-c037d5a0fb0e");katex.render("\\vec \\nabla \\times \\vec E_{ext} =0", element, {displayMode:false,throwOnError:false});$ in this region, we define a scalar field variable $var element = document.getElementById("moose-equation-c06031e6-9fe2-4fcd-82bf-bf5931cf6592");katex.render("v", element, {displayMode:false,throwOnError:false});$ defined on a (nodal) $var element = document.getElementById("moose-equation-a696996f-476b-4da6-acdf-83047545eaf0");katex.render("H^1", element, {displayMode:false,throwOnError:false});$ conforming finite element space on this submesh such that $var element = document.getElementById("moose-equation-a0b584d8-fedd-4f37-baba-bbd901171db2");katex.render("\\vec E_{ext} = -\\vec \\nabla v", element, {displayMode:false,throwOnError:false});$ :

$var element = document.getElementById("moose-equation-0d7100b2-27ba-4e54-913a-5816d0ed51bf");katex.render(" v = \\sum_{i \\in \\Omega_t} V_i v_i", element, {displayMode:true,throwOnError:false});$

so $var element = document.getElementById("moose-equation-f2f04c94-de05-4e80-beaf-b604ce34eced");katex.render("{\\vec c_i = \\vec \\nabla v_i}", element, {displayMode:false,throwOnError:false});$ .

Finally, to impose the global voltage constraint $var element = document.getElementById("moose-equation-bb73d519-7435-4311-86c8-86b2d47db2ed");katex.render("\\mathcal V = \\oint_{\\mathcal C} \\vec E_{ext} \\cdot d\\vec l", element, {displayMode:false,throwOnError:false});$ , we identify that

$var element = document.getElementById("moose-equation-49dcb67e-5cb6-44cb-8c03-d7078cede0da");katex.render(" \\mathcal V = \\oint_{\\mathcal C} \\vec E_{ext} \\cdot d\\vec l = v|_{\\Gamma_c} - v|_{\\partial \\Omega_t \\backslash \\Gamma_c}", element, {displayMode:true,throwOnError:false});$

since $var element = document.getElementById("moose-equation-b98f8bc4-3611-4899-a8e0-3842a3bcfbd8");katex.render("\\vec c_i = 0", element, {displayMode:false,throwOnError:false});$ outside $var element = document.getElementById("moose-equation-8d66fc07-a845-47cd-9c06-09408768d8e2");katex.render("\\Omega_t", element, {displayMode:false,throwOnError:false});$ . The loop voltage constraint is thus strongly imposed by Dirichlet conditions on $var element = document.getElementById("moose-equation-9a440abe-834b-4d07-ab03-49dcce1049c5");katex.render("v", element, {displayMode:false,throwOnError:false});$ , setting $var element = document.getElementById("moose-equation-4b9547c6-c1d4-4723-a74f-c7eb1f620df8");katex.render("v=0", element, {displayMode:false,throwOnError:false});$ on $var element = document.getElementById("moose-equation-dd0b3323-be62-40cd-bd38-baee73fdb4fb");katex.render("\\partial \\Omega_t \\backslash \\Gamma_c", element, {displayMode:false,throwOnError:false});$ and $var element = document.getElementById("moose-equation-dd85ad5e-b0f9-4b9d-a021-08797b37a10e");katex.render("v =\\mathcal V", element, {displayMode:false,throwOnError:false});$ on $var element = document.getElementById("moose-equation-53e8590c-c628-4dc1-87bd-03658dab0c63");katex.render("\\Gamma_c", element, {displayMode:false,throwOnError:false});$ . For lowest order $var element = document.getElementById("moose-equation-bb700c27-215c-4379-b0a8-0eb631c4df0e");katex.render("H^1", element, {displayMode:false,throwOnError:false});$ conforming representations of $var element = document.getElementById("moose-equation-fa061591-ddd5-42b0-b113-c76b8fd8b4bb");katex.render("v", element, {displayMode:false,throwOnError:false});$ , this fully constrains $var element = document.getElementById("moose-equation-7dbe3a65-f242-4fb5-9bcb-f9ded608337d");katex.render("\\vec E_{ext}", element, {displayMode:false,throwOnError:false});$ in the transition region.

This is shown explicitly in the example file:

# Solve for the electric field on a closed conductor subject to
# global loop voltage constraint.

initial_coil_domains = 'TorusCore TorusSheath'
coil_cut_surface = 'Cut'
coil_loop_voltage = -1.0
coil_conductivity = 1.0

[Problem<<<{"href": "../Problem/index.html"}>>>]
  type = MFEMProblem
[]

[Mesh<<<{"href": "../Mesh/index.html"}>>>]
    type = MFEMMesh
    file = ../mesh/embedded_concentric_torus.e
[]

[FunctorMaterials<<<{"href": "../FunctorMaterials/index.html"}>>>]
  [Conductor]
    type = MFEMGenericFunctorMaterial
    prop_names = conductivity
    prop_values = ${coil_conductivity}
  []
[]

[ICs<<<{"href": "../ICs/index.html"}>>>]
  [coil_external_potential_ic]
    type = MFEMScalarBoundaryIC
    variable = coil_external_potential
    boundary = ${coil_cut_surface}
    coefficient = ${coil_loop_voltage}
  []
[]

[SubMeshes]
  [cut]
    type = MFEMCutTransitionSubMesh
    cut_boundary = ${coil_cut_surface}
    block = ${initial_coil_domains}
    transition_subdomain = transition_dom
    transition_subdomain_boundary = transition_bdr
    closed_subdomain = coil_dom
  []
  [coil]
    type = MFEMDomainSubMesh
    block = coil_dom
    execution_order_group = 2
  []
[]

[FESpaces]
  [H1FESpace]
    type = MFEMScalarFESpace
    fec_type = H1
    fec_order = FIRST
  []
  [HCurlFESpace]
    type = MFEMVectorFESpace
    fec_type = ND
    fec_order = FIRST
  []
  [CoilH1FESpace]
    type = MFEMScalarFESpace
    fec_type = H1
    fec_order = FIRST
    submesh = coil
  []
  [CoilHCurlFESpace]
    type = MFEMVectorFESpace
    fec_type = ND
    fec_order = FIRST
    submesh = coil
  []
  [TransitionH1FESpace]
    type = MFEMScalarFESpace
    fec_type = H1
    fec_order = FIRST
    submesh = cut
  []
  [TransitionHCurlFESpace]
    type = MFEMVectorFESpace
    fec_type = ND
    fec_order = FIRST
    submesh = cut
  []
[]

[Variables<<<{"href": "../Variables/index.html"}>>>]
  [coil_induced_potential]
    type = MFEMVariable
    fespace = CoilH1FESpace
  []
[]

[AuxVariables<<<{"href": "../AuxVariables/index.html"}>>>]
  [coil_external_potential]
    type = MFEMVariable
    fespace = CoilH1FESpace
  []
  [transition_external_potential]
    type = MFEMVariable
    fespace = TransitionH1FESpace
  []
  [transition_external_e_field]
    type = MFEMVariable
    fespace = TransitionHCurlFESpace
  []
  [induced_potential]
    type = MFEMVariable
    fespace = H1FESpace
  []
  [induced_e_field]
    type = MFEMVariable
    fespace = HCurlFESpace
  []
  [external_e_field]
    type = MFEMVariable
    fespace = HCurlFESpace
  []
  [e_field]
    type = MFEMVariable
    fespace = HCurlFESpace
  []
[]

[AuxKernels<<<{"href": "../AuxKernels/index.html"}>>>]
  [update_induced_e_field]
    type = MFEMGradAux
    variable = induced_e_field
    source = induced_potential
    scale_factor = -1.0
    execute_on = TIMESTEP_END
  []
  [update_external_e_field]
    type = MFEMGradAux
    variable = transition_external_e_field
    source = transition_external_potential
    scale_factor = -1.0
    execute_on = TIMESTEP_END
  []
  [update_total_e_field]
    type = MFEMSumAux
    variable = e_field
    first_source_variable = induced_e_field
    second_source_variable = external_e_field
    execute_on = TIMESTEP_END
    execution_order_group = 3
  []
[]

[Kernels<<<{"href": "../Kernels/index.html"}>>>]
  [diff]
    type = MFEMDiffusionKernel
    variable = coil_induced_potential
    coefficient = conductivity
  []
  [source]
    type = MFEMMixedGradGradKernel
    trial_variable = coil_external_potential
    variable = coil_induced_potential
    coefficient = conductivity
    block = 'transition_dom'
  []
[]

[Solver]
  type = MFEMSuperLU
[]

[Executioner<<<{"href": "../Executioner/index.html"}>>>]
  type = MFEMSteady
[]

[Transfers<<<{"href": "../Transfers/index.html"}>>>]
  [submesh_transfer_from_coil]
    type = MFEMSubMeshTransfer
    from_variable = coil_induced_potential
    to_variable = induced_potential
    execute_on = TIMESTEP_END
  []
  [submesh_transfer_to_transition]
    type = MFEMSubMeshTransfer
    from_variable = coil_external_potential
    to_variable = transition_external_potential
    execute_on = TIMESTEP_END
  []
  [submesh_transfer_from_transition]
    type = MFEMSubMeshTransfer
    from_variable = transition_external_e_field
    to_variable = external_e_field
    execute_on = TIMESTEP_END
    execution_order_group = 2
  []
[]

[Outputs<<<{"href": "../Mastodon/Outputs/index.html"}>>>]
  [GlobalParaViewDataCollection]
    type = MFEMParaViewDataCollection
    file_base = OutputData/WholePotentialCoil
    vtk_format = ASCII
  []
  [TransitionParaViewDataCollection]
    type = MFEMParaViewDataCollection
    file_base = OutputData/CutPotentialCoil
    vtk_format = ASCII
    submesh = cut
  []
  [CoilParaViewDataCollection]
    type = MFEMParaViewDataCollection
    file_base = OutputData/Coil
    vtk_format = ASCII
    submesh = coil
  []
[]

Solving for the Total Electric Field

With $var element = document.getElementById("moose-equation-25704566-5fb7-4b83-9a39-72e72fcdebd2");katex.render("\\vec E_{ext}", element, {displayMode:false,throwOnError:false});$ now constrained, we now must solve for the remaining degrees of freedom in $var element = document.getElementById("moose-equation-7ebabd7f-bb0c-4d25-969f-968948a2a971");katex.render("\\vec E", element, {displayMode:false,throwOnError:false});$ arising due to surface charges in the coil as required for continuity of current in the conductor ( $var element = document.getElementById("moose-equation-8ff38e6b-83ad-4bef-9c0b-b356da8331d8");katex.render("\\vec \\nabla \\cdot \\vec J = 0", element, {displayMode:false,throwOnError:false});$ ), by solving

$var element = document.getElementById("moose-equation-57e90e1f-3de3-42f5-881d-8ba6b3316406");katex.render("(\\sigma \\vec \\nabla \\phi, \\vec \\nabla \\phi')_{\\Omega_c} +(\\sigma \\vec E_{ext}, \\vec \\nabla \\phi')_{\\Omega_t} = 0 \\,\\,\\, \\forall \\phi' \\in V_\\phi", element, {displayMode:true,throwOnError:false});$

for $var element = document.getElementById("moose-equation-6b5b221c-82c7-4bd4-8644-d4057c5e0eb5");katex.render("\\phi", element, {displayMode:false,throwOnError:false});$ . We then can calculate the total electric field from the sum of its induced and external components.

Solving for the Magnetic Flux Density

Finally, with the electric field in the conductor, we can now find the magnetic flux density around the current density flowing in the coil in the magnetostatic limit.

$var element = document.getElementById("moose-equation-b25f17bb-2f1e-4a98-8da0-ab46fee0f23f");katex.render("(\\mu^{-1} \\vec \\nabla \\times \\vec A, \\vec \\nabla \\times \\vec A')_\\Omega - (\\vec H \\times \\vec n, \\vec A')_{\\partial \\Omega} +(\\sigma \\vec \\nabla \\phi, \\vec A')_{\\Omega_c} + (\\sigma \\vec E_{ext}, \\vec A')_\\Omega = 0 \\,\\,\\, \\forall \\vec A' \\in V_A", element, {displayMode:true,throwOnError:false});$

This is typically solved using HypreAMS as a preconditioner. Note that the above magnetostatic problem is singular in non-conductive regions, and thus either singular = true should be passed to the preconditioner or a small mass term $var element = document.getElementById("moose-equation-3cb1c19d-b069-42c3-83b3-17ae0bd7e52a");katex.render("(\\epsilon \\vec A, \\vec A')_\\Omega", element, {displayMode:false,throwOnError:false});$ with $var element = document.getElementById("moose-equation-edc35759-8f4c-4dbd-b496-899b3163ee07");katex.render("\\epsilon << 1", element, {displayMode:false,throwOnError:false});$ should be added for stability.

Example File

The full closed coil magnetostatic example detailed above can be found below:

# Magnetostatic problem solved on a closed conductor subject to
# global loop voltage constraint.

[Mesh<<<{"href": "../Mesh/index.html"}>>>]
  type = MFEMMesh
  file = ../mesh/embedded_concentric_torus.e
[]

[Problem<<<{"href": "../Problem/index.html"}>>>]
  type = MFEMProblem
[]

[FESpaces]
  [HCurlFESpace]
    type = MFEMVectorFESpace
    fec_type = ND
    fec_order = FIRST
  []
  [HDivFESpace]
    type = MFEMVectorFESpace
    fec_type = RT
    fec_order = CONSTANT
  []
[]

[Variables<<<{"href": "../Variables/index.html"}>>>]
  [a_field]
    type = MFEMVariable
    fespace = HCurlFESpace
  []
[]

[AuxVariables<<<{"href": "../AuxVariables/index.html"}>>>]
  [b_field]
    type = MFEMVariable
    fespace = HDivFESpace
  []
  [e_field]
    type = MFEMVariable
    fespace = HCurlFESpace
  []
[]

[AuxKernels<<<{"href": "../AuxKernels/index.html"}>>>]
  [curl]
    type = MFEMCurlAux
    variable = b_field
    source = a_field
    scale_factor = 1.0
    execute_on = TIMESTEP_END
    execution_order_group = 3
  []
[]

[Functions<<<{"href": "../Functions/index.html"}>>>]
  [exact_a_field]
    type = ParsedVectorFunction<<<{"description": "Returns a vector function based on string descriptions for each component.", "href": "../../source/functions/MooseParsedVectorFunction.html"}>>>
    expression_x<<<{"description": "x-component of function."}>>> = '0'
    expression_y<<<{"description": "y-component of function."}>>> = '0'
    expression_z<<<{"description": "z-component of function."}>>> = '0'
  []
[]

[BCs<<<{"href": "../BCs/index.html"}>>>]
  [tangential_a_bdr]
    type = MFEMVectorTangentialDirichletBC
    variable = a_field
    vector_coefficient = exact_a_field
    boundary = 'Exterior'
  []
[]

[FunctorMaterials<<<{"href": "../FunctorMaterials/index.html"}>>>]
  [Vacuum]
    type = MFEMGenericFunctorMaterial
    prop_names = reluctivity
    prop_values = 1.0
  []
  [Conductor]
    type = MFEMGenericFunctorMaterial
    prop_names = conductivity
    prop_values = 1.0
    block = 'TorusCore TorusSheath'
  []
[]

[Kernels<<<{"href": "../Kernels/index.html"}>>>]
  [mass]
    type = MFEMVectorFEMassKernel
    variable = a_field
    coefficient = 1e-10
  []
  [curlcurl]
    type = MFEMCurlCurlKernel
    variable = a_field
    coefficient = reluctivity
  []
  [source]
    type = MFEMMixedVectorMassKernel
    variable = a_field
    trial_variable = e_field
    coefficient = conductivity
    block = 'TorusCore TorusSheath'
  []
[]

[Preconditioner]
  [ams]
    type = MFEMHypreAMS
    fespace = HCurlFESpace
  []
[]

[Solver]
  type = MFEMHyprePCG
  preconditioner = ams
  l_tol = 1e-14
  l_max_its = 1000
[]

[Executioner<<<{"href": "../Executioner/index.html"}>>>]
  type = MFEMSteady
  device = cpu
[]

[MultiApps<<<{"href": "../MultiApps/index.html"}>>>]
  [coil]
    type = FullSolveMultiApp<<<{"description": "Performs a complete simulation during each execution.", "href": "../../source/multiapps/FullSolveMultiApp.html"}>>>
    input_files<<<{"description": "The input file for each App.  If this parameter only contains one input file it will be used for all of the Apps.  When using 'positions_from_file' it is also admissable to provide one input_file per file."}>>> = cut_closed_coil.i
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = INITIAL
  []
[]

[Transfers<<<{"href": "../Transfers/index.html"}>>>]
  [from_coil]
    type = MultiAppMFEMCopyTransfer
    source_variable = e_field
    variable = e_field
    from_multi_app = coil
  []
[]

[Postprocessors<<<{"href": "../Postprocessors/index.html"}>>>]
  [CoilPower]
    type = MFEMVectorFEInnerProductIntegralPostprocessor
    coefficient = conductivity
    dual_variable = e_field
    primal_variable = e_field
    block = 'TorusCore TorusSheath'
  []
[]

[Outputs<<<{"href": "../Mastodon/Outputs/index.html"}>>>]
  [ParaViewDataCollection]
    type = MFEMParaViewDataCollection
    file_base = OutputData/MagnetostaticClosedCoil
    vtk_format = ASCII
  []
  [ReportedPostprocessors]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../source/outputs/CSV.html"}>>>
    file_base<<<{"description": "The desired solution output name without an extension. If not provided, MOOSE sets it with Outputs/file_base when available. Otherwise, MOOSE uses input file name and this object name for a master input or uses master file_base, the subapp name and this object name for a subapp input to set it."}>>> = OutputData/AVMagnetostaticClosedCoilCSV
  []
[]

Global Current Constraints

Strongly constraining the total current through the conductor instead of the loop voltage is also possible, by defining a cut plane in the domain outside the conductor region using the H-phi formulation instead of the A-V formulation detailed above; the details are similar, using a magnetic scalar potential instead of an electric scalar potential. The constraint to be enforced is

$var element = document.getElementById("moose-equation-a776afbe-a570-4132-a2d7-26588ad8f90f");katex.render("\\oint_{\\mathcal C} \\vec H \\cdot d\\vec l = I", element, {displayMode:true,throwOnError:false});$

where $var element = document.getElementById("moose-equation-2e62bbba-dcea-4cdd-9e40-5ace29d5586d");katex.render("\\mathcal C", element, {displayMode:false,throwOnError:false});$ is any path encircling the conductor, and I is the total current passing through the path $var element = document.getElementById("moose-equation-81777081-e180-4103-a39e-8347f0dd9764");katex.render("\\mathcal C", element, {displayMode:false,throwOnError:false});$ . Outside the conductor, the magnetic field $var element = document.getElementById("moose-equation-65af42ab-76ed-4ccd-b70c-b9dc21068f91");katex.render("\\vec H", element, {displayMode:false,throwOnError:false});$ is curl-free, and thus can be expressed by the gradient of a scalar magnetic potential (up to a sign).

An example file for this approach can be found below:

# Solve for the magnetic field around a closed conductor subject to
# global current constraint.

initial_vacuum_domains = 'Exterior'
vacuum_cut_surface = 'Cut'
conductor_current = 1.0
vacuum_permeability = 1.0

[Problem<<<{"href": "../Problem/index.html"}>>>]
  type = MFEMProblem
[]

[Mesh<<<{"href": "../Mesh/index.html"}>>>]
    type = MFEMMesh
    file = ../mesh/split_embedded_concentric_torus.e
[]

[FunctorMaterials<<<{"href": "../FunctorMaterials/index.html"}>>>]
  [Conductor]
    type = MFEMGenericFunctorMaterial
    prop_names = permeability
    prop_values = ${vacuum_permeability}
  []
[]

[ICs<<<{"href": "../ICs/index.html"}>>>]
  [vacuum_cut_potential_ic]
    type = MFEMScalarBoundaryIC
    variable = vacuum_cut_potential
    boundary = ${vacuum_cut_surface}
    coefficient = ${conductor_current}
  []
[]

[SubMeshes]
  [cut]
    type = MFEMCutTransitionSubMesh
    cut_boundary = ${vacuum_cut_surface}
    block = ${initial_vacuum_domains}
    transition_subdomain = transition_dom
    transition_subdomain_boundary = transition_bdr
    closed_subdomain = vacuum_dom
  []
  [vacuum]
    type = MFEMDomainSubMesh
    block = vacuum_dom
    execution_order_group = 2
  []
[]

[FESpaces]
  [VacuumH1FESpace]
    type = MFEMScalarFESpace
    fec_type = H1
    fec_order = FIRST
    submesh = vacuum
  []
  [VacuumHCurlFESpace]
    type = MFEMVectorFESpace
    fec_type = ND
    fec_order = FIRST
    submesh = vacuum
  []
  [TransitionH1FESpace]
    type = MFEMScalarFESpace
    fec_type = H1
    fec_order = FIRST
    submesh = cut
  []
  [TransitionHCurlFESpace]
    type = MFEMVectorFESpace
    fec_type = ND
    fec_order = FIRST
    submesh = cut
  []
  [HCurlFESpace]
    type = MFEMVectorFESpace
    fec_type = ND
    fec_order = FIRST
  []
[]

[Variables<<<{"href": "../Variables/index.html"}>>>]
  [vacuum_magnetic_potential]
    type = MFEMVariable
    fespace = VacuumH1FESpace
  []
[]

[AuxVariables<<<{"href": "../AuxVariables/index.html"}>>>]
  [vacuum_cut_potential]
    type = MFEMVariable
    fespace = VacuumH1FESpace
  []
  [transition_cut_potential]
    type = MFEMVariable
    fespace = TransitionH1FESpace
  []
  [transition_cut_function_field]
    type = MFEMVariable
    fespace = TransitionHCurlFESpace
  []
  [background_h_field]
    type = MFEMVariable
    fespace = VacuumHCurlFESpace
  []
  [cut_function_field]
    type = MFEMVariable
    fespace = VacuumHCurlFESpace
  []
  [vacuum_h_field]
    type = MFEMVariable
    fespace = VacuumHCurlFESpace
  []
  [h_field]
    type = MFEMVariable
    fespace = HCurlFESpace
  []
[]

[AuxKernels<<<{"href": "../AuxKernels/index.html"}>>>]
  [update_background_h_field]
    type = MFEMGradAux
    variable = background_h_field
    source = vacuum_magnetic_potential
    scale_factor = -1.0
    execute_on = TIMESTEP_END
  []
  [update_transition_cut_function_field]
    type = MFEMGradAux
    variable = transition_cut_function_field
    source = transition_cut_potential
    scale_factor = -1.0
    execute_on = TIMESTEP_END
  []
  [update_total_h_field]
    type = MFEMSumAux
    variable = vacuum_h_field
    first_source_variable = background_h_field
    second_source_variable = cut_function_field
    execute_on = TIMESTEP_END
    execution_order_group = 3
  []
[]

[BCs<<<{"href": "../BCs/index.html"}>>>]
  # Set zero of magnetic potential on symmetry plane
  [Exterior]
    type = MFEMScalarDirichletBC
    variable = vacuum_magnetic_potential
    boundary = 'Cut'
    coefficient = 0.0
  []
[]

[Kernels<<<{"href": "../Kernels/index.html"}>>>]
  [diff]
    type = MFEMDiffusionKernel
    variable = vacuum_magnetic_potential
    coefficient = permeability
  []
  [source]
    type = MFEMMixedGradGradKernel
    trial_variable = vacuum_cut_potential
    variable = vacuum_magnetic_potential
    coefficient = permeability
    block = 'transition_dom'
  []
[]

[Preconditioner]
  [boomeramg]
    type = MFEMHypreBoomerAMG
  []
[]

[Solver]
  type = MFEMHypreGMRES
  preconditioner = boomeramg
  l_tol = 1e-8
  l_max_its = 100
[]

[Executioner<<<{"href": "../Executioner/index.html"}>>>]
  type = MFEMSteady
[]

[Transfers<<<{"href": "../Transfers/index.html"}>>>]
  [submesh_transfer_to_transition]
    type = MFEMSubMeshTransfer
    from_variable = vacuum_cut_potential
    to_variable = transition_cut_potential
    execute_on = TIMESTEP_END
  []
  [submesh_transfer_from_transition]
    type = MFEMSubMeshTransfer
    from_variable = transition_cut_function_field
    to_variable = cut_function_field
    execute_on = TIMESTEP_END
    execution_order_group = 2
  []
  [submesh_transfer_from_vacuum]
    type = MFEMSubMeshTransfer
    from_variable = vacuum_h_field
    to_variable = h_field
    execute_on = TIMESTEP_END
    execution_order_group = 4
  []
[]

[Postprocessors<<<{"href": "../Postprocessors/index.html"}>>>]
  [MagneticEnergy]
    type = MFEMVectorFEInnerProductIntegralPostprocessor
    coefficient = ${fparse 0.5*vacuum_permeability}
    dual_variable = vacuum_h_field
    primal_variable = vacuum_h_field
    execution_order_group = 4
    block = 'Exterior'
  []
[]

[Outputs<<<{"href": "../Mastodon/Outputs/index.html"}>>>]
  [ReportedPostprocessors]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../source/outputs/CSV.html"}>>>
    file_base<<<{"description": "The desired solution output name without an extension. If not provided, MOOSE sets it with Outputs/file_base when available. Otherwise, MOOSE uses input file name and this object name for a master input or uses master file_base, the subapp name and this object name for a subapp input to set it."}>>> = OutputData/HPhiMagnetostaticClosedCoilCSV
  []
  [VacuumParaViewDataCollection]
    type = MFEMParaViewDataCollection
    file_base = OutputData/HPhiMagnetostaticClosedCoil
    vtk_format = ASCII
    submesh = vacuum
  []
[]

(moose/test/tests/mfem/submeshes/cut_closed_coil.i)

# Solve for the electric field on a closed conductor subject to
# global loop voltage constraint.

initial_coil_domains = 'TorusCore TorusSheath'
coil_cut_surface = 'Cut'
coil_loop_voltage = -1.0
coil_conductivity = 1.0

[Problem]
  type = MFEMProblem
[]

[Mesh]
    type = MFEMMesh
    file = ../mesh/embedded_concentric_torus.e
[]

[FunctorMaterials]
  [Conductor]
    type = MFEMGenericFunctorMaterial
    prop_names = conductivity
    prop_values = ${coil_conductivity}
  []
[]

[ICs]
  [coil_external_potential_ic]
    type = MFEMScalarBoundaryIC
    variable = coil_external_potential
    boundary = ${coil_cut_surface}
    coefficient = ${coil_loop_voltage}
  []
[]

[SubMeshes]
  [cut]
    type = MFEMCutTransitionSubMesh
    cut_boundary = ${coil_cut_surface}
    block = ${initial_coil_domains}
    transition_subdomain = transition_dom
    transition_subdomain_boundary = transition_bdr
    closed_subdomain = coil_dom
  []
  [coil]
    type = MFEMDomainSubMesh
    block = coil_dom
    execution_order_group = 2
  []
[]

[FESpaces]
  [H1FESpace]
    type = MFEMScalarFESpace
    fec_type = H1
    fec_order = FIRST
  []
  [HCurlFESpace]
    type = MFEMVectorFESpace
    fec_type = ND
    fec_order = FIRST
  []
  [CoilH1FESpace]
    type = MFEMScalarFESpace
    fec_type = H1
    fec_order = FIRST
    submesh = coil
  []
  [CoilHCurlFESpace]
    type = MFEMVectorFESpace
    fec_type = ND
    fec_order = FIRST
    submesh = coil
  []
  [TransitionH1FESpace]
    type = MFEMScalarFESpace
    fec_type = H1
    fec_order = FIRST
    submesh = cut
  []
  [TransitionHCurlFESpace]
    type = MFEMVectorFESpace
    fec_type = ND
    fec_order = FIRST
    submesh = cut
  []
[]

[Variables]
  [coil_induced_potential]
    type = MFEMVariable
    fespace = CoilH1FESpace
  []
[]

[AuxVariables]
  [coil_external_potential]
    type = MFEMVariable
    fespace = CoilH1FESpace
  []
  [transition_external_potential]
    type = MFEMVariable
    fespace = TransitionH1FESpace
  []
  [transition_external_e_field]
    type = MFEMVariable
    fespace = TransitionHCurlFESpace
  []
  [induced_potential]
    type = MFEMVariable
    fespace = H1FESpace
  []
  [induced_e_field]
    type = MFEMVariable
    fespace = HCurlFESpace
  []
  [external_e_field]
    type = MFEMVariable
    fespace = HCurlFESpace
  []
  [e_field]
    type = MFEMVariable
    fespace = HCurlFESpace
  []
[]

[AuxKernels]
  [update_induced_e_field]
    type = MFEMGradAux
    variable = induced_e_field
    source = induced_potential
    scale_factor = -1.0
    execute_on = TIMESTEP_END
  []
  [update_external_e_field]
    type = MFEMGradAux
    variable = transition_external_e_field
    source = transition_external_potential
    scale_factor = -1.0
    execute_on = TIMESTEP_END
  []
  [update_total_e_field]
    type = MFEMSumAux
    variable = e_field
    first_source_variable = induced_e_field
    second_source_variable = external_e_field
    execute_on = TIMESTEP_END
    execution_order_group = 3
  []
[]

[Kernels]
  [diff]
    type = MFEMDiffusionKernel
    variable = coil_induced_potential
    coefficient = conductivity
  []
  [source]
    type = MFEMMixedGradGradKernel
    trial_variable = coil_external_potential
    variable = coil_induced_potential
    coefficient = conductivity
    block = 'transition_dom'
  []
[]

[Solver]
  type = MFEMSuperLU
[]

[Executioner]
  type = MFEMSteady
[]

[Transfers]
  [submesh_transfer_from_coil]
    type = MFEMSubMeshTransfer
    from_variable = coil_induced_potential
    to_variable = induced_potential
    execute_on = TIMESTEP_END
  []
  [submesh_transfer_to_transition]
    type = MFEMSubMeshTransfer
    from_variable = coil_external_potential
    to_variable = transition_external_potential
    execute_on = TIMESTEP_END
  []
  [submesh_transfer_from_transition]
    type = MFEMSubMeshTransfer
    from_variable = transition_external_e_field
    to_variable = external_e_field
    execute_on = TIMESTEP_END
    execution_order_group = 2
  []
[]

[Outputs]
  [GlobalParaViewDataCollection]
    type = MFEMParaViewDataCollection
    file_base = OutputData/WholePotentialCoil
    vtk_format = ASCII
  []
  [TransitionParaViewDataCollection]
    type = MFEMParaViewDataCollection
    file_base = OutputData/CutPotentialCoil
    vtk_format = ASCII
    submesh = cut
  []
  [CoilParaViewDataCollection]
    type = MFEMParaViewDataCollection
    file_base = OutputData/Coil
    vtk_format = ASCII
    submesh = coil
  []
[]

(moose/test/tests/mfem/submeshes/av_magnetostatic.i)

# Magnetostatic problem solved on a closed conductor subject to
# global loop voltage constraint.

[Mesh]
  type = MFEMMesh
  file = ../mesh/embedded_concentric_torus.e
[]

[Problem]
  type = MFEMProblem
[]

[FESpaces]
  [HCurlFESpace]
    type = MFEMVectorFESpace
    fec_type = ND
    fec_order = FIRST
  []
  [HDivFESpace]
    type = MFEMVectorFESpace
    fec_type = RT
    fec_order = CONSTANT
  []
[]

[Variables]
  [a_field]
    type = MFEMVariable
    fespace = HCurlFESpace
  []
[]

[AuxVariables]
  [b_field]
    type = MFEMVariable
    fespace = HDivFESpace
  []
  [e_field]
    type = MFEMVariable
    fespace = HCurlFESpace
  []
[]

[AuxKernels]
  [curl]
    type = MFEMCurlAux
    variable = b_field
    source = a_field
    scale_factor = 1.0
    execute_on = TIMESTEP_END
    execution_order_group = 3
  []
[]

[Functions]
  [exact_a_field]
    type = ParsedVectorFunction
    expression_x = '0'
    expression_y = '0'
    expression_z = '0'
  []
[]

[BCs]
  [tangential_a_bdr]
    type = MFEMVectorTangentialDirichletBC
    variable = a_field
    vector_coefficient = exact_a_field
    boundary = 'Exterior'
  []
[]

[FunctorMaterials]
  [Vacuum]
    type = MFEMGenericFunctorMaterial
    prop_names = reluctivity
    prop_values = 1.0
  []
  [Conductor]
    type = MFEMGenericFunctorMaterial
    prop_names = conductivity
    prop_values = 1.0
    block = 'TorusCore TorusSheath'
  []
[]

[Kernels]
  [mass]
    type = MFEMVectorFEMassKernel
    variable = a_field
    coefficient = 1e-10
  []
  [curlcurl]
    type = MFEMCurlCurlKernel
    variable = a_field
    coefficient = reluctivity
  []
  [source]
    type = MFEMMixedVectorMassKernel
    variable = a_field
    trial_variable = e_field
    coefficient = conductivity
    block = 'TorusCore TorusSheath'
  []
[]

[Preconditioner]
  [ams]
    type = MFEMHypreAMS
    fespace = HCurlFESpace
  []
[]

[Solver]
  type = MFEMHyprePCG
  preconditioner = ams
  l_tol = 1e-14
  l_max_its = 1000
[]

[Executioner]
  type = MFEMSteady
  device = cpu
[]

[MultiApps]
  [coil]
    type = FullSolveMultiApp
    input_files = cut_closed_coil.i
    execute_on = INITIAL
  []
[]

[Transfers]
  [from_coil]
    type = MultiAppMFEMCopyTransfer
    source_variable = e_field
    variable = e_field
    from_multi_app = coil
  []
[]

[Postprocessors]
  [CoilPower]
    type = MFEMVectorFEInnerProductIntegralPostprocessor
    coefficient = conductivity
    dual_variable = e_field
    primal_variable = e_field
    block = 'TorusCore TorusSheath'
  []
[]

[Outputs]
  [ParaViewDataCollection]
    type = MFEMParaViewDataCollection
    file_base = OutputData/MagnetostaticClosedCoil
    vtk_format = ASCII
  []
  [ReportedPostprocessors]
    type = CSV
    file_base = OutputData/AVMagnetostaticClosedCoilCSV
  []
[]

(moose/test/tests/mfem/submeshes/hphi_magnetostatic.i)

# Solve for the magnetic field around a closed conductor subject to
# global current constraint.

initial_vacuum_domains = 'Exterior'
vacuum_cut_surface = 'Cut'
conductor_current = 1.0
vacuum_permeability = 1.0

[Problem]
  type = MFEMProblem
[]

[Mesh]
    type = MFEMMesh
    file = ../mesh/split_embedded_concentric_torus.e
[]

[FunctorMaterials]
  [Conductor]
    type = MFEMGenericFunctorMaterial
    prop_names = permeability
    prop_values = ${vacuum_permeability}
  []
[]

[ICs]
  [vacuum_cut_potential_ic]
    type = MFEMScalarBoundaryIC
    variable = vacuum_cut_potential
    boundary = ${vacuum_cut_surface}
    coefficient = ${conductor_current}
  []
[]

[SubMeshes]
  [cut]
    type = MFEMCutTransitionSubMesh
    cut_boundary = ${vacuum_cut_surface}
    block = ${initial_vacuum_domains}
    transition_subdomain = transition_dom
    transition_subdomain_boundary = transition_bdr
    closed_subdomain = vacuum_dom
  []
  [vacuum]
    type = MFEMDomainSubMesh
    block = vacuum_dom
    execution_order_group = 2
  []
[]

[FESpaces]
  [VacuumH1FESpace]
    type = MFEMScalarFESpace
    fec_type = H1
    fec_order = FIRST
    submesh = vacuum
  []
  [VacuumHCurlFESpace]
    type = MFEMVectorFESpace
    fec_type = ND
    fec_order = FIRST
    submesh = vacuum
  []
  [TransitionH1FESpace]
    type = MFEMScalarFESpace
    fec_type = H1
    fec_order = FIRST
    submesh = cut
  []
  [TransitionHCurlFESpace]
    type = MFEMVectorFESpace
    fec_type = ND
    fec_order = FIRST
    submesh = cut
  []
  [HCurlFESpace]
    type = MFEMVectorFESpace
    fec_type = ND
    fec_order = FIRST
  []
[]

[Variables]
  [vacuum_magnetic_potential]
    type = MFEMVariable
    fespace = VacuumH1FESpace
  []
[]

[AuxVariables]
  [vacuum_cut_potential]
    type = MFEMVariable
    fespace = VacuumH1FESpace
  []
  [transition_cut_potential]
    type = MFEMVariable
    fespace = TransitionH1FESpace
  []
  [transition_cut_function_field]
    type = MFEMVariable
    fespace = TransitionHCurlFESpace
  []
  [background_h_field]
    type = MFEMVariable
    fespace = VacuumHCurlFESpace
  []
  [cut_function_field]
    type = MFEMVariable
    fespace = VacuumHCurlFESpace
  []
  [vacuum_h_field]
    type = MFEMVariable
    fespace = VacuumHCurlFESpace
  []
  [h_field]
    type = MFEMVariable
    fespace = HCurlFESpace
  []
[]

[AuxKernels]
  [update_background_h_field]
    type = MFEMGradAux
    variable = background_h_field
    source = vacuum_magnetic_potential
    scale_factor = -1.0
    execute_on = TIMESTEP_END
  []
  [update_transition_cut_function_field]
    type = MFEMGradAux
    variable = transition_cut_function_field
    source = transition_cut_potential
    scale_factor = -1.0
    execute_on = TIMESTEP_END
  []
  [update_total_h_field]
    type = MFEMSumAux
    variable = vacuum_h_field
    first_source_variable = background_h_field
    second_source_variable = cut_function_field
    execute_on = TIMESTEP_END
    execution_order_group = 3
  []
[]

[BCs]
  # Set zero of magnetic potential on symmetry plane
  [Exterior]
    type = MFEMScalarDirichletBC
    variable = vacuum_magnetic_potential
    boundary = 'Cut'
    coefficient = 0.0
  []
[]

[Kernels]
  [diff]
    type = MFEMDiffusionKernel
    variable = vacuum_magnetic_potential
    coefficient = permeability
  []
  [source]
    type = MFEMMixedGradGradKernel
    trial_variable = vacuum_cut_potential
    variable = vacuum_magnetic_potential
    coefficient = permeability
    block = 'transition_dom'
  []
[]

[Preconditioner]
  [boomeramg]
    type = MFEMHypreBoomerAMG
  []
[]

[Solver]
  type = MFEMHypreGMRES
  preconditioner = boomeramg
  l_tol = 1e-8
  l_max_its = 100
[]

[Executioner]
  type = MFEMSteady
[]

[Transfers]
  [submesh_transfer_to_transition]
    type = MFEMSubMeshTransfer
    from_variable = vacuum_cut_potential
    to_variable = transition_cut_potential
    execute_on = TIMESTEP_END
  []
  [submesh_transfer_from_transition]
    type = MFEMSubMeshTransfer
    from_variable = transition_cut_function_field
    to_variable = cut_function_field
    execute_on = TIMESTEP_END
    execution_order_group = 2
  []
  [submesh_transfer_from_vacuum]
    type = MFEMSubMeshTransfer
    from_variable = vacuum_h_field
    to_variable = h_field
    execute_on = TIMESTEP_END
    execution_order_group = 4
  []
[]

[Postprocessors]
  [MagneticEnergy]
    type = MFEMVectorFEInnerProductIntegralPostprocessor
    coefficient = ${fparse 0.5*vacuum_permeability}
    dual_variable = vacuum_h_field
    primal_variable = vacuum_h_field
    execution_order_group = 4
    block = 'Exterior'
  []
[]

[Outputs]
  [ReportedPostprocessors]
    type = CSV
    file_base = OutputData/HPhiMagnetostaticClosedCoilCSV
  []
  [VacuumParaViewDataCollection]
    type = MFEMParaViewDataCollection
    file_base = OutputData/HPhiMagnetostaticClosedCoil
    vtk_format = ASCII
    submesh = vacuum
  []
[]

Summary
Description
Discretizing the External Electric Field
Solving for the Total Electric Field
Solving for the Magnetic Flux Density
Example File
Global Current Constraints

Usage

Development

Demonstration