EquationSystem.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#ifdef MOOSE_MFEM_ENABLED

#include "EquationSystem.h"
#include "MFEMLinearSolverBase.h"
#include "libmesh/int_range.h"

namespace Moose::MFEM
{

EquationSystem::~EquationSystem()
{
  DeleteHBlocks();
  DeleteJacobianBlocks();
}

void
EquationSystem::DeleteHBlocks()
{
  for (const auto i : make_range(_h_blocks.NumRows()))
    for (const auto j : make_range(_h_blocks.NumCols()))
    {
      if (_jacobian_blocks.NumRows() && _jacobian_blocks(i, j) == _h_blocks(i, j))
        _jacobian_blocks(i, j) = nullptr;
      delete _h_blocks(i, j);
    }
  _h_blocks.DeleteAll();
}

void
EquationSystem::DeleteJacobianBlocks()
{
  for (const auto i : make_range(_jacobian_blocks.NumRows()))
    for (const auto j : make_range(_jacobian_blocks.NumCols()))
      if (!_h_blocks.NumRows() || _jacobian_blocks(i, j) != _h_blocks(i, j))
        delete _jacobian_blocks(i, j);
  _jacobian_blocks.DeleteAll();
}

bool
EquationSystem::VectorContainsName(const std::vector<std::string> & the_vector,
                                   const std::string & name) const
{
  return std::find(the_vector.begin(), the_vector.end(), name) != the_vector.end();
}

void
EquationSystem::AddCoupledVariableNameIfMissing(const std::string & coupled_var_name)
{
  if (!VectorContainsName(_coupled_var_names, coupled_var_name))
    _coupled_var_names.push_back(coupled_var_name);
}

void
EquationSystem::AddEliminatedVariableNameIfMissing(const std::string & eliminated_var_name)
{
  if (!VectorContainsName(_eliminated_var_names, eliminated_var_name))
    _eliminated_var_names.push_back(eliminated_var_name);
}

void
EquationSystem::AddTestVariableNameIfMissing(const std::string & test_var_name)
{
  if (!VectorContainsName(_test_var_names, test_var_name))
    _test_var_names.push_back(test_var_name);
}

void
EquationSystem::SetTrialVariableNames()
{
  // If a coupled variable has an equation associated with it,
  // add it to the set of trial variables.
  for (const auto & test_var_name : _test_var_names)
    if (VectorContainsName(_coupled_var_names, test_var_name))
      _trial_var_names.push_back(test_var_name);

  // Otherwise, add it to the set of eliminated variables.
  for (const auto & coupled_var_name : _coupled_var_names)
    if (!VectorContainsName(_test_var_names, coupled_var_name))
      _eliminated_var_names.push_back(coupled_var_name);
}

void
EquationSystem::AddKernel(std::shared_ptr<MFEMKernel> kernel)
{
  const auto & trial_var_name = kernel->getTrialVariableName();
  const auto & test_var_name = kernel->getTestVariableName();
  AddCoupledVariableNameIfMissing(trial_var_name);
  AddTestVariableNameIfMissing(test_var_name);
  // Register new kernels map if not present for the test variable
  if (!_kernels_map.Has(test_var_name))
  {
    auto kernel_field_map =
        std::make_shared<Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>>();
    _kernels_map.Register(test_var_name, std::move(kernel_field_map));
  }
  // Register new kernels map if not present for the test/trial variable pair
  if (!_kernels_map.Get(test_var_name)->Has(trial_var_name))
  {
    auto kernels = std::make_shared<std::vector<std::shared_ptr<MFEMKernel>>>();
    _kernels_map.Get(test_var_name)->Register(trial_var_name, std::move(kernels));
  }
  _kernels_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(kernel));
}

void
EquationSystem::AddIntegratedBC(std::shared_ptr<MFEMIntegratedBC> bc)
{
  const auto & trial_var_name = bc->getTrialVariableName();
  const auto & test_var_name = bc->getTestVariableName();
  AddCoupledVariableNameIfMissing(trial_var_name);
  AddTestVariableNameIfMissing(test_var_name);
  // Register new integrated bc map if not present for the test variable
  if (!_integrated_bc_map.Has(test_var_name))
  {
    auto integrated_bc_field_map = std::make_shared<
        Moose::MFEM::NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>>();
    _integrated_bc_map.Register(test_var_name, std::move(integrated_bc_field_map));
  }
  // Register new integrated bc map if not present for the test/trial variable pair
  if (!_integrated_bc_map.Get(test_var_name)->Has(trial_var_name))
  {
    auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMIntegratedBC>>>();
    _integrated_bc_map.Get(test_var_name)->Register(trial_var_name, std::move(bcs));
  }
  _integrated_bc_map.GetRef(test_var_name).Get(trial_var_name)->push_back(std::move(bc));
}

void
EquationSystem::AddEssentialBC(std::shared_ptr<MFEMEssentialBC> bc)
{
  const auto & test_var_name = bc->getTestVariableName();
  AddTestVariableNameIfMissing(test_var_name);
  // Register new essential bc map if not present for the test variable
  if (!_essential_bc_map.Has(test_var_name))
  {
    auto bcs = std::make_shared<std::vector<std::shared_ptr<MFEMEssentialBC>>>();
    _essential_bc_map.Register(test_var_name, std::move(bcs));
  }
  _essential_bc_map.GetRef(test_var_name).push_back(std::move(bc));
}

void
EquationSystem::Init(Moose::MFEM::GridFunctions & gridfunctions,
                     Moose::MFEM::ComplexGridFunctions & cmplx_gridfunctions,
                     mfem::AssemblyLevel assembly_level)
{
  _assembly_level = assembly_level;

  if (cmplx_gridfunctions.size())
    mooseError("Complex variables have been created but the executioner numeric type has not been "
               "set to complex. Please set Executioner/numeric_type = complex.");

  // Extract which coupled variables are to be trivially eliminated and which are trial variables
  SetTrialVariableNames();

  for (auto & test_var_name : _test_var_names)
  {
    if (!gridfunctions.Has(test_var_name))
    {
      mooseError("MFEM variable ",
                 test_var_name,
                 " requested by equation system during initialization was "
                 "not found in gridfunctions");
    }
    // Store pointers to test FESpaces
    _test_pfespaces.push_back(gridfunctions.Get(test_var_name)->ParFESpace());
  }

  for (auto & trial_var_name : _trial_var_names)
  {
    if (!gridfunctions.Has(trial_var_name))
    {
      mooseError("MFEM variable ",
                 trial_var_name,
                 " requested by equation system during initialization was "
                 "not found in gridfunctions");
    }
    // Create auxiliary gridfunctions for storing essential constraints from Dirichlet conditions
    _var_ess_constraints.emplace_back(
        std::make_unique<mfem::ParGridFunction>(gridfunctions.Get(trial_var_name)->ParFESpace()));
  }

  // Store pointers to FESpaces of all coupled variables
  for (auto & coupled_var_name : _coupled_var_names)
    _coupled_pfespaces.push_back(gridfunctions.Get(coupled_var_name)->ParFESpace());

  // Store pointers to coupled variable GridFunctions that are to be eliminated prior to forming the
  // jacobian
  for (auto & eliminated_var_name : _eliminated_var_names)
    _eliminated_variables.Register(eliminated_var_name,
                                   gridfunctions.GetShared(eliminated_var_name));

  // Get a reference to the GridFunctions
  _gfuncs = &gridfunctions;
}

void
EquationSystem::ApplyEssentialBC(const std::string & var_name,
                                 mfem::ParGridFunction & trial_gf,
                                 mfem::Array<int> & global_ess_markers)
{
  if (_essential_bc_map.Has(var_name))
  {
    auto & bcs = _essential_bc_map.GetRef(var_name);
    for (auto & bc : bcs)
    {
      // Set constrained DoFs values on essential boundaries
      bc->ApplyBC(trial_gf);
      // Fetch marker array labelling essential boundaries of current BC
      mfem::Array<int> ess_bdrs(bc->getBoundaryMarkers());
      // Add these boundary markers to the set of markers labelling all essential boundaries
      for (const auto i : make_range(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max()))
        global_ess_markers[i] = std::max(global_ess_markers[i], ess_bdrs[i]);
    }
  }
}

void
EquationSystem::ApplyEssentialBCs()
{
  _ess_tdof_lists.resize(_trial_var_names.size());
  for (const auto i : index_range(_trial_var_names))
  {
    const auto & trial_var_name = _trial_var_names.at(i);
    mfem::ParGridFunction & trial_gf = *_var_ess_constraints.at(i);

    // Make sure we update the size, if this mesh has changed recently for instance
    trial_gf.Update();

    // Initial guess for iterative solvers (initial condition or the previous time step solution)
    trial_gf = _gfuncs->GetRef(trial_var_name);

    mfem::Array<int> global_ess_markers(trial_gf.ParFESpace()->GetParMesh()->bdr_attributes.Max());
    global_ess_markers = 0;
    // Set strongly constrained DoFs of trial_gf on essential boundaries and add markers for all
    // essential boundaries to the global_ess_markers array
    ApplyEssentialBC(trial_var_name, trial_gf, global_ess_markers);
    trial_gf.ParFESpace()->GetEssentialTrueDofs(global_ess_markers, _ess_tdof_lists.at(i));
  }
}

void
EquationSystem::EliminateCoupledVariables()
{
  for (const auto & test_var_name : _test_var_names)
    for (const auto & eliminated_var_name : _eliminated_var_names)
      if (_mblfs.Has(test_var_name) && _mblfs.Get(test_var_name)->Has(eliminated_var_name) &&
          !VectorContainsName(_test_var_names, eliminated_var_name))
      {
        auto & mblf = *_mblfs.Get(test_var_name)->Get(eliminated_var_name);
        mblf.AddMult(*_eliminated_variables.Get(eliminated_var_name), *_lfs.Get(test_var_name), -1);
      }
}

void
EquationSystem::FormLinearSystem(mfem::OperatorHandle & op,
                                 mfem::BlockVector & trueX,
                                 mfem::BlockVector & trueRHS)
{
  mooseAssert(_test_var_names.size() == _trial_var_names.size(),
              "Number of test and trial variables must be the same for block matrix assembly.");

  if (_assembly_level == mfem::AssemblyLevel::LEGACY)
    FormSystemMatrix(op, trueX, trueRHS);
  else
  {
    mooseAssert(_test_var_names.size() == 1 && _test_var_names.size() == _trial_var_names.size(),
                "Non-legacy assembly is only supported for single test and trial variable systems");
    FormSystemOperator(op, trueX, trueRHS);
  }
}

void
EquationSystem::FormSystemOperator(mfem::OperatorHandle & op,
                                   mfem::BlockVector & trueX,
                                   mfem::BlockVector & trueRHS)
{
  auto & test_var_name = _test_var_names.at(0);
  mfem::Vector aux_x, aux_rhs;
  mfem::OperatorPtr aux_a;

  auto blf = _blfs.Get(test_var_name);
  blf->FormLinearSystem(_ess_tdof_lists.at(0),
                        *_var_ess_constraints.at(0),
                        *_lfs.Get(test_var_name),
                        aux_a,
                        aux_x,
                        aux_rhs,
                        /*copy_interior=*/true);

  trueX.GetBlock(0) = aux_x;
  trueRHS.GetBlock(0) = aux_rhs;
  trueX.SyncFromBlocks();
  trueRHS.SyncFromBlocks();

  op.Reset(aux_a.Ptr());
  aux_a.SetOperatorOwner(false);
}

void
EquationSystem::FormSystemMatrix(mfem::OperatorHandle & op,
                                 mfem::BlockVector & trueX,
                                 mfem::BlockVector & trueRHS)
{
  // Allocate block operator
  DeleteHBlocks();
  _h_blocks.SetSize(_test_var_names.size(), _trial_var_names.size());
  _h_blocks = nullptr;
  // Zero out RHS and sync memory
  trueRHS = 0.0;
  trueRHS.SyncToBlocks();

  for (const auto i : index_range(_test_var_names))
  {
    auto test_var_name = _test_var_names.at(i);

    for (const auto j : index_range(_trial_var_names))
    {
      auto trial_var_name = _trial_var_names.at(j);

      mfem::Vector aux_x, aux_rhs;
      mfem::ParLinearForm aux_lf(_test_pfespaces.at(i));
      mfem::HypreParMatrix * aux_a = new mfem::HypreParMatrix;

      if (test_var_name == trial_var_name)
      {
        mooseAssert(i == j, "Trial and test variables must have the same ordering.");
        auto blf = _blfs.Get(test_var_name);
        blf->FormLinearSystem(_ess_tdof_lists.at(j),
                              *_var_ess_constraints.at(j),
                              *_lfs.Get(test_var_name),
                              *aux_a,
                              aux_x,
                              aux_rhs,
                              /*copy_interior=*/true);
        trueX.GetBlock(j) = aux_x;
      }
      else if (_mblfs.Has(test_var_name) && _mblfs.Get(test_var_name)->Has(trial_var_name))
      {
        auto mblf = _mblfs.Get(test_var_name)->Get(trial_var_name);
        mblf->FormRectangularLinearSystem(_ess_tdof_lists.at(j),
                                          _ess_tdof_lists.at(i),
                                          *_var_ess_constraints.at(j),
                                          aux_lf = 0,
                                          *aux_a,
                                          aux_x,
                                          aux_rhs);
      }
      else
        continue;

      trueRHS.GetBlock(i) += aux_rhs;
      _h_blocks(i, j) = aux_a;
    }
  }
  // Sync memory
  trueX.SyncFromBlocks();
  trueRHS.SyncFromBlocks();

  // Create monolithic matrix
  op.Reset(mfem::HypreParMatrixFromBlocks(_h_blocks));
}

void
EquationSystem::FormSystem(mfem::BlockVector & trueX, mfem::BlockVector & trueRHS)
{
  height = trueX.Size();
  width = trueRHS.Size();
  // Store block offsets
  _block_true_offsets.SetSize(trueX.NumBlocks() + 1);
  _block_true_offsets[0] = 0;
  for (unsigned i = 0; i < _trial_var_names.size(); i++)
    _block_true_offsets[i + 1] = trueX.BlockSize(i);
  _block_true_offsets.PartialSum();
  FormLinearSystem(_linear_operator, trueX, trueRHS);
}

void
EquationSystem::Mult(const mfem::Vector & sol, mfem::Vector & residual) const
{
  if (_non_linear)
  {
    ComputeNonlinearResidual(sol, residual);
    _linear_operator->AddMult(sol, residual);
  }
  else
  {
    residual = 0.0;
    _linear_operator->Mult(sol, residual);
  }

  sol.HostRead();
  residual.HostRead();
}

void
EquationSystem::ComputeNonlinearResidual(const mfem::Vector & sol, mfem::Vector & residual) const
{
  mooseAssert(_non_linear, "Should not be calling this method if our forms are not nonlinear");
  residual = 0.0;

  const mfem::BlockVector block_solution(const_cast<mfem::Vector &>(sol), _block_true_offsets);
  SetTrialVariablesFromTrueVectors(block_solution);

  mfem::BlockVector block_residual(residual, _block_true_offsets);
  for (unsigned int i = 0; i < _test_var_names.size(); i++)
  {
    auto & test_var_name = _test_var_names.at(i);
    auto nlf = _nlfs.GetShared(test_var_name);
    nlf->Mult(block_solution.GetBlock(i), block_residual.GetBlock(i));
    block_residual.GetBlock(i).SyncAliasMemory(block_residual);
  }
}

void
EquationSystem::FormJacobianMatrix(const mfem::Vector & u)
{
  DeleteJacobianBlocks();
  _jacobian_blocks.SetSize(_test_var_names.size(), _trial_var_names.size());
  _jacobian_blocks = nullptr;

  const mfem::BlockVector update_vector(const_cast<mfem::Vector &>(u), _block_true_offsets);
  for (const auto i : index_range(_test_var_names))
  {
    auto test_var_name = _test_var_names.at(i);
    if (_nlfs.Has(test_var_name))
    {
      auto nlf = _nlfs.Get(test_var_name);
      mfem::HypreParMatrix * nlf_jac =
          dynamic_cast<mfem::HypreParMatrix *>(&nlf->GetGradient(update_vector.GetBlock(i)));
      mooseAssert(nlf_jac,
                  "Jacobian contribution of nonlinear form associated with " + test_var_name +
                      " is not castable into a HypreParMatrix");
      _jacobian_blocks(i, i) = mfem::ParAdd(_h_blocks(i, i), nlf_jac);
    }
    else
      _jacobian_blocks(i, i) = _h_blocks(i, i);
    for (const auto j : index_range(_trial_var_names))
      if (i != j) // nlf->GetGradient only contributes to on-diagonal blocks
        _jacobian_blocks(i, j) = _h_blocks(i, j);
  }
  // Create monolithic matrix
  _jacobian.Reset(mfem::HypreParMatrixFromBlocks(_jacobian_blocks));
}

mfem::Operator &
EquationSystem::GetGradient(const mfem::Vector & u) const
{
  if (_non_linear)
  {
    if (_assembly_level != mfem::AssemblyLevel::LEGACY)
      mooseError("MFEM nonlinear solvers that require GetGradient() currently require legacy "
                 "assembly in EquationSystem.");
    const_cast<EquationSystem *>(this)->FormJacobianMatrix(u);
  }
  else
    _jacobian = _linear_operator;

  return *_jacobian;
}

void
EquationSystem::SetTrialVariablesFromTrueVectors(const mfem::BlockVector & trueX) const
{
  for (const auto i : index_range(_trial_var_names))
  {
    auto & trial_var_name = _trial_var_names.at(i);
    trueX.GetBlock(i).SyncMemory(trueX);
    _gfuncs->Get(trial_var_name)->Distribute(&(trueX.GetBlock(i)));
  }
}

void
EquationSystem::BuildLinearForms()
{
  // Register linear forms
  for (const auto i : index_range(_test_var_names))
  {
    auto test_var_name = _test_var_names.at(i);
    _lfs.Register(test_var_name, std::make_shared<mfem::ParLinearForm>(_test_pfespaces.at(i)));
    _lfs.GetRef(test_var_name) = 0.0;
  }

  for (auto & test_var_name : _test_var_names)
  {
    // Apply kernels
    auto lf = _lfs.GetShared(test_var_name);
    ApplyDomainLFIntegrators(test_var_name, lf, _kernels_map);
    ApplyBoundaryLFIntegrators(test_var_name, lf, _integrated_bc_map);
    lf->Assemble();
  }

  // Apply essential boundary conditions
  ApplyEssentialBCs();

  // Eliminate trivially eliminated variables by subtracting contributions from linear forms
  EliminateCoupledVariables();
}

void
EquationSystem::BuildNonlinearForms()
{
  // Register non-linear Action forms
  for (const auto i : index_range(_test_var_names))
  {
    auto test_var_name = _test_var_names.at(i);
    _nlfs.Register(test_var_name, std::make_shared<mfem::ParNonlinearForm>(_test_pfespaces.at(i)));
    // Apply kernels
    auto nlf = _nlfs.GetShared(test_var_name);
    nlf->SetEssentialTrueDofs(_ess_tdof_lists.at(i));
    ApplyDomainNLFIntegrators(test_var_name, nlf, _kernels_map, std::nullopt);
    ApplyBoundaryNLFIntegrators(test_var_name, nlf, _integrated_bc_map, std::nullopt);
  }
}

void
EquationSystem::BuildBilinearForms()
{
  // Register bilinear forms
  for (const auto i : index_range(_test_var_names))
  {
    auto test_var_name = _test_var_names.at(i);
    _blfs.Register(test_var_name, std::make_shared<mfem::ParBilinearForm>(_test_pfespaces.at(i)));

    // Apply kernels
    auto blf = _blfs.GetShared(test_var_name);
    blf->SetAssemblyLevel(_assembly_level);
    ApplyBoundaryBLFIntegrators<mfem::ParBilinearForm>(
        test_var_name, test_var_name, blf, _integrated_bc_map);
    ApplyDomainBLFIntegrators<mfem::ParBilinearForm>(
        test_var_name, test_var_name, blf, _kernels_map);
    // Assemble
    blf->Assemble();
  }
}

void
EquationSystem::BuildMixedBilinearForms()
{
  // Register mixed bilinear forms. Note that not all combinations may
  // have a kernel.

  // Create mblf for each test/coupled variable pair with an added kernel.
  // Mixed bilinear forms with coupled variables that are not trial variables are
  // associated with contributions from eliminated variables.
  for (const auto i : index_range(_test_var_names))
  {
    auto test_var_name = _test_var_names.at(i);
    auto test_mblfs = std::make_shared<Moose::MFEM::NamedFieldsMap<mfem::ParMixedBilinearForm>>();
    for (const auto j : index_range(_coupled_var_names))
    {
      const auto & coupled_var_name = _coupled_var_names.at(j);
      auto mblf = std::make_shared<mfem::ParMixedBilinearForm>(_coupled_pfespaces.at(j),
                                                               _test_pfespaces.at(i));
      // Register MixedBilinearForm if kernels exist for it, and assemble kernels
      if (_kernels_map.Has(test_var_name) &&
          _kernels_map.Get(test_var_name)->Has(coupled_var_name) &&
          test_var_name != coupled_var_name)
      {
        mblf->SetAssemblyLevel(_assembly_level);
        // Apply all mixed kernels with this test/trial pair
        ApplyDomainBLFIntegrators<mfem::ParMixedBilinearForm>(
            coupled_var_name, test_var_name, mblf, _kernels_map);
        // Assemble mixed bilinear forms
        mblf->Assemble();
        // Register mixed bilinear forms associated with a single trial variable
        // for the current test variable
        test_mblfs->Register(coupled_var_name, mblf);
      }
    }
    // Register all mixed bilinear form sets associated with a single test variable
    _mblfs.Register(test_var_name, test_mblfs);
  }
}

void
EquationSystem::BuildEquationSystem()
{
  BuildBilinearForms();
  BuildMixedBilinearForms();
  BuildLinearForms();
  BuildNonlinearForms();
}

void
EquationSystem::ApplyDomainLFIntegrators(
    const std::string & test_var_name,
    std::shared_ptr<mfem::ParLinearForm> form,
    NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map)
{
  if (kernels_map.Has(test_var_name) && kernels_map.Get(test_var_name)->Has(test_var_name))
  {
    auto kernels = kernels_map.GetRef(test_var_name).GetRef(test_var_name);
    for (auto & kernel : kernels)
    {
      mfem::LinearFormIntegrator * integ = kernel->createLFIntegrator();

      if (integ)
      {
        kernel->isSubdomainRestricted()
            ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
            : form->AddDomainIntegrator(std::move(integ));
      }
    }
  }
}

void
EquationSystem::ApplyDomainNLFIntegrators(
    const std::string & test_var_name,
    std::shared_ptr<mfem::ParNonlinearForm> form,
    NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMKernel>>>> & kernels_map,
    std::optional<mfem::real_t> scale_factor)
{
  if (kernels_map.Has(test_var_name))
    for (const auto & [trial_var_name, kernels] : kernels_map.GetRef(test_var_name))
      for (auto & kernel : *kernels)
        if (auto * integ = kernel->createNLIntegrator())
        {
          if (_solver_requires_gradient && (trial_var_name != test_var_name))
            mooseError("Support for off-diagonal MFEM nonlinear domain integrators in conjunction "
                       "with a nonlinear solver that requires a gradient is not currently "
                       "implemented. Kernel '",
                       kernel->name(),
                       "' contributes to test variable '",
                       test_var_name,
                       "' from trial variable '",
                       trial_var_name,
                       "'.");

          _non_linear = true;
          if (scale_factor.has_value())
            integ = new NLScaleIntegrator(integ, scale_factor.value(), true);
          kernel->isSubdomainRestricted()
              ? form->AddDomainIntegrator(std::move(integ), kernel->getSubdomainMarkers())
              : form->AddDomainIntegrator(std::move(integ));
        }
}

void
EquationSystem::ApplyBoundaryLFIntegrators(
    const std::string & test_var_name,
    std::shared_ptr<mfem::ParLinearForm> form,
    NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
        integrated_bc_map)
{
  if (integrated_bc_map.Has(test_var_name) &&
      integrated_bc_map.Get(test_var_name)->Has(test_var_name))
  {
    auto bcs = integrated_bc_map.GetRef(test_var_name).GetRef(test_var_name);
    for (auto & bc : bcs)
    {
      mfem::LinearFormIntegrator * integ = bc->createLFIntegrator();

      if (integ)
      {
        bc->isBoundaryRestricted()
            ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
            : form->AddBoundaryIntegrator(std::move(integ));
      }
    }
  }
}

void
EquationSystem::ApplyBoundaryNLFIntegrators(
    const std::string & test_var_name,
    std::shared_ptr<mfem::ParNonlinearForm> form,
    NamedFieldsMap<NamedFieldsMap<std::vector<std::shared_ptr<MFEMIntegratedBC>>>> &
        integrated_bc_map,
    std::optional<mfem::real_t> scale_factor)
{
  if (integrated_bc_map.Has(test_var_name))
    for (const auto & [trial_var_name, bcs] : integrated_bc_map.GetRef(test_var_name))
      for (auto & bc : *bcs)
        if (auto * integ = bc->createNLIntegrator())
        {
          if (_solver_requires_gradient && (test_var_name != trial_var_name))
            mooseError(
                "Support for Off-diagonal MFEM nonlinear boundary integrators in conjunction with "
                "a nonlinear solver that requires a gradient is not currently "
                "implemented. Boundary condition '",
                bc->name(),
                "' contributes to test variable '",
                test_var_name,
                "' from trial variable '",
                trial_var_name,
                "'.");

          _non_linear = true;
          if (scale_factor.has_value())
            integ = new NLScaleIntegrator(integ, scale_factor.value(), true);
          bc->isBoundaryRestricted()
              ? form->AddBoundaryIntegrator(std::move(integ), bc->getBoundaryMarkers())
              : form->AddBoundaryIntegrator(std::move(integ));
        }
}

void
EquationSystem::PrepareLinearSolver(LinearSolverBase & solver)
{
  if (solver.IsLOR())
  {
    if (Complex())
      mooseError("LOR solve is not supported for complex equation systems.");
    if (_test_var_names.size() > 1)
      mooseError("LOR solve is only supported for single-variable systems");
    solver.SetupLOR(*_blfs.Get(_test_var_names.at(0)), _ess_tdof_lists.at(0));
  }

  mooseAssert(_linear_operator.Ptr(),
              "If we are preparing a linear solver, we better have a linear operator");
  solver.SetOperator(_linear_operator);
}

} // namespace Moose::MFEM

#endif