From 6b5da70602dbbc755dd79ac8c8d58d2b0376f312 Mon Sep 17 00:00:00 2001 From: Mark Olesen Date: Sat, 7 Dec 2019 16:55:18 +0100 Subject: [PATCH] ENH: add scalarOps with divide-by-zero protection - add functor versions of floor/ceil/round for scalar --- applications/test/scalarOps/Make/files | 3 + applications/test/scalarOps/Make/options | 1 + applications/test/scalarOps/Test-scalarOps.C | 135 +++++++++++++++++++ src/OpenFOAM/primitives/Scalar/Scalar.H | 67 ++++++++- src/OpenFOAM/primitives/ops/scalarOps.H | 92 +++++++++++++ 5 files changed, 296 insertions(+), 2 deletions(-) create mode 100644 applications/test/scalarOps/Make/files create mode 100644 applications/test/scalarOps/Make/options create mode 100644 applications/test/scalarOps/Test-scalarOps.C create mode 100644 src/OpenFOAM/primitives/ops/scalarOps.H diff --git a/applications/test/scalarOps/Make/files b/applications/test/scalarOps/Make/files new file mode 100644 index 0000000000..af52899047 --- /dev/null +++ b/applications/test/scalarOps/Make/files @@ -0,0 +1,3 @@ +Test-scalarOps.C + +EXE = $(FOAM_USER_APPBIN)/Test-scalarOps diff --git a/applications/test/scalarOps/Make/options b/applications/test/scalarOps/Make/options new file mode 100644 index 0000000000..9d6f459ad8 --- /dev/null +++ b/applications/test/scalarOps/Make/options @@ -0,0 +1 @@ +/**/ diff --git a/applications/test/scalarOps/Test-scalarOps.C b/applications/test/scalarOps/Test-scalarOps.C new file mode 100644 index 0000000000..63df213177 --- /dev/null +++ b/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 . + +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 +void testDivide(const List>& list) +{ + const scalarDivideOp 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>& 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> scalars + ({ + {10.0, 15}, + {5.0, 15}, + {5.0, 0}, + }); + + List> vectors + ({ + { {1,2,3}, 15}, + { {4,5,6}, 15}, + { {7,8,9}, 0}, + }); + + + Info<< nl << "Test scalar/scalar division" << nl; + testDivide(scalars); + + Info<< nl << "Test scalar/scalar modulo" << nl; + testModulo(scalars); + + + Info<< nl << "Test vector/scalar division" << nl; + testDivide(vectors); + + + Info<< "\nEnd\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/primitives/Scalar/Scalar.H b/src/OpenFOAM/primitives/Scalar/Scalar.H index b21ca2d9a2..5b4a4805ca 100644 --- a/src/OpenFOAM/primitives/Scalar/Scalar.H +++ b/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 + struct floorOp + { + T operator()(const T& x) const WARNRETURN + { + T ret; + for (direction cmpt=0; cmpt < pTraits::nComponents; ++cmpt) + { + component(ret, cmpt) = std::floor(component(x, cmpt)); + } + return ret; + } + }; + \endcode + SourceFiles Scalar.C @@ -459,7 +479,7 @@ struct equalOp 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 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 struct floorOp; + +//- Round scalar downwards - functor version of std::floor +template<> +struct floorOp +{ + Scalar operator()(const Scalar& x) const + { + return std::floor(x); + } +}; + + +// Default definition in ops.H (future?) +template struct ceilOp; + +//- Round scalar upwards - functor version of std::ceil +template<> +struct ceilOp +{ + Scalar operator()(const Scalar& x) const + { + return std::ceil(x); + } +}; + + +// Default definition in ops.H (future?) +template struct roundOp; + +//- Round scalar - functor version of std::round +template<> +struct roundOp +{ + Scalar operator()(const Scalar& x) const + { + return std::round(x); } }; diff --git a/src/OpenFOAM/primitives/ops/scalarOps.H b/src/OpenFOAM/primitives/ops/scalarOps.H new file mode 100644 index 0000000000..6519a48d61 --- /dev/null +++ b/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 . + +InNamespace + Foam + +Description + Functors that are scalar-specific. + +\*---------------------------------------------------------------------------*/ + +#ifndef scalarOps_H +#define scalarOps_H + +#include "scalar.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +//- Hypot operation (scalar only) +template +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 +struct scalarDivideOp +{ + T operator()(const T& x, const T2& y) const + { + return (x / stabilise(y, pTraits::vsmall)); + } +}; + + +//- Floating point modulo operation with divide-by-zero protection +template +struct scalarModuloOp +{ + T operator()(const T& x, const T2& y) const + { + if (Foam::mag(y) < pTraits::vsmall) + { + return pTraits::zero; + } + return std::fmod(x, y); + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* //