OilfieldToolsNet/IFEM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

changed: use a common implementation for evaluate()

Browse Source 
this transposes tensor input to column-major order
adjust TensorFuncExpr::deriv / TensorFuncExpr::dderiv accordingly
This commit is contained in:
Arne Morten Kvarving 2023-09-20 16:03:36 +02:00
parent de3f5b4930
commit 48172df260
3 changed files with 612 additions and 49 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

57
src/Utility/ExprFunctions.C
View File

@ -391,15 +391,14 @@ void EvalFunctions::addDerivative (const std::string& functions,
}
template<>
Vec3 VecFuncExpr::evaluate (const Vec3& X) const
template <class ParentFunc, class Ret>
Ret EvalMultiFunction<ParentFunc,Ret>::evaluate (const Vec3& X) const
{
Vec3 result;
std::vector<Real> res_array(this->p.size());
for (size_t i = 0; i < this->p.size(); ++i)
res_array[i] = (*this->p[i])(X);
for (size_t i = 0; i < 3 && i < nsd; ++i)
result[i] = (*p[i])(X);
return result;
return Ret(res_array);
}
@ -441,28 +440,14 @@ void TensorFuncExpr::setNoDims ()
}
template<>
Tensor TensorFuncExpr::evaluate (const Vec3& X) const
{
Tensor sigma(nsd);
size_t i, j, k = 0;
for (i = 1; i <= nsd; ++i)
for (j = 1; j <= nsd; ++j)
sigma(i,j) = (*p[k++])(X);
return sigma;
}
template<>
Tensor TensorFuncExpr::deriv (const Vec3& X, int dir) const
{
Tensor sigma(nsd);
size_t i, j, k = 0;
for (i = 1; i <= nsd; ++i)
for (j = 1; j <= nsd; ++j)
size_t k = 0;
for (size_t j = 1; j <= nsd; ++j)
for (size_t i = 1; i <= nsd; ++i)
sigma(i,j) = p[k++]->deriv(X,dir);
return sigma;
@ -474,9 +459,9 @@ Tensor TensorFuncExpr::dderiv (const Vec3& X, int d1, int d2) const
{
Tensor sigma(nsd);
size_t i, j, k = 0;
for (i = 1; i <= nsd; ++i)
for (j = 1; j <= nsd; ++j)
size_t k = 0;
for (size_t j = 1; j <= nsd; ++j)
for (size_t i = 1; i <= nsd; ++i)
sigma(i,j) = p[k++]->dderiv(X,d1,d2);
return sigma;
@ -497,19 +482,6 @@ void STensorFuncExpr::setNoDims ()
}
template<>
SymmTensor STensorFuncExpr::evaluate (const Vec3& X) const
{
SymmTensor sigma(nsd,p.size()==4);
std::vector<Real>& svec = sigma;
for (size_t i = 0; i < svec.size(); i++)
svec[i] = (*p[i])(X);
return sigma;
}
template<>
SymmTensor STensorFuncExpr::deriv (const Vec3& X, int dir) const
{
@ -534,3 +506,8 @@ SymmTensor STensorFuncExpr::dderiv (const Vec3& X, int d1, int d2) const
return sigma;
}
template Vec3 VecFuncExpr::evaluate(const Vec3&) const;
template Tensor TensorFuncExpr::evaluate(const Vec3&) const;
template SymmTensor STensorFuncExpr::evaluate(const Vec3&) const;

7
src/Utility/ExprFunctions.h
View File

@ -229,14 +229,9 @@ using TensorFuncExpr = EvalMultiFunction<TensorFunc,Tensor>;
//! Symmetric tensor-valued function expression
using STensorFuncExpr = EvalMultiFunction<STensorFunc,SymmTensor>;
//! \brief Specialization for vector functions.
template<> Vec3 VecFuncExpr::evaluate(const Vec3& X) const;
//! \brief Specialization for tensor functions.
template<> void TensorFuncExpr::setNoDims();
//! \brief Specialization for tensor functions.
template<> Tensor TensorFuncExpr::evaluate(const Vec3& X) const;
//! \brief Specialization for tensor functions.
template<> Tensor TensorFuncExpr::deriv(const Vec3& X, int dir) const;
//! \brief Specialization for tensor functions.
template<> Tensor TensorFuncExpr::dderiv(const Vec3& X, int d1, int d2) const;
@ -244,8 +239,6 @@ template<> Tensor TensorFuncExpr::dderiv(const Vec3& X, int d1, int d2) const;
//! \brief Specialization for symmetric tensor functions.
template<> void STensorFuncExpr::setNoDims();
//! \brief Specialization for symmetric tensor functions.
template<> SymmTensor STensorFuncExpr::evaluate(const Vec3& X) const;
//! \brief Specialization for symmetric tensor functions.
template<> SymmTensor STensorFuncExpr::deriv(const Vec3& X, int dir) const;
//! \brief Specialization for symmetric tensor functions.
template<> SymmTensor STensorFuncExpr::dderiv(const Vec3& X,

597
src/Utility/Test/TestFunctions.C
View File

@ -12,11 +12,12 @@
#include "ExprFunctions.h"
#include "Functions.h"
#include <cstdlib>
#include <cmath>
#include "matrixnd.h"
#include "gtest/gtest.h"
#include <cstdlib>
#include <cmath>
TEST(TestScalarFunc, ParseDerivative)
{
@ -46,6 +47,598 @@ TEST(TestScalarFunc, ParseDerivative)
}
TEST(TestVecFunc, Evaluate)
{
const char* func = "sin(x) | cos (y) | exp(z)";
VecFuncExpr f(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);
}
}
TEST(TestTensorFunc, Evaluate)
{
const char* func = "sin(x) | cos (y) | exp(z) | sin(x)*cos(y)";
TensorFuncExpr f(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));
}
}
TEST(TestTensorFunction, Gradient2D)
{
const char* g = "sin(x)*sin(y) | x*x*y*y | exp(x)*exp(2*y) | exp(-2*x)*exp(y)";
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);
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);
utl::matrix3d<Real> 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());
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));
}
}
}
TEST(TestTensorFunction, Gradient2DFD)
{
const char* g = "sin(x)*sin(y) | x*x*y*y | exp(x)*exp(2*y) | exp(-2*x)*exp(y)";
const double eps = 1e-6;
TensorFuncExpr 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);
utl::matrix3d<Real> r(2,2,2);
r.fill(std::array{(sin(Xp.x) - sin(Xm.x))*sin(y),
(Xp.x*Xp.x - Xm.x*Xm.x)*y*y,
(exp(Xp.x) - exp(Xm.x))*exp(2*y),
(exp(-2*Xp.x) - exp(-2*Xm.x))*exp(y),
sin(x)*(sin(Xp.y) - sin(Xm.y)),
x*x*(Xp.y*Xp.y - Xm.y*Xm.y),
exp(x)*(exp(2*Xp.y) - exp(2*Xm.y)),
exp(-2*x)*(exp(Xp.y) - exp(Xm.y))}.data());
r.multiply(1.0 / eps);
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_NEAR(dx(i,j), r(i,j,d), 1e-8);
}
}
}
TEST(TestTensorFunction, Gradient3D)
{
const char* g = "sin(x)*sin(y)*sin(z) | x*x*y*y*z | exp(x)*exp(2*y)*z*z |"
"exp(-2*x)*exp(y)*z | x*y*z | x*y*z*z |"
"x | y | z";
const char* g_x = "cos(x)*sin(y)*sin(z) | 2*x*y*y*z | exp(x)*exp(2*y)*z*z |"
"-2*exp(-2*x)*exp(y)*z | y*z | y*z*z |"
"1.0 | 0.0 | 0.0";
const char* g_y = "sin(x)*cos(y)*sin(z) | 2*x*x*y*z | 2*exp(x)*exp(2*y)*z*z |"
"exp(-2*x)*exp(y)*z | x*z | x*z*z |"
"0.0 | 1.0 | 0.0";
const char* g_z = "sin(x)*sin(y)*cos(z) | x*x*y*y | exp(x)*exp(2*y)*2*z |"
"exp(-2*x)*exp(y) | y*z | 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);
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);
utl::matrix3d<Real> r(3,3,3);
r.fill(std::array{cos(x)*sin(y)*sin(z), 2*x*y*y*z, exp(x)*exp(2*y)*z*z,
-2*exp(-2*x)*exp(y)*z, y*z, y*z*z,
1.0, 0.0, 0.0,
sin(x)*cos(y)*sin(z), 2*x*x*y*z, 2*exp(x)*exp(2*y)*z*z,
exp(-2*x)*exp(y)*z, x*z, x*z*z,
0.0, 1.0, 0.0,
sin(x)*sin(y)*cos(z), x*x*y*y, exp(x)*exp(2*y)*2*z,
exp(-2*x)*exp(y), y*z, x*y*2*z,
0.0, 0.0, 1.0}.data());
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));
}
}
}
TEST(TestTensorFunction, Gradient3DFD)
{
const char* g = "sin(x)*sin(y)*sin(z) | x*x*y*y*z | exp(x)*exp(2*y)*z*z |"
"exp(-2*x)*exp(y)*z | x*y*z | x*y*z*z |"
"x | y | z";
const double eps = 1e-6;
TensorFuncExpr 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);
utl::matrix3d<Real> r(3,3,3);
r.fill(std::array{(sin(Xp.x) - sin(Xm.x))*sin(y)*sin(z),
(Xp.x*Xp.x - Xm.x*Xm.x)*y*y*z,
(exp(Xp.x) - exp(Xm.x))*exp(2*y)*z*z,
(exp(-2*Xp.x) - exp(-2*Xm.x))*exp(y)*z,
(Xp.x - Xm.x)*y*z,
(Xp.x - Xm.x)*y*z*z,
Xp.x -Xm.x,
0.0,
0.0,
sin(x)*(sin(Xp.y) - sin(Xm.y))*sin(z),
x*x*(Xp.y*Xp.y - Xm.y*Xm.y)*z,
exp(x)*(exp(2*Xp.y) - exp(2*Xm.y))*z*z,
exp(-2*x)*(exp(Xp.y) - exp(Xm.y))*z,
x*(Xp.y - Xm.y)*z,
x*(Xp.y - Xm.y)*z*z,
0.0,
Xp.y - Xm.y,
0.0,
sin(x)*sin(y)*(sin(Xp.z) - sin(Xm.z)),
x*x*y*y*(Xp.z - Xm.z),
exp(x)*exp(2*y)*(Xp.z*Xp.z - Xm.z*Xm.z),
exp(-2*x)*exp(y)*(Xp.z - Xm.z),
x*y*(Xp.z - Xm.z),
x*y*(Xp.z*Xp.z - Xm.z*Xm.z),
0.0,
0.0,
Xp.z - Xm.z}.data());
r.multiply(1.0 / eps);
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_NEAR(dx(i,j), r(i,j,d), 1e-8);
}
}
}
TEST(TestTensorFunction, Hessian2D)
{
const char* g = "sin(x)*sin(y) | x*x*y*y | exp(x)*exp(2*y) | exp(-2*x)*exp(y)";
const char* g_xx = "-sin(x)*sin(y) | 2*y*y | exp(x)*exp(2*y) | 4*exp(-2*x)*exp(y)";
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);
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);
utl::matrix4d<Real> r(2,2,2,2);
r.fill(std::array{-sin(x)*sin(y), 2.0*y*y,
exp(x)*exp(2*y), 4.0*exp(-2*x)*exp(y),
cos(x)*cos(y), 2*x*2*y,
exp(x)*2*exp(2*y), -2*exp(-2*x)*exp(y),
cos(x)*cos(y), 2*x*2*y,
exp(x)*2*exp(2*y), -2*exp(-2*x)*exp(y),
-sin(x)*sin(y), 2*x*x,
exp(x)*4*exp(2*y), exp(-2*x)*exp(y)}.data());
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));
}
}
}
TEST(TestTensorFunction, Hessian3D)
{
const char* g = "sin(x)*sin(y)*sin(z) | x*x*y*y*z | exp(x)*exp(2*y)*z*z |"
"exp(-2*x)*exp(y)*z | x*y*z | x*y*z*z |"
"x*x | y*y | z*z";
const char* g_xx = "-sin(x)*sin(y)*sin(z) | 2*y*y*z | exp(x)*exp(2*y)*z*z |"
"4*exp(-2*x)*exp(y)*z | 0.0 | 0.0 |"
"2.0 | 0.0 | 0.0";
const char* g_xy = "cos(x)*cos(y)*sin(z) | 2*x*2*y*z | exp(x)*2*exp(2*y)*z*z |"
"-2*exp(-2*x)*exp(y)*z | z | z*z |"
"0.0 | 0.0 | 0.0";
const char* g_xz = "cos(x)*sin(y)*cos(z) | 2*x*y*y | exp(x)*exp(2*y)*2*z |"
"-2*exp(-2*x)*exp(y) | y | y*2*z |"
"0.0 | 0.0 | 0.0";
const char* g_yy = "-sin(x)*sin(y)*sin(z) | x*x*2*z | exp(x)*4*exp(2*y)*z*z |"
"exp(-2*x)*exp(y)*z | 0.0 | 0.0 |"
"0.0 | 2.0 | 0.0";
const char* g_yz = "sin(x)*cos(y)*cos(z) | x*x*2*y | exp(x)*2*exp(2*y)*2*z |"
"exp(-2*x)*exp(y) | x | x*2*z |"
"0.0 | 0.0 | 0.0";
const char* g_zz = "-sin(x)*sin(y)*sin(z) | 0.0 | exp(x)*exp(2*y)*2 |"
"0.0 | 0.0 | 2.0 |"
"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);
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);
utl::matrix4d<Real> r(3,3,3,3);
r.fill(std::array{-sin(x)*sin(y)*sin(z), 2*y*y*z, exp(x)*exp(2*y)*z*z,
4*exp(-2*x)*exp(y)*z, 0.0, 0.0,
2.0, 0.0, 0.0,
cos(x)*cos(y)*sin(z), 2*x*2*y*z, exp(x)*2*exp(2*y)*z*z,
-2*exp(-2*x)*exp(y)*z, z, z*z,
0.0, 0.0, 0.0,
cos(x)*sin(y)*cos(z), 2*x*y*y, exp(x)*exp(2*y)*2*z,
-2*exp(-2*x)*exp(y), y, y*2*z,
0.0, 0.0, 0.0,
cos(x)*cos(y)*sin(z), 2*x*2*y*z, exp(x)*2*exp(2*y)*z*z,
-2*exp(-2*x)*exp(y)*z, z, z*z,
0.0, 0.0, 0.0,
-sin(x)*sin(y)*sin(z), x*x*2*z, exp(x)*4*exp(2*y)*z*z,
exp(-2*x)*exp(y)*z, 0.0, 0.0,
0.0, 2.0, 0.0,
sin(x)*cos(y)*cos(z), x*x*2*y, exp(x)*2*exp(2*y)*2*z,
exp(-2*x)*exp(y), x, x*2*z,
0.0, 0.0, 0.0,
cos(x)*sin(y)*cos(z), 2*x*y*y, exp(x)*exp(2*y)*2*z,
-2*exp(-2*x)*exp(y), y, y*2*z,
0.0, 0.0, 0.0,
sin(x)*cos(y)*cos(z), x*x*2*y, exp(x)*2*exp(2*y)*2*z,
exp(-2*x)*exp(y), x, x*2*z,
0.0, 0.0, 0.0,
-sin(x)*sin(y)*sin(z), 0.0, exp(x)*exp(2*y)*2,
0.0, 0.0, 2.0,
0.0, 0.0, 2.0}.data());
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));
}
}
}
TEST(TestSTensorFunc, Evaluate2D)
{
const char* func = "sin(x) | cos (y) | exp(z)";
STensorFuncExpr f(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), 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));
}
}
TEST(TestSTensorFunc, Evaluate2Dzz)
{
const char* func = "sin(x) | cos (y) | sin(x)*cos(y) | exp(z)";
STensorFuncExpr f(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), 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));
EXPECT_DOUBLE_EQ(fx(3,3), sin(x)*cos(y));
}
}
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);
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));
}
}
TEST(TestSTensorFunction, Gradient2D)
{
const char* g = "sin(x)*sin(y) | exp(x)*exp(2*y) | x*x*y*y";
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);
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);
utl::matrix3d<Real> r(2,2,2);
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());
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));
}
}
}
TEST(TestSTensorFunction, Gradient2Dzz)
{
const char* g = "sin(x)*sin(y) | exp(x)*exp(2*y) | sin(x)*sin(y) | x*x*y*y";
const char* g_x = "cos(x)*sin(y) | exp(x)*exp(2*y) | cos(x)*sin(y) | 2*x*y*y";
const char* g_y = "sin(x)*cos(y) | 2*exp(x)*exp(2*y) | sin(x)*cos(y) | 2*x*x*y";
STensorFuncExpr 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);
utl::matrix3d<Real> r(2,2,2);
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());
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_NEAR(dx(i,j), r(i,j,d), 1e-12);
EXPECT_NEAR(dx(3,3), (d == 1 ? cos(x) : sin(x)) * (d == 2 ? cos(y) : sin(y)), 1e-12);
}
}
}
TEST(TestSTensorFunction, Gradient2DFD)
{
const char* g = "sin(x)*sin(y) | exp(x)*exp(2*y) | x*x*y*y";
const double eps = 1e-6;
STensorFuncExpr 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);
utl::matrix3d<Real> r(2,2,2);
r.fill(std::array{(sin(Xp.x) - sin(Xm.x))*sin(y),
(Xp.x*Xp.x - Xm.x*Xm.x)*y*y,
(Xp.x*Xp.x - Xm.x*Xm.x)*y*y,
(exp(Xp.x) - exp(Xm.x))*exp(2*y),
sin(x)*(sin(Xp.y) - sin(Xm.y)),
x*x*(Xp.y*Xp.y - Xm.y*Xm.y),
x*x*(Xp.y*Xp.y - Xm.y*Xm.y),
exp(x)*(exp(2*Xp.y) - exp(2*Xm.y))}.data());
r.multiply(1.0 / eps);
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_NEAR(dx(i,j), r(i,j,d), 1e-8);
}
}
}
TEST(TestSTensorFunction, Gradient3D)
{
const char* g = "sin(x)*sin(y)*sin(z) | exp(x)*exp(2*y)*exp(z) | x*x*y*y*z*z |"
"x*y*z | x*x*y*z | z";
const char* g_x = "cos(x)*sin(y)*sin(z) | exp(x)*exp(2*y)*exp(z) | 2*x*y*y*z*z |"
"y*z | 2*x*y*z | 0";
const char* g_y = "sin(x)*cos(y)*sin(z) | exp(x)*2*exp(2*y)*exp(z) | x*x*2*y*z*z |"
"x*z | x*x*z | 0";
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);
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);
utl::matrix3d<Real> r(3,3,3);
r.fill(std::array{cos(x)*sin(y)*sin(z), y*z, 0.0,
y*z, exp(x)*exp(2*y)*exp(z), 2*x*y*z,
0.0, 2*x*y*z, 2*x*y*y*z*z,
sin(x)*cos(y)*sin(z), x*z, 0.0,
x*z, exp(x)*2*exp(2*y)*exp(z), x*x*z,
0.0, x*x*z, x*x*2*y*z*z,
sin(x)*sin(y)*cos(z), x*y, 1.0,
x*y, exp(x)*exp(2*y)*exp(z), x*x*y,
1.0, x*x*y, x*x*y*y*2*z}.data());
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));
}
}
}
TEST(TestSTensorFunction, Gradient3DFD)
{
const char* g = "sin(x)*sin(y)*sin(z) | exp(x)*exp(2*y)*exp(z) | x*x*y*y*z*z |"
"x*y*z | x*x*y*z | z";
const double eps = 1e-6;
STensorFuncExpr 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);
utl::matrix3d<Real> r(3,3,3);
r.fill(std::array{(sin(Xp.x) - sin(Xm.x))*sin(y)*sin(z),
(Xp.x - Xm.x)*y*z,
0.0,
(Xp.x - Xm.x)*y*z,
(exp(Xp.x) - exp(Xm.x))*exp(2*y)*exp(z),
(Xp.x*Xp.x - Xm.x*Xm.x)*y*z,
0.0,
(Xp.x*Xp.x - Xm.x*Xm.x)*y*z,
(Xp.x*Xp.x - Xm.x*Xm.x)*y*y*z*z,
sin(x)*(sin(Xp.y) - sin(Xm.y))*sin(z),
x*(Xp.y - Xm.y)*z,
0.0,
x*(Xp.y - Xm.y)*z,
exp(x)*(exp(2*Xp.y) - exp(2*Xm.y))*exp(z),
x*x*(Xp.y - Xm.y)*z,
0.0,
x*x*(Xp.y - Xm.y)*z,
x*x*(Xp.y*Xp.y - Xm.y*Xm.y)*z*z,
sin(x)*sin(y)*(sin(Xp.z) - sin(Xm.z)),
x*y*(Xp.z - Xm.z),
Xp.z - Xm.z,
x*y*(Xp.z - Xm.z),
exp(x)*exp(2*y)*(exp(Xp.z) - exp(Xm.z)),
x*x*y*(Xp.z - Xm.z),
Xp.z - Xm.z,
x*x*y*(Xp.z - Xm.z),
x*x*y*y*(Xp.z*Xp.z - Xm.z*Xm.z)}.data());
r.multiply(1.0 / eps);
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_NEAR(dx(i,j), r(i,j,d), 1e-8);
}
}
}
TEST(TestEvalFunction, ExtraParam)
{
EvalFunction f("x*foo");