CoefficientMap.h

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
 
#pragma once
 
#include <map>
#include <string>
#include <tuple>
#include <utility>
#include <variant>
#include <vector>
 
#include "MooseException.h"
#include "MooseError.h"
 
#include "libmesh/utility.h"
 
namespace Moose::MFEM
{
 
template <class T, class Tpw>
class CoefficientMap
{
public:
  CoefficientMap() = default;
 
  template <class P, class... Args>
  std::shared_ptr<P> make(Args &&... args)
  {
    auto result = std::make_shared<P>(args...);
    this->_iterable_coefficients.push_back(result);
    return result;
  }
 
  void addCoefficient(const std::string & name, std::shared_ptr<T> coeff)
  {
    mooseAssert(std::find(this->_iterable_coefficients.cbegin(),
                          this->_iterable_coefficients.cend(),
                          coeff) != this->_iterable_coefficients.cend(),
                "Coefficient object was not created by this CoefficientMap.");
 
    const auto [_, inserted] = this->_coefficients.emplace(name, std::move(coeff));
    if (!inserted)
      mooseError("Coefficient with name '" + name + "' already present in CoefficientMap object");
  }
 
  void addPiecewiseBlocks(const std::string & name,
                          std::shared_ptr<T> coeff,
                          const std::vector<std::string> & blocks)
  {
    // If list of blocks is empty then treat as a global coefficient
    if (blocks.size() == 0)
    {
      this->addCoefficient(name, coeff);
      return;
    }
 
    // Initialise property with empty coefficients, if it does not already exist
    if (!this->hasCoefficient(name))
      this->_coefficients.insert({name, this->emptyPWData(coeff)});
    PWData * data = std::get_if<PWData>(&this->_coefficients[name]);
    // Throw an exception if the data is not piecewise
    if (!data)
      mooseError("Global coefficient with name '" + name + "' already present in CoefficientMap");
    mooseAssert(std::find(this->_iterable_coefficients.cbegin(),
                          this->_iterable_coefficients.cend(),
                          coeff) != this->_iterable_coefficients.cend(),
                "Coefficient object was not created by the appropriate coefficient manager.");
    auto & [pw_coeff, coeff_map] = *data;
    this->checkPWData(coeff, pw_coeff, name);
 
    for (const auto & block : blocks)
    {
      if (coeff_map.count(block) > 0)
        mooseError("Property with name '" + name + "' already assigned to block " + block +
                   " in CoefficientMap object");
      coeff_map[block] = coeff;
      pw_coeff->UpdateCoefficient(std::stoi(block), *coeff);
    }
  }
 
  T & getCoefficient(const std::string & name) { return *this->getCoefficientPtr(name); }
 
  std::shared_ptr<T> getCoefficientPtr(const std::string & name)
  {
    try
    {
      auto & coeff = this->_coefficients.at(name);
      try
      {
        return std::get<std::shared_ptr<T>>(coeff);
      }
      catch (const std::bad_variant_access &)
      {
        return std::get<0>(std::get<PWData>(coeff));
      }
    }
    catch (const std::out_of_range &)
    {
      mooseError("Property with name '" + name + "' has not been declared.");
    }
  }
 
  bool hasCoefficient(const std::string & name) const
  {
    return this->_coefficients.count(name) > 0;
  }
 
  bool propertyDefinedOnBlock(const std::string & name, const std::string & block) const
  {
    if (!this->hasCoefficient(name))
      return false;
    auto & coeff = libmesh_map_find(this->_coefficients, name);
    if (std::holds_alternative<std::shared_ptr<T>>(coeff))
      return true;
    auto block_map = std::get<1>(std::get<PWData>(coeff));
    return block_map.count(block) > 0;
  }
 
  void setTime(const mfem::real_t time)
  {
    for (auto & coef : this->_iterable_coefficients)
      coef->SetTime(time);
  }
 
  template <class F>
  void apply(F && func)
  {
    for (auto & coef : this->_iterable_coefficients)
      func(*coef);
  }
 
private:
  using PWData = std::tuple<std::shared_ptr<Tpw>, std::map<const std::string, std::shared_ptr<T>>>;
  std::map<const std::string, std::variant<std::shared_ptr<T>, PWData>> _coefficients;
  std::vector<std::shared_ptr<T>> _iterable_coefficients;
 
  PWData emptyPWData(std::shared_ptr<T> /*coeff*/)
  {
    return std::make_tuple(this->template make<Tpw>(),
                           std::map<const std::string, std::shared_ptr<T>>());
  }
  void checkPWData(std::shared_ptr<T> /*coeff*/,
                   std::shared_ptr<Tpw> /* existing_pw*/,
                   const std::string & /*name*/)
  {
  }
};
 
using ScalarMap = CoefficientMap<mfem::Coefficient, mfem::PWCoefficient>;
using VectorMap = CoefficientMap<mfem::VectorCoefficient, mfem::PWVectorCoefficient>;
using MatrixMap = CoefficientMap<mfem::MatrixCoefficient, mfem::PWMatrixCoefficient>;
 
template <>
VectorMap::PWData VectorMap::emptyPWData(std::shared_ptr<mfem::VectorCoefficient> coeff);
 
template <>
MatrixMap::PWData MatrixMap::emptyPWData(std::shared_ptr<mfem::MatrixCoefficient> coeff);
 
template <>
void VectorMap::checkPWData(std::shared_ptr<mfem::VectorCoefficient> coeff,
                            std::shared_ptr<mfem::PWVectorCoefficient> existing_pw,
                            const std::string & name);
 
template <>
void MatrixMap::checkPWData(std::shared_ptr<mfem::MatrixCoefficient> coeff,
                            std::shared_ptr<mfem::PWMatrixCoefficient> existing_pw,
                            const std::string & name);
} // namespace Moose::MFEM
 
#endif