ExprFunctions.h

Go to the documentation of this file.

// $Id$
//==============================================================================
//==============================================================================
 
#ifndef _EXPR_FUNCTIONS_H
#define _EXPR_FUNCTIONS_H
 
#include "Function.h"
#include "TensorFunction.h"
 
#include <array>
#include <algorithm>
#include <memory>
#include <string>
#include <vector>
 
namespace ExprEval {
  template<class ArgType> class Expression;
  template<class ArgType> class FunctionList;
  template<class ArgType> class ValueList;
}
 
 
template<class Scalar>
class EvalFuncScalar : public ScalarFunc
{
  using Expression = ExprEval::Expression<Scalar>;
  using FunctionList = ExprEval::FunctionList<Scalar>;
  using ValueList = ExprEval::ValueList<Scalar>;
  using FuncType = EvalFuncScalar<Scalar>;
 
  std::vector< std::unique_ptr<Expression> > expr;
  std::vector< std::unique_ptr<FunctionList> >  f;
  std::vector< std::unique_ptr<ValueList> >     v;
 
  std::vector<Scalar*> arg; 
 
  std::unique_ptr<FuncType> gradient; 
 
  Real dx; 
 
public:
  static int numError; 
 
  explicit EvalFuncScalar(const char* function, const char* x = "x",
                          Real eps = Real(1.0e-8));
  virtual ~EvalFuncScalar();
 
  void addDerivative(const std::string& function, const char* x = "x");
 
  bool isConstant() const override { return false; }
 
  Real deriv(Real x) const override;
 
protected:
  Real evaluate(const Real& x) const override;
};
 
 
template<class Scalar>
class EvalFuncSpatial : public RealFunc
{
  using Expression = ExprEval::Expression<Scalar>;
  using FunctionList = ExprEval::FunctionList<Scalar>;
  using ValueList = ExprEval::ValueList<Scalar>;
  using FuncType = EvalFuncSpatial<Scalar>;
 
  std::vector< std::unique_ptr<Expression> > expr;
  std::vector< std::unique_ptr<FunctionList> >  f;
  std::vector< std::unique_ptr<ValueList> >     v;
 
  struct Arg
  {
    Scalar* x; 
    Scalar* y; 
    Scalar* z; 
    Scalar* t; 
 
    const Scalar& get(int dir) const
    {
      switch (dir) {
        case  1: return *x;
        case  2: return *y;
        case  3: return *z;
        case  4: return *t;
        default: return *x;
      }
    }
  };
 
  std::vector<Arg> arg; 
 
  std::array<std::unique_ptr<FuncType>,4> derivative1;
  std::array<std::unique_ptr<FuncType>,6> derivative2;
 
  bool IAmConstant; 
 
  Real dx; 
  Real dt; 
 
public:
  explicit EvalFuncSpatial(const char* function,
                           Real epsX = Real(1.0e-8), Real epsT = Real(1.0e-12));
  virtual ~EvalFuncSpatial();
 
  void addDerivative(const std::string& function, const std::string& variables,
                     int d1, int d2 = 0);
 
  bool isConstant() const override { return IAmConstant; }
 
  Real deriv(const Vec3& X, int dir) const override;
  Real dderiv(const Vec3& X, int dir1, int dir2) const override;
 
  void setParam(const std::string& name, Real value) override;
 
  Vec3 gradient(const Vec3& X) const override
  {
    return this->RealFunc::gradient(X);
  }
 
  SymmTensor hessian(const Vec3& X) const override
  {
    return this->RealFunc::hessian(X);
  }
 
protected:
  Real evaluate(const Vec3& X) const override;
};
 
 
template<class Scalar>
class EvalFunctions
{
protected:
  using FuncType = EvalFuncSpatial<Scalar>; 
 
  EvalFunctions(const std::string& functions, const std::string& variables,
                const Real epsX, const Real epsT);
  virtual ~EvalFunctions();
 
public:
  void addDerivative(const std::string& functions,
                     const std::string& variables, int d1, int d2 = 0);
 
  size_t getNoSpaceDim() const { return nsd; }
 
protected:
  std::vector<std::unique_ptr<FuncType>> p; 
  size_t nsd = 0; 
};
 
 
template <class ParentFunc, class Ret, class Scalar>
class EvalMultiFunction : public ParentFunc, public EvalFunctions<Scalar>
{
  using FuncType = typename EvalFunctions<Scalar>::FuncType;
 
public:
  explicit EvalMultiFunction(const std::string& functions,
                             const std::string& variables = "",
                             const Real epsX = 1e-8,
                             const Real epsT = 1e-12)
    : EvalFunctions<Scalar>(functions,variables,epsX,epsT)
  {
    this->setNoDims();
  }
 
  virtual ~EvalMultiFunction() {}
 
  bool isConstant() const override
  {
    return std::all_of(this->p.begin(), this->p.end(),
                       [](const std::unique_ptr<FuncType>& func)
                       { return func->isConstant(); });
  }
 
  unsigned char getType() const override { return 2; }
 
  Ret deriv(const Vec3& X, int dir) const override;
  Ret dderiv(const Vec3& X, int dir1, int dir2) const override;
 
  void setParam(const std::string& name, Real value) override
  {
    for (std::unique_ptr<FuncType>& func : this->p)
      func->setParam(name,value);
  }
 
protected:
  void setNoDims();
 
  Ret evaluate(const Vec3& X) const override;
 
  std::vector<Real> evalGradient(const Vec3& X) const override;
 
  std::vector<Real> evalHessian(const Vec3& X) const override;
 
  std::vector<Real> evalTimeDerivative(const Vec3& X) const override;
};
 
using EvalFunc = EvalFuncScalar<Real>;
using EvalFunction = EvalFuncSpatial<Real>;
using VecFuncExpr = EvalMultiFunction<VecFunc,Vec3,Real>;
using TensorFuncExpr = EvalMultiFunction<TensorFunc,Tensor,Real>;
using STensorFuncExpr = EvalMultiFunction<STensorFunc,SymmTensor,Real>;
 
template<> int EvalFunc::numError;
 
#endif