mooseStreamAll(ss, args...);
  std::string msg = mooseMsgFmt(ss.str(), "*** Warning ***", COLOR_YELLOW);
  if (Moose::_throw_on_warning)
    throw std::runtime_error(msg);

  {
    libMesh::Threads::spin_mutex::scoped_lock lock(moose_stream_lock);

  MFEMSumAux(const InputParameters & parameters);

  virtual ~MFEMSumAux() = default;

  // Computes the auxvariable.
  virtual void execute() override;

  /// Get name of the trial variable (gridfunction) the kernel acts on.
  /// Defaults to the name of the test variable labelling the weak form.
  virtual const VariableName & getTrialVariableName() const override { return _var_dot_name; }

protected:
  /// Name of variable (gridfunction) representing time derivative of variable.

registerMooseObject("MooseApp", MFEMSumAux);

InputParameters
MFEMSumAux::validParams()
{
  InputParameters params = MFEMAuxKernel::validParams();
  params.addClassDescription(
      "Calculates the sum of two variables sharing an FE space, each optionally scaled by a real "
      "constant, and stores the result in a third.");
  params.addRequiredParam<VariableName>("first_source_variable", "First variable to sum.");
  params.addRequiredParam<VariableName>("second_source_variable", "Second variable to sum.");
  params.addParam<mfem::real_t>(
      "first_scale_factor", 1.0, "Factor to scale the first variable by prior to sum.");
  params.addParam<mfem::real_t>(
      "second_scale_factor", 1.0, "Factor to scale the second variable by prior to sum.");
  return params;
}

MFEMSumAux::MFEMSumAux(const InputParameters & parameters)
  : MFEMAuxKernel(parameters),
    _v1_var_name(getParam<VariableName>("first_source_variable")),
    _v2_var_name(getParam<VariableName>("second_source_variable")),
    _v1_var(*getMFEMProblem().getProblemData().gridfunctions.Get(_v1_var_name)),
    _v2_var(*getMFEMProblem().getProblemData().gridfunctions.Get(_v2_var_name)),
    _lambda1(getParam<mfem::real_t>("first_scale_factor")),
    _lambda2(getParam<mfem::real_t>("second_scale_factor"))
{
}

void
MFEMSumAux::execute()
{
  // result = lambda1 * v1 + lambda2 * v2
  add(_lambda1, _v1_var, _lambda2, _v2_var, _result_var);
}

#endif

registerMooseObject("MooseApp", MFEMScalarBoundaryIC);

InputParameters
MFEMScalarBoundaryIC::validParams()
{
  auto params = MFEMInitialCondition::validParams();
  params += MFEMBoundaryRestrictable::validParams();
  params.addClassDescription("Sets the initial values of an MFEM scalar variable from a "
                             "user-specified scalar coefficient.");
  params.addRequiredParam<MFEMScalarCoefficientName>("coefficient", "The scalar coefficient");
  return params;
}

MFEMScalarBoundaryIC::MFEMScalarBoundaryIC(const InputParameters & params)
  : MFEMInitialCondition(params),
    MFEMBoundaryRestrictable(params,
                             *getMFEMProblem()
                                  .getProblemData()
                                  .gridfunctions.GetRef(getParam<VariableName>("variable"))
                                  .ParFESpace()
                                  ->GetParMesh())
{
}

void
MFEMScalarBoundaryIC::execute()
{
  auto & coeff = getScalarCoefficient("coefficient");
  auto grid_function = getMFEMProblem().getGridFunction(getParam<VariableName>("variable"));
  grid_function->ProjectBdrCoefficient(coeff, getBoundaryMarkers());
}

#endif

mfem::Array<int>
MFEMBlockRestrictable::subdomainsToAttributes(const std::vector<SubdomainName> & subdomain_names)
{
  mfem::Array<int> attributes;
  auto & mesh = getMesh();

  for (const SubdomainName & subdomain_name : subdomain_names)
  {
    try
    {
      // Is this a block ID?
      const int attribute_id = std::stoi(subdomain_name);
      attributes.Append(attribute_id);
    }
    catch (...)
    {
      // It was not
      auto & subdomain_ids = mesh.attribute_sets.GetAttributeSet(subdomain_name);
      for (const auto & subdomain_id : subdomain_ids)
        attributes.Append(subdomain_id);
    }
  }
  return attributes;
}

registerMooseObject("MooseApp", MFEMMixedGradGradKernel);

InputParameters
MFEMMixedGradGradKernel::validParams()
{
  InputParameters params = MFEMMixedBilinearFormKernel::validParams();
  params.addClassDescription(
      "Adds the domain integrator to an MFEM problem for the mixed bilinear form "
      "$(k\\vec\\nabla u, \\vec\\nabla v)_\\Omega$.");
  params.addParam<MFEMScalarCoefficientName>("coefficient", "1.", "Name of property k to use.");
  return params;
}

MFEMMixedGradGradKernel::MFEMMixedGradGradKernel(const InputParameters & parameters)
  : MFEMMixedBilinearFormKernel(parameters), _coef(getScalarCoefficient("coefficient"))
// FIXME: The MFEM bilinear form can also handle vector and matrix
// coefficients, so ideally we'd handle all three too.
{
}

mfem::BilinearFormIntegrator *
MFEMMixedGradGradKernel::createMBFIntegrator()
{
  return new mfem::MixedGradGradIntegrator(_coef);
}

#endif

registerMooseObject("MooseApp", MFEMMixedVectorMassKernel);

InputParameters
MFEMMixedVectorMassKernel::validParams()
{
  InputParameters params = MFEMMixedBilinearFormKernel::validParams();
  params.addClassDescription(
      "Adds the domain integrator to an MFEM problem for the mixed bilinear form "
      "$(k\\vec u, \\vec v)_\\Omega$.");
  params.addParam<MFEMScalarCoefficientName>("coefficient", "1.", "Name of property k to use.");
  return params;
}

MFEMMixedVectorMassKernel::MFEMMixedVectorMassKernel(const InputParameters & parameters)
  : MFEMMixedBilinearFormKernel(parameters), _coef(getScalarCoefficient("coefficient"))
// FIXME: The MFEM bilinear form can also handle vector and matrix
// coefficients, so ideally we'd handle all three too.
{
}

mfem::BilinearFormIntegrator *
MFEMMixedVectorMassKernel::createMBFIntegrator()
{
  return new mfem::MixedVectorMassIntegrator(_coef);
}

#endif

registerMooseObject("MooseApp", MFEMTimeDerivativeVectorFEMassKernel);

InputParameters
MFEMTimeDerivativeVectorFEMassKernel::validParams()
{
  InputParameters params = MFEMVectorFEMassKernel::validParams();
  params.addClassDescription("Adds the domain integrator to an MFEM problem for the bilinear form "
                             "$(k \\dot{\\vec u}, \\vec v)_\\Omega$ "
                             "arising from the weak form of the operator "
                             "$k \\dot{\\vec u}$.");
  return params;
}

MFEMTimeDerivativeVectorFEMassKernel::MFEMTimeDerivativeVectorFEMassKernel(
    const InputParameters & parameters)
  : MFEMVectorFEMassKernel(parameters),
    _var_dot_name(Moose::MFEM::GetTimeDerivativeName(_test_var_name))
{
}

#endif

registerMooseObject("MooseApp", MFEMVectorFEInnerProductIntegralPostprocessor);

InputParameters
MFEMVectorFEInnerProductIntegralPostprocessor::validParams()
{
  InputParameters params = MFEMPostprocessor::validParams();
  params += MFEMBlockRestrictable::validParams();
  params.addClassDescription("Calculates the total flux of a vector field through an interface");
  params.addParam<MFEMScalarCoefficientName>(
      "coefficient", "1.", "Name of optional scalar coefficient to scale integrand by.");
  params.addRequiredParam<VariableName>("primal_variable",
                                        "Name of the first vector variable in the inner product.");
  params.addRequiredParam<VariableName>("dual_variable",
                                        "Name of the second vector variable in the inner product.");
  return params;
}

MFEMVectorFEInnerProductIntegralPostprocessor::MFEMVectorFEInnerProductIntegralPostprocessor(
    const InputParameters & parameters)
  : MFEMPostprocessor(parameters),
    MFEMBlockRestrictable(parameters, getMFEMProblem().mesh().getMFEMParMesh()),
    _primal_var(getMFEMProblem().getProblemData().gridfunctions.GetRef(
        getParam<VariableName>("primal_variable"))),
    _dual_var(getMFEMProblem().getProblemData().gridfunctions.GetRef(
        getParam<VariableName>("dual_variable"))),
    _scalar_coef(getScalarCoefficient("coefficient")),
    _dual_var_coef(&_dual_var),
    _scaled_dual_var_coef(_scalar_coef, _dual_var_coef),
    _subdomain_integrator(_primal_var.ParFESpace())
{
  if (isSubdomainRestricted())
  {
    _subdomain_integrator.AddDomainIntegrator(
        new mfem::VectorFEDomainLFIntegrator(_scaled_dual_var_coef), getSubdomainMarkers());
  }
  else
  {
    _subdomain_integrator.AddDomainIntegrator(
        new mfem::VectorFEDomainLFIntegrator(_scaled_dual_var_coef));
  }
}

void
MFEMVectorFEInnerProductIntegralPostprocessor::execute()
{
  _subdomain_integrator.Assemble();
  _integral = _subdomain_integrator(_primal_var);
}

PostprocessorValue
MFEMVectorFEInnerProductIntegralPostprocessor::getValue() const
{
  return _integral;
}

#endif

registerMooseObject("MooseApp", MFEMCutTransitionSubMesh);

InputParameters
MFEMCutTransitionSubMesh::validParams()
{
  InputParameters params = MFEMSubMesh::validParams();
  params += MFEMBlockRestrictable::validParams();
  params.addClassDescription(
      "Class to construct an MFEMSubMesh formed from the set of elements that have least one "
      "vertex on the specified cut plane, that lie on one side of the plane, "
      "and that are restricted to the set of user-specified subdomains.");
  params.addRequiredParam<BoundaryName>("cut_boundary",
                                        "The boundary associated with the mesh cut.");
  params.addRequiredParam<BoundaryName>(
      "transition_subdomain_boundary",
      "Name to assign boundary of transition subdomain not shared with cut surface.");
  params.addRequiredParam<SubdomainName>(
      "transition_subdomain",
      "The name of the subdomain to be created on the mesh comprised of the set of elements "
      "adjacent to the cut surface on one side.");
  params.addRequiredParam<SubdomainName>(
      "closed_subdomain",
      "The name of the subdomain attribute to be created comprised of the set of all elements "
      "of the closed geometry, including the new transition region.");
  return params;
}

MFEMCutTransitionSubMesh::MFEMCutTransitionSubMesh(const InputParameters & parameters)
  : MFEMSubMesh(parameters),
    MFEMBlockRestrictable(parameters, getMFEMProblem().mesh().getMFEMParMesh()),
    _cut_boundary(getParam<BoundaryName>("cut_boundary")),
    _cut_submesh(std::make_shared<mfem::ParSubMesh>(mfem::ParSubMesh::CreateFromBoundary(
        getMFEMProblem().mesh().getMFEMParMesh(),
        getMesh().bdr_attribute_sets.GetAttributeSet(_cut_boundary)))),
    _transition_subdomain_boundary(getParam<BoundaryName>("transition_subdomain_boundary")),
    _transition_subdomain(getParam<SubdomainName>("transition_subdomain")),
    _closed_subdomain(getParam<SubdomainName>("closed_subdomain")),
    _cut_normal(3)
{
}

void
MFEMCutTransitionSubMesh::buildSubMesh()
{
  labelMesh(getMFEMProblem().mesh().getMFEMParMesh());
  _submesh = std::make_shared<mfem::ParSubMesh>(mfem::ParSubMesh::CreateFromDomain(
      getMesh(), getMesh().attribute_sets.GetAttributeSet(_transition_subdomain)));
  _submesh->attribute_sets.attr_sets = getMesh().attribute_sets.attr_sets;
  _submesh->bdr_attribute_sets.attr_sets = getMesh().bdr_attribute_sets.attr_sets;
}

void
MFEMCutTransitionSubMesh::labelMesh(mfem::ParMesh & parent_mesh)
{
  int mpi_comm_rank = 0;
  int mpi_comm_size = 1;
  MPI_Comm_rank(getMFEMProblem().mesh().getMFEMParMesh().GetComm(), &mpi_comm_rank);
  MPI_Comm_size(getMFEMProblem().mesh().getMFEMParMesh().GetComm(), &mpi_comm_size);

  // First determine face normal based on the first boundary element found on the cut
  // to use when determining orientation relative to the cut
  const mfem::Array<int> & parent_cut_element_id_map = _cut_submesh->GetParentElementIDMap();
  int rank_with_submesh = -1;
  if (parent_cut_element_id_map.Size() > 0)
  {
    int reference_face = parent_cut_element_id_map[0];
    _cut_normal = findFaceNormal(parent_mesh, parent_mesh.GetBdrElementFaceIndex(reference_face));
    rank_with_submesh = mpi_comm_rank;
  }
  MPI_Allreduce(MPI_IN_PLACE,
                &rank_with_submesh,
                1,
                MPI_INT,
                MPI_MAX,
                getMFEMProblem().mesh().getMFEMParMesh().GetComm());
  MPI_Bcast(_cut_normal.GetData(),
            _cut_normal.Size(),
            MPI_DOUBLE,
            rank_with_submesh,
            getMFEMProblem().mesh().getMFEMParMesh().GetComm());
  // Iterate over all vertices on cut, find elements with those vertices,
  // and declare them transition elements if they are on the +ve side of the cut
  mfem::Array<int> transition_els;
  std::vector<HYPRE_BigInt> global_cut_vert_ids;
  mfem::Array<HYPRE_BigInt> gi;
  parent_mesh.GetGlobalVertexIndices(gi);
  std::unique_ptr<mfem::Table> vert_to_elem(parent_mesh.GetVertexToElementTable());
  const mfem::Array<int> & cut_to_parent_vertex_id_map = _cut_submesh->GetParentVertexIDMap();
  for (int i = 0; i < _cut_submesh->GetNV(); ++i)
  {
    int cut_vert = cut_to_parent_vertex_id_map[i];
    global_cut_vert_ids.push_back(gi[cut_vert]);
    int ne = vert_to_elem->RowSize(cut_vert); // number of elements touching cut vertex
    const int * els_adj_to_cut = vert_to_elem->GetRow(cut_vert); // elements touching cut vertex
    for (int i = 0; i < ne; i++)
    {
      const int el_adj_to_cut = els_adj_to_cut[i];
      if (isInDomain(el_adj_to_cut, getSubdomainAttributes(), parent_mesh) &&
          sideOfCut(el_adj_to_cut, cut_vert, parent_mesh) == 1)
        transition_els.Append(el_adj_to_cut);
    }
  }

  // share cut verts coords or global ids across all procs
  int n_cut_vertices = global_cut_vert_ids.size();
  std::vector<int> cut_vert_sizes(mpi_comm_size, 0);
  MPI_Allgather(&n_cut_vertices,
                1,
                MPI_INT,
                &cut_vert_sizes[0],

                MPI_INT,
                getMFEMProblem().mesh().getMFEMParMesh().GetComm());
  // Make an offset array and total the sizes.
  std::vector<int> n_vert_offset(mpi_comm_size, 0);
  for (int i = 1; i < mpi_comm_size; i++)
    n_vert_offset[i] = n_vert_offset[i - 1] + cut_vert_sizes[i - 1];
  int global_n_cut_vertices = 0;
  for (int i = 0; i < mpi_comm_size; i++)
    global_n_cut_vertices += cut_vert_sizes[i];

  // Gather the queries to all ranks.
  std::vector<HYPRE_BigInt> all_cut_verts(global_n_cut_vertices, 0);
  MPI_Allgatherv(&global_cut_vert_ids[0],
                 n_cut_vertices,
                 HYPRE_MPI_BIG_INT,
                 &all_cut_verts[0],

                 getMFEMProblem().mesh().getMFEMParMesh().GetComm());

  // Detect shared vertices and add corresponding elements
  for (int g = 1, sv = 0; g < parent_mesh.GetNGroups(); g++)
  {
    for (int gv = 0; gv < parent_mesh.GroupNVertices(g); gv++, sv++)
    {
      // all elements touching this shared vertex should be updated
      int cut_vert = parent_mesh.GroupVertex(g, gv);
      for (std::size_t i = 0; i < all_cut_verts.size(); i += 1)
      {
        if (gi[cut_vert] == all_cut_verts[i]) // check if shared vertex is on the cut plane
        {
          int ne = vert_to_elem->RowSize(cut_vert); // number of elements touching cut vertex
          const int * els_adj_to_cut =
              vert_to_elem->GetRow(cut_vert); // elements touching cut vertex
          for (int i = 0; i < ne; i++)
          {
            const int el_adj_to_cut = els_adj_to_cut[i];
            if (isInDomain(el_adj_to_cut, getSubdomainAttributes(), parent_mesh) &&
                sideOfCut(el_adj_to_cut, cut_vert, parent_mesh) == 1)
              transition_els.Append(el_adj_to_cut);
          }
        }
      }
    }
  }

  transition_els.Sort();
  transition_els.Unique();

  setAttributes(parent_mesh, transition_els);
}

mfem::Vector
MFEMCutTransitionSubMesh::findFaceNormal(const mfem::ParMesh & mesh, const int & face)
{
  MFEM_ASSERT(mesh->SpatialDimension() == 3,
              "MFEMCutTransitionSubMesh only works in 3-dimensional meshes!");
  mfem::Vector normal;
  mfem::Array<int> face_verts;
  std::vector<mfem::Vector> v;
  mesh.GetFaceVertices(face, face_verts);

  // First we get the coordinates of 3 vertices on the face
  for (auto vtx : face_verts)
  {
    mfem::Vector vtx_coords(3);
    for (int j = 0; j < 3; ++j)
      vtx_coords[j] = mesh.GetVertex(vtx)[j];
    v.push_back(vtx_coords);
  }

  // Now we find the unit vector normal to the face
  v[0] -= v[1];
  v[1] -= v[2];
  v[0].cross3D(v[1], normal);
  normal /= normal.Norml2();
  return normal;
}

int
MFEMCutTransitionSubMesh::sideOfCut(const int & el,
                                    const int & el_vertex_on_cut,
                                    mfem::ParMesh & parent_mesh)
{
  const int sdim = parent_mesh.SpaceDimension();
  mfem::Vector el_center(3);
  parent_mesh.GetElementCenter(el, el_center);
  mfem::Vector vertex_coords(parent_mesh.GetVertex(el_vertex_on_cut), sdim);
  mfem::Vector relative_center(sdim);
  for (int j = 0; j < sdim; j++)
    relative_center[j] = el_center[j] - vertex_coords[j];
  double side = _cut_normal * relative_center;
  if (side > 0)
    return 1;
  else
    return -1;
}

void
MFEMCutTransitionSubMesh::setAttributes(mfem::ParMesh & parent_mesh,
                                        mfem::Array<int> & transition_els)
{
  // Generate a set of local new attributes for transition region elements
  const int old_max_attr = parent_mesh.attributes.Max();
  mfem::Array<int> new_attrs(old_max_attr);
  for (int i = 0; i < new_attrs.Size(); ++i)
    new_attrs[i] = i + 1;
  for (int i = 0; i < transition_els.Size(); ++i)
  {
    const auto & transition_el = transition_els[i];
    const int old_attr = parent_mesh.GetAttribute(transition_el);
    new_attrs[old_attr - 1] = old_max_attr + old_attr % old_max_attr + 1;
  }

  // Distribute local attribute IDs to other MPI ranks to create set of new
  // global attribute IDs for the transition region.
  mfem::Array<int> global_new_attrs(old_max_attr);
  global_new_attrs = -1;
  MPI_Allreduce(new_attrs,
                global_new_attrs,
                old_max_attr,
                MPI_INT,

                getMFEMProblem().mesh().getMFEMParMesh().GetComm());

  // Set the new attribute IDs for transition region elements
  for (int i = 0; i < transition_els.Size(); ++i)
  {
    const auto & transition_el = transition_els[i];
    const int old_attr = parent_mesh.GetAttribute(transition_el);
    if (old_attr < old_max_attr + 1)
      parent_mesh.SetAttribute(transition_el, global_new_attrs[old_attr - 1]);
  }

  mfem::AttributeSets & attr_sets = parent_mesh.attribute_sets;
  mfem::AttributeSets & bdr_attr_sets = parent_mesh.bdr_attribute_sets;
  // Create attribute set labelling the newly created transition region
  // on one side of the cut
  attr_sets.CreateAttributeSet(_transition_subdomain);
  // Create attribute set labelling the entire closed geometry
  attr_sets.SetAttributeSet(_closed_subdomain, getSubdomainAttributes());
  // Add the new domain attributes to new attribute sets
  const std::set<std::string> attr_set_names = attr_sets.GetAttributeSetNames();
  for (int old_attr = 1; old_attr < global_new_attrs.Size() + 1; ++old_attr)
  {
    int new_attr = global_new_attrs[old_attr - 1];
    // Add attributes only if they're transition region attributes
    if (new_attr > old_max_attr)
    {
      attr_sets.AddToAttributeSet(_transition_subdomain, new_attr);
      for (auto const & attr_set_name : attr_set_names)
      {
        const mfem::Array<int> & attr_set = attr_sets.GetAttributeSet(attr_set_name);
        // If attribute set has the old attribute of the transition region, add the new one
        if (attr_set.Find(old_attr) != -1)
          attr_sets.AddToAttributeSet(attr_set_name, new_attr);
      }
    }
  }

  // Create boundary attribute set labelling the exterior of the newly created
  // transition region, excluding the cut
  bdr_attr_sets.CreateAttributeSet(_transition_subdomain_boundary);
  bdr_attr_sets.AddToAttributeSet(_transition_subdomain_boundary,
                                  parent_mesh.bdr_attributes.Max() + 1);

  parent_mesh.SetAttributes();
}

bool
MFEMCutTransitionSubMesh::isInDomain(const int & element,
                                     const mfem::Array<int> & subdomains,
                                     const mfem::ParMesh & mesh)
{
  // element<0 for ghost elements
  if (element < 0)
    return true;

  for (const auto & subdomain : subdomains)
    if (mesh.GetAttribute(element) == subdomain)
      return true;
  return false;
}

#endif

  params += MFEMBlockRestrictable::validParams();
  params.addClassDescription("Class to construct an MFEMSubMesh formed from the subspace of the "
                             "parent mesh restricted to the set of user-specified subdomains.");
  params.addParam<BoundaryName>("submesh_boundary", "Name to assign submesh boundary.");
  return params;
}

  _submesh->attribute_sets.attr_sets = getMesh().attribute_sets.attr_sets;
  _submesh->bdr_attribute_sets.attr_sets = getMesh().bdr_attribute_sets.attr_sets;

  if (isParamSetByUser("submesh_boundary"))
  {
    const BoundaryName & submesh_boundary = getParam<BoundaryName>("submesh_boundary");
    _submesh->bdr_attribute_sets.CreateAttributeSet(submesh_boundary);
    _submesh->bdr_attribute_sets.AddToAttributeSet(submesh_boundary,
                                                   getMesh().bdr_attributes.Max() + 1);
  }
}