parsed_function.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
#include <kernel/analytic/function.hpp>
#include <kernel/util/exception.hpp>
 
// The contents of this file require the 'fparser' third-party library.
#if defined(FEAT_HAVE_FPARSER) || defined(DOXYGEN)
 
#include <fparser.hh>
 
#include <array>
#include <vector>
#include <map>
 
namespace FEAT
{
  namespace Analytic
  {
    class ParsedFunctionParseError : public ParseError
    {
    public:
      explicit ParsedFunctionParseError(const String& msg) : ParseError(msg) {}
    };
 
    class ParsedFunctionEvalError : public Exception
    {
    public:
      explicit ParsedFunctionEvalError(const String& msg) : Exception(msg) {}
    };
 
    template<int dim_>
    class ParsedScalarFunction :
      public Analytic::Function
    {
    public:
      static_assert((dim_ >= 1) && (dim_ <= 3), "unsupported function dimension");
 
      static constexpr int domain_dim = dim_;
 
      typedef Analytic::Image::Scalar ImageType;
 
      static constexpr bool can_value = true;
 
      static constexpr bool can_grad = false;
 
      static constexpr bool can_hess = false;
 
    private:
      ::FunctionParser _parser;
      std::map<String, double> _extra_vars;
 
    public:
      explicit ParsedScalarFunction() :
        _parser(),
        _extra_vars()
      {
        // add constants to our parser
        _parser.AddConstant("pi", Math::pi<double>());
        _parser.AddConstant("eps", Math::eps<double>());
      }
 
      explicit ParsedScalarFunction(const String& function) :
        ParsedScalarFunction()
      {
        parse(function);
      }
 
      void add_constant(const String& name, double value)
      {
        _parser.AddConstant(name.c_str(), value);
      }
 
      void add_variable(const String& name, double value = 0.0)
      {
        XASSERTM(_extra_vars.find(name) == _extra_vars.end(), "extra variable '" + name + "' already added");
        XASSERTM(name.compare_no_case("x") != 0, "variable name 'x' is reserved and cannot be used as extra variable");
        XASSERTM(name.compare_no_case("y") != 0, "variable name 'y' is reserved and cannot be used as extra variable");
        XASSERTM(name.compare_no_case("z") != 0, "variable name 'z' is reserved and cannot be used as extra variable");
        _extra_vars.emplace(name, value);
      }
 
      void set_variable(const String& name, double value)
      {
        auto it = _extra_vars.find(name);
        XASSERTM(it != _extra_vars.end(), "no variable named '" + name + "' found");
        it->second = value;
      }
 
      void parse(const String& function)
      {
          // add variables to our parser
        String vars("x");
        if(dim_ > 1) vars += ",y";
        if(dim_ > 2) vars += ",z";
 
        // add extra variables
        for(auto& v : _extra_vars)
          (vars += ",") += v.first;
 
        // try to parse the function
        const int ret = _parser.Parse(function.c_str(), vars.c_str());
        if(ret >= 0)
        {
          String msg(_parser.ErrorMsg());
          msg.append("\n>>> '");
          msg.append(function);
          msg.append("'");
          if(ret < int(function.size()))
          {
            // ret contains the index of the first invalid character in the input string
            // append an additional line to mark the faulty character
            msg.append("\n>>>");
            msg.append(String(std::size_t(ret+2), '-'));
            msg.append("^");
          }
          throw ParsedFunctionParseError(msg);
        }
 
        // optimize the parsed function
        _parser.Optimize();
      }
 
      template<typename Traits_>
      class Evaluator :
        public Analytic::Function::Evaluator<Traits_>
      {
      public:
        typedef typename Traits_::DataType DataType;
        typedef typename Traits_::PointType PointType;
        typedef typename Traits_::ValueType ValueType;
 
      private:
        // Note: We need to create a copy of the parser rather than just a reference to it,
        // as the 'Eval' function of the parser is not thread-safe, whereas our Evaluators
        // are always meant to be so!
        ::FunctionParser _parser;
        std::vector<double> _vars;
 
      public:
        explicit Evaluator(const ParsedScalarFunction& function) :
          _parser(function._parser),
          _vars(std::size_t(dim_), 0.0)
        {
          // force deep copy to avoid race conditions in case of multi-threading
          _parser.ForceDeepCopy();
          // set extra vars
          for(const auto& v : function._extra_vars)
            _vars.push_back(v.second);
        }
 
        ValueType value(const PointType& point)
        {
          // copy coordinates
          for(int i(0); i < dim_; ++i)
            _vars[std::size_t(i)] = double(point[i]);
 
          // evaluate the parser
          const double val = _parser.Eval(_vars.data());
 
          // check for errors
          switch(_parser.EvalError())
          {
          case 0: // no error
            break;
 
          case 1: // division by zero
            throw ParsedFunctionEvalError("Error in ParsedScalarFunction evaluation: division by zero");
 
          case 2: // sqrt error
          case 3: // log error
          case 4: // trigonometric error
            throw ParsedFunctionEvalError("Error in ParsedScalarFunction evaluation: illegal input value");
 
          case 5: // recursion error
            throw ParsedFunctionEvalError("Error in ParsedScalarFunction evaluation: maximum recursion depth reached");
 
          default: // ???
            throw ParsedFunctionEvalError("Error in ParsedScalarFunction evaluation: unknown error");
          }
 
          // return value
          return ValueType(val);
        }
      }; // class ParsedScalarFunction::Evaluator<...>
    }; // class ParsedScalarFunction
 
    template<int dim_>
    using ParsedFunction = ParsedScalarFunction<dim_>;
 
 
    template<int dom_dim_, int img_dim_ = dom_dim_>
    class ParsedVectorFunction :
      public Analytic::Function
    {
    public:
      static_assert((dom_dim_ >= 1) && (dom_dim_ <= 3), "unsupported domain dimension");
      static_assert(img_dim_ >= 1, "unsupported image dimension");
 
      static constexpr int domain_dim = dom_dim_;
 
      typedef Analytic::Image::Vector<img_dim_> ImageType;
 
      static constexpr bool can_value = true;
 
      static constexpr bool can_grad = false;
 
      static constexpr bool can_hess = false;
 
    private:
      std::map<String, double> _extra_vars;
 
      std::array<::FunctionParser, img_dim_> _parsers;
 
    public:
      explicit ParsedVectorFunction()
      {
        // add constants to our parsers
        add_constant("pi", Math::pi<double>());
        add_constant("eps", Math::eps<double>());
      }
 
      explicit ParsedVectorFunction(const String& function) :
        ParsedVectorFunction()
      {
        parse(function);
      }
 
      void add_constant(const String& name, double value)
      {
        for(auto& p : _parsers)
          p.AddConstant(name.c_str(), value);
      }
 
      void add_variable(const String& name, double value = 0.0)
      {
        XASSERTM(_extra_vars.find(name) == _extra_vars.end(), "extra variable '" + name + "' already added");
        XASSERTM(name.compare_no_case("x") != 0, "variable name 'x' is reserved and cannot be used as extra variable");
        XASSERTM(name.compare_no_case("y") != 0, "variable name 'y' is reserved and cannot be used as extra variable");
        XASSERTM(name.compare_no_case("z") != 0, "variable name 'z' is reserved and cannot be used as extra variable");
        _extra_vars.emplace(name, value);
      }
 
      void set_variable(const String& name, double value)
      {
        auto it = _extra_vars.find(name);
        XASSERTM(it != _extra_vars.end(), "no variable named '" + name + "' found");
        it->second = value;
      }
 
      void parse(const String& function)
      {
        // add variables to our parser
        String vars("x");
        if(dom_dim_ > 1) vars += ",y";
        if(dom_dim_ > 2) vars += ",z";
 
        // add extra variables
        for(auto& v : _extra_vars)
          (vars += ",") += v.first;
 
        String sfunc(function);
        if(sfunc.starts_with('['))
          sfunc.pop_front();
        if(sfunc.ends_with(']'))
          sfunc.pop_back();
        sfunc.trim_me();
 
        // split function string
        std::deque<String> sfuncs = sfunc.split_by_charset("'");
        if(sfuncs.size() != _parsers.size())
        {
          String msg("Invalid number of formulae: expected ");
          msg.append(stringify(_parsers.size()));
          msg.append(" but got ");
          msg.append(stringify(sfuncs.size()));
          throw ParsedFunctionParseError(msg);
        }
 
        // parse all functions
        for(std::size_t i(0); i < _parsers.size(); ++i)
        {
          // try to parse the function
          const int ret = _parsers.at(i).Parse(sfuncs.at(i).c_str(), vars.c_str());
          if(ret >= 0)
          {
            String msg(_parsers.at(i).ErrorMsg());
            msg.append("\n>>> '");
            msg.append(sfuncs.at(i));
            msg.append("'");
            if(ret < int(sfuncs.at(i).size()))
            {
              // ret contains the index of the first invalid character in the input string
              // append an additional line to mark the faulty character
              msg.append("\n>>>");
              msg.append(String(std::size_t(ret+2), '-'));
              msg.append("^");
            }
            throw ParsedFunctionParseError(msg);
          }
 
          // optimize the parsed function
          _parsers.at(i).Optimize();
        }
      }
 
      template<typename Traits_>
      class Evaluator :
        public Analytic::Function::Evaluator<Traits_>
      {
      public:
        typedef typename Traits_::DataType DataType;
        typedef typename Traits_::PointType PointType;
        typedef typename Traits_::ValueType ValueType;
 
      private:
        std::array< ::FunctionParser, img_dim_> _parsers;
        std::vector<double> _vars;
 
      public:
        explicit Evaluator(const ParsedVectorFunction& function) :
          _parsers(function._parsers),
          _vars(std::size_t(dom_dim_), 0.0)
        {
          // force deep copy to avoid race conditions in case of multi-threading
          for(auto& p : _parsers)
            p.ForceDeepCopy();
 
          // set extra vars
          for(const auto& v : function._extra_vars)
            _vars.push_back(v.second);
        }
 
        ValueType value(const PointType& point)
        {
          // copy coordinates
          for(int i(0); i < dom_dim_; ++i)
            _vars[std::size_t(i)] = double(point[i]);
 
          // evaluate the parser
          ValueType val;
          for(std::size_t i(0); i < _parsers.size(); ++i)
          {
            // evaluate the parser
            val[int(i)] = DataType(_parsers[i].Eval(_vars.data()));
 
            // check for errors
            switch(_parsers[i] .EvalError())
            {
            case 0: // no error
              break;
 
            case 1: // division by zero
              throw ParsedFunctionEvalError("Error in ParsedScalarFunction evaluation: division by zero");
 
            case 2: // sqrt error
            case 3: // log error
            case 4: // trigonometric error
              throw ParsedFunctionEvalError("Error in ParsedScalarFunction evaluation: illegal input value");
 
            case 5: // recursion error
              throw ParsedFunctionEvalError("Error in ParsedScalarFunction evaluation: maximum recursion depth reached");
 
            default: // ???
              throw ParsedFunctionEvalError("Error in ParsedScalarFunction evaluation: unknown error");
            }
          }
 
          return val;
        }
      }; // class ParsedVectorFunction::Evaluator<...>
    }; // class ParsedVectorFunction
  } // namespace Analytic
} // namespace FEAT
 
#endif // FEAT_HAVE_FPARSER