Added: Time-derivative of time-dependent spatial functions

src/Utility/ExprFunctions.C

@@ -15,6 +15,9 @@
 #include "Vec3.h"
 #include "Tensor.h"
 #include "expreval.h"
+#ifdef USE_OPENMP
+#include <omp.h>
+#endif
 
 
 /*!
@@ -181,7 +184,8 @@ Real EvalFunc::deriv (Real x) const
 }
 
 
-EvalFunction::EvalFunction (const char* function) : gradient{}, dgradient{}
+EvalFunction::EvalFunction (const char* function, Real epsX, Real epsT)
+  : gradient{}, dgradient{}, dx(epsX), dt(epsT)
 {
   try {
 #ifdef USE_OPENMP
@@ -318,10 +322,29 @@ Real EvalFunction::evaluate (const Vec3& X) const
 
 Real EvalFunction::deriv (const Vec3& X, int dir) const
 {
-  if (dir < 1 || dir > 3 || !gradient[--dir])
+  if (dir < 1)
     return Real(0);
+  else if (dir < 4)
+  {
+    if (gradient[--dir])
+      return gradient[dir]->evaluate(X);
 
-  return gradient[dir]->evaluate(X);
+    // Evaluate spatial derivative using central difference
+    Vec4 X0, X1;
+    X0.assign(X); X0[dir] -= 0.5*dx;
+    X1.assign(X); X1[dir] += 0.5*dx;
+    return (this->evaluate(X1) - this->evaluate(X0)) / dx;
+  }
+  else if (!IAmConstant)
+  {
+    // Evaluate time-derivative using central difference
+    Vec4 X0, X1;
+    X0.assign(X); X0.t -= 0.5*dt;
+    X1.assign(X); X1.t += 0.5*dt;
+    return (this->evaluate(X1) - this->evaluate(X0)) / dt;
+  }
+  else
+    return Real(0);
 }
 
 

src/Utility/ExprFunctions.h

@@ -19,9 +19,6 @@
 #include <string>
 #include <vector>
 #include <array>
-#ifdef USE_OPENMP
-#include <omp.h>
-#endif
 
 namespace ExprEval {
   class Expression;
@@ -50,7 +47,8 @@ public:
   static int numError; //!< Error counter - set by the exception handler
 
   //! \brief The constructor parses the expression string.
-  explicit EvalFunc(const char* function, const char* x = "x", Real eps = Real(1.0e-8));
+  explicit EvalFunc(const char* function, const char* x = "x",
+                    Real eps = Real(1.0e-8));
   //! \brief The destructor frees the dynamically allocated objects.
   virtual ~EvalFunc();
 
@@ -102,9 +100,13 @@ class EvalFunction : public RealFunc
 
   bool IAmConstant; //!< Indicates whether the time coordinate is given or not
 
+  Real dx; //!< Domain increment for calculation of numerical derivative
+  Real dt; //!< Domain increment for calculation of numerical time-derivative
+
 public:
   //! \brief The constructor parses the expression string.
-  explicit EvalFunction(const char* function);
+  explicit EvalFunction(const char* function,
+                        Real epsX = Real(1.0e-8), Real epsT = Real(1.0e-12));
   //! \brief The destructor frees the dynamically allocated objects.
   virtual ~EvalFunction();
 

src/Utility/Function.h

@@ -241,6 +241,9 @@ class TractionFunc : public utl::Function2<Vec3,Vec3>
 public:
   //! \brief Returns whether the traction is always normal to the face or not.
   virtual bool isNormalPressure() const { return false; }
+
+  //! \brief Returns the time-derivative of the function.
+  virtual Vec3 deriv(const Vec3&, const Vec3&) const { return Vec3(); }
 };
 
 #endif

src/Utility/Functions.C

@@ -271,10 +271,16 @@ Real SineFunc::deriv (Real x) const
 Real ConstTimeFunc::evaluate (const Vec3& X) const
 {
   const Vec4* Xt = dynamic_cast<const Vec4*>(&X);
-  if (Xt)
-    return (*tfunc)(Xt->t);
-  else
-    return (*tfunc)(Real(0));
+  return (*tfunc)(Xt ? Xt->t : Real(0));
+}
+
+
+Real ConstTimeFunc::deriv (const Vec3& X, int dir) const
+{
+  if (dir < 4) return Real(0);
+
+  const Vec4* Xt = dynamic_cast<const Vec4*>(&X);
+  return tfunc->deriv(Xt ? Xt->t : Real(0));
 }
 
 
@@ -288,14 +294,24 @@ Real SpaceTimeFunc::evaluate (const Vec3& X) const
 Real SpaceTimeFunc::deriv (const Vec3& X, int dir) const
 {
   const Vec4* Xt = dynamic_cast<const Vec4*>(&X);
-  return sfunc->deriv(X,dir) * (*tfunc)(Xt ? Xt->t : Real(0));
+  if (dir < 4)
+    return sfunc->deriv(X,dir) * (*tfunc)(Xt ? Xt->t : Real(0));
+  else
+    return (*sfunc)(X) * tfunc->deriv(Xt ? Xt->t : Real(0));
 }
 
 
 Real SpaceTimeFunc::dderiv (const Vec3& X, int dir1, int dir2) const
 {
   const Vec4* Xt = dynamic_cast<const Vec4*>(&X);
-  return sfunc->dderiv(X,dir1,dir2) * (*tfunc)(Xt ? Xt->t : Real(0));
+  if (dir1 < 4 && dir2 < 4)
+    return sfunc->dderiv(X,dir1,dir2) * (*tfunc)(Xt ? Xt->t : Real(0));
+  else if (dir1 < 4)
+    return sfunc->deriv(X,dir1) * tfunc->deriv(Xt ? Xt->t : Real(0));
+  else if (dir2 < 4)
+    return sfunc->deriv(X,dir2) * tfunc->deriv(Xt ? Xt->t : Real(0));
+  else
+    return Real(0);
 }
 
 
@@ -465,6 +481,26 @@ Real Interpolate1D::evaluate (const Vec3& X) const
 }
 
 
+Real Interpolate1D::deriv (const Vec3& X, int ddir) const
+{
+  Real res = Real(0);
+  if (ddir == dir+1)
+    res = lfunc.deriv(X[dir]);
+  else if (ddir > 3)
+    res = lfunc(X[dir]);
+  else
+    return res;
+
+  const Vec4* Xt = dynamic_cast<const Vec4*>(&X);
+  if (Xt && time > Real(0) && Xt->t < time)
+    res *= (ddir < 4 ? Xt->t : 1.0)/time;
+  else if (ddir > 3)
+    res = Real(0);
+
+  return res;
+}
+
+
 /*!
   The functions are assumed on the general form
   \f[ f({\bf X},t) = A * g({\bf X}) * h(t) \f]

src/Utility/Functions.h

@@ -235,7 +235,7 @@ protected:
 
 class ConstTimeFunc : public RealFunc
 {
-  const ScalarFunc* tfunc; //!< The time dependent function value
+  const ScalarFunc* tfunc; //!< The time-dependent function value
 
 public:
   //! \brief Constructor initializing the function value.
@@ -248,6 +248,9 @@ public:
   //! \brief Returns whether the function is time-independent or not.
   virtual bool isConstant() const { return tfunc->isZero(); }
 
+  //! \brief Returns first-derivative of the function.
+  virtual Real deriv(const Vec3& X, int dir) const;
+
 protected:
   //! \brief Evaluates the time-varying function.
   virtual Real evaluate(const Vec3& X) const;
@@ -497,8 +500,8 @@ class StepXFunc : public RealFunc
 
 public:
   //! \brief Constructor initializing the function parameters.
-  explicit StepXFunc(Real v, Real X0 = Real(0), Real X1 = Real(1), char dir = 'X')
-    : fv(v), x0(X0), x1(X1), d(dir) {}
+  explicit StepXFunc(Real v, Real a = Real(0), Real b = Real(1), char dir = 'X')
+    : fv(v), x0(a), x1(b), d(dir) {}
 
   //! \brief Returns whether the function is identically zero or not.
   virtual bool isZero() const { return fv == Real(0); }
@@ -561,6 +564,9 @@ public:
   //! \brief Returns whether the function is time-independent or not.
   virtual bool isConstant() const { return time <= Real(0); }
 
+  //! \brief Returns first-derivative of the function.
+  virtual Real deriv(const Vec3&, int ddir) const;
+
 protected:
   //! \brief Evaluates the function by interpolation.
   virtual Real evaluate(const Vec3& X) const;

src/Utility/TractionField.C

@@ -41,6 +41,17 @@ Vec3 TractionField::evaluate (const Vec3& x, const Vec3& n) const
 }
 
 
+Vec3 TractionField::deriv (const Vec3& x, const Vec3& n) const
+{
+  if (sigma) // symmetric tensor field
+    return sigma->deriv(x,4) * n;
+  else if (sigmaN) // non-symmetric tensor field
+    return sigmaN->deriv(x,4) * n;
+  else // zero tensor field
+    return Vec3();
+}
+
+
 bool TractionField::isZero () const
 {
   if (sigma)
@@ -79,6 +90,31 @@ Vec3 PressureField::evaluate (const Vec3& x, const Vec3& n) const
 }
 
 
+Vec3 PressureField::deriv (const Vec3& x, const Vec3& n) const
+{
+  if (!pressure) // zero pressure field
+    return Vec3();
+
+  if (pdir < 1 && !pdfn) // normal pressure
+    return pressure->deriv(x,4) * n;
+
+  Vec3 t;
+  if (pdfn) // pressure direction specified as a function
+  {
+    t = (*pdfn)(x);
+    t.normalize();
+    t *= pressure->deriv(x,4);
+  }
+  else if (pdir > 0) // pressure acting in global pdir direction
+    t[(pdir-1)%3] = pressure->deriv(x,4);
+
+  if (pdir > 3) // normal pressure in global pdir direction
+    t = (t*n) * n;
+
+  return t;
+}
+
+
 ForceDirField::~ForceDirField ()
 {
   delete force;

src/Utility/TractionField.h

@@ -47,6 +47,11 @@ public:
   //! \brief Returns whether the function is identically zero or not.
   virtual bool isZero() const;
 
+  //! \brief Returns the time-derivative of the function.
+  //! \param[in] x Global coordinates of evaluation point
+  //! \param[in] n Outward-directed unit normal vector at evaluation point
+  virtual Vec3 deriv(const Vec3& x, const Vec3& n) const;
+
 protected:
   //! \brief Evaluates the traction field function at the specified point.
   //! \param[in] x Global coordinates of evaluation point
@@ -95,6 +100,11 @@ public:
   //! \brief Returns whether the function is identically zero or not.
   virtual bool isZero() const { return pressure ? pressure->isZero() : true; }
 
+  //! \brief Returns the time-derivative of the function.
+  //! \param[in] x Global coordinates of evaluation point
+  //! \param[in] n Outward-directed unit normal vector at evaluation point
+  virtual Vec3 deriv(const Vec3& x, const Vec3& n) const;
+
 protected:
   //! \brief Evaluates the traction field function at the specified point.
   //! \param[in] x Global coordinates of evaluation point