FunctionDictionary.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/FunctionDictionary.hpp"
#include "pde/common/readers/CoeffReader.hpp"
 
FunctionDictionary::FunctionDictionary()
{
  addEntry("dicofunc_time_selector",new DicoFuncTimeSelector());
}
 
FunctionDictionary::~FunctionDictionary()
{
  for (auto& e: m_scalars)
    if (not m_doNotDelete.contains(e.first))
      ceps::destroyObject(e.second);
  for (auto& e: m_vectors)
    if (not m_doNotDelete.contains(e.first))
      ceps::destroyObject(e.second);
  for (auto& e: m_tensors)
    if (not m_doNotDelete.contains(e.first))
      ceps::destroyObject(e.second);
}
 
void
FunctionDictionary::add(const CepsString &label, const CepsString &options, Geometry* geom)
{
 
  CepsVector<CepsString> splitted = ceps::splitAtNthDelimiter(options,1);
  CepsString fieldType = ceps::toUpper(splitted[0]);
  CepsString params    = splitted[1];
 
  if (fieldType == "CONSTANT")
    return addConstant(label,params);
  if (fieldType == "FILE")
    return addFromFile(label,params,geom);
  if (fieldType == "FUNCTION")
    return addFunction(label,params);
  if (fieldType == "CPIECEWISE")
    return addCstPiecewise(label,params);
  if (fieldType == "FPIECEWISE")
    return addFctPiecewise(label,params);
  if (fieldType == "OPERATOR")
    return addOperator(label,params);
  if (fieldType == "ALIAS")
    return addAlias(label,params);
 
  CEPS_ABORT(
    "When adding object with label " << label << std::endl <<
    "   unrecognized type: " << fieldType << std::endl << 
    "   Valid types are : CONSTANT, FILE, FUNCTION, CPIECEWISE, FPIECEWISE, OPERATOR, ALIAS.\n" <<
    m_address
  );
}
 
void
FunctionDictionary::addConstant(const CepsString &label, const CepsString &params)
{
  const CepsVector<CepsString> &v = ceps::split(params);
  switch (v.size())
  {
    case 1: 
      return addConstant(label,ceps::toMathScalar(v));
    case 3: 
      return addConstant(label,ceps::toMathVertex(v));
    case 9: 
      return addConstant(label,ceps::toMathTensor(v));
    default: break;
  }
  CEPS_ABORT(
    "When adding constant field "<< label <<" to the dictionary, incorrect number of values provided\n" <<
    "   Either 1, 3 or 9 values must be given, for scalars, vectors or tensors, respectively.\n" <<
    m_address
  );
}
 
void
FunctionDictionary::addCstPiecewise(const CepsString &label, const CepsString &params)
{
  CepsVector<CepsString> attrValues = ceps::split(params,",");
 
  // Detect the dimensionality of the field
  CepsVector<CepsString> av = ceps::split(attrValues[0]);
  CEPS_ABORT_IF(av.size()!=2 and av.size()!=4 and av.size()!=10,
    "When initializing piecewise constant field " << label << ",\n" <<
    "   wrong number of words in description <attribute> <value>.\n" << m_address
  );
  CepsUInt dim = av.size()-1;  
 
  if (dim==1)
  {
    CepsMap<CepsAttribute,CepsMathScalar> map;
    for (CepsString attrValue : attrValues)
    {
      CepsVector<CepsString> avv = ceps::splitUntilNthDelimiter(attrValue,1);
      CepsVector<CepsString> val = ceps::split(avv[1]);
      map.insert(std::make_pair(ceps::toInt(avv[0]),ceps::toMathScalar(val)));
    }
    return addEntry(label,newCstPiecewiseSAFunc(map));
  }
  else if (dim==3)
  {
    CepsMap<CepsAttribute,CepsMathVertex> map;
    for (CepsString attrValue : attrValues)
    {
      CepsVector<CepsString> avv = ceps::splitUntilNthDelimiter(attrValue,1);
      CepsVector<CepsString> val = ceps::split(avv[1]);
      map.insert(std::make_pair(ceps::toInt(avv[0]),ceps::toMathVertex(val)));
    }
    return addEntry(label,newCstPiecewiseSAFunc(map));
  }
  // dim==9
   
    CepsMap<CepsAttribute,CepsMathTensor> map;
    for (CepsString attrValue : attrValues)
    {
      CepsVector<CepsString> avv = ceps::splitUntilNthDelimiter(attrValue,1);
      CepsVector<CepsString> val = ceps::split(avv[1]);
      map.insert(std::make_pair(ceps::toInt(avv[0]),ceps::toMathTensor(val)));
    }
    return addEntry(label,newCstPiecewiseSAFunc(map));
   
}
 
void
FunctionDictionary::addFctPiecewise(const CepsString &label, const CepsString &params)
{
  CepsVector<CepsString> attrValues = ceps::split(params,",");
  CepsMap<CepsAttribute,SAFunc<CepsMathScalar>*> maps;
  CepsMap<CepsAttribute,SAFunc<CepsMathVertex>*> mapv;
  CepsMap<CepsAttribute,SAFunc<CepsMathTensor>*> mapt;
  for (CepsString attrValue : attrValues)
  {
    CepsVector<CepsString> av = ceps::splitUntilNthDelimiter(attrValue,1);
    CEPS_ABORT_IF(av.size()<2,
      "When adding piecewise function " << label << " to dictionary, not enough parameters\n" << 
      "   given after attribute " << attrValue
    );
 
    // Check if function is already defined in dictionary, may be scalar, vector tensor
    ScalarEntry fptrs = getScalar(av[1],nullptr);
    VectorEntry fptrv = getVector(av[1],nullptr);
    TensorEntry fptrt = getTensor(av[1],nullptr);
    if (not ceps::isNullPtr(fptrs))
      maps.insert(std::make_pair(ceps::toInt(av[0]),ceps::runtimeCast<SAFunc<CepsMathScalar>*>(fptrs)));
    else if (not ceps::isNullPtr(fptrv))
      mapv.insert(std::make_pair(ceps::toInt(av[0]),ceps::runtimeCast<SAFunc<CepsMathVertex>*>(fptrv)));
    else if (not ceps::isNullPtr(fptrt))
      mapt.insert(std::make_pair(ceps::toInt(av[0]),ceps::runtimeCast<SAFunc<CepsMathTensor>*>(fptrt)));
    // If not, build from parameters, scalar only
    else
    {
      fptrs = parseFunction(av[1]);
      maps.insert(std::make_pair(ceps::toInt(av[0]),ceps::runtimeCast<SAFunc<CepsMathScalar>*>(fptrs)));
    }
  }
 
  if (not maps.empty())
    return addEntry(label,newFctPiecewiseSAFunc(maps));
  else if (not mapv.empty())
    return addEntry(label,newFctPiecewiseSAFunc(mapv));
  return addEntry(label,newFctPiecewiseSAFunc(mapt));
  
}
 
void
FunctionDictionary::addFromFile(
  const CepsString& label, 
  const CepsString& params, 
  Geometry*         geom,
  CepsBool          neverMind
)
{
  const CepsVector<CepsString> &v = ceps::split(params);
  CepsUInt nWords  = v.size();
 
  CEPS_ABORT_IF(nWords<1,
    "When initializing dictionary entry " << label << " with FILE labelword, \n" <<
    "   no file name was provided !"
  );
  // Check if file extension is compatible with coeff reader, if not, create interpolation
  CepsString ext = ceps::getExtension(v[0]);
 
  if (ext=="vtk" or ext=="msh")
  {
    CEPS_ABORT_IF(nWords<4,
      "When initializing dictionary entry " << label << " with FILE labelword, missing either\n" <<
      "   - fileName\n" <<
      "   - SCALAR, VECTOR, or TENSOR labelword\n" <<
      "   - name of data array in file\n" <<
      "   - CELL or NODE location labelword\n" <<
      m_address
    );
  
    CepsString       locWord = ceps::toKey(v[3]);
    CepsFunctionFlag loc     = locWord == "NODE" ? CepsFunctionFlag::NodeIndex :
                              (locWord == "CELL" ? CepsFunctionFlag::CellIndex :
                                                   CepsFunctionFlag::None);
    
    CEPS_ABORT_IF(loc==CepsFunctionFlag::None,
      "When initializing dictionary entry " << label << " with FILE labelword,\n" <<
      "   the last parameter should be either NODE or CELL.\n" <<
      m_address
    );
  
    CepsString arrayName = v[2];
    CepsString dim       = ceps::toKey(v[1]);
    
    if (neverMind)
    {
      CoeffReader reader(
        v[0],loc==CepsFunctionFlag::NodeIndex ? CepsLocationFlag::Point : CepsLocationFlag::Cell,geom
      );
      if (not reader.hasArray(arrayName))
        return;
    }
    
    if (dim=="SCALAR")
      return addEntry(label,new ArraySAFunc<CepsMathScalar>(v[0],arrayName,loc,geom));
    else if (dim=="VECTOR")
      return addEntry(label,new ArraySAFunc<CepsMathVertex>(v[0],arrayName,loc,geom));
    else if (dim=="TENSOR")
      return addEntry(label,new ArraySAFunc<CepsMathTensor>(v[0],arrayName,loc,geom));
    else  
      CEPS_ABORT(
        "When initializing dictionary entry " << label << " with FILE labelword,\n" <<
        "   the last parameter should be either NODE or CELL\n" <<
        m_address
      );
 
  }
  else
  {
    std::ifstream ifs(v[0]);
    CEPS_ABORT_IF(not ifs.good(),
      "Impossible to open file " << v[0] << " to create interpolation table"
    );
 
    CepsString line="";
    while (line.empty())
      std::getline(ifs,line);
    // Determine which type we are reading
    CepsUInt nElems = ceps::split(line).size()-1;
    CepsBool time   = v.size()>1 and ceps::toUpper(v[1])=="TIME";
    if (nElems>8)
      return addInterpolationTable<9,CepsMathTensor>(label,ifs,v[0],time);
    else if (nElems>2)
      return addInterpolationTable<3,CepsMathVertex>(label,ifs,v[0],time);
    return addInterpolationTable<1,CepsMathScalar>(label,ifs,v[0],time);
 
  }
 
  return;
 
}
 
 
void
FunctionDictionary::addFunction (const CepsString &label, const CepsString &params)
{
  return addEntry(label,parseFunction(ceps::toUpper(params)));
}
 
void
FunctionDictionary::addOperator (const CepsString &label, const CepsString &params)
{
  auto v = ceps::split (params);
  CEPS_ABORT_IF (v.size()!=3,
    "When adding function " << label << " to the dictionary with keyword OPERATOR,\n" <<
    "   operator definition requires three arguments (func1 operator func2)\n" <<
    m_address
  );
 
  const CepsChar &op        = v[1].at(0);   // operator
  CepsString      str   [2] = {ceps::toKey(v[0]), ceps::toKey(v[2])};  // left and right side of the operator
  CepsString      lr    [2] = { "left", "right"};
  ScalarEntry     scalar[2] = {nullptr, nullptr};
  VectorEntry     vector[2] = {nullptr, nullptr};
  TensorEntry     tensor[2] = {nullptr, nullptr};
 
  // Loop on left and right
  for (CepsUInt i = 0; i < 2; ++i)
  {
    // Check if operand is a real number
    // If it is the case, a constant entry is added to the dictionary
    if (not hasScalar(str[i]) and not hasVector(str[i]) and not hasTensor(str[i]))
    {
      CEPS_ABORT_IF(not ceps::isNumber(str[i]),
        "When adding function " << label << " to the dictionary with keyword OPERATOR,\n" <<
        "   " << lr[i] << "operand is neither an existing function of the dictionary\n" <<
        "   or a real constant. Got: " << str[i] << "\n" << m_address
      );
      addConstant(label + "_" + lr[i], ceps::toReal(str[i]));
      str[i] = label + "_" + lr[i];
    }
 
    scalar[i] = getScalar(str[i], nullptr);
    vector[i] = getVector(str[i], nullptr);
    tensor[i] = getTensor(str[i], nullptr);
  }
 
  auto tryToAdd = [&]<class _SharedPtr> (_SharedPtr obj) -> void {
    if (ceps::isValidPtr(scalar[1]))
      return addOperator(label,obj,op,scalar[1]);
    if (ceps::isValidPtr(vector[1]))
      return addOperator(label,obj,op,vector[1]);
    return addOperator(label,obj,op,tensor[1]);
  };
 
  if (ceps::isValidPtr(scalar[0]))  // left = scalar
    return tryToAdd(scalar[0]);
  if (ceps::isValidPtr(vector[0]))  // left = vector
    return tryToAdd(vector[0]);
  // if (ceps::isValidPtr(tensor[0]))  // left = tensor
  return tryToAdd(tensor[0]);
}
 
void
FunctionDictionary::addAlias(const CepsString& label, const CepsString& params)
{
  CepsString LABEL  = ceps::toKey(label);
  CepsString PARAMS = ceps::toKey(params);
  if (hasScalar(params))
    m_scalars[LABEL] = m_scalars[PARAMS];
  else if (hasVector(params))
    m_vectors[LABEL] = m_vectors[PARAMS];
  else if (hasTensor(params))
    m_tensors[LABEL] = m_tensors[PARAMS];
  else
    CEPS_ABORT(
      "Cannot create alias to function " << params << " that is not in the dictionary"
    );
 
  m_doNotDelete.insert(LABEL);
 
  return;
}
 
void
FunctionDictionary::addEntry(const CepsString &label, ScalarEntry func)
{
  return addInternal(label,func,m_scalars,"Scalar");
}
 
void
FunctionDictionary::addEntry(const CepsString &label, VectorEntry func)
{
  return addInternal(label,func,m_vectors,"Vector");
}
 
void
FunctionDictionary::addEntry(const CepsString &label, TensorEntry func)
{
  return addInternal(label,func,m_tensors,"Tensor");
}
 
// -------------------------------------------------------------------------
// Has, get
 
CepsBool
FunctionDictionary::hasScalar (const CepsString &label) const
{
  return m_scalars.contains(ceps::toKey(label));
}
 
CepsBool
FunctionDictionary::hasVector (const CepsString &label) const
{
  return m_vectors.contains(ceps::toKey(label));
}
 
CepsBool
FunctionDictionary::hasTensor (const CepsString &label) const
{
  return m_tensors.contains(ceps::toKey(label));
}
 
// -----------------------------------------------------------------
// Get functions
// -----------------------------------------------------------------
 
const typename FunctionDictionary::ScalarMap &
FunctionDictionary::getScalars() const
{
  return m_scalars;
}
 
const typename FunctionDictionary::VectorMap &
FunctionDictionary::getVectors() const
{
  return m_vectors;
}
 
const typename FunctionDictionary::TensorMap &
FunctionDictionary::getTensors() const
{
  return m_tensors;
}
 
const typename FunctionDictionary::ScalarEntry
FunctionDictionary::getScalar (const CepsString &label) const
{
  CepsString LABEL = ceps::toKey(label);
  CEPS_ABORT_IF(not hasScalar(LABEL),
    "FunctionDictionary does not have any scalar entry " << std::quoted(LABEL)
  );
  return m_scalars.at(LABEL);
}
 
const typename FunctionDictionary::VectorEntry
FunctionDictionary::getVector (const CepsString &label) const
{
  CepsString LABEL = ceps::toKey(label);
  CEPS_ABORT_IF(not hasVector(LABEL),
    "FunctionDictionary does not have any vector entry " << std::quoted(LABEL)
  );
  return m_vectors.at(LABEL);
}
 
const typename FunctionDictionary::TensorEntry
FunctionDictionary::getTensor (const CepsString &label) const
{
  CepsString LABEL = ceps::toKey(label);
  CEPS_ABORT_IF(not hasTensor(LABEL),
    "FunctionDictionary does not have any tensor entry " << std::quoted(LABEL)
  );
  return m_tensors.at(LABEL);
}
 
const typename FunctionDictionary::ScalarEntry
FunctionDictionary::getScalar (const CepsString &label, const ScalarEntry def) const
{
  CepsString LABEL = ceps::toKey(label);
  if (hasScalar(LABEL))
    return m_scalars.at(LABEL);
  return def;
}
 
const typename FunctionDictionary::VectorEntry
FunctionDictionary::getVector (const CepsString &label, const VectorEntry def) const
{
  CepsString LABEL = ceps::toKey(label);
  if (hasVector(LABEL))
    return m_vectors.at(LABEL);
  return def;
}
 
const typename FunctionDictionary::TensorEntry
FunctionDictionary::getTensor (const CepsString &label, const TensorEntry def) const
{
  CepsString LABEL = ceps::toKey(label);
  if (hasTensor(LABEL))
    return m_tensors.at(LABEL);
  return def;
}
 
// -----------------------------------------------------------------
// Delete functions
// -----------------------------------------------------------------
 
void
FunctionDictionary::deleteScalar(const CepsString& label)
{
  CepsString LABEL = ceps::toKey(label);
  ceps::destroyObject(getScalar(LABEL));
  m_scalars.erase(LABEL);
}
 
void
FunctionDictionary::deleteVector(const CepsString& label)
{
  CepsString LABEL = ceps::toKey(label);
  ceps::destroyObject(getVector(LABEL));
  m_vectors.erase(LABEL);
}
 
void
FunctionDictionary::deleteTensor(const CepsString& label)
{
  CepsString LABEL = ceps::toKey(label);
  ceps::destroyObject(getTensor(LABEL));
  m_tensors.erase(LABEL);
}
 
// -----------------------------------------------------------------
// Private function
// -----------------------------------------------------------------
 
 
typename FunctionDictionary::ScalarEntry
FunctionDictionary::parseFunction (const CepsString &params)
{
  const CepsString &TIME  = "TIME";
  const CepsString &SPACE = "SPACE";
  const CepsString &UFUNC = "UFUNC";
 
  CepsSet<CepsString>                        tags = {TIME, SPACE, UFUNC};
  CepsMap<CepsString,CepsVector<CepsString>> args = ceps::splitByKeyword(ceps::toUpper(params),tags);
  CepsMap<CepsString,CepsBool>               is   = {};
  CepsMap<CepsString,ScalarEntry>            fun  = {};
 
  // Detect which kind of function (TIME SPACE UFUNC) is being added
  // Use previously defined functions of the dictionary if available
  for (CepsString tag : tags)
  {
    is  [tag] = ceps::argVector(ceps::split(ceps::toUpper(params)),tag)!=-1;
    fun [tag] = nullptr;
 
    if (is[tag] and args[tag].size()>0)
      fun[tag] = getScalar(args[tag][0],nullptr);
  }
 
  CEPS_ABORT_IF(not is[TIME] and not is[SPACE] and not is[UFUNC],
    "When adding function to the dictionary, no keyword TIME, SPACE or UFUNC found"
  );
 
  // TIME : parse string arguments
  if (is[TIME] and ceps::isNullPtr(fun[TIME]))
    fun[TIME] = parseTimeFunction(args[TIME]);
 
  // SPACE : parse string arguments
  if (is[SPACE] and ceps::isNullPtr(fun[SPACE]))
    fun[SPACE] = parseSpatialFunction(args[SPACE]);
 
  // UFUNC : no parsing
 
  // In the case of 2 multiple valid functions, return the product
  // if only 1 function is begin given, return it.
  ScalarEntry res = nullptr;
  for (auto [tag1,obj1]: fun)
    if (ceps::isValidPtr(obj1))
    {
      for (auto [tag2,obj2] : fun)
        if (tag1 != tag2 and ceps::isValidPtr(obj2))
        {
          auto res = newProductSAFunc(obj1,obj2);
          res->setOwned(true,true);
          return res; 
        }
      res = obj1;
    }
  return res; // never reached
}
 
typename FunctionDictionary::ScalarEntry
FunctionDictionary::parseTimeFunction (const CepsVector<CepsString> &words)
{
  checkParam(words,{"START","END","DURATION","PERIOD","AMPLITUDE"});
  
  CepsString profile = ceps::toUpper(ceps::readParams (words, "PROFILE", "CONSTANT", 1));
  if (profile == "SRECT" or profile == "SMOOTH_RECTANGLE")
    profile += " " + ceps::readParams (words, "PROFILE", "CONSTANT", 2);
 
  return ceps::getNew<ExcitationTimeFunction>(
    profile,                                                     //
    ceps::toReal(ceps::readParams(words,"START"    ,    "0.0")), //
    ceps::toReal(ceps::readParams(words,"END"      , "5000.0")), //
    ceps::toReal(ceps::readParams(words,"DURATION" ,    "1.0")), //
    ceps::toReal(ceps::readParams(words,"PERIOD"   ,   "-1.0")), //
    ceps::toReal(ceps::readParams(words,"AMPLITUDE",    "1.0"))  //
  );
}
 
typename FunctionDictionary::ScalarEntry
FunctionDictionary::parseSpatialFunction (const CepsVector<CepsString> &words)
{
  checkParam(words,{"CENTER","DIAMETER","AMPLITUDE"},{3,1,1});
 
  CepsString profile = ceps::toUpper(ceps::readParams (words, "PROFILE", "CONSTANT", 1));
  if (profile == "SRECT" or profile == "SMOOTH_RECTANGLE")
    profile += " " + ceps::readParams (words, "PROFILE", "CONSTANT", 2);
 
  return ceps::getNew<ExcitationSpatialFunction>(
    profile,                                                        //
    CepsReal3D({                                                    //
      ceps::toReal (ceps::readParams (words, "CENTER", "0.0", 1)),  //
      ceps::toReal (ceps::readParams (words, "CENTER", "0.0", 2)),  //
      ceps::toReal (ceps::readParams (words, "CENTER", "0.0", 3))   //
    }),                                                             //
    ceps::toReal (ceps::readParams (words, "DIAMETER", "1.0")),     //
    ceps::toReal (ceps::readParams (words, "AMPLITUDE","1.0"))      //
  );
}
  
void
FunctionDictionary::checkParam(const CepsVector<CepsString>& words, 
                       const CepsVector<CepsString>& keys,
                       const CepsVector<CepsUInt>& nToRead) const
{
  for(CepsUInt i=0; i<keys.size(); i++)
  {
    CepsInt pos = ceps::argVector(words,keys[i]);
    if (pos!=-1)
    {
      CepsUInt n = nToRead.size()!=0 ? nToRead[i] : 1;
      CEPS_ABORT_IF(pos>=CepsInt(words.size()-n),
        "When adding function into the dictionary, missing parameter(s)\n" <<
        "   after keyword " << keys[i] << "\n" << m_address
      );
      for(CepsUInt j=0;j<n;j++)
        CEPS_ABORT_IF(not ceps::isNumber(words[pos+1+j]),
          "When adding function into the dictionary, wrong parameter\n" <<
          "   after keyword " << keys[i] << ", got: " << words[pos+1+j] << "\n" << m_address
        );
    }
  } 
}
 
CepsEnum
DicoFuncTimeSelector::getFlags() const
{
  return CepsEnum(CepsFunctionFlag::Time);
}
 
CepsReal
DicoFuncTimeSelector::eval(CepsStandardArgs args)
{
  return args.t;
}