Changed: ScalarFunc is now a sub-class of utl::Function<Real,Real>.

Added: Virtual method for calculation of derivative of scalar functions.
For expression functions the derivative can be evaluated through a separate
user-specified expression, or numerically by a finite difference approach.

src/Utility/ExprFunctions.C

@@ -83,10 +83,12 @@ static void ExprException (const ExprEval::Exception& exc, const char* task,
   EvalFunc::numError++;
 }
 
+
 int EvalFunc::numError = 0;
 
 
-EvalFunc::EvalFunc (const char* function, const char* x)
+EvalFunc::EvalFunc (const char* function, const char* x, Real eps)
+  : gradient(nullptr), dx(eps)
 {
   try {
     size_t nalloc = 1;
@@ -104,7 +106,7 @@ EvalFunc::EvalFunc (const char* function, const char* x)
       v[i] = new ExprEval::ValueList;
       f[i]->AddDefaultFunctions();
       v[i]->AddDefaultValues();
-      v[i]->Add(x,0,false);
+      v[i]->Add(x,0.0,false);
       expr[i]->SetFunctionList(f[i]);
       expr[i]->SetValueList(v[i]);
       expr[i]->Parse(function);
@@ -112,7 +114,7 @@ EvalFunc::EvalFunc (const char* function, const char* x)
     }
   }
   catch (ExprEval::Exception e) {
-    cleanup();
+    this->cleanup();
     ExprException(e,"parsing",function);
   }
 }
@@ -120,34 +122,42 @@ EvalFunc::EvalFunc (const char* function, const char* x)
 
 EvalFunc::~EvalFunc ()
 {
-  cleanup();
+  this->cleanup();
 }
 
 
-void EvalFunc::cleanup()
+void EvalFunc::cleanup ()
 {
-  for (auto& it : expr)
+  for (ExprEval::Expression* it : expr)
     delete it;
-  for (auto& it : f)
+  for (ExprEval::FunctionList* it : f)
     delete it;
-  for (auto& it : v)
+  for (ExprEval::ValueList* it : v)
     delete it;
-  arg.clear();
+  delete gradient;
   expr.clear();
   f.clear();
   v.clear();
+  arg.clear();
+}
+
+
+void EvalFunc::derivative (const std::string& function, const char* x)
+{
+  if (!gradient)
+    gradient = new EvalFunc(function.c_str(),x);
 }
 
 
 Real EvalFunc::evaluate (const Real& x) const
 {
+  Real result = Real(0);
   size_t i = 0;
 #ifdef USE_OPENMP
   i = omp_get_thread_num();
 #endif
   if (i >= arg.size())
-    return 0;
+    return result;
-  Real result = Real(0);
   try {
     *arg[i] = x;
     result = expr[i]->Evaluate();
@@ -160,6 +170,16 @@ Real EvalFunc::evaluate (const Real& x) const
 }
 
 
+Real EvalFunc::deriv (Real x) const
+{
+  if (gradient)
+    return gradient->evaluate(x);
+
+  // Evaluate derivative using central difference
+  return (this->evaluate(x+0.5*dx) - this->evaluate(x-0.5*dx)) / dx;
+}
+
+
 EvalFunction::EvalFunction (const char* function) : gradient{}, dgradient{}
 {
   try {
@@ -178,10 +198,10 @@ EvalFunction::EvalFunction (const char* function) : gradient{}, dgradient{}
       v[i] = new ExprEval::ValueList;
       f[i]->AddDefaultFunctions();
       v[i]->AddDefaultValues();
-      v[i]->Add("x",0,false);
+      v[i]->Add("x",0.0,false);
-      v[i]->Add("y",0,false);
+      v[i]->Add("y",0.0,false);
-      v[i]->Add("z",0,false);
+      v[i]->Add("z",0.0,false);
-      v[i]->Add("t",0,false);
+      v[i]->Add("t",0.0,false);
       expr[i]->SetFunctionList(f[i]);
       expr[i]->SetValueList(v[i]);
       expr[i]->Parse(function);
@@ -192,7 +212,7 @@ EvalFunction::EvalFunction (const char* function) : gradient{}, dgradient{}
     }
   }
   catch (ExprEval::Exception e) {
-    cleanup();
+    this->cleanup();
     ExprException(e,"parsing",function);
   }
 
@@ -205,26 +225,28 @@ EvalFunction::EvalFunction (const char* function) : gradient{}, dgradient{}
 
 EvalFunction::~EvalFunction ()
 {
-  cleanup();
+  this->cleanup();
 }
 
 
-void EvalFunction::cleanup()
+void EvalFunction::cleanup ()
 {
-  for (auto& it : expr)
+  for (ExprEval::Expression* it : expr)
     delete it;
-  for (auto& it : f)
+  for (ExprEval::FunctionList* it : f)
     delete it;
-  for (auto& it : v)
+  for (ExprEval::ValueList* it : v)
     delete it;
-  for (auto& it : gradient)
+  for (EvalFunction* it : gradient)
     delete it;
-  for (auto& it : dgradient)
+  for (EvalFunction* it : dgradient)
     delete it;
-  arg.clear();
   expr.clear();
   f.clear();
   v.clear();
+  arg.clear();
+  gradient.fill(nullptr);
+  dgradient.fill(nullptr);
 }
 
 
@@ -276,7 +298,7 @@ Real EvalFunction::evaluate (const Vec3& X) const
     i = omp_get_thread_num();
 #endif
     if (i >= arg.size())
-      return 0;
+      return result;
 
     *arg[i].x = X.x;
     *arg[i].y = X.y;
@@ -353,14 +375,14 @@ EvalFunctions::EvalFunctions (const std::string& functions,
                               const std::string& variables)
 {
   std::vector<std::string> components = splitComps(functions,variables);
-  for (const auto& comp : components)
+  for (const std::string& comp : components)
     p.push_back(new EvalFunction(comp.c_str()));
 }
 
 
 EvalFunctions::~EvalFunctions ()
 {
-  for (auto& it : p)
+  for (EvalFunction* it : p)
     delete it;
 }
 

src/Utility/ExprFunctions.h

@@ -42,23 +42,36 @@ class EvalFunc : public ScalarFunc
 
   std::vector<Real*> arg; //!< Function argument values
 
+  EvalFunc* gradient; //!< First derivative expression
+
+  Real dx; //!< Domain increment for calculation of numerical derivative
+
 public:
+  static int numError; //!< Error counter - set by the exception handler
+
   //! \brief The constructor parses the expression string.
-  EvalFunc(const char* function, const char* x = "x" );
+  EvalFunc(const char* function, const char* x = "x", Real eps = Real(1.0e-8));
   //! \brief The destructor frees the dynamically allocated objects.
   virtual ~EvalFunc();
 
-  static int numError; //!< Error counter - set by the exception handler
+  //! \brief Adds an expression function for a first derivative.
+  void derivative(const std::string& function, const char* x = "x");
+
+  //! \brief Returns whether the function is time-independent or not.
+  virtual bool isConstant() const { return false; }
+
+  //! \brief Returns the first-derivative of the function.
+  virtual Real deriv(Real x) const;
 
 protected:
   //! \brief Non-implemented copy constructor to disallow copying.
-  EvalFunc(const EvalFunc&);
+  EvalFunc(const EvalFunc&) = delete;
-  //! \brief Non-implemented assigment operator to disallow copying.
+  //! \brief Non-implemented assignment operator to disallow copying.
-  EvalFunc& operator=(const EvalFunc&);
+  EvalFunc& operator=(const EvalFunc&) = delete;
   //! \brief Evaluates the function expression.
   virtual Real evaluate(const Real& x) const;
 
-  //! \brief Cleanup allocated data.
+  //! \brief Cleans up the allocated data.
   void cleanup();
 };
 
@@ -108,14 +121,14 @@ public:
   virtual Real dderiv(const Vec3& X, int dir1, int dir2) const;
 
 protected:
+  //! \brief Non-implemented copy constructor to disallow copying.
+  EvalFunction(const EvalFunction&) = delete;
+  //! \brief Non-implemented assignment operator to disallow copying.
+  EvalFunction& operator=(const EvalFunction&) = delete;
   //! \brief Evaluates the function expression.
   virtual Real evaluate(const Vec3& X) const;
-  //! \brief Non-implemented copy constructor to disallow copying.
-  EvalFunction(const EvalFunction&);
-  //! \brief Non-implemented assignment operator to disallow copying.
-  EvalFunction& operator=(const EvalFunction&);
 
-  //! \brief Cleanup allocated data.
+  //! \brief Cleans up the allocated data.
   void cleanup();
 };
 

src/Utility/Function.h

@@ -117,8 +117,23 @@ namespace utl
 }
 
 
-//! \brief Scalar-valued unary function of a scalar value.
+/*!
-typedef utl::Function<Real,Real> ScalarFunc;
+  \brief Scalar-valued unary function of a scalar value.
+*/
+
+class ScalarFunc : public utl::Function<Real,Real>
+{
+protected:
+  //! \brief The constructor is protected to allow sub-class instances only.
+  ScalarFunc() {}
+
+public:
+  //! \brief Empty destructor.
+  virtual ~ScalarFunc() {}
+
+  //! \brief Returns the first-derivative of the function.
+  virtual Real deriv(Real x) const { return Real(0); }
+};
 
 
 /*!

src/Utility/Functions.C

@@ -34,6 +34,12 @@ Real RampFunc::evaluate (const Real& x) const
 }
 
 
+Real RampFunc::deriv (Real x) const
+{
+  return x < xmax ? fval/xmax : Real(0);
+}
+
+
 Real DiracFunc::evaluate (const Real& x) const
 {
   return fabs(x-xmax) < 1.0e-4 ? amp : Real(0);
@@ -52,6 +58,12 @@ Real SineFunc::evaluate (const Real& x) const
 }
 
 
+Real SineFunc::deriv (Real x) const
+{
+  return freq*scale*cos(freq*x+phase);
+}
+
+
 Real ConstTimeFunc::evaluate (const Vec3& X) const
 {
   const Vec4* Xt = dynamic_cast<const Vec4*>(&X);
@@ -467,14 +479,21 @@ const ScalarFunc* utl::parseTimeFunction (const char* type, char* cline, Real C)
 {
   if (strncasecmp(type,"expr",4) == 0 && cline != nullptr)
   {
+    cline = strtok(cline,":");
     IFEM::cout << cline;
     EvalFunc::numError = 0;
-    ScalarFunc* sf = new EvalFunc(cline,"t");
+    ScalarFunc* sf = new EvalFunc(cline,"t",C);
     if (EvalFunc::numError > 0)
     {
       delete sf;
       sf = nullptr;
     }
+    // The derivative can be specified as a second expression after the colon
+    if ((cline = strtok(nullptr,":")))
+    {
+      IFEM::cout <<" (derivative: "<< cline <<")";
+      static_cast<EvalFunc*>(sf)->derivative(cline,"t");
+    }
     return sf;
   }
   else if (strncasecmp(type,"Ramp",4) == 0 || strcmp(type,"Tinit") == 0)
@@ -519,7 +538,8 @@ const ScalarFunc* utl::parseTimeFunction (const char* type, char* cline, Real C)
 }
 
 
-ScalarFunc* utl::parseTimeFunc (const char* func, const std::string& type)
+ScalarFunc* utl::parseTimeFunc (const char* func, const std::string& type,
+                                Real eps)
 {
   char* cstr = nullptr;
   const ScalarFunc* sf = nullptr;
@@ -527,7 +547,7 @@ ScalarFunc* utl::parseTimeFunc (const char* func, const std::string& type)
   {
     IFEM::cout <<"(expression) ";
     if (func) cstr = strdup(func);
-    sf = parseTimeFunction("expression",cstr);
+    sf = parseTimeFunction("expression",cstr,eps);
   }
   else if (type == "linear")
     sf = parseTimeFunction(func,cstr);

src/Utility/Functions.h

@@ -53,6 +53,9 @@ public:
   //! \brief Returns whether the function is identically zero or not.
   virtual bool isZero() const { return scale == Real(0); }
 
+  //! \brief Returns the first-derivative of the function.
+  virtual Real deriv(Real x) const { return scale; }
+
 protected:
   //! \brief Evaluates the function at \a x.
   virtual Real evaluate(const Real& x) const { return scale*x; }
@@ -75,6 +78,9 @@ public:
   //! \brief Returns whether the function is identically zero or not.
   virtual bool isZero() const { return fval == Real(0); }
 
+  //! \brief Returns the first-derivative of the function.
+  virtual Real deriv(Real x) const;
+
 protected:
   //! \brief Evaluates the function at \a x.
   virtual Real evaluate(const Real& x) const;
@@ -143,6 +149,9 @@ public:
   //! \brief Returns whether the function is identically zero or not.
   virtual bool isZero() const { return scale == Real(0); }
 
+  //! \brief Returns the first-derivative of the function.
+  virtual Real deriv(Real x) const;
+
 protected:
   //! \brief Evaluates the function at \a x.
   virtual Real evaluate(const Real& x) const;
@@ -540,8 +549,10 @@ namespace utl
   //! \brief Creates a time function by parsing a character string.
   //! \param[in] func Character string to parse function definition from
   //! \param[in] type Function definition type flag
+  //! \param[in] eps Domain increment for calculation of numerical derivative
   ScalarFunc* parseTimeFunc(const char* func,
-                            const std::string& type = "expression");
+                            const std::string& type = "expression",
+                            Real eps = Real(1.0e-8));
 
   //! \brief Creates a scalar-valued function by parsing a character string.
   //! \param[in] func Character string to parse function definition from

src/Utility/Test/TestFunctions.C

@@ -0,0 +1,42 @@
+//==============================================================================
+//!
+//! \file TestFunctions.C
+//!
+//! \date Apr 28 2018
+//!
+//! \author Knut Morten Okstad / SINTEF
+//!
+//! \brief Tests for parsing of functions.
+//!
+//==============================================================================
+
+#include "Functions.h"
+#include <cstdlib>
+#include <cmath>
+
+#include "gtest/gtest.h"
+
+
+TEST(TestScalarFunc, ParseDerivative)
+{
+  const char* func1 = "sin(1.5*t)*t";
+  const char* func2 = "sin(1.5*t)*t:1.5*cos(1.5*t)*t+sin(1.5*t)";
+
+  ScalarFunc* f1 = utl::parseTimeFunc(func1,"expression");
+  ScalarFunc* f2 = utl::parseTimeFunc(func2,"expression");
+
+  ASSERT_TRUE(f1 != nullptr);
+  ASSERT_TRUE(f2 != nullptr);
+
+  double t = 0.0;
+  for (int i = 0; i < 20; i++)
+  {
+    t += 0.314*(double)random()/(double)RAND_MAX;
+    std::cout <<"f("<< t <<") = "<< (*f1)(t)
+	      <<"  f'("<< t <<") = "<< f1->deriv(t) << std::endl;
+    EXPECT_FLOAT_EQ((*f1)(t),sin(1.5*t)*t);
+    EXPECT_FLOAT_EQ((*f2)(t),sin(1.5*t)*t);
+    EXPECT_FLOAT_EQ(f1->deriv(t),1.5*cos(1.5*t)*t+sin(1.5*t));
+    EXPECT_FLOAT_EQ(f2->deriv(t),1.5*cos(1.5*t)*t+sin(1.5*t));
+  }
+}