BUG: incorrect scalar/complex division (#1331)

ENH: define addition/subtraction operations for scalar and complex

- required since construct complex from scalar is explicit
- additional tests in Test-complex

applications/test/complex/Test-complex.C

@@ -118,75 +118,135 @@ int main(int argc, char *argv[])
     // Info<< "sqrt   : " << sqrt(fld1) << nl;
     // Info<< "pow(2) : " << pow(fld1, 2) << nl;
 
-    Info<< nl << "Elementary complex arithmetic operations:" << nl;
+    #if 1
+    Info<< nl << "## Elementary complex-complex arithmetic operations:" << nl;
     {
-        complex a(6, 1);
+        const complex a(6, 1);
         complex b = a;
 
-        Info << "Compound assignment operations:" << nl;
+        Info << "# Compound assignment operations:" << nl;
 
-        // Multiplication
-        b *= a;
-        Info<< "b *= a:" << tab << "b=" << b << nl;
+        Info<< "a = " << a << ", b = " << b << nl;
 
         // Addition
         b += a;
-        Info<< "b += a:" << tab << "b=" << b << nl;
+        Info<< "b += a:" << tab << "b =" << b << nl;
 
         // Subtraction
         b -= a;
-        Info<< "b -= a:" << tab << "b=" << b << nl;
+        Info<< "b -= a:" << tab << "b =" << b << nl;
+
+        // Multiplication
+        b *= a;
+        Info<< "b *= a:" << tab << "b =" << b << nl;
 
         // Division
         b /= a;
-        Info<< "b /= a:" << tab << "b=" << b << nl;
+        Info<< "b /= a:" << tab << "b =" << b << nl;
+    }
+    #endif
+
+
+    #if 1
+    Info<< nl << "## Elementary complex-scalar arithmetic operations:" << nl;
+    {
+        const scalar a = 5;
+        complex b(6, 1);
+
+        Info << "# Non-assignment operations:" << nl;
+
+        Info<< "(scalar) a = " << a << ", b = " << b << nl;
+
+        // Addition
+        b = a + b;
+        Info<< "b = a + b: " << tab << b << nl;
 
-        Info << "Operations with scalars:" << nl;
+        b = b + a;
+        Info<< "b = b + a: " << tab << b << nl;
 
-        Info<< "b=" << b << nl;
+        // Subtraction
+        b = a - b;
+        Info<< "b = a - b: " << tab << b << nl;
 
-        // Scalar multiplication
-        b *= 2.0;
-        Info<< "b*2 (elementwise multiplication):" << tab << b << nl;
+        b = b - a;
+        Info<< "b = b - a: " << tab << b << nl;
 
-        // Scalar addition
-        b += 1.0;
-        Info<< "b + 1 (only real part):" << tab << b << nl;
+        // Multiplication
+        b = a*b;
+        Info<< "b = a*b: " << tab << b << nl;
 
-        // Scalar subtraction
-        b -= 1.0;
-        Info<< "b - 1 (only real part):" << tab << b << nl;
+        b = b*a;
+        Info<< "b = b*a: " << tab << b << nl;
+
+        // Division
+        b = a/b;
+        Info<< "b = a/b = scalar(a)/b = complex(a)/b:" << tab << b << nl;
+
+        b = b/a;
+        Info<< "b = b/a: " << tab << b << nl;
+
+
+        Info << "# Compound assignment operations:" << nl;
+
+        Info<< "(scalar) a = " << a << ", b = " << b << nl;
+
+        // Addition: complex+scalar
+        b += a;
+        Info<< "b += a (only real part):" << tab << b << nl;
+
+        // Subtraction: complex-scalar
+        b -= a;
+        Info<< "b -= a (only real part):" << tab << b << nl;
+
+        // Multiplication: complex*scalar
+        b *= a;
+        Info<< "b *= a (real and imag parts):" << tab << b << nl;
+
+        // Division: complex/scalar
+        b /= a;
+        Info<< "b /= a (real and imag parts):" << tab << b << nl;
 
-        // Scalar division
-        b = 1.0/b;
-        Info<< "1/b (elementwise division):" << tab << b << nl;
     }
+    #endif
+
 
-    Info<< nl << "Other mathematical expressions:" << nl;
+    #if 1
+    Info<< nl << "## Other mathematical expressions:" << nl;
     {
-        complex a(4.3, -3.14);
-        complex b(0, -4.3);
+        const complex a(4.3, -3.14);
+        const complex b(0, -4.3);
+        const complex c(-4.3, 0);
 
-        Info<< "a=" << a << tab << "b=" << b << nl;
+        Info<< "a = " << a << ", b = " << b << ", c = " << c << nl;
 
         // Square-root
-        //Info<< "sqrt(a)=" << sqrt(a) << tab << "sqrt(b)=" << sqrt(b) << nl;
+        Info<< "sqrt(a) = " << Foam::sqrt(a) << ", "
+            << "sqrt(b) = " << Foam::sqrt(b) << ", "
+            << "sqrt(c) = " << Foam::sqrt(c) << nl;
 
         // Square
-        Info<< "sqr(a)=" << sqr(a) << tab << "sqr(b)=" << sqr(b) << nl;
+        Info<< "sqr(a) = " << sqr(a) << ", "
+            << "sqr(b) = " << sqr(b) << ", "
+            << "sqr(c) = " << sqr(c) << nl;
 
         // n^th power
-        //Info<< "pow(a,-1)=" << pow(a,-1) << tab
-        //      << "pow(b,-1)=" << pow(b,-1) << nl;
+        Info<< "pow(a, -1) = " << pow(a, -1) << ", "
+            << "pow(b, -1) = " << pow(b, -1) << ", "
+            << "pow(c, -1) = " << pow(c, -1) << nl;
 
         // Exponential
-        //Info<< "exp(a)=" << exp(a) << tab << "exp(b)=" << exp(b) << nl;
+        Info<< "exp(a) = " << exp(a) << ", "
+            << "exp(b) = " << exp(b) << ", "
+            << "exp(c) = " << exp(c) << nl;
 
         // Natural logarithm
-        //Info<< "log(a)=" << log(a) << tab << "log(b)=" << log(b) << nl;
+        Info<< "log(a) = " << log(a) << ", "
+            << "log(b) = " << log(b) << ", "
+            << "log(c) = " << log(c) << nl;
+
     }
+    #endif
 
-    Info<< nl << "End" << nl;
 
     // Make some changes
     {

src/OpenFOAM/primitives/complex/complex.H

@@ -60,13 +60,17 @@ inline const complex& min(const complex&, const complex&);
 inline const complex& max(const complex&, const complex&);
 inline complex limit(const complex&, const complex&);
 inline const complex& sum(const complex&);
-inline complex operator+(const complex&, const complex&);
 inline complex operator-(const complex&);
+inline complex operator+(const complex&, const complex&);
+inline complex operator+(const complex&, const scalar);
+inline complex operator+(const scalar, const complex&);
 inline complex operator-(const complex&, const complex&);
+inline complex operator-(const complex&, const scalar);
+inline complex operator-(const scalar, const complex&);
 inline complex operator*(const complex&, const complex&);
-inline complex operator/(const complex&, const complex&);
-inline complex operator*(const scalar, const complex&);
 inline complex operator*(const complex&, const scalar);
+inline complex operator*(const scalar, const complex&);
+inline complex operator/(const complex&, const complex&);
 inline complex operator/(const complex&, const scalar);
 inline complex operator/(const scalar, const complex&);
 
@@ -174,15 +178,19 @@ public:
         //- Assign zero
         inline void operator=(const Foam::zero);
 
-        inline void operator+=(const complex& c);
-        inline void operator-=(const complex& c);
-        inline void operator*=(const complex& c);
-        inline void operator/=(const complex& c);
-
+        //- Assign scalar (imag = zero)
         inline void operator=(const scalar s);
+
+        inline void operator+=(const complex& c);
         inline void operator+=(const scalar s);
+
+        inline void operator-=(const complex& c);
         inline void operator-=(const scalar s);
+
+        inline void operator*=(const complex& c);
         inline void operator*=(const scalar s);
+
+        inline void operator/=(const complex& c);
         inline void operator/=(const scalar s);
 
         inline bool operator==(const complex& c) const;
@@ -209,14 +217,21 @@ public:
 
     // Friend Operators
 
-        friend complex operator+(const complex& c1, const complex& c2);
         friend complex operator-(const complex& c);
+
+        friend complex operator+(const complex& c1, const complex& c2);
+        friend complex operator+(const complex& c, const scalar s);
+        friend complex operator+(const scalar s, const complex& c);
+
         friend complex operator-(const complex& c1, const complex& c2);
-        friend complex operator*(const complex& c1, const complex& c2);
-        friend complex operator/(const complex& c1, const complex& c2);
+        friend complex operator-(const complex& c, const scalar s);
+        friend complex operator-(const scalar s, const complex& c);
 
-        friend complex operator*(const scalar s, const complex& c);
+        friend complex operator*(const complex& c1, const complex& c2);
         friend complex operator*(const complex& c, const scalar s);
+        friend complex operator*(const scalar s, const complex& c);
+
+        friend complex operator/(const complex& c1, const complex& c2);
         friend complex operator/(const complex& c, const scalar s);
         friend complex operator/(const scalar s, const complex& c);
 };

src/OpenFOAM/primitives/complex/complexI.H

@@ -129,55 +129,55 @@ inline void Foam::complex::operator=(const Foam::zero)
 }
 
 
-inline void Foam::complex::operator+=(const complex& c)
+inline void Foam::complex::operator=(const scalar s)
 {
-    re += c.re;
-    im += c.im;
+    re = s;
+    im = 0;
 }
 
 
-inline void Foam::complex::operator-=(const complex& c)
+inline void Foam::complex::operator+=(const complex& c)
 {
-    re -= c.re;
-    im -= c.im;
+    re += c.re;
+    im += c.im;
 }
 
 
-inline void Foam::complex::operator*=(const complex& c)
+inline void Foam::complex::operator+=(const scalar s)
 {
-    *this = (*this)*c;
+    re += s;
 }
 
 
-inline void Foam::complex::operator/=(const complex& c)
+inline void Foam::complex::operator-=(const complex& c)
 {
-    *this = *this/c;
+    re -= c.re;
+    im -= c.im;
 }
 
 
-inline void Foam::complex::operator=(const scalar s)
+inline void Foam::complex::operator-=(const scalar s)
 {
-    re = s;
-    im = 0;
+    re -= s;
 }
 
 
-inline void Foam::complex::operator+=(const scalar s)
+inline void Foam::complex::operator*=(const complex& c)
 {
-    re += s;
+    *this = (*this)*c;
 }
 
 
-inline void Foam::complex::operator-=(const scalar s)
+inline void Foam::complex::operator*=(const scalar s)
 {
-    re -= s;
+    re *= s;
+    im *= s;
 }
 
 
-inline void Foam::complex::operator*=(const scalar s)
+inline void Foam::complex::operator/=(const complex& c)
 {
-    re *= s;
-    im *= s;
+    *this = *this/c;
 }
 
 
@@ -288,6 +288,12 @@ inline complex transform(const Tensor<scalar>&, const complex c)
 
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
+inline complex operator-(const complex& c)
+{
+    return complex(-c.re, -c.im);
+}
+
+
 inline complex operator+(const complex& c1, const complex& c2)
 {
     return complex
@@ -298,13 +304,15 @@ inline complex operator+(const complex& c1, const complex& c2)
 }
 
 
-inline complex operator-(const complex& c)
+inline complex operator+(const complex& c, const scalar s)
 {
-    return complex
-    (
-        -c.re,
-        -c.im
-    );
+    return complex(c.re + s, c.im);
+}
+
+
+inline complex operator+(const scalar s, const complex& c)
+{
+    return complex(c.re + s, c.im);
 }
 
 
@@ -318,6 +326,18 @@ inline complex operator-(const complex& c1, const complex& c2)
 }
 
 
+inline complex operator-(const complex& c, const scalar s)
+{
+    return complex(c.re - s, c.im);
+}
+
+
+inline complex operator-(const scalar s, const complex& c)
+{
+    return complex(s - c.re, -c.im);
+}
+
+
 inline complex operator*(const complex& c1, const complex& c2)
 {
     return complex
@@ -328,15 +348,9 @@ inline complex operator*(const complex& c1, const complex& c2)
 }
 
 
-inline complex operator/(const complex& c1, const complex& c2)
+inline complex operator*(const complex& c, const scalar s)
 {
-    const scalar sqrC2 = magSqr(c2);
-
-    return complex
-    (
-        (c1.re*c2.re + c1.im*c2.im)/sqrC2,
-        (c1.im*c2.re - c1.re*c2.im)/sqrC2
-    );
+    return complex(s*c.re, s*c.im);
 }
 
 
@@ -346,9 +360,15 @@ inline complex operator*(const scalar s, const complex& c)
 }
 
 
-inline complex operator*(const complex& c, const scalar s)
+inline complex operator/(const complex& c1, const complex& c2)
 {
-    return complex(s*c.re, s*c.im);
+    const scalar sqrC2 = magSqr(c2);
+
+    return complex
+    (
+        (c1.re*c2.re + c1.im*c2.im)/sqrC2,
+        (c1.im*c2.re - c1.re*c2.im)/sqrC2
+    );
 }
 
 
@@ -360,7 +380,7 @@ inline complex operator/(const complex& c, const scalar s)
 
 inline complex operator/(const scalar s, const complex& c)
 {
-    return complex(s/c.re, s/c.im);
+    return complex(s)/c;
 }
 
 

[Kutalmis Bercin]

A note for the next generation: Please note that due to the squashing the diffs were messed up. Please manually compare the commits in comparison to the prev files, if necessary.