eba7a485ba
COMP: define labelSphericalTensor::I - remove spurious 'labelI' global constant (labelSphericalTensor::I) STYLE: replace use of deprecated Tensor vectorComponent STYLE: avoid bit-wise assignment of bool (VectorSpace compare ops)
1177 lines
32 KiB
C
1177 lines
32 KiB
C
/*---------------------------------------------------------------------------*\
|
|
========= |
|
|
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
|
\\ / O peration |
|
|
\\ / A nd | www.openfoam.com
|
|
\\/ M anipulation |
|
|
-------------------------------------------------------------------------------
|
|
Copyright (C) 2018 OpenFOAM Foundation
|
|
Copyright (C) 2019-2020 OpenCFD Ltd.
|
|
-------------------------------------------------------------------------------
|
|
License
|
|
This file is part of OpenFOAM.
|
|
|
|
OpenFOAM is free software: you can redistribute it and/or modify it
|
|
under the terms of the GNU General Public License as published by
|
|
the Free Software Foundation, either version 3 of the License, or
|
|
(at your option) any later version.
|
|
|
|
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
|
|
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
for more details.
|
|
|
|
You should have received a copy of the GNU General Public License
|
|
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
Application
|
|
Test-Tensor
|
|
|
|
Description
|
|
Tests for \c Tensor constructors, member functions and operators
|
|
using \c floatScalar, \c doubleScalar, and \c complex base types.
|
|
|
|
Eigen decomposition tests for \c tensor, i.e. Tensor<scalar>.
|
|
|
|
Cross-checks were obtained from 'NumPy 1.15.1' and 'SciPy 1.1.0' if no
|
|
theoretical cross-check exists (like eigendecomposition relations), and
|
|
were hard-coded for elementwise comparisons.
|
|
|
|
For \c complex base type, the cross-checks do only involve zero imag part.
|
|
|
|
\*---------------------------------------------------------------------------*/
|
|
|
|
#include "tensor.H"
|
|
#include "transform.H"
|
|
#include "Random.H"
|
|
#include "scalar.H"
|
|
#include "complex.H"
|
|
|
|
using namespace Foam;
|
|
|
|
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
|
|
|
// Total number of unit tests
|
|
unsigned nTest_ = 0;
|
|
|
|
|
|
// Total number of failed unit tests
|
|
unsigned nFail_ = 0;
|
|
|
|
|
|
// Create a random tensor
|
|
tensor makeRandomContainer(Random& rnd)
|
|
{
|
|
tensor A(Zero);
|
|
std::generate(A.begin(), A.end(), [&]{ return rnd.GaussNormal<scalar>(); });
|
|
return A;
|
|
}
|
|
|
|
|
|
// Compare two floating point types, and print output.
|
|
// Do ++nFail_ if values of two objects are not equal within a given tolerance.
|
|
// The function is converted from PEP-485.
|
|
template<class Type>
|
|
typename std::enable_if<pTraits<Type>::rank == 0, void>::type
|
|
cmp
|
|
(
|
|
const word& msg,
|
|
const Type& x,
|
|
const Type& y,
|
|
const scalar absTol = 0, //<! useful for cmps near zero
|
|
const scalar relTol = 1e-8 //<! are values the same within 8 decimals
|
|
)
|
|
{
|
|
Info<< msg << x << "?=" << y << endl;
|
|
|
|
unsigned nFail = 0;
|
|
|
|
if (max(absTol, relTol*max(mag(x), mag(y))) < mag(x - y))
|
|
{
|
|
++nFail;
|
|
}
|
|
|
|
if (nFail)
|
|
{
|
|
Info<< nl
|
|
<< " #### Fail in " << nFail << " comps ####" << nl << endl;
|
|
++nFail_;
|
|
}
|
|
++nTest_;
|
|
}
|
|
|
|
|
|
// Compare two containers elementwise, and print output.
|
|
// Do ++nFail_ if two components are not equal within a given tolerance.
|
|
// The function is converted from PEP-485
|
|
template<class Type>
|
|
typename std::enable_if<pTraits<Type>::rank != 0, void>::type
|
|
cmp
|
|
(
|
|
const word& msg,
|
|
const Type& x,
|
|
const Type& y,
|
|
const scalar absTol = 0,
|
|
const scalar relTol = 1e-8
|
|
)
|
|
{
|
|
Info<< msg << x << "?=" << y << endl;
|
|
|
|
unsigned nFail = 0;
|
|
|
|
for (direction i = 0; i < pTraits<Type>::nComponents; ++i)
|
|
{
|
|
if (max(absTol, relTol*max(mag(x[i]), mag(y[i]))) < mag(x[i] - y[i]))
|
|
{
|
|
++nFail;
|
|
}
|
|
}
|
|
|
|
if (nFail)
|
|
{
|
|
Info<< nl
|
|
<< " #### Fail in " << nFail << " comps ####" << nl << endl;
|
|
++nFail_;
|
|
}
|
|
++nTest_;
|
|
}
|
|
|
|
|
|
// Create each constructor of Tensor<Type>, and print output
|
|
template<class Type>
|
|
void test_constructors(Type)
|
|
{
|
|
{
|
|
Info<< "# Construct initialized to zero:" << nl;
|
|
const Tensor<Type> T(Zero);
|
|
Info<< T << endl;
|
|
}
|
|
{
|
|
Info<< "# Construct given MatrixSpace of the same rank:" << nl;
|
|
const MatrixSpace<Tensor<Type>, Type, 3, 3> M(Zero);
|
|
const Tensor<Type> T(M);
|
|
Info<< T << endl;
|
|
}
|
|
{
|
|
Info<< "# Construct given VectorSpace of the same rank:" << nl;
|
|
const VectorSpace<Tensor<Type>, Type, 9> V(Zero);
|
|
const Tensor<Type> T(V);
|
|
Info<< T << endl;
|
|
}
|
|
{
|
|
Info<< "# Construct given SphericalTensor<Type>:" << nl;
|
|
const SphericalTensor<Type> Sp(Type(5));
|
|
const Tensor<Type> T(Sp);
|
|
Info<< T << endl;
|
|
}
|
|
{
|
|
Info<< "# Construct given SymmTensor<Type>:" << nl;
|
|
const SymmTensor<Type> S
|
|
(
|
|
Type(1), Type(2), Type(3),
|
|
Type(5), Type(6),
|
|
Type(9)
|
|
);
|
|
const Tensor<Type> T(S);
|
|
Info<< T << endl;
|
|
}
|
|
{
|
|
Info<< "# Construct given triad of row vectors,"
|
|
<< " optionally treated as transposed (ie, column vectors)" << nl;
|
|
const Vector<Vector<Type>> vecs
|
|
(
|
|
Vector<Type>(Type(1), Type(2), Type(3)),
|
|
Vector<Type>(Type(4), Type(5), Type(6)),
|
|
Vector<Type>(Type(7), Type(8), Type(9))
|
|
);
|
|
const Tensor<Type> T(vecs);
|
|
Info<< T << nl;
|
|
|
|
const Tensor<Type> transposedT(vecs, true);
|
|
Info<< transposedT << endl;
|
|
}
|
|
{
|
|
Info<< "# Construct given the three row vectors,"
|
|
<< " optionally treated as transposed (ie, column vectors)" << nl;
|
|
const Vector<Type> a(Type(1), Type(2), Type(3));
|
|
const Vector<Type> b(Type(4), Type(5), Type(6));
|
|
const Vector<Type> c(Type(7), Type(8), Type(9));
|
|
const Tensor<Type> T(a, b, c);
|
|
Info<< T << nl;
|
|
|
|
const Tensor<Type> transposedT(a, b, c, true);
|
|
Info<< transposedT << endl;
|
|
}
|
|
{
|
|
Info<< "# Construct given the nine components:" << nl;
|
|
const Tensor<Type> T
|
|
(
|
|
Type(1), Type(2), Type(-3),
|
|
Type(4), Type(5), Type(-6),
|
|
Type(7), Type(8), Type(-9)
|
|
);
|
|
Info<< T << endl;
|
|
}
|
|
{
|
|
Info<< "# Copy construct:" << nl;
|
|
const Tensor<Type> T(Zero);
|
|
const Tensor<Type> copyT(T);
|
|
Info<< T << tab << copyT << endl;
|
|
}
|
|
}
|
|
|
|
|
|
// Execute each member function of Tensor<Type>, and print output
|
|
template<class Type>
|
|
void test_member_funcs(Type)
|
|
{
|
|
Tensor<Type> T
|
|
(
|
|
Type(1), Type(2), Type(-3),
|
|
Type(4), Type(5), Type(-6),
|
|
Type(7), Type(8), Type(-9)
|
|
);
|
|
Tensor<Type> Tbak = T;
|
|
const Tensor<Type> cT
|
|
(
|
|
Type(-9), Type(8), Type(7),
|
|
Type(-6), Type(5), Type(4),
|
|
Type(-3), Type(2), Type(1)
|
|
);
|
|
|
|
Info<< "# Operand: " << nl
|
|
<< " Tensor = " << T << endl;
|
|
|
|
|
|
{
|
|
Info<< "# Component access:" << nl;
|
|
|
|
Tensor<Type> cpT
|
|
(
|
|
T.xx(), T.xy(), T.xz(),
|
|
T.yx(), T.yy(), T.yz(),
|
|
T.zx(), T.zy(), T.zz()
|
|
);
|
|
cmp(" 'Tensor' access:", T, cpT);
|
|
|
|
const Tensor<Type> cpcT
|
|
(
|
|
cT.xx(), cT.xy(), cT.xz(),
|
|
cT.yx(), cT.yy(), cT.yz(),
|
|
cT.zx(), cT.zy(), cT.zz()
|
|
);
|
|
cmp(" 'const Tensor' access:", cT, cpcT);
|
|
}
|
|
{
|
|
Info<< "# Column-vector access:" << nl;
|
|
cmp(" cx():", T.cx(), Vector<Type>(Type(1), Type(4), Type(7)));
|
|
cmp(" cy():", T.cy(), Vector<Type>(Type(2), Type(5), Type(8)));
|
|
cmp(" cz():", T.cz(), Vector<Type>(Type(-3), Type(-6), Type(-9)));
|
|
cmp(" col(0):", T.col(0), Vector<Type>(Type(1), Type(4), Type(7)));
|
|
cmp(" col(1):", T.col(1), Vector<Type>(Type(2), Type(5), Type(8)));
|
|
cmp(" col(2):", T.col(2), Vector<Type>(Type(-3), Type(-6), Type(-9)));
|
|
cmp
|
|
(
|
|
" col<0>:",
|
|
T.template col<0>(),
|
|
Vector<Type>(Type(1), Type(4), Type(7))
|
|
);
|
|
cmp
|
|
(
|
|
" col<1>:",
|
|
T.template col<1>(),
|
|
Vector<Type>(Type(2), Type(5), Type(8))
|
|
);
|
|
cmp
|
|
(
|
|
" col<2>:",
|
|
T.template col<2>(),
|
|
Vector<Type>(Type(-3), Type(-6), Type(-9))
|
|
);
|
|
// Compilation error: Info << " col<3> = " << T.col<3>() << nl;
|
|
|
|
|
|
Info<< "# Column-vector manipulation:" << nl;
|
|
T.col(1, Vector<Type>(Type(0), Type(1), Type(99)));
|
|
cmp
|
|
(
|
|
" col(1, Vector):",
|
|
T.col(1),
|
|
Vector<Type>(Type(0), Type(1), Type(99))
|
|
);
|
|
|
|
T.cols
|
|
(
|
|
Vector<Type>(Type(1), Type(1), Type(1)),
|
|
Vector<Type>(Type(-1), Type(1), Type(2)),
|
|
Vector<Type>(Type(1), Type(1), Type(3))
|
|
);
|
|
cmp
|
|
(
|
|
" cols(Vectors):",
|
|
T,
|
|
Tensor<Type>
|
|
(
|
|
Type(1), Type(-1), Type(1),
|
|
Type(1), Type(1), Type(1),
|
|
Type(1), Type(2), Type(3)
|
|
)
|
|
);
|
|
}
|
|
{
|
|
Info<< "# Row-vector access:" << nl;
|
|
T = Tbak;
|
|
cmp(" x():", T.x(), Vector<Type>(Type(1), Type(2), Type(-3)));
|
|
cmp(" y():", T.y(), Vector<Type>(Type(4), Type(5), Type(-6)));
|
|
cmp(" z():", T.z(), Vector<Type>(Type(7), Type(8), Type(-9)));
|
|
cmp(" row(0):", T.row(0), Vector<Type>(Type(1), Type(2), Type(-3)));
|
|
cmp(" row(1):", T.row(1), Vector<Type>(Type(4), Type(5), Type(-6)));
|
|
cmp(" row(2):", T.row(2), Vector<Type>(Type(7), Type(8), Type(-9)));
|
|
cmp
|
|
(
|
|
" row<0>:",
|
|
T.template row<0>(),
|
|
Vector<Type>(Type(1), Type(2), Type(-3))
|
|
);
|
|
cmp
|
|
(
|
|
" row<1>:",
|
|
T.template row<1>(),
|
|
Vector<Type>(Type(4), Type(5), Type(-6))
|
|
);
|
|
cmp
|
|
(
|
|
" row<2>:",
|
|
T.template row<2>(),
|
|
Vector<Type>(Type(7), Type(8), Type(-9))
|
|
);
|
|
// Compilation error: Info << " row<3> = " << T.row<3>() << nl;
|
|
|
|
|
|
Info<< "# Row-vector manipulation:" << nl;
|
|
T.row(1, Vector<Type>(Type(0), Type(1), Type(99)));
|
|
cmp
|
|
(
|
|
" row(1, Vector):",
|
|
T.row(1),
|
|
Vector<Type>(Type(0), Type(1), Type(99))
|
|
);
|
|
|
|
T.rows
|
|
(
|
|
Vector<Type>(Type(1), Type(1), Type(1)),
|
|
Vector<Type>(Type(-1), Type(1), Type(2)),
|
|
Vector<Type>(Type(1), Type(1), Type(3))
|
|
);
|
|
cmp
|
|
(
|
|
" rows(Vectors):",
|
|
T,
|
|
Tensor<Type>
|
|
(
|
|
Type(1), Type(1), Type(1),
|
|
Type(-1), Type(1), Type(2),
|
|
Type(1), Type(1), Type(3)
|
|
)
|
|
);
|
|
}
|
|
{
|
|
Info<< "# Diagonal access:" << nl;
|
|
|
|
T = Tbak;
|
|
cmp
|
|
(
|
|
" 'Tensor'.diag():",
|
|
T.diag(),
|
|
Vector<Type>(Type(1), Type(5), Type(-9))
|
|
);
|
|
cmp
|
|
(
|
|
" 'const Tensor'.diag():",
|
|
cT.diag(),
|
|
Vector<Type>(Type(-9), Type(5), Type(1))
|
|
);
|
|
|
|
|
|
Info<< "# Diagonal manipulation:" << nl;
|
|
|
|
T.diag(Vector<Type>(Type(-10), Type(-15), Type(-20)));
|
|
cmp
|
|
(
|
|
" 'Tensor'.diag('Vector'):",
|
|
T.diag(),
|
|
Vector<Type>(Type(-10), Type(-15), Type(-20))
|
|
);
|
|
}
|
|
{
|
|
Info<< "# Tensor operations:" << nl;
|
|
|
|
T = Tbak;
|
|
cmp(" Transpose:", T, (T.T()).T());
|
|
cmp // Singular matrix
|
|
(
|
|
" Inverse:",
|
|
T.inv(),
|
|
Tensor<Type>
|
|
(
|
|
Type(-4.50359963e+15), Type(9.00719925e+15), Type(-4.50359963e+15),
|
|
Type(9.00719925e+15), Type(-1.80143985e+16), Type(9.00719925e+15),
|
|
Type(4.50359963e+15), Type(-9.00719925e+15), Type(4.50359963e+15)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Inner-product:",
|
|
T.inner(T),
|
|
Tensor<Type>
|
|
(
|
|
Type(-12), Type(-12), Type(12),
|
|
Type(-18), Type(-15), Type(12),
|
|
Type(-24), Type(-18), Type(12)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Schur-product:",
|
|
T.schur(T),
|
|
Tensor<Type>
|
|
(
|
|
Type(1), Type(4), Type(9),
|
|
Type(16), Type(25), Type(36),
|
|
Type(49), Type(64), Type(81)
|
|
)
|
|
);
|
|
}
|
|
{
|
|
Info<< "# Member operators:" << nl;
|
|
|
|
T = Tbak;
|
|
T &= T;
|
|
cmp
|
|
(
|
|
" Assign inner-product of this with another Tensor:",
|
|
T,
|
|
(Tbak & Tbak)
|
|
);
|
|
|
|
T = VectorSpace<Tensor<Type>, Type, 9>(Zero);
|
|
cmp
|
|
(
|
|
" Assign to an equivalent vector space:",
|
|
T,
|
|
Tensor<Type>(Zero)
|
|
);
|
|
|
|
T = SphericalTensor<Type>(Type(5));
|
|
cmp
|
|
(
|
|
" Assign to a SphericalTensor:",
|
|
T,
|
|
Tensor<Type>
|
|
(
|
|
Type(5), Zero, Zero,
|
|
Zero, Type(5), Zero,
|
|
Zero, Zero, Type(5)
|
|
)
|
|
);
|
|
|
|
T = SymmTensor<Type>
|
|
(
|
|
Type(1), Type(2), Type(-3),
|
|
Type(5), Type(-6),
|
|
Type(-9)
|
|
);
|
|
cmp
|
|
(
|
|
" Assign to a SymmTensor:",
|
|
T,
|
|
Tensor<Type>
|
|
(
|
|
Type(1), Type(2), Type(-3),
|
|
Type(2), Type(5), Type(-6),
|
|
Type(-3), Type(-6), Type(-9)
|
|
)
|
|
);
|
|
|
|
T = Vector<Vector<Type>>
|
|
(
|
|
Vector<Type>(Type(-1), Type(2), Type(3)),
|
|
Vector<Type>(Type(-4), Type(5), Type(6)),
|
|
Vector<Type>(Type(4), Type(-5), Type(6))
|
|
);
|
|
cmp
|
|
(
|
|
" Assign to a triad of row vectors:",
|
|
T,
|
|
Tensor<Type>
|
|
(
|
|
Type(-1), Type(2), Type(3),
|
|
Type(-4), Type(5), Type(6),
|
|
Type(4), Type(-5), Type(6)
|
|
)
|
|
);
|
|
}
|
|
}
|
|
|
|
|
|
// Execute each global function of Tensor<Type>, and print output
|
|
template<class Type>
|
|
void test_global_funcs(Type)
|
|
{
|
|
const Tensor<Type> T
|
|
(
|
|
Type(-1), Type(2), Type(-3),
|
|
Type(4), Type(5), Type(-6),
|
|
Type(7), Type(8), Type(-9)
|
|
);
|
|
const SymmTensor<Type> sT
|
|
(
|
|
Type(-1), Type(2), Type(-3),
|
|
Type(5), Type(-6),
|
|
Type(-9)
|
|
);
|
|
|
|
Info<< "# Operands: " << nl
|
|
<< " Tensor = " << T << nl
|
|
<< " SymmTensor = " << sT << endl;
|
|
|
|
|
|
cmp(" Trace = ", tr(T), Type(-5));
|
|
cmp(" Spherical part = ", sph(T), SphericalTensor<Type>(tr(T)/Type(3)));
|
|
cmp
|
|
(
|
|
" Symmetric part = ",
|
|
symm(T),
|
|
SymmTensor<Type>
|
|
(
|
|
Type(-1), Type(3), Type(2),
|
|
Type(5), Type(1),
|
|
Type(-9)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Twice the symmetric part = ",
|
|
twoSymm(T),
|
|
SymmTensor<Type>
|
|
(
|
|
Type(-2), Type(6), Type(4),
|
|
Type(10), Type(2),
|
|
Type(-18)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Skew-symmetric part = ",
|
|
skew(T),
|
|
Tensor<Type>
|
|
(
|
|
Type(0), Type(-1), Type(-5),
|
|
Type(1), Type(0), Type(-7),
|
|
Type(5), Type(7), Type(0)
|
|
)
|
|
);
|
|
/*
|
|
// Complex-type is not supported for this function.
|
|
cmp
|
|
(
|
|
" Skew-symmetric part of a SymmTensor = ",
|
|
skew(sT),
|
|
Tensor<Type>(Zero)
|
|
);
|
|
*/
|
|
cmp
|
|
(
|
|
" Deviatoric part = ",
|
|
dev(T),
|
|
Tensor<Type>
|
|
(
|
|
Type(0.66666667), Type(2), Type(-3),
|
|
Type(4), Type(6.66666667), Type(-6),
|
|
Type(7), Type(8), Type(-7.33333333)
|
|
),
|
|
1e-7
|
|
);
|
|
cmp(" Two-third deviatoric part = ", dev2(T), T - 2*sph(T));
|
|
cmp(" Determinant = ", det(T), Type(-6.000000000000005));
|
|
cmp
|
|
(
|
|
" Cofactor tensor = ",
|
|
cof(T),
|
|
Tensor<Type>
|
|
(
|
|
Type(3), Type(-6), Type(-3),
|
|
Type(-6), Type(30), Type(22),
|
|
Type(3), Type(-18), Type(-13)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Inverse = ",
|
|
inv(T, det(T)),
|
|
Tensor<Type>
|
|
(
|
|
Type(-0.5), Type(1), Type(-0.5),
|
|
Type(1), Type(-5), Type(3),
|
|
Type(0.5), Type(-3.66666667), Type(2.16666667)
|
|
),
|
|
1e-8
|
|
);
|
|
cmp
|
|
(
|
|
" Inverse (another) = ",
|
|
inv(T),
|
|
Tensor<Type>
|
|
(
|
|
Type(-0.5), Type(1), Type(-0.5),
|
|
Type(1), Type(-5), Type(3),
|
|
Type(0.5), Type(-3.66666667), Type(2.16666667)
|
|
),
|
|
1e-8
|
|
);
|
|
cmp
|
|
(
|
|
" Inverse (another) = ",
|
|
T.inv(),
|
|
Tensor<Type>
|
|
(
|
|
Type(-0.5), Type(1), Type(-0.5),
|
|
Type(1), Type(-5), Type(3),
|
|
Type(0.5), Type(-3.66666667), Type(2.16666667)
|
|
),
|
|
1e-8
|
|
);
|
|
cmp(" First invariant = ", invariantI(T), Type(-5));
|
|
cmp(" Second invariant = ", invariantII(T), Type(20));
|
|
cmp(" Third invariant = ", invariantIII(T), Type(-6.000000000000005));
|
|
}
|
|
|
|
|
|
// Execute each global operator of Tensor<Type>, and print output
|
|
template<class Type>
|
|
void test_global_opers(Type)
|
|
{
|
|
const Tensor<Type> T
|
|
(
|
|
Type(-1), Type(2), Type(-3),
|
|
Type(4), Type(5), Type(-6),
|
|
Type(7), Type(8), Type(-9)
|
|
);
|
|
const SymmTensor<Type> sT
|
|
(
|
|
Type(-1), Type(2), Type(-3),
|
|
Type(5), Type(-6),
|
|
Type(-9)
|
|
);
|
|
const SphericalTensor<Type> spT(Type(1));
|
|
const Vector<Type> v(Type(3), Type(2), Type(1));
|
|
const Type x(4);
|
|
|
|
Info<< "# Operands:" << nl
|
|
<< " Tensor = " << T << nl
|
|
<< " SymmTensor = " << sT << nl
|
|
<< " SphericalTensor = " << spT << nl
|
|
<< " Vector = " << v << nl
|
|
<< " Type = " << x << endl;
|
|
|
|
|
|
cmp
|
|
(
|
|
" Sum of SpTensor-Tensor = ",
|
|
(spT + T),
|
|
Tensor<Type>
|
|
(
|
|
Type(0), Type(2), Type(-3),
|
|
Type(4), Type(6), Type(-6),
|
|
Type(7), Type(8), Type(-8)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Sum of Tensor-SpTensor = ",
|
|
(T + spT),
|
|
Tensor<Type>
|
|
(
|
|
Type(0), Type(2), Type(-3),
|
|
Type(4), Type(6), Type(-6),
|
|
Type(7), Type(8), Type(-8)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Sum of SymmTensor-Tensor = ",
|
|
(sT + T),
|
|
Tensor<Type>
|
|
(
|
|
Type(-2), Type(4), Type(-6),
|
|
Type(6), Type(10), Type(-12),
|
|
Type(4), Type(2), Type(-18)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Sum of Tensor-SymmTensor = ",
|
|
(T + sT),
|
|
Tensor<Type>
|
|
(
|
|
Type(-2), Type(4), Type(-6),
|
|
Type(6), Type(10), Type(-12),
|
|
Type(4), Type(2), Type(-18)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Subtract Tensor from SpTensor = ",
|
|
(spT - T),
|
|
Tensor<Type>
|
|
(
|
|
Type(2), Type(-2), Type(3),
|
|
Type(-4), Type(-4), Type(6),
|
|
Type(-7), Type(-8), Type(10)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Subtract SpTensor from Tensor = ",
|
|
(T - spT),
|
|
Tensor<Type>
|
|
(
|
|
Type(-2), Type(2), Type(-3),
|
|
Type(4), Type(4), Type(-6),
|
|
Type(7), Type(8), Type(-10)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Subtract Tensor from SymmTensor = ",
|
|
(sT - T),
|
|
Tensor<Type>
|
|
(
|
|
Type(0), Type(0), Type(0),
|
|
Type(-2), Type(0), Type(0),
|
|
Type(-10), Type(-14), Type(0)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Subtract SymmTensor from Tensor = ",
|
|
(T - sT),
|
|
Tensor<Type>
|
|
(
|
|
Type(0), Type(0), Type(0),
|
|
Type(2), Type(0), Type(0),
|
|
Type(10), Type(14), Type(0)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Hodge dual of a Tensor = ",
|
|
*T,
|
|
Vector<Type>(T.yz(), -T.xz(), T.xy())
|
|
);
|
|
cmp
|
|
(
|
|
" Hodge dual of a Vector = ",
|
|
*v,
|
|
Tensor<Type>
|
|
(
|
|
Zero, -v.z(), v.y(),
|
|
v.z(), Zero, -v.x(),
|
|
-v.y(), v.x(), Zero
|
|
)
|
|
);
|
|
/*cmp
|
|
(
|
|
" Division of Vector by Tensor = ",
|
|
(v/T),
|
|
Tensor<Type>
|
|
(
|
|
Type(-3), Type(1), Type(-0.33333333),
|
|
Type(0.75), Type(0.4), Type(-0.16666667),
|
|
Type(0.42857143), Type(0.25), Type(-0.11111111)
|
|
)
|
|
);*/
|
|
cmp
|
|
(
|
|
" Division of Tensor by Type = ",
|
|
(T/x),
|
|
Tensor<Type>
|
|
(
|
|
Type(-0.25), Type(0.5), Type(-0.75),
|
|
Type(1), Type(1.25), Type(-1.5),
|
|
Type(1.75), Type(2), Type(-2.25)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Inner-product of Tensor-Tensor = ",
|
|
(T & T),
|
|
Tensor<Type>
|
|
(
|
|
Type(-12), Type(-16), Type(18),
|
|
Type(-26), Type(-15), Type(12),
|
|
Type(-38), Type(-18), Type(12)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Inner-product of SpTensor-Tensor = ",
|
|
(spT & T),
|
|
Tensor<Type>
|
|
(
|
|
Type(-1), Type(2), Type(-3),
|
|
Type(4), Type(5), Type(-6),
|
|
Type(7), Type(8), Type(-9)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Inner-product of Tensor-SpTensor = ",
|
|
(T & spT),
|
|
Tensor<Type>
|
|
(
|
|
Type(-1), Type(2), Type(-3),
|
|
Type(4), Type(5), Type(-6),
|
|
Type(7), Type(8), Type(-9)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Inner-product of SymmTensor-Tensor = ",
|
|
(sT & T),
|
|
Tensor<Type>
|
|
(
|
|
Type(-12), Type(-16), Type(18),
|
|
Type(-24), Type(-19), Type(18),
|
|
Type(-84), Type(-108), Type(126)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Inner-product of Tensor-SymmTensor = ",
|
|
(T & sT),
|
|
Tensor<Type>
|
|
(
|
|
Type(14), Type(26), Type(18),
|
|
Type(24), Type(69), Type(12),
|
|
Type(36), Type(108), Type(12)
|
|
)
|
|
);
|
|
cmp
|
|
(
|
|
" Inner-product of Tensor-Vector = ",
|
|
(T & v),
|
|
Vector<Type>(Type(-2), Type(16), Type(28)) // Column-vector
|
|
);
|
|
cmp
|
|
(
|
|
" Inner-product of Vector-Tensor = ",
|
|
(v & T),
|
|
Vector<Type>(Type(12), Type(24), Type(-30)) // Row-vector
|
|
);
|
|
cmp(" D-inner-product of SpTensor-Tensor = ", (spT && T), Type(-5));
|
|
cmp(" D-inner-product of Tensor-SpTensor = ", (T && spT), Type(-5));
|
|
cmp(" D-inner-product of SymmTensor-Tensor = ", (sT && T), Type(95));
|
|
cmp(" D-inner-product of Tensor-SymmTensor = ", (T && sT), Type(95));
|
|
cmp
|
|
(
|
|
" Outer-product of Vector-Vector = ",
|
|
(v*v),
|
|
Tensor<Type>
|
|
(
|
|
Type(9), Type(6), Type(3),
|
|
Type(6), Type(4), Type(2),
|
|
Type(3), Type(2), Type(1)
|
|
)
|
|
);
|
|
}
|
|
|
|
|
|
// Return false if given eigenvalues fail to satisy eigenvalue relations
|
|
// Relations: (Beauregard & Fraleigh (1973), ISBN 0-395-14017-X, p. 307)
|
|
void test_eigenvalues
|
|
(
|
|
const tensor& T,
|
|
const Vector<complex>& EVals,
|
|
const bool prod = true
|
|
)
|
|
{
|
|
if (prod)
|
|
{
|
|
const scalar determinant = det(T);
|
|
// In case of complex EVals, the production is effectively scalar
|
|
// due to the (complex*complex conjugate) results in zero imag part
|
|
const scalar EValsProd = ((EVals.x()*EVals.y()*EVals.z()).real());
|
|
cmp("# Product of eigenvalues = det(T):", EValsProd, determinant, 1e-8);
|
|
}
|
|
|
|
{
|
|
const scalar trace = tr(T);
|
|
scalar EValsSum = 0.0;
|
|
// In case of complex EVals, the summation is effectively scalar
|
|
// due to the (complex+complex conjugate) results in zero imag part
|
|
for (const auto& val : EVals)
|
|
{
|
|
EValsSum += val.real();
|
|
}
|
|
cmp("# Sum of eigenvalues = trace(T):", EValsSum, trace);
|
|
}
|
|
}
|
|
|
|
|
|
// Return false if a given eigenvalue-eigenvector pair
|
|
// fails to satisfy the characteristic equation
|
|
void test_characteristic_equation
|
|
(
|
|
const tensor& T,
|
|
const Vector<complex>& EVals,
|
|
const Tensor<complex>& EVecs
|
|
)
|
|
{
|
|
Info<< "# Characteristic equation:" << nl;
|
|
|
|
Tensor<complex> Tc(Zero);
|
|
forAll(T, i)
|
|
{
|
|
Tc[i] = complex(T[i], 0);
|
|
}
|
|
|
|
for (direction dir = 0; dir < pTraits<vector>::nComponents; ++dir)
|
|
{
|
|
const Vector<complex> leftSide(Tc & EVecs.row(dir));
|
|
const Vector<complex> rightSide(EVals[dir]*EVecs.row(dir));
|
|
const Vector<complex> X(leftSide - rightSide);
|
|
|
|
for (const auto x : X)
|
|
{
|
|
cmp(" (sT & EVec - EVal*EVec) = 0:", mag(x), 0.0, 1e-5);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
// Return false if the eigen functions fail to satisfy relations
|
|
void test_eigen_funcs(const tensor& T, const bool prod = true)
|
|
{
|
|
Info<< "# Operand:" << nl
|
|
<< " tensor = " << T << nl;
|
|
|
|
|
|
Info<< "# Return eigenvalues of a given tensor:" << nl;
|
|
const Vector<complex> EVals(eigenValues(T));
|
|
Info<< EVals << endl;
|
|
test_eigenvalues(T, EVals, prod);
|
|
|
|
Info<< "# Return an eigenvector of a given tensor in a given direction"
|
|
<< " corresponding to a given eigenvalue:" << nl;
|
|
const Vector<complex> standardBasis1(Zero, pTraits<complex>::one, Zero);
|
|
const Vector<complex> standardBasis2(Zero, Zero, pTraits<complex>::one);
|
|
const Vector<complex> EVec
|
|
(
|
|
eigenVector(T, EVals.x(), standardBasis1, standardBasis2)
|
|
);
|
|
Info<< EVec << endl;
|
|
|
|
Info<< "# Return eigenvectors of a given tensor corresponding to"
|
|
<< " given eigenvalues:" << nl;
|
|
const Tensor<complex> EVecs0(eigenVectors(T, EVals));
|
|
Info<< EVecs0 << endl;
|
|
test_characteristic_equation(T, EVals, EVecs0);
|
|
|
|
Info<< "# Return eigenvectors of a given tensor by computing"
|
|
<< " the eigenvalues of the tensor in the background:" << nl;
|
|
const Tensor<complex> EVecs1(eigenVectors(T));
|
|
Info<< EVecs1 << endl;
|
|
}
|
|
|
|
|
|
// Do compile-time recursion over the given types
|
|
template<std::size_t I = 0, typename... Tp>
|
|
inline typename std::enable_if<I == sizeof...(Tp), void>::type
|
|
run_tests(const std::tuple<Tp...>& types, const List<word>& typeID){}
|
|
|
|
|
|
template<std::size_t I = 0, typename... Tp>
|
|
inline typename std::enable_if<I < sizeof...(Tp), void>::type
|
|
run_tests(const std::tuple<Tp...>& types, const List<word>& typeID)
|
|
{
|
|
Info<< nl << " ## Test constructors: "<< typeID[I] <<" ##" << nl;
|
|
test_constructors(std::get<I>(types));
|
|
|
|
Info<< nl << " ## Test member functions: "<< typeID[I] <<" ##" << nl;
|
|
test_member_funcs(std::get<I>(types));
|
|
|
|
Info<< nl << " ## Test global functions: "<< typeID[I] << " ##" << nl;
|
|
test_global_funcs(std::get<I>(types));
|
|
|
|
Info<< nl << " ## Test global operators: "<< typeID[I] <<" ##" << nl;
|
|
test_global_opers(std::get<I>(types));
|
|
|
|
run_tests<I + 1, Tp...>(types, typeID);
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * //
|
|
|
|
int main()
|
|
{
|
|
const std::tuple<floatScalar, doubleScalar, complex> types
|
|
(
|
|
std::make_tuple(Zero, Zero, Zero)
|
|
);
|
|
|
|
const List<word> typeID
|
|
({
|
|
"Tensor<floatScalar>",
|
|
"Tensor<doubleScalar>",
|
|
"Tensor<complex>"
|
|
});
|
|
|
|
run_tests(types, typeID);
|
|
|
|
|
|
Info<< nl << " ## Test tensor eigen functions: ##" << nl;
|
|
const label numberOfTests = 10000;
|
|
Random rndGen(1234);
|
|
|
|
for (label i = 0; i < numberOfTests; ++i)
|
|
{
|
|
const tensor T(makeRandomContainer(rndGen));
|
|
test_eigen_funcs(T);
|
|
}
|
|
|
|
{
|
|
Info<< nl << " ## A zero tensor: ##"<< nl;
|
|
const tensor zeroT(Zero);
|
|
test_eigen_funcs(zeroT);
|
|
}
|
|
{
|
|
Info<< nl
|
|
<< " ## A skew-symmetric tensor with no-real eigenvalues: ##"
|
|
<< nl;
|
|
const tensor T
|
|
(
|
|
0, 1, 1,
|
|
-1, 0, 1,
|
|
-1, -1, 0
|
|
);
|
|
test_eigen_funcs(T);
|
|
}
|
|
{
|
|
Info<< nl
|
|
<< " ## A stiff tensor: ##"
|
|
<< nl;
|
|
const tensor stiff
|
|
(
|
|
pow(10.0, 10), pow(10.0, 8), pow(10.0, 7),
|
|
pow(10.0, -8), pow(10.0, 9), pow(10.0, -8),
|
|
pow(10.0, 10), pow(10.0, 8), pow(10.0, 7)
|
|
);
|
|
// Although eigendecomposition is successful for the stiff tensors,
|
|
// cross-check between prod(eigenvalues) ?= det(stiff) is inherently
|
|
// problematic; therefore, eigenvalues of the stiff tensors are
|
|
// cross-checked by only sum(eigenvalues) ?= trace(stiff)
|
|
const bool testProd = false;
|
|
test_eigen_funcs(stiff, testProd);
|
|
}
|
|
{
|
|
Info<< nl
|
|
<< " ## Random tensor with tiny off-diag elements: ##"
|
|
<< nl;
|
|
|
|
const List<scalar> epsilons
|
|
({
|
|
0, SMALL, Foam::sqrt(SMALL), sqr(SMALL), Foam::cbrt(SMALL),
|
|
-SMALL, -Foam::sqrt(SMALL), -sqr(SMALL), -Foam::cbrt(SMALL)
|
|
});
|
|
|
|
for (label i = 0; i < numberOfTests; ++i)
|
|
{
|
|
for (const auto& eps : epsilons)
|
|
{
|
|
{
|
|
tensor T(makeRandomContainer(rndGen));
|
|
T.xy() = eps*rndGen.GaussNormal<scalar>();
|
|
test_eigen_funcs(T);
|
|
}
|
|
{
|
|
tensor T(makeRandomContainer(rndGen));
|
|
T.xy() = eps*rndGen.GaussNormal<scalar>();
|
|
T.xz() = eps*rndGen.GaussNormal<scalar>();
|
|
test_eigen_funcs(T);
|
|
}
|
|
{
|
|
tensor T(makeRandomContainer(rndGen));
|
|
T.xy() = eps*rndGen.GaussNormal<scalar>();
|
|
T.xz() = eps*rndGen.GaussNormal<scalar>();
|
|
T.yz() = eps*rndGen.GaussNormal<scalar>();
|
|
test_eigen_funcs(T);
|
|
}
|
|
{
|
|
tensor T(makeRandomContainer(rndGen));
|
|
T.xy() = eps*rndGen.GaussNormal<scalar>();
|
|
T.xz() = eps*rndGen.GaussNormal<scalar>();
|
|
T.yz() = eps*rndGen.GaussNormal<scalar>();
|
|
T.yx() = eps*rndGen.GaussNormal<scalar>();
|
|
test_eigen_funcs(T);
|
|
}
|
|
{
|
|
tensor T(makeRandomContainer(rndGen));
|
|
T.xy() = eps*rndGen.GaussNormal<scalar>();
|
|
T.xz() = eps*rndGen.GaussNormal<scalar>();
|
|
T.yz() = eps*rndGen.GaussNormal<scalar>();
|
|
T.yx() = eps*rndGen.GaussNormal<scalar>();
|
|
T.zx() = eps*rndGen.GaussNormal<scalar>();
|
|
test_eigen_funcs(T);
|
|
}
|
|
{
|
|
tensor T(makeRandomContainer(rndGen));
|
|
T.xy() = eps*rndGen.GaussNormal<scalar>();
|
|
T.xz() = eps*rndGen.GaussNormal<scalar>();
|
|
T.yz() = eps*rndGen.GaussNormal<scalar>();
|
|
T.yx() = eps*rndGen.GaussNormal<scalar>();
|
|
T.zx() = eps*rndGen.GaussNormal<scalar>();
|
|
T.zy() = eps*rndGen.GaussNormal<scalar>();
|
|
test_eigen_funcs(T);
|
|
}
|
|
{
|
|
tensor T(makeRandomContainer(rndGen));
|
|
T.xy() = 0;
|
|
T.xz() = eps*rndGen.GaussNormal<scalar>();
|
|
T.yz() = 0;
|
|
T.yx() = eps*rndGen.GaussNormal<scalar>();
|
|
T.zx() = eps*rndGen.GaussNormal<scalar>();
|
|
T.zy() = 0;
|
|
test_eigen_funcs(T);
|
|
}
|
|
{
|
|
tensor T(makeRandomContainer(rndGen));
|
|
T.xy() = eps;
|
|
T.xz() = eps;
|
|
T.yz() = eps;
|
|
T.yx() = eps;
|
|
T.zx() = eps;
|
|
T.zy() = eps;
|
|
test_eigen_funcs(T);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
if (nFail_)
|
|
{
|
|
Info<< nl << " #### "
|
|
<< "Failed in " << nFail_ << " tests "
|
|
<< "out of total " << nTest_ << " tests "
|
|
<< "####\n" << endl;
|
|
return 1;
|
|
}
|
|
|
|
Info<< nl << " #### Passed all " << nTest_ <<" tests ####\n" << endl;
|
|
return 0;
|
|
}
|
|
|
|
|
|
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|