BoundaryConditionManager.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/BoundaryConditionManager.hpp"
 
#include "geometry/selectors/BoundarySelector.hpp"
#include "geometry/selectors/AttributesSelector.hpp"
#include "geometry/selectors/CustomSelector.hpp"
 
BoundaryConditionManager::BoundaryConditionManager(FunctionDictionary* dico) :
    CepsObject(),
    ceps::HoldsDimension(),
    m_dofs           (nullptr),
    m_rowsAlreadyUsed({}),
    m_dictionary     (dico)
{
  m_dirichlet = ceps::getNew<Manager>();
  m_neumann   = ceps::getNew<Manager>();
  m_robin     = ceps::getNew<Manager>();
}
 
BoundaryConditionManager::~BoundaryConditionManager()
{
  this->reset();
}
 
void 
BoundaryConditionManager::setDofsInfos(DistributedInfos<DegreeOfFreedom*>* dofs)
{
  m_dofs = dofs;
  return;
}
 
void
BoundaryConditionManager::add(const CepsString& params)
{
  auto       kr   = ceps::splitAtNthDelimiter(params,1);
  CepsString key  = ceps::toKey(kr[0]);
  CepsString rest = kr[1];
 
  CEPS_ABORT_IF(rest.empty(),
    "When adding boundary condition, no parameters found after key " << key
  );
 
  add(key,rest);
  return;
}
 
void
BoundaryConditionManager::add(const CepsString& key, const CepsString& params)
{
 
  auto       kr   = ceps::splitAtNthDelimiter(params,1);
  CepsString type = ceps::toKey(kr[0]);
  CepsString rest = kr[1];
 
  CEPS_ABORT_IF(rest.empty(),
    "When adding boundary condition, no parameters found after keyword " << type
  );
 
  if (type == "DIRICHLET")
    add(key,CepsBoundaryConditionFlag::Dirichlet,rest);
  else if (type == "NEUMANN")
    add(key,CepsBoundaryConditionFlag::Neumann  ,rest);
  else if (type == "ROBIN")
    add(key,CepsBoundaryConditionFlag::Robin    ,rest);
  else
  {
    CEPS_ABORT(
      "When adding boundary condition with key " << key << ", unknown type " << type << "\n" <<
      "   Valid types are DIRICHLET, NEUMANN and ROBIN"
    );
  }
  return;
}
 
void
BoundaryConditionManager::add(
  const CepsString&                key,
  const CepsBoundaryConditionFlag& type,
  const CepsString&                params
)
{
 
  auto       kr      = ceps::splitAtNthDelimiter(params,1);
  CepsString dicoKey = ceps::toKey(kr[0]);
  CepsString rest    = kr[1];
 
  CEPS_ABORT_IF(rest.empty(),
    "When adding boundary condition, no parameters found after dictionary entry " << dicoKey
  );
 
  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));
 
 
  return add(key,type,rest,func);
}
 
void
BoundaryConditionManager::add(
  const CepsString&                key,
  const CepsBoundaryConditionFlag& type,
  const CepsString&                params,
  ScalarSAFunc*                    functor
)
{
  CEPS_ABORT_IF((type==CepsBoundaryConditionFlag::Neumann or
                 type==CepsBoundaryConditionFlag::Robin)
             and (not (functor->hasOption(CepsFunctionFlag::Space)     or
                       functor->hasOption(CepsFunctionFlag::CellIndex) or
                       functor->hasOption(CepsFunctionFlag::Attribute) or
                       functor->hasOption(CepsFunctionFlag::None))
                   or  functor->hasOption(CepsFunctionFlag::Solution)),
     "Function passed for Neumann/Robin BC must be defined either as\n"                    <<
     "     - a constant function (in space, may be time dependant)\n"                      <<
     "     - a function of any point x,\n"                                                 <<
     "     - a piecewise function using region attributes to determine the value,\n"       <<
     "     - a piecewise function using cell ID to determine the value,\n"                 <<
     "   This function may be evaluated at points of the boundary which do not coincide\n" <<
     "   with degrees of freedom or geometrical nodes"
  );
 
  // Set base BC internal multipliers to activate or deactivate terms in alpha u + beta dnu = g
  if (not m_dictionary->hasScalar("0.0"))  
    m_dictionary->addConstant("0.0",0.0);
  if (not m_dictionary->hasScalar("1.0"))  
    m_dictionary->addConstant("1.0",1.0);
 
  // -------------------------------------------------- //
  // Read parameters. currently supported keywords :
  //  - UNKNOWN
  //  - ATTRIBUTE(S)
  //  - SCALE
  //  - VOLUMIC    [only for dirichlet...]
  //  - ROBINCOEFF   [for Robin BC...]
  // -------------------------------------------------- //
 
  const CepsString &UNKNOWN    = "UNKNOWN";
  const CepsString &ATTRIBUTE  = "ATTRIBUTE";
  const CepsString &ATTRIBUTES = "ATTRIBUTES";
  const CepsString &SCALE      = "SCALE";
  const CepsString &VOLUMIC    = "VOLUMIC";
  const CepsString &ROBINCOEFF = "ROBINCOEFF";
 
  auto options = ceps::splitByKeyword(params,{UNKNOWN,ATTRIBUTE,ATTRIBUTES,SCALE,VOLUMIC,ROBINCOEFF});
 
  // Unknown
  CepsUnknownIndex unknown = 0u;
  if (options.contains(UNKNOWN))
  {
    CEPS_ABORT_IF(options[UNKNOWN].empty(),
      "When adding boundary condition, missing unknown identifier after "<< UNKNOWN << " keyword"
    );
    unknown = CepsUnknownIndex(ceps::toInt(options[UNKNOWN][0]));
  }
 
  // Attributes
  CepsSet<CepsUInt> attributes = {};
 
  for (auto label : {ATTRIBUTE, ATTRIBUTES})
    if (options.contains(label))
    {
      CEPS_ABORT_IF(options[label].empty(),
        "When adding boundary condition, 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 boundary condition, could not find real number after SCALE keyword"
    );
    scale = ceps::toReal(options[SCALE].at(0));
  }
 
  // Robin coefficient beta in dnu + beta u = g (alpha will be 1)
  typename FunctionDictionary::ScalarEntry beta = m_dictionary->getScalar("1.0");
  if (options.contains(ROBINCOEFF))
  {
    CepsString opt = options[ROBINCOEFF][0];
    CEPS_ABORT_IF(not m_dictionary->hasScalar(opt) and not ceps::isNumber(opt),
      "When adding boundary condition, unable to find "<< std::quoted(opt) << " in dictionary\n"
      "   or convert it to a real constant"
    );
    if (m_dictionary->hasScalar(opt))
      beta = m_dictionary->getScalar(opt);
    else
    {
      m_dictionary->addConstant("robin_coeff_"+key,ceps::toReal(opt));
      beta = m_dictionary->getScalar("robin_coeff_"+key);
    }
  }
 
  // "Volumic" bc for "dirichlet"
  CepsBool volumic = options.contains(VOLUMIC);
 
  // Now create the boundary condition
  ScalarBoundaryCondition *bc = ceps::getNew<ScalarBoundaryCondition>(functor,&(m_dofs->getOwned()),unknown);
  bc->setBoundaryConditionType(type);
  bc->setScaleFactor          (scale);
  
  // needed during the loop on elements, support contains halo
  if (type != CepsBoundaryConditionFlag::Dirichlet)
    bc->addDomain(&(m_dofs->getHalo()));
 
  // Compute support of boundary condition, mandatory for BCs
  // First, boundary filter
  if (not volumic)
  {
    BoundarySelector<CepsVector<DegreeOfFreedom*>::const_iterator> selecb;
    selecb.onlyOnThisProc(type == CepsBoundaryConditionFlag::Dirichlet);
    bc->computeSupport(&selecb);
  }
 
  // Attributes ?
  if (not attributes.empty())
  {
    AttributesSelector<CepsVector<DegreeOfFreedom*>::const_iterator> selec;
    selec.setAttributes(attributes.begin(),attributes.end());
    selec.onlyOnThisProc(type == CepsBoundaryConditionFlag::Dirichlet);
    selec.onlyOnBoundary(not volumic); 
    bc->computeSupport(&selec,not volumic);
    bc->setAttributes(attributes.begin(), attributes.end());
  }
 
  // Filter by unknown
  std::function<CepsBool(DegreeOfFreedom*)> unknownSelector = [&] (DegreeOfFreedom* d) -> CepsBool {
    return d->getUnknown()->getIdentifier() == unknown;
  };
  CustomSelector<CepsVector<DegreeOfFreedom*>::const_iterator> selec2(&unknownSelector);
  selec2.onlyOnThisProc(type == CepsBoundaryConditionFlag::Dirichlet);
  selec2.onlyOnBoundary(not volumic); 
  bc->computeSupport(&selec2,not volumic);
 
  // Store dofIds that are added, to avoid multiple BCs on a single dof
  if (type == CepsBoundaryConditionFlag::Dirichlet)
    for (DegreeOfFreedom* dof: bc->getSupport())
    {
      auto [_,added] = m_rowsAlreadyUsed.insert(dof);
      CEPS_ABORT_IF(not added,
        "When adding boundary condition " << key << ", one of the degrees of freedom\n" <<
        "   selected from the given options already has a boundary condition assigned to it" 
      );
    }
 
  typename FunctionDictionary::ScalarEntry Zero = m_dictionary->getScalar("0.0");
  typename FunctionDictionary::ScalarEntry One  = m_dictionary->getScalar("1.0");
 
  switch (type)
  {
    case CepsBoundaryConditionFlag::Dirichlet:
      bc->setAlpha(Zero);
      bc->setBeta (One);
      break;
    case CepsBoundaryConditionFlag::Neumann:
      bc->setAlpha(One);
      bc->setBeta (Zero);
      break;
    case CepsBoundaryConditionFlag::Robin:
      bc->setAlpha(One);
      bc->setBeta (beta);
  }
 
  add(key,bc);
  ceps::barrier();
  return;
}
 
void
BoundaryConditionManager::add(const CepsString &key, ScalarBoundaryCondition *bc)
{
  CepsString nkey    = ceps::toKey(key);
  Manager   *manager = getBCs(bc->getBoundaryConditionType());
 
  if (manager->contains(nkey))
  {
    CEPS_WARNS("already contains a boundary condition with the key " << nkey << ". Ignoring.");
    ceps::destroyObject(bc);
    return;
  }
 
  manager->insert(std::make_pair(nkey,bc));
 
  CepsString bctype = "not set";
  switch (bc->getBoundaryConditionType())
  { 
    case CepsBoundaryConditionFlag::Dirichlet: 
      bctype = "Dirichlet"; 
      break;
    case CepsBoundaryConditionFlag::Neumann: 
      bctype = "Neumann"; 
      break;
    case CepsBoundaryConditionFlag::Robin: 
      bctype = "Robin"; 
      break;
  }
                       
  CEPS_SAYS ("  Added new boundary condition : " << key << " (type: " + bctype + ")");
 
  return;
 
}
 
ScalarBoundaryCondition*
BoundaryConditionManager::getBC(
  CepsString                       key,
  const CepsBoundaryConditionFlag &type
) const
{
  return getBCs(type)->at(ceps::toKey(key));
}
 
BoundaryConditionManager::Manager*
BoundaryConditionManager::getBCs(const CepsBoundaryConditionFlag &type) const
{
  if (type == CepsBoundaryConditionFlag::Dirichlet)
    return m_dirichlet;
  else if (type == CepsBoundaryConditionFlag::Neumann)
    return m_neumann;
  return m_robin;
}
 
BoundaryConditionManager::Manager*
BoundaryConditionManager::getDirichletBCs() const
{
  return m_dirichlet;
}
 
BoundaryConditionManager::Manager*
BoundaryConditionManager::getNeumannBCs() const
{
  return m_neumann;
}
 
BoundaryConditionManager::Manager*
BoundaryConditionManager::getRobinBCs() const
{
  return m_robin;
}
 
void
BoundaryConditionManager::reset()
{
  for (auto toParse : {m_dirichlet,m_neumann,m_robin})
  {
    for (auto [name,bc]: *toParse)
      ceps::destroyObject(bc);
    ceps::destroyObject(toParse);
  }
  return;
}
 
void
BoundaryConditionManager::actualizeAll(CepsReal time)
{
  actualizeAllDirichlet(time);
  actualizeAllNeumann  (time);
  actualizeAllRobin    (time);
 
  return;
}
 
void
BoundaryConditionManager::actualizeAllDirichlet(CepsReal time)
{
  for (Manager::value_type obj : *m_dirichlet)
    obj.second->actualize(time);
  return;
}
 
void
BoundaryConditionManager::actualizeAllNeumann(CepsReal time)
{
  for (Manager::value_type obj : *m_neumann)
    obj.second->actualize(time);
 
  return;
}
 
void
BoundaryConditionManager::actualizeAllRobin(CepsReal time)
{
  for (Manager::value_type obj : *m_robin)
    obj.second->actualize(time);
 
  return;
}
 
// std::ostream &
// operator<< (std::ostream &os, const BoundaryConditionManager &manager)
// {
//   if (manager.m_dirichlet->size () > 0)
//     os << "--> Dirichlet " << *manager.m_dirichlet << std::endl;
 
//   for (CepsUInt i = 0; i < 3; ++i)
//     if (manager.m_neumann[i]->size () > 0)
//       os << "--> Neumann " + std::to_string (i) + "D " << *manager.m_neumann[i] << std::endl;
 
//   for (CepsUInt i = 0; i < 3; ++i)
//     if (manager.m_robin[i]->size () > 0)
//       os << "--> Robin " + std::to_string (i) + "D " << *manager.m_robin[i] << std::endl;
 
//   return os;
// }