ENH: correct sumProd return type (#1086)

- previously returned scalar, but now return pTraits cmptType
  which is the same as scalarProduct / outerProduct type.

applications/test/complex/Test-complex.C

@@ -97,6 +97,7 @@ int main(int argc, char *argv[])
         c.Im() *= 5;
     }
 
+    Info<< "sumProd: " << sumProd(fld1, fld2) << nl;
 
     fld1 *= 10;
     Info<< "scalar multiply: " << flatOutput(fld1) << nl;

src/OpenFOAM/fields/Fields/Field/FieldFunctions.C

@@ -423,14 +423,17 @@ Type minMagSqr(const UList<Type>& f)
 TMP_UNARY_FUNCTION(Type, minMagSqr)
 
 template<class Type>
-scalar sumProd(const UList<Type>& f1, const UList<Type>& f2)
+typename pTraits<Type>::cmptType
+sumProd(const UList<Type>& f1, const UList<Type>& f2)
 {
-    scalar SumProd = 0;
+    typedef typename pTraits<Type>::cmptType outType;
+
+    outType result = Zero;
     if (f1.size() && (f1.size() == f2.size()))
     {
-        TFOR_ALL_S_OP_F_OP_F(scalar, SumProd, +=, Type, f1, &&, Type, f2)
+        TFOR_ALL_S_OP_F_OP_F(outType, result, +=, Type, f1, &&, Type, f2)
     }
-    return SumProd;
+    return result;
 }
 
 
@@ -544,16 +547,18 @@ G_UNARY_FUNCTION(scalarMinMax, gMinMaxMag, minMaxMag, minMaxMag)
 
 
 template<class Type>
-scalar gSumProd
+typename pTraits<Type>::cmptType gSumProd
 (
     const UList<Type>& f1,
     const UList<Type>& f2,
     const label comm
 )
 {
-    scalar SumProd = sumProd(f1, f2);
+    typedef typename pTraits<Type>::cmptType outType;
-    reduce(SumProd, sumOp<scalar>(), Pstream::msgType(), comm);
+
-    return SumProd;
+    outType result = sumProd(f1, f2);
+    reduce(result, sumOp<outType>(), Pstream::msgType(), comm);
+    return result;
 }
 
 template<class Type>
@@ -658,7 +663,7 @@ tmp<Field<typename product<Type1, Type2>::type>>                               \
 operator Op(const UList<Type1>& f1, const tmp<Field<Type2>>& tf2)              \
 {                                                                              \
     typedef typename product<Type1, Type2>::type productType;                  \
-    auto tres = reuseTmp<productType, Type2>::New(tf2);     \
+    auto tres = reuseTmp<productType, Type2>::New(tf2);                        \
     OpFunc(tres.ref(), f1, tf2());                                             \
     tf2.clear();                                                               \
     return tres;                                                               \
@@ -669,7 +674,7 @@ tmp<Field<typename product<Type1, Type2>::type>>                               \
 operator Op(const tmp<Field<Type1>>& tf1, const UList<Type2>& f2)              \
 {                                                                              \
     typedef typename product<Type1, Type2>::type productType;                  \
-    auto tres = reuseTmp<productType, Type1>::New(tf1);     \
+    auto tres = reuseTmp<productType, Type1>::New(tf1);                        \
     OpFunc(tres.ref(), tf1(), f2);                                             \
     tf1.clear();                                                               \
     return tres;                                                               \

src/OpenFOAM/fields/Fields/Field/FieldFunctions.H

@@ -197,7 +197,11 @@ TMP_UNARY_FUNCTION(Type, minMagSqr)
 
 
 template<class Type>
-scalar sumProd(const UList<Type>& f1, const UList<Type>& f2);
+typename pTraits<Type>::cmptType sumProd
+(
+    const UList<Type>& f1,
+    const UList<Type>& f2
+);
 
 template<class Type>
 Type sumCmptProd(const UList<Type>& f1, const UList<Type>& f2);
@@ -244,7 +248,7 @@ G_UNARY_FUNCTION(scalarMinMax, gMinMaxMag, minMaxMag, minMaxMag)
 #undef G_UNARY_FUNCTION
 
 template<class Type>
-scalar gSumProd
+typename pTraits<Type>::cmptType gSumProd
 (
     const UList<Type>& f1,
     const UList<Type>& f2,

src/OpenFOAM/fields/Fields/complex/complexField.C

@@ -132,6 +132,18 @@ Foam::scalarField Foam::Im(const UList<complex>& cf)
 namespace Foam
 {
 
+template<>
+complex sumProd(const UList<complex>& f1, const UList<complex>& f2)
+{
+    complex result = Zero;
+    if (f1.size() && (f1.size() == f2.size()))
+    {
+        TFOR_ALL_S_OP_F_OP_F(complex, result, +=, complex, f1, *, complex, f2)
+    }
+    return result;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 UNARY_FUNCTION(complex, complex, pow3)

src/OpenFOAM/fields/Fields/complex/complexField.H

@@ -74,6 +74,12 @@ scalarField Im(const UList<complex>& cf);
 //- Sum real and imag components
 scalarField ReImSum(const UList<complex>& cf);
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Sum product
+template<>
+complex sumProd(const UList<complex>& f1, const UList<complex>& f2);
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/fields/Fields/scalarField/scalarField.C

@@ -92,14 +92,26 @@ tmp<scalarField> stabilise(const tmp<scalarField>& tsf, const scalar s)
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<>
-scalar sumProd(const UList<scalar>& f1, const UList<scalar>& f2)
+float sumProd(const UList<float>& f1, const UList<float>& f2)
 {
-    scalar SumProd = 0.0;
+    float result = 0.0;
     if (f1.size() && (f1.size() == f2.size()))
     {
-        TFOR_ALL_S_OP_F_OP_F(scalar, SumProd, +=, scalar, f1, *, scalar, f2)
+        TFOR_ALL_S_OP_F_OP_F(float, result, +=, float, f1, *, float, f2)
     }
-    return SumProd;
+    return result;
+}
+
+
+template<>
+double sumProd(const UList<double>& f1, const UList<double>& f2)
+{
+    double result = 0.0;
+    if (f1.size() && (f1.size() == f2.size()))
+    {
+        TFOR_ALL_S_OP_F_OP_F(double, result, +=, double, f1, *, double, f2)
+    }
+    return result;
 }
 
 

src/OpenFOAM/fields/Fields/scalarField/scalarField.H

@@ -76,8 +76,13 @@ tmp<scalarField> stabilise(const tmp<scalarField>&, const scalar s);
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+//- Sum product for float
 template<>
-scalar sumProd(const UList<scalar>& f1, const UList<scalar>& f2);
+float sumProd(const UList<float>& f1, const UList<float>& f2);
+
+//- Sum product for double
+template<>
+double sumProd(const UList<double>& f1, const UList<double>& f2);
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //