diff --git a/contrib/hephaestus b/contrib/hephaestus
index f6105d01..6368185b 160000
--- a/contrib/hephaestus
+++ b/contrib/hephaestus
@@ -1 +1 @@
-Subproject commit f6105d01171420dc64d9f2d4f8a9a3f0159e4a00
+Subproject commit 6368185b7fa86f581684718eb9396e7aadbdca25
diff --git a/examples/MFEM/coil/AFormCoil.i b/examples/MFEM/coil/AFormCoil.i
index 8681c9e1..3dfcf1cc 100644
--- a/examples/MFEM/coil/AFormCoil.i
+++ b/examples/MFEM/coil/AFormCoil.i
@@ -6,6 +6,7 @@
 
 [Problem]
   type = MFEMProblem
+  use_glvis = true
 []
 
 [Formulation]
@@ -14,6 +15,7 @@
   magnetic_reluctivity_name = magnetic_reluctivity
   magnetic_permeability_name = magnetic_permeability
   electric_conductivity_name = electrical_conductivity
+  magnetic_flux_density_name = magnetic_flux_density
 []
 
 [FESpaces]
@@ -43,6 +45,10 @@
     type = MFEMVariable
     fespace = HDivFESpace
   []
+  [source_current_density]
+    type = MFEMVariable
+    fespace = HCurlFESpace
+  []
   [electric_potential]
     type = MFEMVariable
     fespace = H1FESpace
@@ -50,18 +56,12 @@
 []
 
 [Functions]
-  [potential_high]
+  [current_magnitude]
     type = ParsedFunction
     value = cos(2.0*pi*freq*t)
     vars = 'freq'
     vals = '0.01666667'
   []
-  [potential_low]
-    type = ParsedFunction
-    value = -cos(2.0*pi*freq*t)
-    vars = 'freq'
-    vals = '0.01666667'
-  []
 []
 
 [BCs]
@@ -71,18 +71,6 @@
     vector_coefficient = TangentialECoef
     boundary = '1 2 3 4'
   []
-  [high_terminal]
-    type = MFEMScalarDirichletBC
-    variable = electric_potential
-    boundary = '1'
-    function = potential_high
-  []
-  [low_terminal]
-    type = MFEMScalarDirichletBC
-    variable = electric_potential
-    boundary = '2'
-    function = potential_low
-  []
 []
 
 [Materials]
@@ -158,20 +146,20 @@
     value = 0.0
   []
 
-  [OneCoef]
-    type = MFEMConstantCoefficient
-    value = 1.0e12
+  [CurrentCoef]
+    type = MFEMFunctionCoefficient
+    function = current_magnitude
   []
 []
 
 [Sources]
   [SourcePotential]
-    type = MFEMScalarPotentialSource
-    potential = electric_potential
-    conductivity = electrical_conductivity
-    hcurl_fespace = HCurlFESpace
-    h1_fespace = H1FESpace
-    solver_max_its = 1000
+    type = MFEMOpenCoilSource
+    total_current_coef = CurrentCoef
+    source_current_density_gridfunction = source_current_density
+    source_potential_gridfunction = electric_potential
+    coil_in_boundary = 1
+    coil_out_boundary = 2
     block = 1
   []
 []
diff --git a/include/sources/MFEMOpenCoilSource.h b/include/sources/MFEMOpenCoilSource.h
new file mode 100644
index 00000000..4673ec72
--- /dev/null
+++ b/include/sources/MFEMOpenCoilSource.h
@@ -0,0 +1,29 @@
+#pragma once
+#include "MFEMSource.h"
+#include "MFEMVariable.h"
+#include "MFEMCoefficient.h"
+
+class MFEMOpenCoilSource : public MFEMSource
+{
+public:
+  static InputParameters validParams();
+
+  MFEMOpenCoilSource(const InputParameters & parameters);
+  virtual ~MFEMOpenCoilSource();
+
+  virtual void execute() override {}
+  virtual void initialize() override {}
+  virtual void finalize() override {}
+
+  virtual hephaestus::Source * getSource() override;
+  virtual void storeCoefficients(hephaestus::Coefficients & coefficients) override;
+
+protected:
+  const MFEMVariable & _source_current_density_gridfunction;
+  const MFEMVariable & _source_potential_gridfunction;  
+  const MFEMCoefficient & _total_current_coef;
+  const hephaestus::InputParameters _open_coil_params;
+  mfem::Array<int> _coil_domains;
+  const int _coil_in_id, _coil_out_id;
+  const std::pair<int,int> _electrodes;
+};
diff --git a/src/sources/MFEMOpenCoilSource.C b/src/sources/MFEMOpenCoilSource.C
new file mode 100644
index 00000000..878df240
--- /dev/null
+++ b/src/sources/MFEMOpenCoilSource.C
@@ -0,0 +1,62 @@
+#include "MFEMOpenCoilSource.h"
+
+registerMooseObject("ApolloApp", MFEMOpenCoilSource);
+
+InputParameters
+MFEMOpenCoilSource::validParams()
+{
+  InputParameters params = MFEMSource::validParams();
+
+  params.addParam<UserObjectName>(
+      "total_current_coef",
+      "The total current ($I$) flowing through the coil. May be time dependent.");
+
+  params.addParam<UserObjectName>("source_current_density_gridfunction",
+                                  "The gridfunction to store the source current density.");
+
+  params.addParam<UserObjectName>("source_potential_gridfunction",
+                                  "The gridfunction to store the scalar potential in the coil.");
+
+  params.addParam<BoundaryName>(
+      "coil_in_boundary",
+      "The boundary (id or name) from the mesh determinining the coil input terminal");
+
+  params.addParam<BoundaryName>(
+      "coil_out_boundary",
+      "The boundary (id or name) from the mesh determinining the coil output terminal");
+  return params;
+}
+
+MFEMOpenCoilSource::MFEMOpenCoilSource(const InputParameters & parameters)
+  : MFEMSource(parameters),
+    _source_current_density_gridfunction(
+        getUserObject<MFEMVariable>("source_current_density_gridfunction")),
+    _source_potential_gridfunction(getUserObject<MFEMVariable>("source_potential_gridfunction")),
+    _total_current_coef(getUserObject<MFEMCoefficient>("total_current_coef")),
+    _open_coil_params({{"SourceName", _source_current_density_gridfunction.name()},
+                       {"PotentialName", _source_potential_gridfunction.name()},
+                       {"IFuncCoefName", _total_current_coef.name()}}),
+    _coil_domains(blocks.size()),
+    _coil_in_id(std::stoi(getParam<BoundaryName>("coil_in_boundary"))),
+    _coil_out_id(std::stoi(getParam<BoundaryName>("coil_out_boundary"))),
+    _electrodes(std::make_pair(_coil_in_id, _coil_out_id))
+{
+  for (unsigned int bid = 0; bid < blocks.size(); ++bid)
+  {
+    _coil_domains[bid] = blocks[bid];
+  }
+  _source = new hephaestus::OpenCoilSolver(_open_coil_params, _coil_domains, _electrodes);
+}
+
+hephaestus::Source *
+MFEMOpenCoilSource::getSource()
+{
+  return _source;
+}
+
+void
+MFEMOpenCoilSource::storeCoefficients(hephaestus::Coefficients & coefficients)
+{
+}
+
+MFEMOpenCoilSource::~MFEMOpenCoilSource() {}