diff --git a/src/Utility/ExprFunctions.C b/src/Utility/ExprFunctions.C index 93222731..3d3c91b9 100644 --- a/src/Utility/ExprFunctions.C +++ b/src/Utility/ExprFunctions.C @@ -305,7 +305,9 @@ Real EvalFuncScalar::deriv (Real x) const } -EvalFunction::EvalFunction (const char* function, Real epsX, Real epsT) +template +EvalFuncSpatial:: +EvalFuncSpatial (const char* function, Real epsX, Real epsT) : dx(epsX), dt(epsT) { try { @@ -354,27 +356,31 @@ EvalFunction::EvalFunction (const char* function, Real epsX, Real epsT) } -EvalFunction::~EvalFunction () = default; +template +EvalFuncSpatial::~EvalFuncSpatial () = default; -void EvalFunction::addDerivative (const std::string& function, - const std::string& variables, - int d1, int d2) +template +void EvalFuncSpatial:: +addDerivative (const std::string& function, + const std::string& variables, + int d1, int d2) { if (d1 > 0 && d1 <= 4 && d2 < 1) // A first order derivative is specified { if (!derivative1[--d1]) - derivative1[d1] = std::make_unique((variables+function).c_str()); + derivative1[d1] = std::make_unique((variables+function).c_str()); } else if ((d1 = voigtIdx(d1,d2)) >= 0) // A second order derivative is specified { if (!derivative2[d1]) - derivative2[d1] = std::make_unique((variables+function).c_str()); + derivative2[d1] = std::make_unique((variables+function).c_str()); } } -Real EvalFunction::evaluate (const Vec3& X) const +template +Real EvalFuncSpatial::evaluate (const Vec3& X) const { const Vec4* Xt = dynamic_cast(&X); Real result = Real(0); @@ -390,7 +396,10 @@ Real EvalFunction::evaluate (const Vec3& X) const *arg[i].y = X.y; *arg[i].z = X.z; *arg[i].t = Xt ? Xt->t : Real(0); - result = expr[i]->Evaluate(); + if constexpr (std::is_same_v) + result = expr[i]->Evaluate(); + else + result = expr[i]->Evaluate().expr->val; } catch (ExprEval::Exception& e) { ExprException(e,"evaluating expression"); @@ -400,7 +409,8 @@ Real EvalFunction::evaluate (const Vec3& X) const } -Real EvalFunction::deriv (const Vec3& X, int dir) const +template<> +Real EvalFuncSpatial::deriv (const Vec3& X, int dir) const { if (dir < 1) return Real(0); @@ -431,7 +441,31 @@ Real EvalFunction::deriv (const Vec3& X, int dir) const } -Real EvalFunction::dderiv (const Vec3& X, int d1, int d2) const +template<> +Real EvalFuncSpatial::deriv (const Vec3& X, int dir) const +{ + if (dir < 1 || dir > 4) + return Real(0); + + size_t i = 0; +#ifdef USE_OPENMP + i = omp_get_thread_num(); +#endif + + const Vec4* Xt = dynamic_cast(&X); + *arg[i].x = X.x; + *arg[i].y = X.y; + *arg[i].z = X.z; + *arg[i].t = Xt ? Xt->t : Real(0); + + // Evaluate spatial derivative using auto-diff + return derivativesx(expr[i]->Evaluate(), + autodiff::wrt(arg[i].get(dir)))[0].expr->val; +} + + +template<> +Real EvalFuncSpatial::dderiv (const Vec3& X, int d1, int d2) const { if ((d1 = voigtIdx(d1,d2)) < 0) return Real(0); @@ -440,10 +474,88 @@ Real EvalFunction::dderiv (const Vec3& X, int d1, int d2) const } -void EvalFunction::setParam (const std::string& name, double value) +template<> +Real EvalFuncSpatial::dderiv (const Vec3& X, int d1, int d2) const +{ + if (d1 < 1 || d1 > 3 || + d2 < 1 || d2 > 3) + return Real(0); + + size_t i = 0; +#ifdef USE_OPENMP + i = omp_get_thread_num(); +#endif + + const Vec4* Xt = dynamic_cast(&X); + *arg[i].x = X.x; + *arg[i].y = X.y; + *arg[i].z = X.z; + *arg[i].t = Xt ? Xt->t : Real(0); + + return derivativesx(derivativesx(expr[i]->Evaluate(), + autodiff::wrt(arg[i].get(d1)))[0], + autodiff::wrt(arg[i].get(d2)))[0].expr->val; +} + + +template<> +Vec3 EvalFuncSpatial::gradient (const Vec3& X) const +{ + size_t i = 0; +#ifdef USE_OPENMP + i = omp_get_thread_num(); +#endif + + const Vec4* Xt = dynamic_cast(&X); + double t; + *arg[i].x = X.x; + *arg[i].y = X.y; + *arg[i].z = X.z; + *arg[i].t = t = Xt ? Xt->t : Real(0); + + const auto dx = derivativesx(expr[i]->Evaluate(), + autodiff::wrt(*arg[i].x, *arg[i].y, *arg[i].z)); + + return Vec4(dx[0].expr->val, dx[1].expr->val, dx[2].expr->val, t); +} + + +template<> +SymmTensor EvalFuncSpatial::hessian (const Vec3& X) const +{ + size_t i = 0; +#ifdef USE_OPENMP + i = omp_get_thread_num(); +#endif + + const Vec4* Xt = dynamic_cast(&X); + *arg[i].x = X.x; + *arg[i].y = X.y; + *arg[i].z = X.z; + *arg[i].t = Xt ? Xt->t : Real(0); + + const auto dx = + derivativesx(expr[i]->Evaluate(), autodiff::wrt(*arg[i].x, *arg[i].y, *arg[i].z)); + + const auto [uxx, uxy, uxz] = + derivativesx(dx[0], autodiff::wrt(*arg[i].x, *arg[i].y, *arg[i].z)); + + const auto [uyy, uyz] = + derivativesx(dx[1], autodiff::wrt(*arg[i].y, *arg[i].z)); + + const auto [uzz] = + derivativesx(dx[2], autodiff::wrt(*arg[i].z)); + + return SymmTensor({uxx.expr->val, uyy.expr->val, uzz.expr->val, + uxy.expr->val, uyz.expr->val, uxz.expr->val}); +} + + +template +void EvalFuncSpatial::setParam (const std::string& name, double value) { for (std::unique_ptr& v1 : v) { - double* address = v1->GetAddress(name); + Scalar* address = v1->GetAddress(name); if (!address) v1->Add(name,value,false); else @@ -452,21 +564,24 @@ void EvalFunction::setParam (const std::string& name, double value) } -EvalFunctions::EvalFunctions (const std::string& functions, - const std::string& variables, - const Real epsX, const Real epsT) +template +EvalFunctions::EvalFunctions (const std::string& functions, + const std::string& variables, + const Real epsX, const Real epsT) { std::vector components = splitComps(functions,variables); for (const std::string& comp : components) - p.emplace_back(std::make_unique(comp.c_str(),epsX,epsT)); + p.emplace_back(std::make_unique(comp.c_str(),epsX,epsT)); } -EvalFunctions::~EvalFunctions () = default; +template +EvalFunctions::~EvalFunctions () = default; -void EvalFunctions::addDerivative (const std::string& functions, - const std::string& variables, int d1, int d2) +template +void EvalFunctions::addDerivative (const std::string& functions, + const std::string& variables, int d1, int d2) { std::vector components = splitComps(functions,variables); for (size_t i = 0; i < p.size() && i < components.size(); i++) @@ -474,8 +589,9 @@ void EvalFunctions::addDerivative (const std::string& functions, } -template -Ret EvalMultiFunction::evaluate (const Vec3& X) const +template +Ret EvalMultiFunction:: +evaluate (const Vec3& X) const { std::vector res_array(this->p.size()); for (size_t i = 0; i < this->p.size(); ++i) @@ -485,44 +601,47 @@ Ret EvalMultiFunction::evaluate (const Vec3& X) const } -template -void EvalMultiFunction::setNoDims () +template +void EvalMultiFunction::setNoDims () { std::tie(nsd, this->ncmp) = getNoDims(this->p.size()); } -template -Ret EvalMultiFunction::deriv (const Vec3& X, int dir) const +template +Ret EvalMultiFunction:: +deriv (const Vec3& X, int dir) const { - std::vector tmp(p.size()); - for (size_t i = 0; i < p.size(); ++i) - tmp[i] = p[i]->deriv(X,dir); + std::vector tmp(this->p.size()); + for (size_t i = 0; i < this->p.size(); ++i) + tmp[i] = this->p[i]->deriv(X,dir); return Ret(tmp); } -template -Ret EvalMultiFunction::dderiv (const Vec3& X, int d1, int d2) const +template +Ret EvalMultiFunction:: +dderiv (const Vec3& X, int d1, int d2) const { - std::vector tmp(p.size()); - for (size_t i = 0; i < p.size(); ++i) - tmp[i] = p[i]->dderiv(X,d1,d2); + std::vector tmp(this->p.size()); + for (size_t i = 0; i < this->p.size(); ++i) + tmp[i] = this->p[i]->dderiv(X,d1,d2); return Ret(tmp); } -template +template std::vector -EvalMultiFunction::evalGradient (const Vec3& X) const +EvalMultiFunction:: +evalGradient (const Vec3& X) const { std::vector result; result.reserve(this->ncmp*this->nsd); std::vector dx; dx.reserve(this->p.size()); - for (const std::unique_ptr& f : this->p) + for (const std::unique_ptr& f : this->p) dx.push_back(f->gradient(X)); for (size_t d = 1; d <= this->nsd; ++d) @@ -533,15 +652,16 @@ EvalMultiFunction::evalGradient (const Vec3& X) const } -template +template std::vector -EvalMultiFunction::evalHessian (const Vec3& X) const +EvalMultiFunction:: +evalHessian (const Vec3& X) const { std::vector result; result.reserve(this->p.size()*this->nsd*this->nsd); std::vector dx; dx.reserve(this->p.size()); - for (const std::unique_ptr& f : this->p) + for (const std::unique_ptr& f : this->p) dx.push_back(f->hessian(X)); for (size_t d2 = 1; d2 <= this->nsd; ++d2) @@ -553,14 +673,14 @@ EvalMultiFunction::evalHessian (const Vec3& X) const } -template +template std::vector -EvalMultiFunction::evalTimeDerivative (const Vec3& X) const +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)); + for (const std::unique_ptr& f : this->p) + result.push_back(f->timeDerivative(X)); return result; } @@ -568,6 +688,13 @@ EvalMultiFunction::evalTimeDerivative (const Vec3& X) const template class EvalFuncScalar; template class EvalFuncScalar; -template class EvalMultiFunction; -template class EvalMultiFunction; -template class EvalMultiFunction; +template class EvalFuncSpatial; +template class EvalFuncSpatial; +template class EvalFunctions; +template class EvalFunctions; +template class EvalMultiFunction; +template class EvalMultiFunction; +template class EvalMultiFunction; +template class EvalMultiFunction; +template class EvalMultiFunction; +template class EvalMultiFunction; diff --git a/src/Utility/ExprFunctions.h b/src/Utility/ExprFunctions.h index 31d1d30f..857b720e 100644 --- a/src/Utility/ExprFunctions.h +++ b/src/Utility/ExprFunctions.h @@ -80,11 +80,13 @@ protected: \brief A scalar-valued spatial function, general function expression. */ -class EvalFunction : public RealFunc +template +class EvalFuncSpatial : public RealFunc { - using Expression = ExprEval::Expression; //!< Type alias for expression tree - using FunctionList = ExprEval::FunctionList; //!< Type alias for function list - using ValueList = ExprEval::ValueList; //!< Type alias for value list + using Expression = ExprEval::Expression; //!< Type alias for expression tree + using FunctionList = ExprEval::FunctionList; //!< Type alias for function list + using ValueList = ExprEval::ValueList; //!< Type alias for value list + using FuncType = EvalFuncSpatial; //!< Type alias for function std::vector> expr; //!< Roots of the expression tree std::vector> f; //!< Lists of functions std::vector> v; //!< Lists of variables and constants @@ -92,16 +94,27 @@ class EvalFunction : public RealFunc //! \brief A struct representing a spatial function argument. struct Arg { - Real* x; //!< X-coordinate - Real* y; //!< Y-coordinate - Real* z; //!< Z-coordinate - Real* t; //!< Time + Scalar* x; //!< X-coordinate + Scalar* y; //!< Y-coordinate + Scalar* z; //!< Z-coordinate + Scalar* t; //!< Time + + const Scalar& get(int dir) const + { + switch (dir) { + default: + case 1: return *x; + case 2: return *y; + case 3: return *z; + case 4: return *t; + } + } }; std::vector arg; //!< Function argument values - std::array,4> derivative1; //!< First order derivative expressions - std::array,6> derivative2; //!< Second order derivative expressions + std::array,4> derivative1; //!< First order derivative expressions + std::array,6> derivative2; //!< Second order derivative expressions bool IAmConstant; //!< Indicates whether the time coordinate is given or not @@ -110,12 +123,12 @@ class EvalFunction : public RealFunc public: //! \brief The constructor parses the expression string. - explicit EvalFunction(const char* function, - Real epsX = Real(1.0e-8), Real epsT = Real(1.0e-12)); + explicit EvalFuncSpatial(const char* function, + Real epsX = Real(1.0e-8), Real epsT = Real(1.0e-12)); //! \brief Defaulted destructor. //! \details The implementation needs to be in compile unit so we have the //! definition for the types of the unique_ptr's. - virtual ~EvalFunction(); + virtual ~EvalFuncSpatial(); //! \brief Adds an expression function for a first or second derivative. void addDerivative(const std::string& function, const std::string& variables, @@ -132,11 +145,19 @@ public: //! \brief Set an additional parameter in the function. void setParam(const std::string& name, double value); + //! \brief Evaluates first derivatives of the function. + Vec3 gradient(const Vec3& X) const override + { + return this->RealFunc::gradient(X); + } + + //! \brief Evaluates first derivatives of the function. + SymmTensor hessian(const Vec3& X) const override + { + return this->RealFunc::hessian(X); + } + 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. Real evaluate(const Vec3& X) const override; }; @@ -146,9 +167,12 @@ protected: \brief A base class for multi-component expression functions. */ +template class EvalFunctions { protected: + using FuncType = EvalFuncSpatial; //!< Type alias for function + //! \brief The constructor parses the expression string for each component. EvalFunctions(const std::string& functions, const std::string& variables, const Real epsX, const Real epsT); @@ -163,7 +187,7 @@ public: const std::string& variables, int d1, int d2 = 0); protected: - std::vector> p; //!< Array of component expressions + std::vector> p; //!< Array of component expressions }; @@ -172,10 +196,11 @@ protected: \details The function is implemented as an array of EvalFunction objects. */ -template -class EvalMultiFunction : public ParentFunc, public EvalFunctions +template +class EvalMultiFunction : public ParentFunc, public EvalFunctions { size_t nsd; //!< Number of spatial dimensions + using FuncType = typename EvalFunctions::FuncType; //!< Type alias for function public: //! \brief The constructor parses the expression string for each component. @@ -183,7 +208,7 @@ public: const std::string& variables = "", const Real epsX = 1e-8, const Real epsT = 1e-12) - : EvalFunctions(functions,variables,epsX,epsT), nsd(0) { this->setNoDims(); } + : EvalFunctions(functions,variables,epsX,epsT), nsd(0) { this->setNoDims(); } //! \brief Empty destructor. virtual ~EvalMultiFunction() {} @@ -191,7 +216,7 @@ public: //! \brief Returns whether the function is time-independent or not. bool isConstant() const override { - for (const std::unique_ptr& func : p) + for (const std::unique_ptr& func : this->p) if (!func->isConstant()) return false; return true; @@ -208,7 +233,7 @@ public: //! \brief Set an additional parameter in the function. void setParam(const std::string& name, double value) { - for (std::unique_ptr& func : p) + for (std::unique_ptr& func : this->p) func->setParam(name, value); } @@ -234,11 +259,13 @@ protected: //! Scalar-valued function expression using EvalFunc = EvalFuncScalar; +//! Scalar-valued spatial function expression +using EvalFunction = EvalFuncSpatial; //! Vector-valued function expression -using VecFuncExpr = EvalMultiFunction; +using VecFuncExpr = EvalMultiFunction; //! Tensor-valued function expression -using TensorFuncExpr = EvalMultiFunction; +using TensorFuncExpr = EvalMultiFunction; //! Symmetric tensor-valued function expression -using STensorFuncExpr = EvalMultiFunction; +using STensorFuncExpr = EvalMultiFunction; #endif diff --git a/src/Utility/Test/TestFunctions.C b/src/Utility/Test/TestFunctions.C index 90279aa4..6c355b0b 100644 --- a/src/Utility/Test/TestFunctions.C +++ b/src/Utility/Test/TestFunctions.C @@ -23,6 +23,13 @@ #include +using EvalFuncAd = EvalFuncScalar; +using EvalFunctionAd = EvalFuncSpatial; +using VecFuncExprAd = EvalMultiFunction; +using TensorFuncExprAd = EvalMultiFunction; +using STensorFuncExprAd = EvalMultiFunction; + + TEST(TestScalarFunc, ParseDerivative) { const char* func1 = "sin(1.5*t)*t"; @@ -31,7 +38,7 @@ TEST(TestScalarFunc, ParseDerivative) ScalarFunc* f1 = utl::parseTimeFunc(func1,"expression"); ScalarFunc* f2 = utl::parseTimeFunc(func2,"expression"); - EvalFuncScalar f3(func1,"t"); + EvalFuncAd f3(func1,"t"); ASSERT_TRUE(f1 != nullptr); ASSERT_TRUE(f2 != nullptr); @@ -56,17 +63,20 @@ 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)"; + const char* g = "sin(x)*sin(y)*sin(z)"; + const char* g_x = "cos(x)*sin(y)*sin(z)"; + const char* g_y = "sin(x)*cos(y)*sin(z)"; + const char* g_z = "sin(x)*sin(y)*cos(z)"; - EvalFunction f(f1); - f.addDerivative(d1, "", 1); - f.addDerivative(d2, "", 2); - f.addDerivative(d3, "", 3); + EvalFunction f1(g); + f1.addDerivative(g_x, "", 1); + f1.addDerivative(g_y, "", 2); + f1.addDerivative(g_z, "", 3); - EXPECT_TRUE(f.isConstant()); + EvalFunctionAd f2(g); + + EXPECT_TRUE(f1.isConstant()); + EXPECT_TRUE(f2.isConstant()); for (double x : {0.1, 0.2, 0.3}) for (double y : {0.5, 0.6, 0.7}) @@ -76,10 +86,13 @@ TEST(TestRealFunc, Gradient) 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]); + for (const RealFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const Vec3 grad = fp->gradient(X); + for (size_t i = 1; i <= 3; ++i) { + EXPECT_DOUBLE_EQ(fp->deriv(X, i), r[i-1]); + EXPECT_DOUBLE_EQ(grad[i-1], r[i-1]); + } } } } @@ -116,23 +129,26 @@ TEST(TestRealFunc, GradientFD) TEST(TestRealFunc, Hessian) { - const char* f1 = "sin(x)*sin(y)*sin(z)"; - const char* d11 = "-sin(x)*sin(y)*sin(z)"; - const char* d22 = "-sin(x)*sin(y)*sin(z)"; - const char* d33 = "-sin(x)*sin(y)*sin(z)"; - const char* d12 = "cos(x)*cos(y)*sin(z)"; - const char* d13 = "cos(x)*sin(y)*cos(z)"; - const char* d23 = "sin(x)*cos(y)*cos(z)"; + const char* g = "sin(x)*sin(y)*sin(z)"; + const char* g_xx = "-sin(x)*sin(y)*sin(z)"; + const char* g_xy = "cos(x)*cos(y)*sin(z)"; + const char* g_xz = "cos(x)*sin(y)*cos(z)"; + const char* g_yy = "-sin(x)*sin(y)*sin(z)"; + const char* g_zz = "-sin(x)*sin(y)*sin(z)"; + const char* g_yz = "sin(x)*cos(y)*cos(z)"; - EvalFunction f(f1); - f.addDerivative(d11,"",1,1); - f.addDerivative(d12,"",1,2); - f.addDerivative(d13,"",1,3); - f.addDerivative(d22,"",2,2); - f.addDerivative(d23,"",2,3); - f.addDerivative(d33,"",3,3); + EvalFunction f1(g); + f1.addDerivative(g_xx,"",1,1); + f1.addDerivative(g_xy,"",1,2); + f1.addDerivative(g_xz,"",1,3); + f1.addDerivative(g_yy,"",2,2); + f1.addDerivative(g_yz,"",2,3); + f1.addDerivative(g_zz,"",3,3); - EXPECT_TRUE(f.isConstant()); + EvalFunctionAd f2(g); + + EXPECT_TRUE(f1.isConstant()); + EXPECT_TRUE(f2.isConstant()); for (double x : {0.1, 0.2, 0.3}) for (double y : {0.5, 0.6, 0.7}) @@ -145,11 +161,14 @@ TEST(TestRealFunc, Hessian) sin(x)*cos(y)*cos(z), cos(x)*sin(y)*cos(z)}); - const SymmTensor hess = f.hessian(X); - for (size_t i = 1; i <= 3; ++i) - for (size_t j = 1; j <= 3; ++j) { - EXPECT_DOUBLE_EQ(f.dderiv(X,i,j), r(i,j)); - EXPECT_DOUBLE_EQ(hess(i,j), r(i,j)); + for (const RealFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const SymmTensor hess = fp->hessian(X); + for (size_t i = 1; i <= 3; ++i) + for (size_t j = 1; j <= 3; ++j) { + EXPECT_DOUBLE_EQ(fp->dderiv(X,i,j), r(i,j)); + EXPECT_DOUBLE_EQ(hess(i,j), r(i,j)); + } } } } @@ -159,17 +178,21 @@ TEST(TestVecFunc, Evaluate) { const char* func = "sin(x) | cos (y) | exp(z)"; - VecFuncExpr f(func); + VecFuncExpr f1(func); + VecFuncExprAd f2(func); 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 fx = f(X); const Vec3 r(sin(x), cos(y), exp(z)); - EXPECT_DOUBLE_EQ(fx.x, r.x); - EXPECT_DOUBLE_EQ(fx.y, r.y); - EXPECT_DOUBLE_EQ(fx.z, r.z); + for (const VecFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const Vec3 fx = (*fp)(X); + EXPECT_DOUBLE_EQ(fx.x, r.x); + EXPECT_DOUBLE_EQ(fx.y, r.y); + EXPECT_DOUBLE_EQ(fx.z, r.z); + } } } @@ -180,11 +203,14 @@ TEST(TestVecFunction, Gradient2D) const char* g_x = "cos(x)*sin(y) | 2*x*y*y"; const char* g_y = "sin(x)*cos(y) | 2*x*x*y"; - VecFuncExpr f(g); - f.addDerivative(g_x,"",1); - f.addDerivative(g_y,"",2); + VecFuncExpr f1(g); + f1.addDerivative(g_x,"",1); + f1.addDerivative(g_y,"",2); - EXPECT_TRUE(f.isConstant()); + VecFuncExprAd f2(g); + + EXPECT_TRUE(f1.isConstant()); + EXPECT_TRUE(f2.isConstant()); for (double x : {0.1, 0.2, 0.3}) for (double y : {0.5, 0.6, 0.7}) { @@ -192,12 +218,15 @@ TEST(TestVecFunction, Gradient2D) const Tensor r({cos(x)*sin(y), 2*x*y*y, sin(x)*cos(y), 2*x*x*y}); - const Tensor grad = f.gradient(X); - for (size_t d = 1; d <= 2; ++d) { - const Vec3 dx = f.deriv(X,d); - for (size_t i = 1; i <= 2; ++i) { - EXPECT_DOUBLE_EQ(dx[i-1], r(i,d)); - EXPECT_DOUBLE_EQ(grad(i,d), r(i,d)); + for (const VecFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const Tensor grad = fp->gradient(X); + for (size_t d = 1; d <= 2; ++d) { + const Vec3 dx = fp->deriv(X,d); + for (size_t i = 1; i <= 2; ++i) { + EXPECT_DOUBLE_EQ(dx[i-1], r(i,d)); + EXPECT_DOUBLE_EQ(grad(i,d), r(i,d)); + } } } } @@ -243,12 +272,15 @@ TEST(TestVecFunction, Gradient3D) const char* g_y = "sin(x)*cos(y)*sin(z) | 2*x*x*y*z*z*z | 2*exp(-x)*exp(2*y)*exp(z*z)"; const char* g_z = "sin(x)*sin(y)*cos(z) | 3*x*x*y*y*z*z | 2*z*exp(-x)*exp(2*y)*exp(z*z)"; - VecFuncExpr f(g); - f.addDerivative(g_x,"",1); - f.addDerivative(g_y,"",2); - f.addDerivative(g_z,"",3); + VecFuncExpr f1(g); + f1.addDerivative(g_x,"",1); + f1.addDerivative(g_y,"",2); + f1.addDerivative(g_z,"",3); - EXPECT_TRUE(f.isConstant()); + VecFuncExprAd f2(g); + + EXPECT_TRUE(f1.isConstant()); + EXPECT_TRUE(f2.isConstant()); for (double x : {0.1, 0.2, 0.3}) for (double y : {0.5, 0.6, 0.7}) @@ -258,12 +290,15 @@ TEST(TestVecFunction, Gradient3D) sin(x)*cos(y)*sin(z), 2*x*x*y*z*z*z, 2*exp(-x)*exp(2*y)*exp(z*z), sin(x)*sin(y)*cos(z), 3*x*x*y*y*z*z, 2*z*exp(-x)*exp(2*y)*exp(z*z)}); - const Tensor grad = f.gradient(X); - for (size_t d = 1; d <= 3; ++d) { - const Vec3 dx = f.deriv(X,d); - for (size_t i = 1; i <= 3; ++i) { - EXPECT_DOUBLE_EQ(dx[i-1], r(i,d)); - EXPECT_DOUBLE_EQ(grad(i,d), r(i,d)); + for (const VecFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const Tensor grad = fp->gradient(X); + for (size_t d = 1; d <= 3; ++d) { + const Vec3 dx = fp->deriv(X,d); + for (size_t i = 1; i <= 3; ++i) { + EXPECT_DOUBLE_EQ(dx[i-1], r(i,d)); + EXPECT_DOUBLE_EQ(grad(i,d), r(i,d)); + } } } } @@ -314,12 +349,15 @@ TEST(TestVecFunction, Hessian2D) const char* g_xy = "cos(x)*cos(y) | 2*x*2*y"; const char* g_yy = "-sin(x)*sin(y) | 2*x*x"; - VecFuncExpr f(g); - f.addDerivative(g_xx,"",1,1); - f.addDerivative(g_yy,"",2,2); - f.addDerivative(g_xy,"",1,2); + VecFuncExpr f1(g); + f1.addDerivative(g_xx,"",1,1); + f1.addDerivative(g_yy,"",2,2); + f1.addDerivative(g_xy,"",1,2); - EXPECT_TRUE(f.isConstant()); + VecFuncExprAd f2(g); + + EXPECT_TRUE(f1.isConstant()); + EXPECT_TRUE(f2.isConstant()); for (double x : {0.1, 0.2, 0.3}) for (double y : {0.5, 0.6, 0.7}) { @@ -330,15 +368,18 @@ TEST(TestVecFunction, Hessian2D) cos(x)*cos(y), 2*x*2*y, -sin(x)*sin(y), 2*x*x}.data()); - const utl::matrix3d hess = f.hessian(X); - for (size_t d1 = 1; d1 <= 2; ++d1) - for (size_t d2 = 1; d2 <= 2; ++d2) { - const Vec3 d2x = f.dderiv(X,d1,d2); - for (size_t i = 1; i <= 2; ++i) { - EXPECT_DOUBLE_EQ(hess(i,d1,d2), r(i,d1,d2)); - EXPECT_DOUBLE_EQ(d2x[i-1], r(i,d1,d2)); + for (const VecFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const utl::matrix3d hess = fp->hessian(X); + for (size_t d1 = 1; d1 <= 2; ++d1) + for (size_t d2 = 1; d2 <= 2; ++d2) { + const Vec3 d2x = fp->dderiv(X,d1,d2); + for (size_t i = 1; i <= 2; ++i) { + EXPECT_DOUBLE_EQ(hess(i,d1,d2), r(i,d1,d2)); + EXPECT_DOUBLE_EQ(d2x[i-1], r(i,d1,d2)); + } } - } + } } } @@ -353,15 +394,18 @@ TEST(TestVecFunction, Hessian3D) const char* g_yz = "sin(x)*cos(y)*cos(z) | x*x*2*y*2*z | exp(x)*2*exp(2*y)*3*exp(3*z)"; const char* g_zz = "-sin(x)*sin(y)*sin(z) | x*x*y*y*2 | exp(x)*exp(2*y)*9*exp(3*z)"; - VecFuncExpr f(g); - f.addDerivative(g_xx,"",1,1); - f.addDerivative(g_yy,"",2,2); - f.addDerivative(g_zz,"",3,3); - f.addDerivative(g_xy,"",1,2); - f.addDerivative(g_xz,"",1,3); - f.addDerivative(g_yz,"",2,3); + VecFuncExpr f1(g); + f1.addDerivative(g_xx,"",1,1); + f1.addDerivative(g_yy,"",2,2); + f1.addDerivative(g_zz,"",3,3); + f1.addDerivative(g_xy,"",1,2); + f1.addDerivative(g_xz,"",1,3); + f1.addDerivative(g_yz,"",2,3); - EXPECT_TRUE(f.isConstant()); + VecFuncExprAd f2(g); + + EXPECT_TRUE(f1.isConstant()); + EXPECT_TRUE(f2.isConstant()); for (double x : {0.1, 0.2, 0.3}) for (double y : {0.5, 0.6, 0.7}) @@ -380,15 +424,19 @@ TEST(TestVecFunction, Hessian3D) sin(x)*cos(y)*cos(z), x*x*2*y*2*z, exp(x)*2*exp(2*y)*3*exp(3*z), -sin(x)*sin(y)*sin(z), x*x*y*y*2, exp(x)*exp(2*y)*9*exp(3*z)}.data()); - const utl::matrix3d hess = f.hessian(X); - for (size_t d1 = 1; d1 <= 3; ++d1) - for (size_t d2 = 1; d2 <= 3; ++d2) { - const Vec3 d2x = f.dderiv(X,d1,d2); - for (size_t i = 1; i <= 3; ++i) { - EXPECT_DOUBLE_EQ(d2x[i-1], r(i,d1,d2)); - EXPECT_DOUBLE_EQ(hess(i,d1,d2), r(i,d1,d2)); + + for (const VecFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const utl::matrix3d hess = fp->hessian(X); + for (size_t d1 = 1; d1 <= 3; ++d1) + for (size_t d2 = 1; d2 <= 3; ++d2) { + const Vec3 d2x = fp->dderiv(X,d1,d2); + for (size_t i = 1; i <= 3; ++i) { + EXPECT_DOUBLE_EQ(d2x[i-1], r(i,d1,d2)); + EXPECT_DOUBLE_EQ(hess(i,d1,d2), r(i,d1,d2)); + } } - } + } } } @@ -470,17 +518,22 @@ TEST(TestTensorFunc, Evaluate) { const char* func = "sin(x) | cos (y) | exp(z) | sin(x)*cos(y)"; - TensorFuncExpr f(func); + TensorFuncExpr f1(func); + TensorFuncExprAd f2(func); 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 Tensor fx = f(X); const Tensor r({sin(x), cos(y), exp(z), sin(x)*cos(y)}); - for (size_t i = 1; i <= 2; ++i) - for (size_t j = 1; j <= 2; ++j) - EXPECT_DOUBLE_EQ(fx(i,j), r(i,j)); + + for (const TensorFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const Tensor fx = (*fp)(X); + for (size_t i = 1; i <= 2; ++i) + for (size_t j = 1; j <= 2; ++j) + EXPECT_DOUBLE_EQ(fx(i,j), r(i,j)); + } } } @@ -491,11 +544,13 @@ TEST(TestTensorFunction, Gradient2D) const char* g_x = "cos(x)*sin(y) | 2*x*y*y | exp(x)*exp(2*y) | -2*exp(-2*x)*exp(y)"; const char* g_y = "sin(x)*cos(y) | 2*x*x*y | 2*exp(x)*exp(2*y) | exp(-2*x)*exp(y)"; - TensorFuncExpr f(g); - f.addDerivative(g_x,"",1); - f.addDerivative(g_y,"",2); + TensorFuncExpr f1(g); + f1.addDerivative(g_x,"",1); + f1.addDerivative(g_y,"",2); - EXPECT_TRUE(f.isConstant()); + TensorFuncExprAd f2(g); + + EXPECT_TRUE(f1.isConstant()); for (double x : {0.1, 0.2, 0.3}) for (double y : {0.5, 0.6, 0.7}) { @@ -503,15 +558,17 @@ TEST(TestTensorFunction, Gradient2D) utl::matrix3d r(2,2,2); r.fill(std::array{cos(x)*sin(y), 2*x*y*y, exp(x)*exp(2*y), -2*exp(-2*x)*exp(y), sin(x)*cos(y), 2*x*x*y, 2*exp(x)*exp(2*y), exp(-2*x)*exp(y)}.data()); - - const utl::matrix3d grad = f.gradient(X); - for (size_t d = 1; d <= 2; ++d) { - const Tensor dx = f.deriv(X,d); - for (size_t i = 1; i <= 2; ++i) - for (size_t j = 1; j <= 2; ++j) { - EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d)); - EXPECT_DOUBLE_EQ(grad(i,j,d), r(i,j,d)); - } + for (const TensorFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const utl::matrix3d grad = fp->gradient(X); + for (size_t d = 1; d <= 2; ++d) { + const Tensor dx = fp->deriv(X,d); + for (size_t i = 1; i <= 2; ++i) + for (size_t j = 1; j <= 2; ++j) { + EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d)); + EXPECT_DOUBLE_EQ(grad(i,j,d), r(i,j,d)); + } + } } } } @@ -570,12 +627,15 @@ TEST(TestTensorFunction, Gradient3D) "exp(-2*x)*exp(y) | x*y | x*y*2*z |" "0.0 | 0.0 | 1.0"; - TensorFuncExpr f(g); - f.addDerivative(g_x,"",1); - f.addDerivative(g_y,"",2); - f.addDerivative(g_z,"",3); + TensorFuncExpr f1(g); + f1.addDerivative(g_x,"",1); + f1.addDerivative(g_y,"",2); + f1.addDerivative(g_z,"",3); - EXPECT_TRUE(f.isConstant()); + TensorFuncExprAd f2(g); + + EXPECT_TRUE(f1.isConstant()); + EXPECT_TRUE(f2.isConstant()); for (double x : {0.1, 0.2, 0.3}) for (double y : {0.5, 0.6, 0.7}) @@ -594,14 +654,17 @@ TEST(TestTensorFunction, Gradient3D) exp(-2*x)*exp(y), x*y, x*y*2*z, 0.0, 0.0, 1.0}.data()); - const utl::matrix3d grad = f.gradient(X); - for (size_t d = 1; d <= 3; ++d) { - const Tensor dx = f.deriv(X,d); - for (size_t i = 1; i <= 3; ++i) - for (size_t j = 1; j <= 3; ++j) { - EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d)); - EXPECT_DOUBLE_EQ(grad(i,j,d), r(i,j,d)); - } + for (const TensorFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const utl::matrix3d grad = fp->gradient(X); + for (size_t d = 1; d <= 3; ++d) { + const Tensor dx = fp->deriv(X,d); + for (size_t i = 1; i <= 3; ++i) + for (size_t j = 1; j <= 3; ++j) { + EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d)); + EXPECT_DOUBLE_EQ(grad(i,j,d), r(i,j,d)); + } + } } } } @@ -674,12 +737,15 @@ TEST(TestTensorFunction, Hessian2D) const char* g_xy = "cos(x)*cos(y) | 2*x*2*y | exp(x)*2*exp(2*y) | -2*exp(-2*x)*exp(y)"; const char* g_yy = "-sin(x)*sin(y) | 2*x*x | exp(x)*4*exp(2*y) | exp(-2*x)*exp(y)"; - TensorFuncExpr f(g); - f.addDerivative(g_xx,"",1,1); - f.addDerivative(g_yy,"",2,2); - f.addDerivative(g_xy,"",1,2); + TensorFuncExpr f1(g); + f1.addDerivative(g_xx,"",1,1); + f1.addDerivative(g_yy,"",2,2); + f1.addDerivative(g_xy,"",1,2); - EXPECT_TRUE(f.isConstant()); + TensorFuncExprAd f2(g); + + EXPECT_TRUE(f1.isConstant()); + EXPECT_TRUE(f2.isConstant()); for (double x : {0.1, 0.2, 0.3}) for (double y : {0.5, 0.6, 0.7}) { @@ -698,16 +764,19 @@ TEST(TestTensorFunction, Hessian2D) -sin(x)*sin(y), 2*x*x, exp(x)*4*exp(2*y), exp(-2*x)*exp(y)}.data()); - const utl::matrix4d hess = f.hessian(X); - for (size_t d1 = 1; d1 <= 2; ++d1) - for (size_t d2 = 1; d2 <= 2; ++d2) { - const Tensor dx = f.dderiv(X,d1,d2); - for (size_t i = 1; i <= 2; ++i) - for (size_t j = 1; j <= 2; ++j) { - EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d1,d2)); - EXPECT_DOUBLE_EQ(hess(i,j,d1,d2), r(i,j,d1,d2)); - } - } + for (const TensorFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const utl::matrix4d hess = fp->hessian(X); + for (size_t d1 = 1; d1 <= 2; ++d1) + for (size_t d2 = 1; d2 <= 2; ++d2) { + const Tensor dx = fp->dderiv(X,d1,d2); + for (size_t i = 1; i <= 2; ++i) + for (size_t j = 1; j <= 2; ++j) { + EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d1,d2)); + EXPECT_DOUBLE_EQ(hess(i,j,d1,d2), r(i,j,d1,d2)); + } + } + } } } @@ -736,15 +805,18 @@ TEST(TestTensorFunction, Hessian3D) "0.0 | 0.0 | 2.0*x*y |" "0.0 | 0.0 | 2.0"; - TensorFuncExpr f(g); - f.addDerivative(g_xx,"",1,1); - f.addDerivative(g_yy,"",2,2); - f.addDerivative(g_zz,"",3,3); - f.addDerivative(g_xy,"",1,2); - f.addDerivative(g_xz,"",1,3); - f.addDerivative(g_yz,"",2,3); + TensorFuncExpr f1(g); + f1.addDerivative(g_xx,"",1,1); + f1.addDerivative(g_yy,"",2,2); + f1.addDerivative(g_zz,"",3,3); + f1.addDerivative(g_xy,"",1,2); + f1.addDerivative(g_xz,"",1,3); + f1.addDerivative(g_yz,"",2,3); - EXPECT_TRUE(f.isConstant()); + TensorFuncExprAd f2(g); + + EXPECT_TRUE(f1.isConstant()); + EXPECT_TRUE(f2.isConstant()); for (double x : {0.1, 0.2, 0.3}) for (double y : {0.5, 0.6, 0.7}) @@ -787,16 +859,19 @@ TEST(TestTensorFunction, Hessian3D) 0.0, 0.0, 2.0*x*y, 0.0, 0.0, 2.0}.data()); - const utl::matrix4d hess = f.hessian(X); - for (size_t d1 = 1; d1 <= 3; ++d1) - for (size_t d2 = 1; d2 <= 3; ++d2) { - const Tensor dx = f.dderiv(X,d1,d2); - for (size_t i = 1; i <= 3; ++i) - for (size_t j = 1; j <= 3; ++j) { - EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d1,d2)); - EXPECT_DOUBLE_EQ(hess(i,j,d1,d2), r(i,j,d1,d2)); - } - } + for (const TensorFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const utl::matrix4d hess = fp->hessian(X); + for (size_t d1 = 1; d1 <= 3; ++d1) + for (size_t d2 = 1; d2 <= 3; ++d2) { + const Tensor dx = fp->dderiv(X,d1,d2); + for (size_t i = 1; i <= 3; ++i) + for (size_t j = 1; j <= 3; ++j) { + EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d1,d2)); + EXPECT_DOUBLE_EQ(hess(i,j,d1,d2), r(i,j,d1,d2)); + } + } + } } } @@ -854,17 +929,24 @@ TEST(TestSTensorFunc, Evaluate2D) { const char* func = "sin(x) | cos (y) | exp(z)"; - STensorFuncExpr f(func); + STensorFuncExpr f1(func); + STensorFuncExprAd f2(func); + + EXPECT_TRUE(f1.isConstant()); + EXPECT_TRUE(f2.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 SymmTensor fx = f(X); const Tensor r({sin(x), exp(z), exp(z), cos(y)}); - for (size_t i = 1; i <= 2; ++i) - for (size_t j = 1; j <= 2; ++j) - EXPECT_DOUBLE_EQ(fx(i,j), r(i,j)); + for (const STensorFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const SymmTensor fx = (*fp)(X); + for (size_t i = 1; i <= 2; ++i) + for (size_t j = 1; j <= 2; ++j) + EXPECT_DOUBLE_EQ(fx(i,j), r(i,j)); + } } } @@ -894,19 +976,23 @@ TEST(TestSTensorFunc, Evaluate3D) const char* func = "sin(x) | cos (y) | exp(z) |" "sin(x)*sin(y) | sin(x)*cos(y) | exp(x)*exp(y)"; - STensorFuncExpr f(func); + STensorFuncExpr f1(func); + STensorFuncExprAd f2(func); 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 SymmTensor fx = f(X); const Tensor r({sin(x), sin(x)*sin(y), exp(x)*exp(y), sin(x)*sin(y), cos(y), sin(x)*cos(y), exp(x)*exp(y), sin(x)*cos(y), exp(z)}); - for (size_t i = 1; i <= 3; ++i) - for (size_t j = 1; j <= 3; ++j) - EXPECT_DOUBLE_EQ(fx(i,j), r(i,j)); + for (const STensorFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const SymmTensor fx = (*fp)(X); + for (size_t i = 1; i <= 3; ++i) + for (size_t j = 1; j <= 3; ++j) + EXPECT_DOUBLE_EQ(fx(i,j), r(i,j)); + } } } @@ -917,11 +1003,14 @@ TEST(TestSTensorFunction, Gradient2D) const char* g_x = "cos(x)*sin(y) | exp(x)*exp(2*y) | 2*x*y*y"; const char* g_y = "sin(x)*cos(y) | 2*exp(x)*exp(2*y) | 2*x*x*y"; - STensorFuncExpr f(g); - f.addDerivative(g_x,"",1); - f.addDerivative(g_y,"",2); + STensorFuncExpr f1(g); + f1.addDerivative(g_x,"",1); + f1.addDerivative(g_y,"",2); - EXPECT_TRUE(f.isConstant()); + STensorFuncExprAd f2(g); + + EXPECT_TRUE(f1.isConstant()); + EXPECT_TRUE(f2.isConstant()); for (double x : {0.1, 0.2, 0.3}) for (double y : {0.5, 0.6, 0.7}) { @@ -930,14 +1019,17 @@ TEST(TestSTensorFunction, Gradient2D) r.fill(std::array{cos(x)*sin(y), 2*x*y*y, 2*x*y*y, exp(x)*exp(2*y), sin(x)*cos(y), 2*x*x*y, 2*x*x*y, 2*exp(x)*exp(2*y)}.data()); - const utl::matrix3d grad = f.gradient(X); - for (size_t d = 1; d <= 2; ++d) { - const SymmTensor dx = f.deriv(X,d); - for (size_t i = 1; i <= 2; ++i) - for (size_t j = 1; j <= 2; ++j) { - EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d)); - EXPECT_DOUBLE_EQ(grad(i,j,d), r(i,j,d)); - } + for (const STensorFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const utl::matrix3d grad = fp->gradient(X); + for (size_t d = 1; d <= 2; ++d) { + const SymmTensor dx = fp->deriv(X,d); + for (size_t i = 1; i <= 2; ++i) + for (size_t j = 1; j <= 2; ++j) { + EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d)); + EXPECT_DOUBLE_EQ(grad(i,j,d), r(i,j,d)); + } + } } } } @@ -1022,12 +1114,15 @@ TEST(TestSTensorFunction, Gradient3D) const char* g_z = "sin(x)*sin(y)*cos(z) | exp(x)*exp(2*y)*exp(z) | x*x*y*y*2*z |" "x*y | x*x*y | 1"; - STensorFuncExpr f(g); - f.addDerivative(g_x,"",1); - f.addDerivative(g_y,"",2); - f.addDerivative(g_z,"",3); + STensorFuncExpr f1(g); + f1.addDerivative(g_x,"",1); + f1.addDerivative(g_y,"",2); + f1.addDerivative(g_z,"",3); - EXPECT_TRUE(f.isConstant()); + STensorFuncExprAd f2(g); + + EXPECT_TRUE(f1.isConstant()); + EXPECT_TRUE(f2.isConstant()); for (double x : {0.1, 0.2, 0.3}) for (double y : {0.5, 0.6, 0.7}) @@ -1046,14 +1141,17 @@ TEST(TestSTensorFunction, Gradient3D) x*y, exp(x)*exp(2*y)*exp(z), x*x*y, 1.0, x*x*y, x*x*y*y*2*z}.data()); - const utl::matrix3d grad = f.gradient(X); - for (size_t d = 1; d <= 3; ++d) { - const SymmTensor dx = f.deriv(X,d); - for (size_t i = 1; i <= 3; ++i) - for (size_t j = 1; j <= 3; ++j) { - EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d)); - EXPECT_DOUBLE_EQ(grad(i,j,d), r(i,j,d)); - } + for (const STensorFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const utl::matrix3d grad = fp->gradient(X); + for (size_t d = 1; d <= 3; ++d) { + const SymmTensor dx = fp->deriv(X,d); + for (size_t i = 1; i <= 3; ++i) + for (size_t j = 1; j <= 3; ++j) { + EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d)); + EXPECT_DOUBLE_EQ(grad(i,j,d), r(i,j,d)); + } + } } } } @@ -1145,12 +1243,15 @@ TEST(TestSTensorFunction, Hessian2D) const char* g_xy = "cos(x)*cos(y) | exp(x)*2*exp(2*y) | 2*x*2*y"; const char* g_yy = "-sin(x)*sin(y) | exp(x)*4*exp(2*y) | x*x*2"; - STensorFuncExpr f(g); - f.addDerivative(g_xx,"",1,1); - f.addDerivative(g_xy,"",1,2); - f.addDerivative(g_yy,"",2,2); + STensorFuncExpr f1(g); + f1.addDerivative(g_xx,"",1,1); + f1.addDerivative(g_xy,"",1,2); + f1.addDerivative(g_yy,"",2,2); - EXPECT_TRUE(f.isConstant()); + STensorFuncExprAd f2(g); + + EXPECT_TRUE(f1.isConstant()); + EXPECT_TRUE(f2.isConstant()); for (double x : {0.1, 0.2, 0.3}) for (double y : {0.5, 0.6, 0.7}) { @@ -1168,16 +1269,19 @@ TEST(TestSTensorFunction, Hessian2D) -sin(x)*sin(y), x*x*2, x*x*2, exp(x)*4*exp(2*y)}.data()); - const utl::matrix4d hess = f.hessian(X); - for (size_t d1 = 1; d1 <= 2; ++d1) - for (size_t d2 = 1; d2 <= 2; ++d2) { - const SymmTensor dx = f.dderiv(X,d1,d2); - for (size_t i = 1; i <= 2; ++i) - for (size_t j = 1; j <= 2; ++j) { - EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d1,d2)); - EXPECT_DOUBLE_EQ(hess(i,j,d1,d2), r(i,j,d1,d2)); - } - } + for (const STensorFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const utl::matrix4d hess = fp->hessian(X); + for (size_t d1 = 1; d1 <= 2; ++d1) + for (size_t d2 = 1; d2 <= 2; ++d2) { + const SymmTensor dx = fp->dderiv(X,d1,d2); + for (size_t i = 1; i <= 2; ++i) + for (size_t j = 1; j <= 2; ++j) { + EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d1,d2)); + EXPECT_DOUBLE_EQ(hess(i,j,d1,d2), r(i,j,d1,d2)); + } + } + } } } @@ -1200,15 +1304,18 @@ TEST(TestSTensorFunction, Hessian3D) "0.0 | 0.0 | 0.0"; - STensorFuncExpr f(g); - f.addDerivative(g_xx,"",1,1); - f.addDerivative(g_xy,"",1,2); - f.addDerivative(g_xz,"",1,3); - f.addDerivative(g_yy,"",2,2); - f.addDerivative(g_yz,"",2,3); - f.addDerivative(g_zz,"",3,3); + STensorFuncExpr f1(g); + f1.addDerivative(g_xx,"",1,1); + f1.addDerivative(g_xy,"",1,2); + f1.addDerivative(g_xz,"",1,3); + f1.addDerivative(g_yy,"",2,2); + f1.addDerivative(g_yz,"",2,3); + f1.addDerivative(g_zz,"",3,3); - EXPECT_TRUE(f.isConstant()); + STensorFuncExprAd f2(g); + + EXPECT_TRUE(f1.isConstant()); + EXPECT_TRUE(f2.isConstant()); for (double x : {0.1, 0.2, 0.3}) for (double y : {0.5, 0.6, 0.7}) @@ -1251,17 +1358,20 @@ TEST(TestSTensorFunction, Hessian3D) 0.0, exp(x)*exp(2*y)*exp(z), 0.0, 0.0, 0.0, x*x*y*y*2}.data()); - const utl::matrix4d hess = f.hessian(X); - for (size_t d1 = 1; d1 <= 3; ++d1) - for (size_t d2 = 1; d2 <= 3; ++d2) { - const SymmTensor dx = f.dderiv(X,d1,d2); - for (size_t i = 1; i <= 3; ++i) - for (size_t j = 1; j <= 3; ++j) { - EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d1,d2)); - EXPECT_DOUBLE_EQ(hess(i,j,d1,d2), r(i,j,d1,d2)); - } - } - } + for (const STensorFunc* fp : {static_cast(&f1), + static_cast(&f2)}) { + const utl::matrix4d hess = fp->hessian(X); + for (size_t d1 = 1; d1 <= 3; ++d1) + for (size_t d2 = 1; d2 <= 3; ++d2) { + const SymmTensor dx = fp->dderiv(X,d1,d2); + for (size_t i = 1; i <= 3; ++i) + for (size_t j = 1; j <= 3; ++j) { + EXPECT_DOUBLE_EQ(dx(i,j), r(i,j,d1,d2)); + EXPECT_DOUBLE_EQ(hess(i,j,d1,d2), r(i,j,d1,d2)); + } + } + } + } } @@ -1384,5 +1494,5 @@ TEST(TestEvalFunction, Derivatives) EXPECT_DOUBLE_EQ(f.dderiv(X,3,1), cos(x)*sin(y)*cos(z)*sin(t)); EXPECT_DOUBLE_EQ(f.dderiv(X,3,2), sin(x)*cos(y)*cos(z)*sin(t)); EXPECT_DOUBLE_EQ(f.dderiv(X,3,3), -sin(x)*sin(y)*sin(z)*sin(t)); - } + } }