From 23e2d0bdff0bbfb4777da7c4960cb7cb4d8f06da Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Mon, 23 May 2016 12:04:02 +0100
Subject: [PATCH] DimensionedScalarField, GeometricScalarField: Added more
 rigorous dimension-checking for 'pow' functions

---
 .../DimensionedScalarField.C                  | 130 +++++++++++++++---
 .../GeometricScalarField.C                    | 126 ++++++++++++++---
 2 files changed, 212 insertions(+), 44 deletions(-)

diff --git a/src/OpenFOAM/fields/DimensionedFields/DimensionedScalarField/DimensionedScalarField.C b/src/OpenFOAM/fields/DimensionedFields/DimensionedScalarField/DimensionedScalarField.C
index 9eb3fcf53f..d3d20ce0e2 100644
--- a/src/OpenFOAM/fields/DimensionedFields/DimensionedScalarField/DimensionedScalarField.C
+++ b/src/OpenFOAM/fields/DimensionedFields/DimensionedScalarField/DimensionedScalarField.C
@@ -106,6 +106,20 @@ tmp<DimensionedField<scalar, GeoMesh>> pow
     const DimensionedField<scalar, GeoMesh>& dsf2
 )
 {
+    if (!dsf1.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Base field is not dimensionless: " << dsf1.dimensions()
+            << exit(FatalError);
+    }
+
+    if (!dsf2.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Exponent field is not dimensionless: " << dsf2.dimensions()
+            << exit(FatalError);
+    }
+
     tmp<DimensionedField<scalar, GeoMesh>> tPow
     (
         new DimensionedField<scalar, GeoMesh>
@@ -117,11 +131,7 @@ tmp<DimensionedField<scalar, GeoMesh>> pow
                 dsf1.db()
             ),
             dsf1.mesh(),
-            pow
-            (
-                dsf1.dimensions(),
-                dimensionedScalar("1", 1.0, dsf2.dimensions())
-            )
+            dimless
         )
     );
 
@@ -140,15 +150,25 @@ tmp<DimensionedField<scalar, GeoMesh>> pow
 {
     const DimensionedField<scalar, GeoMesh>& dsf1 = tdsf1();
 
+    if (!dsf1.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Base field is not dimensionless: " << dsf1.dimensions()
+            << exit(FatalError);
+    }
+
+    if (!dsf2.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Exponent field is not dimensionless: " << dsf2.dimensions()
+            << exit(FatalError);
+    }
+
     tmp<DimensionedField<scalar, GeoMesh>> tPow = New
     (
         tdsf1,
         "pow(" + dsf1.name() + ',' + dsf2.name() + ')',
-        pow
-        (
-            dsf1.dimensions(),
-            dimensionedScalar("1", 1.0, dsf2.dimensions())
-        )
+        dimless
     );
 
     pow(tPow.ref().field(), dsf1.field(), dsf2.field());
@@ -168,15 +188,25 @@ tmp<DimensionedField<scalar, GeoMesh>> pow
 {
     const DimensionedField<scalar, GeoMesh>& dsf2 = tdsf2();
 
+    if (!dsf1.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Base field is not dimensionless: " << dsf1.dimensions()
+            << exit(FatalError);
+    }
+
+    if (!dsf2.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Exponent field is not dimensionless: " << dsf2.dimensions()
+            << exit(FatalError);
+    }
+
     tmp<DimensionedField<scalar, GeoMesh>> tPow = New
     (
         tdsf2,
         "pow(" + dsf1.name() + ',' + dsf2.name() + ')',
-        pow
-        (
-            dsf1.dimensions(),
-            dimensionedScalar("1", 1.0, dsf2.dimensions())
-        )
+        dimless
     );
 
     pow(tPow.ref().field(), dsf1.field(), dsf2.field());
@@ -186,6 +216,7 @@ tmp<DimensionedField<scalar, GeoMesh>> pow
     return tPow;
 }
 
+
 template<class GeoMesh>
 tmp<DimensionedField<scalar, GeoMesh>> pow
 (
@@ -196,6 +227,20 @@ tmp<DimensionedField<scalar, GeoMesh>> pow
     const DimensionedField<scalar, GeoMesh>& dsf1 = tdsf1();
     const DimensionedField<scalar, GeoMesh>& dsf2 = tdsf2();
 
+    if (!dsf1.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Base field is not dimensionless: " << dsf1.dimensions()
+            << exit(FatalError);
+    }
+
+    if (!dsf2.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Exponent field is not dimensionless: " << dsf2.dimensions()
+            << exit(FatalError);
+    }
+
     tmp<DimensionedField<scalar, GeoMesh>> tPow =
         reuseTmpTmpDimensionedField<scalar, scalar, scalar, scalar, GeoMesh>::
         New
@@ -203,11 +248,7 @@ tmp<DimensionedField<scalar, GeoMesh>> pow
             tdsf1,
             tdsf2,
             "pow(" + dsf1.name() + ',' + dsf2.name() + ')',
-            pow
-            (
-                dsf1.dimensions(),
-                dimensionedScalar("1", 1.0, dsf2.dimensions())
-            )
+            dimless
         );
 
     pow(tPow.ref().field(), dsf1.field(), dsf2.field());
@@ -226,6 +267,13 @@ tmp<DimensionedField<scalar, GeoMesh>> pow
     const dimensionedScalar& ds
 )
 {
+    if (!ds.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Exponent is not dimensionless: " << ds.dimensions()
+            << exit(FatalError);
+    }
+
     tmp<DimensionedField<scalar, GeoMesh>> tPow
     (
         new DimensionedField<scalar, GeoMesh>
@@ -246,6 +294,7 @@ tmp<DimensionedField<scalar, GeoMesh>> pow
     return tPow;
 }
 
+
 template<class GeoMesh>
 tmp<DimensionedField<scalar, GeoMesh>> pow
 (
@@ -253,6 +302,13 @@ tmp<DimensionedField<scalar, GeoMesh>> pow
     const dimensionedScalar& ds
 )
 {
+    if (!ds.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Exponent is not dimensionless: " << ds.dimensions()
+            << exit(FatalError);
+    }
+
     const DimensionedField<scalar, GeoMesh>& dsf = tdsf();
 
     tmp<DimensionedField<scalar, GeoMesh>> tPow = New
@@ -269,6 +325,7 @@ tmp<DimensionedField<scalar, GeoMesh>> pow
     return tPow;
 }
 
+
 template<class GeoMesh>
 tmp<DimensionedField<scalar, GeoMesh>> pow
 (
@@ -279,6 +336,7 @@ tmp<DimensionedField<scalar, GeoMesh>> pow
     return pow(dsf, dimensionedScalar(s));
 }
 
+
 template<class GeoMesh>
 tmp<DimensionedField<scalar, GeoMesh>> pow
 (
@@ -297,6 +355,20 @@ tmp<DimensionedField<scalar, GeoMesh>> pow
     const DimensionedField<scalar, GeoMesh>& dsf
 )
 {
+    if (!ds.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Base scalar is not dimensionless: " << ds.dimensions()
+            << exit(FatalError);
+    }
+
+    if (!dsf.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Exponent field is not dimensionless: " << dsf.dimensions()
+            << exit(FatalError);
+    }
+
     tmp<DimensionedField<scalar, GeoMesh>> tPow
     (
         new DimensionedField<scalar, GeoMesh>
@@ -308,7 +380,7 @@ tmp<DimensionedField<scalar, GeoMesh>> pow
                 dsf.db()
             ),
             dsf.mesh(),
-            pow(ds, dsf.dimensions())
+            dimless
         )
     );
 
@@ -327,11 +399,25 @@ tmp<DimensionedField<scalar, GeoMesh>> pow
 {
     const DimensionedField<scalar, GeoMesh>& dsf = tdsf();
 
+    if (!ds.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Base scalar is not dimensionless: " << ds.dimensions()
+            << exit(FatalError);
+    }
+
+    if (!dsf.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Exponent field is not dimensionless: " << dsf.dimensions()
+            << exit(FatalError);
+    }
+
     tmp<DimensionedField<scalar, GeoMesh>> tPow = New
     (
         tdsf,
         "pow(" + ds.name() + ',' + dsf.name() + ')',
-        pow(ds, dsf.dimensions())
+        dimless
     );
 
     pow(tPow.ref().field(), ds.value(), dsf.field());
diff --git a/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C b/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C
index f05d3c5cb0..0362e15dbe 100644
--- a/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C
+++ b/src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.C
@@ -138,6 +138,20 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> pow
     const GeometricField<scalar, PatchField, GeoMesh>& gsf2
 )
 {
+    if (!gsf1.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Base field is not dimensionless: " << gsf1.dimensions()
+            << exit(FatalError);
+    }
+
+    if (!gsf2.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Exponent field is not dimensionless: " << gsf2.dimensions()
+            << exit(FatalError);
+    }
+
     tmp<GeometricField<scalar, PatchField, GeoMesh>> tPow
     (
         new GeometricField<scalar, PatchField, GeoMesh>
@@ -151,11 +165,7 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> pow
                 IOobject::NO_WRITE
             ),
             gsf1.mesh(),
-            pow
-            (
-                gsf1.dimensions(),
-                dimensionedScalar("1", gsf2.dimensions(), 1.0)
-            )
+            dimless
         )
     );
 
@@ -174,17 +184,27 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> pow
 {
     const GeometricField<scalar, PatchField, GeoMesh>& gsf1 = tgsf1();
 
+    if (!gsf1.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Base field is not dimensionless: " << gsf1.dimensions()
+            << exit(FatalError);
+    }
+
+    if (!gsf2.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Exponent field is not dimensionless: " << gsf2.dimensions()
+            << exit(FatalError);
+    }
+
     tmp<GeometricField<scalar, PatchField, GeoMesh>> tPow
     (
         New
         (
             tgsf1,
             "pow(" + gsf1.name() + ',' + gsf2.name() + ')',
-            pow
-            (
-                gsf1.dimensions(),
-                dimensionedScalar("1", gsf2.dimensions(), 1.0)
-            )
+            dimless
         )
     );
 
@@ -205,17 +225,27 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> pow
 {
     const GeometricField<scalar, PatchField, GeoMesh>& gsf2 = tgsf2();
 
+    if (!gsf1.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Base field is not dimensionless: " << gsf1.dimensions()
+            << exit(FatalError);
+    }
+
+    if (!gsf2.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Exponent field is not dimensionless: " << gsf2.dimensions()
+            << exit(FatalError);
+    }
+
     tmp<GeometricField<scalar, PatchField, GeoMesh>> tPow
     (
         New
         (
             tgsf2,
             "pow(" + gsf1.name() + ',' + gsf2.name() + ')',
-            pow
-            (
-                gsf1.dimensions(),
-                dimensionedScalar("1", gsf2.dimensions(), 1.0)
-            )
+            dimless
         )
     );
 
@@ -236,6 +266,20 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> pow
     const GeometricField<scalar, PatchField, GeoMesh>& gsf1 = tgsf1();
     const GeometricField<scalar, PatchField, GeoMesh>& gsf2 = tgsf2();
 
+    if (!gsf1.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Base field is not dimensionless: " << gsf1.dimensions()
+            << exit(FatalError);
+    }
+
+    if (!gsf2.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Exponent field is not dimensionless: " << gsf2.dimensions()
+            << exit(FatalError);
+    }
+
     tmp<GeometricField<scalar, PatchField, GeoMesh>> tPow
     (
         reuseTmpTmpGeometricField
@@ -244,11 +288,7 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> pow
             tgsf1,
             tgsf2,
             "pow(" + gsf1.name() + ',' + gsf2.name() + ')',
-            pow
-            (
-                gsf1.dimensions(),
-                dimensionedScalar("1", gsf2.dimensions(), 1.0)
-            )
+            dimless
         )
     );
 
@@ -281,6 +321,13 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> pow
     const dimensionedScalar& ds
 )
 {
+    if (!ds.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Exponent is not dimensionless: " << ds.dimensions()
+            << exit(FatalError);
+    }
+
     tmp<GeometricField<scalar, PatchField, GeoMesh>> tPow
     (
         new GeometricField<scalar, PatchField, GeoMesh>
@@ -310,6 +357,13 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> pow
     const dimensionedScalar& ds
 )
 {
+    if (!ds.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Exponent is not dimensionless: " << ds.dimensions()
+            << exit(FatalError);
+    }
+
     const GeometricField<scalar, PatchField, GeoMesh>& gsf = tgsf();
 
     tmp<GeometricField<scalar, PatchField, GeoMesh>> tPow
@@ -370,6 +424,20 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> pow
     const GeometricField<scalar, PatchField, GeoMesh>& gsf
 )
 {
+    if (!ds.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Base scalar is not dimensionless: " << ds.dimensions()
+            << exit(FatalError);
+    }
+
+    if (!gsf.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Exponent field is not dimensionless: " << gsf.dimensions()
+            << exit(FatalError);
+    }
+
     tmp<GeometricField<scalar, PatchField, GeoMesh>> tPow
     (
         new GeometricField<scalar, PatchField, GeoMesh>
@@ -383,7 +451,7 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> pow
                 IOobject::NO_WRITE
             ),
             gsf.mesh(),
-            pow(ds, gsf.dimensions())
+            dimless
         )
     );
 
@@ -402,13 +470,27 @@ tmp<GeometricField<scalar, PatchField, GeoMesh>> pow
 {
     const GeometricField<scalar, PatchField, GeoMesh>& gsf = tgsf();
 
+    if (!ds.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Base scalar is not dimensionless: " << ds.dimensions()
+            << exit(FatalError);
+    }
+
+    if (!gsf.dimensions().dimensionless())
+    {
+        FatalErrorInFunction
+            << "Exponent field is not dimensionless: " << gsf.dimensions()
+            << exit(FatalError);
+    }
+
     tmp<GeometricField<scalar, PatchField, GeoMesh>> tPow
     (
         New
         (
             tgsf,
             "pow(" + ds.name() + ',' + gsf.name() + ')',
-            pow(ds, gsf.dimensions())
+            dimless
         )
     );