From dbcbede16ce3562957111e168325311b3c87ac66 Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Thu, 21 Sep 2023 12:26:52 +0200 Subject: [PATCH] added: (Tensor|Vec)Func::Gradient returns the gradient of the function --- src/Utility/ExprFunctions.C | 18 ++++ src/Utility/ExprFunctions.h | 3 + src/Utility/Function.h | 11 ++ src/Utility/TensorFunction.C | 52 +++++++++ src/Utility/TensorFunction.h | 10 ++ src/Utility/Test/TestFunctions.C | 175 +++++++++++++++++++++++++++++-- 6 files changed, 261 insertions(+), 8 deletions(-) create mode 100644 src/Utility/TensorFunction.C diff --git a/src/Utility/ExprFunctions.C b/src/Utility/ExprFunctions.C index c3c3e28e..99962b66 100644 --- a/src/Utility/ExprFunctions.C +++ b/src/Utility/ExprFunctions.C @@ -481,6 +481,24 @@ Ret EvalMultiFunction::dderiv (const Vec3& X, int d1, int d2) co } +template +std::vector +EvalMultiFunction::evalGradient (const Vec3& X) const +{ + std::vector result(this->ncmp*this->nsd); + std::vector dx; + for (const std::unique_ptr& f : this->p) + dx.push_back(f->gradient(X)); + + size_t k = 0; + for (size_t d = 1; d <= this->nsd; ++d) + for (size_t i = 1; i <= this->ncmp; ++i) + result[k++] = dx[i-1][d-1]; + + return result; +} + + template class EvalMultiFunction; template class EvalMultiFunction; template class EvalMultiFunction; diff --git a/src/Utility/ExprFunctions.h b/src/Utility/ExprFunctions.h index 9ebed7a6..edcdcf1a 100644 --- a/src/Utility/ExprFunctions.h +++ b/src/Utility/ExprFunctions.h @@ -223,6 +223,9 @@ protected: //! \brief Evaluates the function expressions. Ret evaluate(const Vec3& X) const override; + + //! \brief Returns the gradient of the function as a 1D array. + std::vector evalGradient(const Vec3& X) const override; }; //! Vector-valued function expression diff --git a/src/Utility/Function.h b/src/Utility/Function.h index 1174b06b..b3c9e112 100644 --- a/src/Utility/Function.h +++ b/src/Utility/Function.h @@ -110,6 +110,9 @@ namespace utl virtual Result dderiv(const Vec3&, int, int) const { return zero; } protected: + //! \brief Returns the gradient of the function as a 1D array. + virtual std::vector evalGradient(const Vec3&) const { return {}; } + Result zero; //!< Return value for default implementations of derivatives }; } @@ -251,6 +254,14 @@ public: { return this->evaluate(X).length(); } + + //! \brief Evaluates first derivatives of the function. + Tensor gradient(const Vec3& X) const + { + Tensor result(ncmp); + result = this->evalGradient(X); + return result; + } }; diff --git a/src/Utility/TensorFunction.C b/src/Utility/TensorFunction.C new file mode 100644 index 00000000..24e7acac --- /dev/null +++ b/src/Utility/TensorFunction.C @@ -0,0 +1,52 @@ +// $Id$ +//============================================================================== +//! +//! \file TensorFunction.C +//! +//! \date May 10 2017 +//! +//! \author Knut Morten Okstad / SINTEF +//! +//! \brief Spatial tensor-valued functions. +//! +//============================================================================== + +#include "TensorFunction.h" + + +utl::matrix3d TensorFunc::gradient(const Vec3& X) const +{ + const size_t nsd = sqrt(ncmp); + utl::matrix3d result(nsd, nsd, nsd); + result.fill(this->evalGradient(X).data()); + return result; +} + + +size_t STensorFunc::index(size_t nsd, size_t i, size_t j) const +{ + if (i == j) + return i-1; // diagonal term + else if (nsd == 2) + return ncmp-1; // off-diagonal term (2D) + + if (i == j+1 || i+2 == j) + std::swap(i,j); + + return i+2; // upper triangular term (3D) +} + + +utl::matrix3d STensorFunc::gradient(const Vec3& X) const +{ + const size_t nsd = ncmp > 5 ? 3 : (ncmp > 2 ? 2 : 1); + utl::matrix3d result(nsd,nsd,nsd); + const std::vector temp = this->evalGradient(X); + + for (size_t d = 1; d <= nsd; ++d) + for (size_t i = 1; i <= nsd; ++i) + for (size_t j = 1; j <= nsd; ++j) + result(i,j,d) = temp[index(nsd,i,j) + (d-1)*ncmp]; + + return result; +} diff --git a/src/Utility/TensorFunction.h b/src/Utility/TensorFunction.h index af1c2edd..adc4dcd6 100644 --- a/src/Utility/TensorFunction.h +++ b/src/Utility/TensorFunction.h @@ -15,6 +15,7 @@ #define _TENSOR_FUNCTION_H #include "Function.h" +#include "matrixnd.h" #include "Tensor.h" @@ -49,6 +50,9 @@ public: { return this->evaluate(X).trace(); } + + //! \brief Evaluates first derivatives of the function. + utl::matrix3d gradient(const Vec3& X) const; }; @@ -59,6 +63,9 @@ public: class STensorFunc : public utl::SpatialFunction, public FunctionBase { + //! \brief Returns the flat indices of the symmetric tensor. + size_t index(size_t nsd, size_t i, size_t j) const; + protected: //! \brief The constructor is protected to allow sub-class instances only. STensorFunc(size_t n = 0, bool with33 = false) @@ -85,6 +92,9 @@ public: { return this->evaluate(X).trace(); } + + //! \brief Evaluates first derivatives of the function. + utl::matrix3d gradient(const Vec3& X) const; }; #endif diff --git a/src/Utility/Test/TestFunctions.C b/src/Utility/Test/TestFunctions.C index 859965b6..9539190f 100644 --- a/src/Utility/Test/TestFunctions.C +++ b/src/Utility/Test/TestFunctions.C @@ -11,6 +11,7 @@ //============================================================================== #include "ExprFunctions.h" +#include "TensorFunction.h" #include "Functions.h" #include "matrixnd.h" @@ -19,6 +20,7 @@ #include #include + TEST(TestScalarFunc, ParseDerivative) { const char* func1 = "sin(1.5*t)*t"; @@ -167,6 +169,139 @@ TEST(TestVecFunc, Evaluate) } +TEST(TestVecFunction, Gradient2D) +{ + const char* g = "sin(x)*sin(y) | x*x*y*y"; + 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); + + EXPECT_TRUE(f.isConstant()); + + for (double x : {0.1, 0.2, 0.3}) + for (double y : {0.5, 0.6, 0.7}) { + const Vec3 X(x,y); + 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)); + } + } + } +} + + +TEST(TestVecFunction, Gradient2DFD) +{ + const char* g = "sin(x)*sin(y) | x*x*y*y"; + + const double eps = 1e-6; + + VecFuncExpr f(g,"",eps); + + EXPECT_TRUE(f.isConstant()); + + for (double x : {0.1, 0.2, 0.3}) + for (double y : {0.5, 0.6, 0.7}) { + const Vec3 X(x,y); + const Vec3 Xp(x + 0.5*eps, y + 0.5*eps); + const Vec3 Xm(x - 0.5*eps, y - 0.5*eps); + const Tensor r({(sin(Xp.x) - sin(Xm.x))*sin(y) / eps, + (Xp.x*Xp.x - Xm.x*Xm.x)*y*y / eps, + sin(x)*(sin(Xp.y) - sin(Xm.y)) / eps, + x*x*(Xp.y*Xp.y - Xm.y*Xm.y) / eps}); + + 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_NEAR(dx[i-1], r(i,d), 1e-8); + EXPECT_NEAR(grad(i,d), r(i,d), 1e-8); + } + } + } +} + + +TEST(TestVecFunction, Gradient3D) +{ + const char* g = "sin(x)*sin(y)*sin(z) | x*x*y*y*z*z*z | exp(-x)*exp(2*y)*exp(z*z)"; + const char* g_x = "cos(x)*sin(y)*sin(z) | 2*x*y*y*z*z*z | -exp(-x)*exp(2*y)*exp(z*z)"; + 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); + + 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 Tensor r({cos(x)*sin(y)*sin(z), 2*x*y*y*z*z*z, -exp(-x)*exp(2*y)*exp(z*z), + 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)); + } + } + } +} + + +TEST(TestVecFunction, Gradient3DFD) +{ + const char* g = "sin(x)*sin(y)*sin(z) | x*x*y*y*z*z*z | exp(-x)*exp(2*y)*exp(z*z)"; + + const double eps = 1e-6; + VecFuncExpr f(g,"",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 Tensor r({(sin(Xp.x) - sin(Xm.x))*sin(y)*sin(z) / eps, + (Xp.x*Xp.x - Xm.x*Xm.x)*X.y*X.y*X.z*X.z*X.z / eps, + (exp(-Xp.x) - exp(-Xm.x))*exp(2*y)*exp(z*z) / eps, + sin(x)*(sin(Xp.y) - sin(Xm.y))*sin(z) / eps, + x*x*(Xp.y*Xp.y - Xm.y*Xm.y)*z*z*z / eps, + exp(-x)*(exp(2*Xp.y) - exp(2*Xm.y))*exp(z*z) / eps, + sin(x)*sin(y)*(sin(Xp.z) - sin(Xm.z)) / eps, + x*x*y*y*(Xp.z*Xp.z*Xp.z - Xm.z*Xm.z*Xm.z) / eps, + exp(-x)*exp(2*y)*(exp(Xp.z*Xp.z) - exp(Xm.z*Xm.z)) / eps}); + + 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_NEAR(dx[i-1], r(i,d), 1e-8); + EXPECT_NEAR(grad(i,d), r(i,d), 1e-8); + } + } + } +} + + TEST(TestVecFuncExpr, NumDimensions) { const char* func1 = "x"; @@ -223,11 +358,14 @@ TEST(TestTensorFunction, Gradient2D) 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) + 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)); + } } } } @@ -258,11 +396,14 @@ TEST(TestTensorFunction, Gradient2DFD) exp(-2*x)*(exp(Xp.y) - exp(Xm.y))}.data()); r.multiply(1.0 / eps); + 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) + for (size_t j = 1; j <= 2; ++j) { EXPECT_NEAR(dx(i,j), r(i,j,d), 1e-8); + EXPECT_NEAR(grad(i,j,d), r(i,j,d), 1e-8); + } } } } @@ -305,11 +446,14 @@ TEST(TestTensorFunction, Gradient3D) exp(-2*x)*exp(y), y*z, 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) + 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)); + } } } } @@ -362,11 +506,14 @@ TEST(TestTensorFunction, Gradient3DFD) Xp.z - Xm.z}.data()); r.multiply(1.0 / eps); + 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) + for (size_t j = 1; j <= 3; ++j) { EXPECT_NEAR(dx(i,j), r(i,j,d), 1e-8); + EXPECT_NEAR(grad(i,j,d), r(i,j,d), 1e-8); + } } } } @@ -598,11 +745,14 @@ 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) + 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)); + } } } } @@ -663,11 +813,14 @@ TEST(TestSTensorFunction, Gradient2DFD) exp(x)*(exp(2*Xp.y) - exp(2*Xm.y))}.data()); r.multiply(1.0 / eps); + 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) + for (size_t j = 1; j <= 2; ++j) { EXPECT_NEAR(dx(i,j), r(i,j,d), 1e-8); + EXPECT_NEAR(grad(i,j,d), r(i,j,d), 1e-8); + } } } } @@ -708,11 +861,14 @@ 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) + 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)); + } } } } @@ -766,11 +922,14 @@ TEST(TestSTensorFunction, Gradient3DFD) x*x*y*y*(Xp.z*Xp.z - Xm.z*Xm.z)}.data()); r.multiply(1.0 / eps); + 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) + for (size_t j = 1; j <= 3; ++j) { EXPECT_NEAR(dx(i,j), r(i,j,d), 1e-8); + EXPECT_NEAR(grad(i,j,d), r(i,j,d), 1e-8); + } } } }