SourceTermManager.cpp

Go to the documentation of this file.

/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
   This file is part of CEPS.
 
   CEPS is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.
 
   CEPS is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License
   along with CEPS (see file LICENSE at root of project).
   If not, see <https://www.gnu.org/licenses/>.
 
 
   Copyright 2019-2024 Inria, Universite de Bordeaux
 
   Authors, in alphabetical order:
 
     Pierre-Elliott BECUE, Florian CARO, Yves COUDIERE(*), Andjela DAVIDOVIC, 
     Charlie DOUANLA-LONTSI, Marc FUENTES, Mehdi JUHOOR, Michael LEGUEBE(*), 
     Pauline MIGERDITICHAN, Valentin PANNETIER(*), Nejib ZEMZEMI.
     * : currently active authors
 
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
#include "pde/common/functions/SourceTermManager.hpp"
#include "geometry/selectors/AttributesSelector.hpp"
#include "geometry/selectors/CustomSelector.hpp"
 
SourceTermManager::SourceTermManager(FunctionDictionary *dico) :
    CepsObject(),
    m_dofs(nullptr),
    m_manager(ceps::getNew<Manager>()),
    m_dictionary(dico)
{}
 
SourceTermManager::~SourceTermManager()
{
  for (auto [name,st]: *m_manager)
    ceps::destroyObject(st);
  ceps::destroyObject(m_manager);
}
 
void 
SourceTermManager::setDofsInfos(DistributedInfos<DegreeOfFreedom*>* dofs)
{
  m_dofs = dofs;
  return;
}
 
void
SourceTermManager::add(const CepsString &params,CepsBool computeSupport)
{
  // -------------------------------------------------- //
  // simple parse of params to deduce key
  // -------------------------------------------------- //
 
  CepsString::size_type pos  = params.find(" ");
  CepsString            key  = ceps::toKey(params.substr(0,pos));
  CepsString            rest = pos == std::string::npos ? "" : params.substr(pos+1);
 
  CEPS_ABORT_IF(rest.empty(),
    "When adding source term, no parameters found after key " << key
  );
 
  add(key,rest,computeSupport);
  return;
}
 
void
SourceTermManager::add(const CepsString &key, const CepsString &params, CepsBool computeSupport)
{
  // -------------------------------------------------- //
  // deduce the dico_key or the value
  // -------------------------------------------------- //
 
  CepsString::size_type pos     = params.find(" ");
  CepsString            dicokey = ceps::toKey(params.substr(0,pos));
  CepsString            rest    = pos == std::string::npos ? "" : params.substr(pos+1);
 
  CEPS_ABORT_IF(not ceps::isNumber(dicokey) and rest.empty(),
    "When adding source term, no parameters found after function name " << key << " " << dicokey
  );
 
  add(key,dicokey,rest,computeSupport);
  return;
}
 
void
SourceTermManager::add(const CepsString& key, CepsString& dicokey, 
                       const CepsString& params, CepsBool computeSupport)
{
  // -------------------------------------------------- //
  // extract dico_tag or value
  // -------------------------------------------------- //
 
  ScalarSAFunc* func;
 
  if (not m_dictionary->hasScalar(dicokey))
  {
    CepsString ndicokey = key + "_CONSTANT";
    m_dictionary->addConstant(ndicokey,dicokey);
    func = ceps::runtimeCast<ScalarSAFunc*>(m_dictionary->getScalar(ndicokey));
  }
  else
    func = ceps::runtimeCast<ScalarSAFunc*>(m_dictionary->getScalar(dicokey));
 
  add(key,params,func,computeSupport);
  return;
}
 
void
SourceTermManager::add(
  const CepsString& key,
  const CepsString& params,
  ScalarSAFunc*     functor,
  CepsBool          computeSupport
)
{
  // -------------------------------------------------- //
  // Read parameters. Currently supported keywords :
  //  - COMPONENT
  //  - ATTRIBUTE(S)
  //  - SCALE
  //  - TYPE [IAPP, KTERM, DEFAULT]
  //  - UNIT
  // -------------------------------------------------- //
 
  const CepsString& COMPONENT  = "UNKNOWN";
  const CepsString& ATTRIBUTE  = "ATTRIBUTE";
  const CepsString& ATTRIBUTES = "ATTRIBUTES";
  const CepsString& SCALE      = "SCALE";
  const CepsString& UNIT       = "UNIT";
  const CepsString& TYPE       = "TYPE";
 
  auto options = ceps::splitByKeyword(params,{COMPONENT,ATTRIBUTE,ATTRIBUTES,SCALE,UNIT,TYPE});
 
  // Unknown (previously called component)
  CepsUnknownIndex component = 0;
  if (options.contains(COMPONENT))
  {
    CEPS_ABORT_IF(options[COMPONENT].empty(),
      "When adding source term, missing unknown identifier after "<< COMPONENT << " keyword"
    );
    component = CepsUnknownIndex(ceps::toInt(options[COMPONENT][0]));
  }
    
  // Attributes
  CepsSet<CepsUInt> attributes = {};
 
  for (auto label : {ATTRIBUTE, ATTRIBUTES})
    if (options.contains(label))
    {
      CEPS_ABORT_IF(options[label].empty(),
        "When adding source term, could not find any attribute after " << label << " keyword"      
      );
      for (CepsString s : options[label])
        attributes.insert(CepsAttribute(ceps::toInt(s)));
    }
 
  if (attributes.empty())
    attributes.insert(CepsUniversal);
 
  // Scaling factor
  CepsReal scale = 1.;
  if (options.contains(SCALE))
  {
    CEPS_ABORT_IF(options[SCALE].empty() or not ceps::isNumber(options[SCALE][0]),
      "When adding source term, could not find real number after SCALE keyword"
    );
    scale = ceps::toReal(options[SCALE].at(0));
  }
 
  // Type of source term
  CepsSourceTermFlag flag = CepsSourceTermFlag::Default;
 
  if (options.contains(TYPE))
  {
    CEPS_ABORT_IF(options[TYPE].empty(),
      "When adding source term, could not find any entry after TYPE keyword"
    );
    if (options[TYPE].at(0) == "KTERM")
      flag = CepsSourceTermFlag::Laplace;
    else if (options[TYPE].at(0) == "IAPP")
      flag = CepsSourceTermFlag::Stimulation;
    else
    {
      CEPS_ABORT_IF(options[TYPE].at(0)!= "DEFAULT", CepsString(
        "When adding source term, invalid entry after TYPE keyword. Valid types are\n") +
        "   KTERM (for source terms to be mulitplied XXX FIXME XXX )\n" + 
        "   IAPP (volumic source term specific to cardiac PDEs)\n" + 
        "   DEFAULT (regular source term)"
      );
    }
  }
 
  // Now create the field 
  ScalarSourceTerm* st = ceps::getNew<ScalarSourceTerm>(functor,&(m_dofs->getOwned()),component);
  st->setScaleFactor(scale);
  st->setFlag       (flag);
 
  if (computeSupport)
  {
    // Compute support of source term
 
    // If there is a defined region
    if (not attributes.empty())
    {
      AttributesSelector<CepsVector<DegreeOfFreedom*>::const_iterator> selec;
      selec.setAttributes(attributes.begin(),attributes.end());
      selec.onlyOnThisProc(true);
      selec.onlyOnBoundary(false); 
      st->computeSupport(&selec);
    }
 
    // Then remove dofs for other unknowns
    std::function<CepsBool(DegreeOfFreedom*)> unknownSelector = [&] (DegreeOfFreedom* d) -> CepsBool {
      return d->getUnknown()->getIdentifier() == component;
    };
    CustomSelector<CepsVector<DegreeOfFreedom*>::const_iterator> selec2(&unknownSelector);
    selec2.onlyOnThisProc(true);
    selec2.onlyOnBoundary(false); 
    st->computeSupport(&selec2,true);
 
  }
 
  add(key,st);
  return;
}
 
void
SourceTermManager::add(const CepsString &key, ScalarSourceTerm *st)
{
  CepsString nkey = ceps::toKey(key);
 
  if (m_manager->contains(nkey))
  {
    CEPS_WARNS ("already contains a source term with the key " << nkey << ". Ignoring.");
    ceps::destroyObject(st);
    return;
  }
 
  m_manager->insert(std::make_pair(nkey,st));
 
  CEPS_SAYS ("  Added new source term : " << key);
  return;
}
 
void
SourceTermManager::actualizeAll(CepsReal time)
{
  for (auto obj : *m_manager)
    obj.second->actualize(time);
 
  return;
}
 
ScalarSourceTerm*
SourceTermManager::getSourceTerm(CepsString label) const
{
  CepsString KEY = ceps::toKey(label);
  CEPS_ABORT_IF(not m_manager->contains(KEY),
    "Source term manager has no source term with label " + label
  );
  return m_manager->at(KEY);
}
 
CepsVector<ScalarSourceTerm*>
SourceTermManager::asVector() const
{
  return ceps::valuesOf(*m_manager);
}
 
typename SourceTermManager::Manager* const
SourceTermManager::getManager() const
{
  return m_manager;
}
 
CepsUInt
SourceTermManager::getNbOfSourceTermsOf(CepsSourceTermFlag flag) const
{
  return getNbOfSourceTermsOf(ceps::toIntegral(flag));
}
 
 
CepsUInt
SourceTermManager::getNbOfSourceTermsOf(CepsEnum flag) const
{
  CepsUInt res = 0;
  for (auto it: *m_manager)
  {
    if (isEnabledOption(flag, it.second->getFlag()))
      res++;
  }
  return res;
}