diff --git a/applications/test/complex/Test-complex.C b/applications/test/complex/Test-complex.C index 938f63fb8eb3db5d501368b56bf9b59bbb54fa9d..8a4442335b32b12b0da5ddeadb1b66488e97f700 100644 --- a/applications/test/complex/Test-complex.C +++ b/applications/test/complex/Test-complex.C @@ -24,11 +24,12 @@ License Application Description - Some tests for complex numbers + Tests for complex numbers \*---------------------------------------------------------------------------*/ #include "argList.H" +#include "complex.H" #include "complexFields.H" #include "ops.H" #include "ListOps.H" @@ -41,8 +42,7 @@ void print1(const complex& z) } -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -// Main program: +// * * * * * * * * * * * * * * * Main Program * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { @@ -113,11 +113,81 @@ int main(int argc, char *argv[]) // Info<< "operator * : " << (fld1 * fld2) << nl; // Info<< "operator / : " << (fld1 / fld2) << nl; - Info<< "operator / : " << (fld1 / 2) << nl; + // Info<< "operator / : " << (fld1 / 2) << nl; // Info<< "operator / : " << (fld1 / fld2) << nl; - Info<< "sqrt : " << sqrt(fld1) << nl; + // Info<< "sqrt : " << sqrt(fld1) << nl; // Info<< "pow(2) : " << pow(fld1, 2) << nl; + Info<< nl << "Elementary complex arithmetic operations:" << nl; + { + complex a (6, 1); + complex b = a; + + Info << "Compound assignment operations:" << nl; + + // Multiplication + b *= a; + Info<< "b *= a:" << tab << "b=" << b << nl; + + // Addition + b += a; + Info<< "b += a:" << tab << "b=" << b << nl; + + // Subtraction + b -= a; + Info<< "b -= a:" << tab << "b=" << b << nl; + + // Division + b /= a; + Info<< "b /= a:" << tab << "b=" << b << nl; + + Info << "Operations with scalars:" << nl; + + Info<< "b=" << b << nl; + + // Scalar multiplication + b *= 2.0; + Info<< "b*2 (elementwise multiplication):" << tab << b << nl; + + // Scalar addition + b += 1.0; + Info<< "b + 1 (only real part):" << tab << b << nl; + + // Scalar subtraction + b -= 1.0; + Info<< "b - 1 (only real part):" << tab << b << nl; + + // Scalar division + b = 1.0/b; + Info<< "1/b (elementwise division):" << tab << b << nl; + + } + + Info<< nl << "Other mathematical expressions:" << nl; + { + complex a (4.3, -3.14); + complex b (0, -4.3); + + Info<< "a=" << a << tab << "b=" << b << nl; + + // Square-root + //Info<< "sqrt(a)=" << sqrt(a) << tab << "sqrt(b)=" << sqrt(b) << nl; + + // Square + Info<< "sqr(a)=" << sqr(a) << tab << "sqr(b)=" << sqr(b) << nl; + + // n^th power + //Info<< "pow(a,-1)=" << pow(a,-1) << tab + // << "pow(b,-1)=" << pow(b,-1) << nl; + + // Exponential + //Info<< "exp(a)=" << exp(a) << tab << "exp(b)=" << exp(b) << nl; + + // Natural logarithm + //Info<< "log(a)=" << log(a) << tab << "log(b)=" << log(b) << nl; + } + + Info<< nl << "End" << nl; // Make some changes { @@ -150,6 +220,7 @@ int main(int argc, char *argv[]) // Info<< "min/max = " << MinMax<complex>(fld1) << nl; Info<< "\nEnd\n" << endl; + return 0; } diff --git a/src/OpenFOAM/primitives/complex/complex.C b/src/OpenFOAM/primitives/complex/complex.C index 3926ff29d20686e5cc98313233dc0b8577a07a6d..df3c0ca872623a7e60357395c4f86d2ab6d09a36 100644 --- a/src/OpenFOAM/primitives/complex/complex.C +++ b/src/OpenFOAM/primitives/complex/complex.C @@ -30,7 +30,6 @@ License // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // -const char* const Foam::complex::typeName = "complex"; const Foam::complex Foam::complex::zero(0, 0); const Foam::complex Foam::complex::one(1, 0); @@ -69,7 +68,6 @@ Foam::pTraits<Foam::complex>::pTraits(Istream& is) } - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::complex::complex(Istream& is) @@ -78,7 +76,7 @@ Foam::complex::complex(Istream& is) } -// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * // Foam::word Foam::name(const complex& c) { diff --git a/src/OpenFOAM/primitives/complex/complex.H b/src/OpenFOAM/primitives/complex/complex.H index c67eafb6ccfab1ce81c3cb6eb39fbe553b4e4b6b..c765813a1ce6d73924dfb2b9c798a5a8a0486cc2 100644 --- a/src/OpenFOAM/primitives/complex/complex.H +++ b/src/OpenFOAM/primitives/complex/complex.H @@ -39,6 +39,7 @@ SourceFiles #define complex_H #include <complex> +#include <type_traits> #include "scalar.H" #include "word.H" #include "zero.H" @@ -91,7 +92,7 @@ public: // Static Data Members //- The type name is "complex" - static const char* const typeName; + static constexpr const char* const typeName = "complex"; //- A complex zero (0,0) static const complex zero; @@ -209,12 +210,17 @@ public: friend scalar magSqr(const complex& c); friend scalar mag(const complex& c); friend complex sqr(const complex& c); - friend const complex& min(const complex&, const complex&); - friend const complex& max(const complex&, const complex&); - friend complex limit(const complex&, const complex&); + //- sgn() https://en.wikipedia.org/wiki/Sign_function#Complex_signum + friend complex sign(const complex& c); + + //- csgn() https://en.wikipedia.org/wiki/Sign_function#Complex_signum + friend scalar csign(const complex& c); - friend const complex& sum(const complex&); + friend const complex& min(const complex& c1, const complex& c2); + friend const complex& max(const complex& c1, const complex& c2); + friend complex limit(const complex& c1, const complex& c2); + friend const complex& sum(const complex& c); // Friend Operators @@ -232,6 +238,10 @@ public: }; +/*---------------------------------------------------------------------------*\ + Class pTraits<complex> Declaration +\*---------------------------------------------------------------------------*/ + // Template specialisation for pTraits<complex> template<> class pTraits<complex> @@ -247,7 +257,7 @@ public: typedef label labelType; - // Member constants + // Member Constants //- Dimensionality of space static constexpr direction dim = 3; @@ -296,6 +306,41 @@ public: }; +/*---------------------------------------------------------------------------*\ + Namespace Detail +\*---------------------------------------------------------------------------*/ + +namespace Detail +{ + // Helper functions for complex, in Detail namespace to avoid possible + // name collisions (could change in the future) + + //- The 'conjugate' of non-complex returns itself (pass-through) + //- it does not return a complex! + template<class T> + typename std::enable_if + < + !std::is_same<complex, T>::value, + const T& + >::type conj(const T& val) + { + return val; + } + + //- The conjugate of a complex number + template<class T> + typename std::enable_if + < + std::is_same<complex, T>::value, + complex + >::type conj(const T& val) + { + return val.conjugate(); + } + +} // End namespace Detail + + // * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * // Istream& operator>>(Istream& is, complex& c); diff --git a/src/OpenFOAM/primitives/complex/complexI.H b/src/OpenFOAM/primitives/complex/complexI.H index 7038d12da33bba227f1dbd9c3c14a8b3388ec8c6..7780282dee0868aa2e2621a7d9e8c78afad1865c 100644 --- a/src/OpenFOAM/primitives/complex/complexI.H +++ b/src/OpenFOAM/primitives/complex/complexI.H @@ -221,13 +221,26 @@ inline scalar magSqr(const complex& c) inline scalar mag(const complex& c) { - return sqrt(magSqr(c)); // OR std::hypot(c.re, c.im); + return std::hypot(c.re, c.im); } inline complex sqr(const complex& c) { - return c * c; + return c*c; +} + + +inline complex sign(const complex& c) +{ + const scalar s(mag(c)); + return s < ROOTVSMALL ? Zero : c/s; +} + + +inline scalar csign(const complex& c) +{ + return equal(c.Re(), 0) ? sign(c.Im()) : sign(c.Re()); }