ENH: add scalarOps with divide-by-zero protection

- add functor versions of floor/ceil/round for scalar

applications/test/scalarOps/Make/files

@@ -0,0 +1,3 @@
+Test-scalarOps.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-scalarOps

applications/test/scalarOps/Make/options

@@ -0,0 +1 @@
+/**/

applications/test/scalarOps/Test-scalarOps.C

@@ -0,0 +1,135 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 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-scalarOps
+
+Description
+    Test scalar-only ops
+
+\*---------------------------------------------------------------------------*/
+
+#include "IOstreams.H"
+#include "labelList.H"
+#include "scalarList.H"
+#include "FlatOutput.H"
+#include "Tuple2.H"
+#include "ops.H"
+#include "scalarOps.H"
+#include "vector.H"
+#include "Tuple2.H"
+
+using namespace Foam;
+
+
+template<class T>
+void testDivide(const List<Tuple2<T, scalar>>& list)
+{
+    const scalarDivideOp<T> bop;
+
+    for (const auto& pair : list)
+    {
+        Info<< "num=" << pair.first()
+            << " den=" << pair.second() << flush;
+
+        Info<< " = " << bop(pair.first(), pair.second())
+            << endl;
+    }
+
+    Info<< "----" << nl;
+
+    for (const auto& pair : list)
+    {
+        Info<< "num=" << pair.first()
+            << " den=" << pair.second() << flush;
+
+        Info<< " = " << (pair.first() / pair.second()) << endl;
+    }
+}
+
+
+void testModulo(const List<Tuple2<scalar, scalar>>& list)
+{
+    const scalarModuloOp<> bop;
+
+    for (const auto& pair : list)
+    {
+        Info<< "num=" << pair.first()
+            << " den=" << pair.second() << flush;
+
+        Info<< " = " << bop(pair.first(), pair.second())
+            << endl;
+    }
+
+    Info<< "----" << nl;
+
+    for (const auto& pair : list)
+    {
+        Info<< "num=" << pair.first()
+            << " den=" << pair.second() << flush;
+
+        Info<< " = " << std::fmod(pair.first(), pair.second()) << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+//  Main program:
+
+int main(int argc, char *argv[])
+{
+    List<Tuple2<scalar, scalar>> scalars
+    ({
+        {10.0, 15},
+        {5.0, 15},
+        {5.0, 0},
+    });
+
+    List<Tuple2<vector, scalar>> vectors
+    ({
+        { {1,2,3}, 15},
+        { {4,5,6}, 15},
+        { {7,8,9}, 0},
+    });
+
+
+    Info<< nl << "Test scalar/scalar division" << nl;
+    testDivide<scalar>(scalars);
+
+    Info<< nl << "Test scalar/scalar modulo" << nl;
+    testModulo(scalars);
+
+
+    Info<< nl << "Test vector/scalar division" << nl;
+    testDivide<vector>(vectors);
+
+
+    Info<< "\nEnd\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/primitives/Scalar/Scalar.H

@@ -30,6 +30,26 @@ Typedef
 Description
     Floating-point number (float or double)
 
+Note
+    The floor/ceil/round operations could easily be extended to handle
+    VectorSpace when the need arises. For example,
+
+    \code
+    template<class T>
+    struct floorOp
+    {
+        T operator()(const T& x) const WARNRETURN
+        {
+            T ret;
+            for (direction cmpt=0; cmpt < pTraits<T>::nComponents; ++cmpt)
+            {
+                component(ret, cmpt) = std::floor(component(x, cmpt));
+            }
+            return ret;
+        }
+    };
+    \endcode
+
 SourceFiles
     Scalar.C
 
@@ -459,7 +479,7 @@ struct equalOp<Scalar>
 
     bool operator()(const Scalar& a, const Scalar& b) const
     {
-        return Foam::mag(a - b) <= tolerance;
+        return mag(a - b) <= tolerance;
     }
 };
 
@@ -481,7 +501,50 @@ struct notEqualOp<Scalar>
 
     bool operator()(const Scalar& a, const Scalar& b) const
     {
-        return Foam::mag(a - b) > tolerance;
+        return mag(a - b) > tolerance;
+    }
+};
+
+
+
+// Default definition in ops.H (future?)
+template<class T> struct floorOp;
+
+//- Round scalar downwards - functor version of std::floor
+template<>
+struct floorOp<Scalar>
+{
+    Scalar operator()(const Scalar& x) const
+    {
+        return std::floor(x);
+    }
+};
+
+
+// Default definition in ops.H (future?)
+template<class T> struct ceilOp;
+
+//- Round scalar upwards - functor version of std::ceil
+template<>
+struct ceilOp<Scalar>
+{
+    Scalar operator()(const Scalar& x) const
+    {
+        return std::ceil(x);
+    }
+};
+
+
+// Default definition in ops.H (future?)
+template<class T> struct roundOp;
+
+//- Round scalar - functor version of std::round
+template<>
+struct roundOp<Scalar>
+{
+    Scalar operator()(const Scalar& x) const
+    {
+        return std::round(x);
     }
 };
 

src/OpenFOAM/primitives/ops/scalarOps.H

@@ -0,0 +1,92 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2019 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/>.
+
+InNamespace
+    Foam
+
+Description
+    Functors that are scalar-specific.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef scalarOps_H
+#define scalarOps_H
+
+#include "scalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//- Hypot operation (scalar only)
+template<class T>
+struct hypotOp
+{
+    T operator()(const T& x, const T& y) const
+    {
+        return std::hypot(x, y);
+    }
+};
+
+
+//- Scalar division with divide-by-zero protection
+//  Uses stabilise, but could also handle as per modulo and return zero
+template<class T, class T2=Foam::scalar>
+struct scalarDivideOp
+{
+    T operator()(const T& x, const T2& y) const
+    {
+        return (x / stabilise(y, pTraits<T2>::vsmall));
+    }
+};
+
+
+//- Floating point modulo operation with divide-by-zero protection
+template<class T, class T2=Foam::scalar>
+struct scalarModuloOp
+{
+    T operator()(const T& x, const T2& y) const
+    {
+        if (Foam::mag(y) < pTraits<T2>::vsmall)
+        {
+            return pTraits<T>::zero;
+        }
+        return std::fmod(x, y);
+    }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //