From e2c76eaa085d8bd134485057d1d89595336f392a Mon Sep 17 00:00:00 2001 From: Knut Morten Okstad Date: Sun, 30 Jan 2022 13:19:05 +0100 Subject: [PATCH] Added: Time-derivative of time-dependent spatial functions --- src/Utility/ExprFunctions.C | 29 +++++++++++++++++++--- src/Utility/ExprFunctions.h | 12 ++++++---- src/Utility/Function.h | 3 +++ src/Utility/Functions.C | 48 ++++++++++++++++++++++++++++++++----- src/Utility/Functions.h | 12 +++++++--- src/Utility/TractionField.C | 36 ++++++++++++++++++++++++++++ src/Utility/TractionField.h | 10 ++++++++ 7 files changed, 133 insertions(+), 17 deletions(-) diff --git a/src/Utility/ExprFunctions.C b/src/Utility/ExprFunctions.C index 02d1e79e..cad159bb 100644 --- a/src/Utility/ExprFunctions.C +++ b/src/Utility/ExprFunctions.C @@ -15,6 +15,9 @@ #include "Vec3.h" #include "Tensor.h" #include "expreval.h" +#ifdef USE_OPENMP +#include +#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); } diff --git a/src/Utility/ExprFunctions.h b/src/Utility/ExprFunctions.h index 5757e680..182e1e72 100644 --- a/src/Utility/ExprFunctions.h +++ b/src/Utility/ExprFunctions.h @@ -19,9 +19,6 @@ #include #include #include -#ifdef USE_OPENMP -#include -#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(); diff --git a/src/Utility/Function.h b/src/Utility/Function.h index 5e9e49d5..df8bbf27 100644 --- a/src/Utility/Function.h +++ b/src/Utility/Function.h @@ -241,6 +241,9 @@ class TractionFunc : public utl::Function2 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 diff --git a/src/Utility/Functions.C b/src/Utility/Functions.C index 5b2214a8..cc623c55 100644 --- a/src/Utility/Functions.C +++ b/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(&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(&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(&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(&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(&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] diff --git a/src/Utility/Functions.h b/src/Utility/Functions.h index 585b4075..eb03caad 100644 --- a/src/Utility/Functions.h +++ b/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; diff --git a/src/Utility/TractionField.C b/src/Utility/TractionField.C index 7ca381ec..1987db19 100644 --- a/src/Utility/TractionField.C +++ b/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; diff --git a/src/Utility/TractionField.h b/src/Utility/TractionField.h index 223d11a5..87f041d1 100644 --- a/src/Utility/TractionField.h +++ b/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