added: RealFunc::gradient

returns the gradient of the function as a Vec3

src/Utility/ExprFunctions.C

@@ -324,15 +324,15 @@ void EvalFunction::addDerivative (const std::string& function,
       return 5;
   };
 
-  if (d1 > 0 && d1 <= 4 && d2 < 1) // A first derivative is specified
+  if (d1 > 0 && d1 <= 4 && d2 < 1) // A first order derivative is specified
   {
-    if (!gradient[--d1])
-      gradient[d1] = std::make_unique<EvalFunction>((variables+function).c_str());
+    if (!derivative1[--d1])
+      derivative1[d1] = std::make_unique<EvalFunction>((variables+function).c_str());
   }
-  else if ((d1 = voigtIdx(d1,d2)) >= 0) // A second derivative is specified
+  else if ((d1 = voigtIdx(d1,d2)) >= 0) // A second order derivative is specified
   {
-    if (!dgradient[d1])
-      dgradient[d1] = std::make_unique<EvalFunction>((variables+function).c_str());
+    if (!derivative2[d1])
+      derivative2[d1] = std::make_unique<EvalFunction>((variables+function).c_str());
   }
 }
 
@@ -369,8 +369,8 @@ Real EvalFunction::deriv (const Vec3& X, int dir) const
     return Real(0);
   else if (dir < 4)
   {
-    if (gradient[--dir])
-      return gradient[dir]->evaluate(X);
+    if (derivative1[--dir])
+      return derivative1[dir]->evaluate(X);
 
     // Evaluate spatial derivative using central difference
     Vec4 X0, X1;
@@ -380,8 +380,8 @@ Real EvalFunction::deriv (const Vec3& X, int dir) const
   }
   else if (!IAmConstant)
   {
-    if (gradient[3])
-      return gradient[3]->evaluate(X);
+    if (derivative1[3])
+      return derivative1[3]->evaluate(X);
 
     // Evaluate time-derivative using central difference
     Vec4 X0, X1;
@@ -410,7 +410,7 @@ Real EvalFunction::dderiv (const Vec3& X, int d1, int d2) const
   else // off-diagonal term, 13
     d1 = 5;
 
-  return dgradient[d1] ? dgradient[d1]->evaluate(X) : Real(0);
+  return derivative2[d1] ? derivative2[d1]->evaluate(X) : Real(0);
 }
 
 

src/Utility/ExprFunctions.h

@@ -102,8 +102,8 @@ class EvalFunction : public RealFunc
 
   std::vector<Arg> arg; //!< Function argument values
 
-  std::array<std::unique_ptr<EvalFunction>,4> gradient;  //!< First derivative expressions
-  std::array<std::unique_ptr<EvalFunction>,6> dgradient; //!< Second derivative expressions
+  std::array<std::unique_ptr<EvalFunction>,4> derivative1; //!< First order derivative expressions
+  std::array<std::unique_ptr<EvalFunction>,6> derivative2; //!< Second order derivative expressions
 
   bool IAmConstant; //!< Indicates whether the time coordinate is given or not
 

src/Utility/Function.h

@@ -193,6 +193,16 @@ public:
     return std::vector<Real>(1,this->evaluate(X));
   }
 
+  //! \brief Evaluates first derivatives of the function.
+  virtual Vec3 gradient(const Vec3& X) const
+  {
+    Vec3 result;
+    for (size_t d = 1; d <= 3; ++d)
+      result[d-1] = this->deriv(X,d);
+
+    return result;
+  }
+
   //! \brief Returns a representative scalar equivalent of the function value.
   virtual Real getScalarValue(const Vec3& X) const { return this->evaluate(X); }
 };

src/Utility/Test/TestFunctions.C

@@ -47,6 +47,66 @@ TEST(TestScalarFunc, ParseDerivative)
 }
 
 
+TEST(TestRealFunc, Gradient)
+{
+  const char* f1 = "sin(x)*sin(y)*sin(z)";
+  const char* d1 = "cos(x)*sin(y)*sin(z)";
+  const char* d2 = "sin(x)*cos(y)*sin(z)";
+  const char* d3 = "sin(x)*sin(y)*cos(z)";
+
+  EvalFunction f(f1);
+  f.addDerivative(d1, "", 1);
+  f.addDerivative(d2, "", 2);
+  f.addDerivative(d3, "", 3);
+
+  EXPECT_TRUE(f.isConstant());
+
+  for (double x : {0.1, 0.2, 0.3})
+    for (double y : {0.5, 0.6, 0.7})
+      for (double z : {0.8, 0.9, 1.0}) {
+        const Vec3 X(x,y,z);
+        const Vec3 r(cos(x)*sin(y)*sin(z),
+                     sin(x)*cos(y)*sin(z),
+                     sin(x)*sin(y)*cos(z));
+
+        const Vec3 grad = f.gradient(X);
+        for (size_t i = 1; i <= 3; ++i) {
+          EXPECT_DOUBLE_EQ(f.deriv(X, i), r[i-1]);
+          EXPECT_DOUBLE_EQ(grad[i-1], r[i-1]);
+        }
+      }
+}
+
+
+TEST(TestRealFunc, GradientFD)
+{
+  const char* f1 = "sin(x)*sin(y)*sin(z)";
+
+  const double eps = 1e-6;
+
+  EvalFunction f(f1, eps);
+
+  EXPECT_TRUE(f.isConstant());
+
+  for (double x : {0.1, 0.2, 0.3})
+    for (double y : {0.5, 0.6, 0.7})
+      for (double z : {0.8, 0.9, 1.0}) {
+        const Vec3 X(x,y,z);
+        const Vec3 Xp(x+0.5*eps, y+0.5*eps, z+0.5*eps);
+        const Vec3 Xm(x-0.5*eps, y-0.5*eps, z-0.5*eps);
+        const Vec3 r((sin(Xp.x) - sin(Xm.x))*sin(y)*sin(z) / eps,
+                     sin(x)*(sin(Xp.y) - sin(Xm.y))*sin(z) / eps,
+                     sin(x)*sin(y)*(sin(Xp.z) - sin(Xm.z)) / eps);
+
+        const Vec3 grad = f.gradient(X);
+        for (size_t i = 1; i <= 3; ++i) {
+          EXPECT_NEAR(f.deriv(X, i), r[i-1], 1e-8);
+          EXPECT_NEAR(grad[i-1], r[i-1], 1e-8);
+        }
+      }
+}
+
+
 TEST(TestVecFunc, Evaluate)
 {
   const char* func = "sin(x) | cos (y) | exp(z)";