diff --git a/src/Utility/ExprFunctions.C b/src/Utility/ExprFunctions.C index 99962b66..fc5932f1 100644 --- a/src/Utility/ExprFunctions.C +++ b/src/Utility/ExprFunctions.C @@ -499,6 +499,19 @@ EvalMultiFunction::evalGradient (const Vec3& X) const } +template +std::vector +EvalMultiFunction::evalTimeDerivative (const Vec3& X) const +{ + std::vector result; + result.reserve(this->ncmp); + for (const std::unique_ptr& f : this->p) + result.push_back(f->deriv(X,4)); + + return result; +} + + template class EvalMultiFunction; template class EvalMultiFunction; template class EvalMultiFunction; diff --git a/src/Utility/ExprFunctions.h b/src/Utility/ExprFunctions.h index edcdcf1a..4d84743b 100644 --- a/src/Utility/ExprFunctions.h +++ b/src/Utility/ExprFunctions.h @@ -226,6 +226,9 @@ protected: //! \brief Returns the gradient of the function as a 1D array. std::vector evalGradient(const Vec3& X) const override; + + //! \brief Returns the time derivatives of the function as a 1D array. + std::vector evalTimeDerivative(const Vec3& X) const override; }; //! Vector-valued function expression diff --git a/src/Utility/Function.h b/src/Utility/Function.h index b3c9e112..9f02ca58 100644 --- a/src/Utility/Function.h +++ b/src/Utility/Function.h @@ -112,6 +112,8 @@ namespace utl protected: //! \brief Returns the gradient of the function as a 1D array. virtual std::vector evalGradient(const Vec3&) const { return {}; } + //! \brief Returns the time derivatives of the function as a 1D array. + virtual std::vector evalTimeDerivative(const Vec3&) const { return {}; } Result zero; //!< Return value for default implementations of derivatives }; @@ -220,6 +222,9 @@ public: //! \brief Returns a representative scalar equivalent of the function value. virtual Real getScalarValue(const Vec3& X) const { return this->evaluate(X); } + + //! \brief Returns the time derivative of the function. + Real timeDerivative(const Vec3& X) const { return this->deriv(X,4); } }; @@ -262,6 +267,9 @@ public: result = this->evalGradient(X); return result; } + + //! \brief Evaluates time derivatives of the function. + Vec3 timeDerivative(const Vec3& X) const { return this->evalTimeDerivative(X); } }; @@ -276,7 +284,7 @@ public: virtual bool isNormalPressure() const { return false; } //! \brief Returns the time-derivative of the function. - virtual Vec3 deriv(const Vec3&, const Vec3&) const { return Vec3(); } + virtual Vec3 timeDerivative(const Vec3&, const Vec3&) const { return Vec3(); } }; #endif diff --git a/src/Utility/TensorFunction.C b/src/Utility/TensorFunction.C index 24e7acac..29c39c14 100644 --- a/src/Utility/TensorFunction.C +++ b/src/Utility/TensorFunction.C @@ -23,6 +23,15 @@ utl::matrix3d TensorFunc::gradient(const Vec3& X) const } +Tensor TensorFunc::timeDerivative(const Vec3& X) const +{ + const size_t nsd = sqrt(ncmp); + Tensor result(nsd); + result = this->evalTimeDerivative(X); + return result; +} + + size_t STensorFunc::index(size_t nsd, size_t i, size_t j) const { if (i == j) @@ -50,3 +59,12 @@ utl::matrix3d STensorFunc::gradient(const Vec3& X) const return result; } + + +SymmTensor STensorFunc::timeDerivative(const Vec3& X) const +{ + const size_t nsd = ncmp > 5 ? 3 : (ncmp > 2 ? 2 : 1); + SymmTensor result(nsd, ncmp == 4); + result = this->evalTimeDerivative(X); + return result; +} diff --git a/src/Utility/TensorFunction.h b/src/Utility/TensorFunction.h index adc4dcd6..0be131fc 100644 --- a/src/Utility/TensorFunction.h +++ b/src/Utility/TensorFunction.h @@ -53,6 +53,9 @@ public: //! \brief Evaluates first derivatives of the function. utl::matrix3d gradient(const Vec3& X) const; + + //! \brief Evaluates time derivatives of the function. + Tensor timeDerivative(const Vec3& X) const; }; @@ -95,6 +98,9 @@ public: //! \brief Evaluates first derivatives of the function. utl::matrix3d gradient(const Vec3& X) const; + + //! \brief Evaluates time derivatives of the function. + SymmTensor timeDerivative(const Vec3& X) const; }; #endif diff --git a/src/Utility/Test/TestFunctions.C b/src/Utility/Test/TestFunctions.C index 9539190f..0b451a5e 100644 --- a/src/Utility/Test/TestFunctions.C +++ b/src/Utility/Test/TestFunctions.C @@ -320,6 +320,61 @@ TEST(TestVecFuncExpr, NumDimensions) } +TEST(TestVecFuncExpr, TimeDerivative) +{ + const char* g = "sin(x)*sin(y)*sin(z)*sin(t) | x*x*y*y*z*t*t | exp(-2*t)"; + const char* g_t = "sin(x)*sin(y)*sin(z)*cos(t) | x*x*y*y*z*2*t | -2*exp(-2*t)"; + + VecFuncExpr f(g); + f.addDerivative(g_t,"",4); + + EXPECT_FALSE(f.isConstant()); + + for (double t : {0.1, 0.2, 0.3}) + 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 Vec4 X(x,y,z,t); + const Vec3 r(sin(x)*sin(y)*sin(z)*cos(t), + x*x*y*y*z*2*t, + -2*exp(-2*t)); + const Vec3 dt = f.timeDerivative(X); + EXPECT_DOUBLE_EQ(dt[0], r[0]); + EXPECT_DOUBLE_EQ(dt[1], r[1]); + EXPECT_DOUBLE_EQ(dt[2], r[2]); + } +} + + +TEST(TestVecFuncExpr, TimeDerivativeFD) +{ + const char* g = "sin(x)*sin(y)*sin(z)*sin(t) | x*x*y*y*z*t*t | exp(-2*t)"; + + const double eps = 1e-6; + VecFuncExpr f(g,"",1e-8,eps); + + EXPECT_FALSE(f.isConstant()); + + for (double t : {0.1, 0.2, 0.3}) + 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 Vec4 X(x,y,z,t); + const Vec4 Xp(x,y,z,t + 0.5*eps); + const Vec4 Xm(x,y,z,t - 0.5*eps); + Vec3 r(sin(x)*sin(y)*sin(z)*(sin(Xp.t) - sin(Xm.t)), + x*x*y*y*z*(Xp.t*Xp.t - Xm.t*Xm.t), + exp(-2*Xp.t) - exp(-2*Xm.t)); + r *= 1.0 / eps; + + const Vec3 dt = f.timeDerivative(X); + EXPECT_NEAR(dt[0], r[0], 1e-8); + EXPECT_NEAR(dt[1], r[1], 1e-8); + EXPECT_NEAR(dt[2], r[2], 1e-8); + } +} + + TEST(TestTensorFunc, Evaluate) { const char* func = "sin(x) | cos (y) | exp(z) | sin(x)*cos(y)"; @@ -515,7 +570,7 @@ TEST(TestTensorFunction, Gradient3DFD) EXPECT_NEAR(grad(i,j,d), r(i,j,d), 1e-8); } } - } + } } @@ -647,6 +702,39 @@ TEST(TestTensorFunction, Hessian3D) } +TEST(TestTensorFunction, TimeDerivative) +{ + const char* g = "sin(x)*sin(y)*sin(z)*sin(t) | x*x*y*y*z*t*t | exp(x)*exp(2*y)*z*z*t |" + "exp(-2*x)*exp(y)*z*sin(t) | x*y*z*t*t | x*y*z*z*t |" + "x*sin(t) | y*t*t | z*t"; + const char* g_t = "sin(x)*sin(y)*sin(z)*cos(t) | x*x*y*y*z*2*t | exp(x)*exp(2*y)*z*z |" + "exp(-2*x)*exp(y)*z*cos(t) | x*y*z*2*t | x*y*z*z |" + "x*cos(t) | y*2*t | z"; + + TensorFuncExpr f(g); + f.addDerivative(g_t,"",4); + + EXPECT_FALSE(f.isConstant()); + + for (double t : {0.1, 0.2, 0.3}) + 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 Vec4 X(x,y,z,t); + const Tensor r({sin(x)*sin(y)*sin(z)*cos(t), x*x*y*y*z*2*t, exp(x)*exp(2*y)*z*z, + exp(-2*x)*exp(y)*z*cos(t), x*y*z*2*t, x*y*z*z, + x*cos(t), y*2*t, z}); + + const Tensor dt = f.timeDerivative(X); + for (size_t i = 1; i <= 3; ++i) + for (size_t j = 1; j <= 3; ++j) + EXPECT_DOUBLE_EQ(dt(i,j), r(i,j)); + } +} + + + + TEST(TestTensorFuncExpr, NumDimensions) { const char* func1 = "x"; @@ -953,6 +1041,58 @@ TEST(TestSTensorFuncExpr, NumDimensions) } +TEST(TestSTensorFunction, TimeDerivative) +{ + const char* g = "sin(x)*sin(y)*sin(t) | exp(x)*exp(2*y)*exp(-4*t) | x*x*y*y*t*t"; + const char* g_t = "sin(x)*sin(y)*cos(t) | exp(x)*exp(2*y)*-4*exp(-4*t) | x*x*y*y*2*t"; + + STensorFuncExpr f(g); + f.addDerivative(g_t,"",4); + + EXPECT_FALSE(f.isConstant()); + + for (double t : {0.1, 0.2, 0.3}) + for (double x : {0.1, 0.2, 0.3}) + for (double y : {0.5, 0.6, 0.7}) { + const Vec4 X(x,y,0,t); + const SymmTensor r({sin(x)*sin(y)*cos(t), exp(x)*exp(2*y)*-4*exp(-4*t), x*x*y*y*2*t}); + + const SymmTensor dt = f.timeDerivative(X); + for (size_t i = 1; i <= 2; ++i) + for (size_t j = 1; j <= 2; ++j) + EXPECT_DOUBLE_EQ(dt(i,j), r(i,j)); + } +} + + +TEST(TestSTensorFunction, TimeDerivativeFD) +{ + const char* g = "sin(x)*sin(y)*sin(t) | exp(x)*exp(2*y)*exp(-4*t) | x*x*y*y*t*t"; + + const double eps = 1e-6; + STensorFuncExpr f(g,"",1e-8,eps); + + EXPECT_FALSE(f.isConstant()); + + for (double t : {0.1, 0.2, 0.3}) + for (double x : {0.1, 0.2, 0.3}) + for (double y : {0.5, 0.6, 0.7}) { + const Vec4 X(x,y,0,t); + const Vec4 Xp(x,y,0,t + 0.5*eps); + const Vec4 Xm(x,y,0,t - 0.5*eps); + SymmTensor r({sin(x)*sin(y)*(sin(Xp.t) - sin(Xm.t)), + exp(x)*exp(2*y)*(exp(-4*Xp.t) - exp(-4*Xm.t)), + x*x*y*y*(Xp.t*Xp.t - Xm.t*Xm.t)}); + r *= 1.0 / eps; + + const SymmTensor dt = f.timeDerivative(X); + for (size_t i = 1; i <= 2; ++i) + for (size_t j = 1; j <= 2; ++j) + EXPECT_NEAR(dt(i,j), r(i,j), 1e-8); + } +} + + TEST(TestEvalFunction, ExtraParam) { EvalFunction f("x*foo");