BUG: integratedNonUniformTable: correct offsets (fixes #2614)

- with geometryOrder=1, calculate the edge normals from the adjacent
  faces (area-weighted, inverse distance squared) and also
  use that for the Le() calculation.

  Includes the contributions from processor edge neighbours, so it
  should be consistent on both sides.

  This new method (consider as 'beta') contrasts with the current
  standard method that first calculates area-weighted point normals
  and uses the average of them for the edge normal.

  Enable for testing either with a controlDict OptimisationSwitch entry
  "fa:geometryOrder", or on the command-line:

      solverName -opt-switch=fa:geometryOrder=1

src/finiteArea/faMesh/faMesh.C

@@ -40,12 +40,24 @@ License
 #include "processorFaPatch.H"
 #include "wedgeFaPatch.H"
 #include "faPatchData.H"
+#include "registerSwitch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
     defineTypeNameAndDebug(faMesh, 0);
+
+    int faMesh::geometryOrder_
+    (
+        debug::optimisationSwitch("fa:geometryOrder", 2)
+    );
+    registerOptSwitch
+    (
+        "fa:geometryOrder",
+        int,
+        faMesh::geometryOrder_
+    );
 }
 
 
@@ -53,8 +65,6 @@ const Foam::word Foam::faMesh::prefix("finite-area");
 
 Foam::word Foam::faMesh::meshSubDir = "faMesh";
 
-int Foam::faMesh::geometryOrder_ = 1;      // 1: Standard treatment
-
 const int Foam::faMesh::quadricsFit_ = 0;  // Tuning (experimental)
 
 
@@ -228,7 +238,7 @@ void Foam::faMesh::clearGeomNotAreas() const
     deleteDemandDrivenData(patchStartsPtr_);
     deleteDemandDrivenData(LePtr_);
     deleteDemandDrivenData(magLePtr_);
-    deleteDemandDrivenData(centresPtr_);
+    deleteDemandDrivenData(faceCentresPtr_);
     deleteDemandDrivenData(edgeCentresPtr_);
     deleteDemandDrivenData(faceAreaNormalsPtr_);
     deleteDemandDrivenData(edgeAreaNormalsPtr_);
@@ -373,7 +383,7 @@ Foam::faMesh::faMesh
     patchStartsPtr_(nullptr),
     LePtr_(nullptr),
     magLePtr_(nullptr),
-    centresPtr_(nullptr),
+    faceCentresPtr_(nullptr),
     edgeCentresPtr_(nullptr),
     faceAreaNormalsPtr_(nullptr),
     edgeAreaNormalsPtr_(nullptr),
@@ -479,7 +489,7 @@ Foam::faMesh::faMesh
     patchStartsPtr_(nullptr),
     LePtr_(nullptr),
     magLePtr_(nullptr),
-    centresPtr_(nullptr),
+    faceCentresPtr_(nullptr),
     edgeCentresPtr_(nullptr),
     faceAreaNormalsPtr_(nullptr),
     edgeAreaNormalsPtr_(nullptr),
@@ -560,7 +570,7 @@ Foam::faMesh::faMesh
     patchStartsPtr_(nullptr),
     LePtr_(nullptr),
     magLePtr_(nullptr),
-    centresPtr_(nullptr),
+    faceCentresPtr_(nullptr),
     edgeCentresPtr_(nullptr),
     faceAreaNormalsPtr_(nullptr),
     edgeAreaNormalsPtr_(nullptr),
@@ -780,12 +790,12 @@ const Foam::edgeScalarField& Foam::faMesh::magLe() const
 
 const Foam::areaVectorField& Foam::faMesh::areaCentres() const
 {
-    if (!centresPtr_)
+    if (!faceCentresPtr_)
     {
-        calcAreaCentres();
+        calcFaceCentres();
     }
 
-    return *centresPtr_;
+    return *faceCentresPtr_;
 }
 
 
@@ -999,14 +1009,11 @@ bool Foam::faMesh::movePoints()
 
 bool Foam::faMesh::correctPatchPointNormals(const label patchID) const
 {
-    if ((patchID < 0) || (patchID >= boundary().size()))
+    if
-    {
+    (
-        FatalErrorInFunction
+        bool(correctPatchPointNormalsPtr_)
-            << "patchID is not in the valid range"
+     && patchID >= 0 && patchID < boundary().size()
-            << abort(FatalError);
+    )
-    }
-
-    if (correctPatchPointNormalsPtr_)
     {
         return (*correctPatchPointNormalsPtr_)[patchID];
     }

src/finiteArea/faMesh/faMesh.H

@@ -287,7 +287,7 @@ class faMesh
         mutable edgeScalarField* magLePtr_;
 
         //- Face centres
-        mutable areaVectorField* centresPtr_;
+        mutable areaVectorField* faceCentresPtr_;
 
         //- Edge centres
         mutable edgeVectorField* edgeCentresPtr_;
@@ -328,9 +328,6 @@ class faMesh
 
     // Static Private Data
 
-        //- Geometry treatment (0: primitive, 1: standard)
-        static int geometryOrder_;
-
         //- Quadrics fit for pointAreaNormals (experimental)
         static const int quadricsFit_;
 
@@ -378,6 +375,11 @@ class faMesh
             //- are related to the areaMesh
             void calcWhichPatchFaces() const;
 
+            //- Calculate the edge normals (direct calculation)
+            //- with flat boundary addressing
+            //  Possible communication
+            tmp<vectorField> calcRawEdgeNormals(int calcType) const;
+
             //- Calculate edge lengths
             //  Triggers communication via calcEdgeAreaNormals
             void calcLe() const;
@@ -388,7 +390,7 @@ class faMesh
 
             //- Calculate face centres
             //  No communication
-            void calcAreaCentres() const;
+            void calcFaceCentres() const;
 
             //- Calculate edge centres
             //  No communication
@@ -505,6 +507,12 @@ public:
         typedef faBoundaryMesh BoundaryMesh;
 
 
+    // Tuning switches
+
+        //- Geometry treatment
+        static int geometryOrder_;
+
+
     //- Runtime type information
     TypeName("faMesh");
 
@@ -568,8 +576,8 @@ public:
 
     // Static Functions
 
-        //- Return the current geometry treatment (0: primitive, 1: standard)
+        //- Return the current geometry treatment
-        //  A zero level is with restricted neighbour information
+        //  A zero level or negative is with restricted neighbour information
         static int geometryOrder() noexcept
         {
             return geometryOrder_;
@@ -577,7 +585,7 @@ public:
 
         //- Set the preferred geometry treatment
         //  \return the previous value
-        static int geometryOrder(const int order) noexcept
+        static int geometryOrder(int order) noexcept
         {
             int old(geometryOrder_);
             geometryOrder_ = order;
@@ -770,13 +778,13 @@ public:
             //- Face centres of boundary halo neighbours
             const pointField& haloFaceCentres() const;
 
-            //- Face normals of boundary halo neighbours
+            //- Face unit-normals of boundary halo neighbours
             const vectorField& haloFaceNormals() const;
 
             //- Face centres of boundary halo neighbours for specified patch
             tmp<pointField> haloFaceCentres(const label patchi) const;
 
-            //- Face normals of boundary halo neighbours for specified patch
+            //- Face unit-normals of boundary halo neighbours for specified patch
             tmp<vectorField> haloFaceNormals(const label patchi) const;
 
 
@@ -786,7 +794,7 @@ public:
         labelList faceCells() const;
 
 
-    // Storage management
+    // Storage Management
 
         //- Remove all files from mesh instance
         void removeFiles(const fileName& instanceDir) const;
@@ -795,6 +803,46 @@ public:
         void removeFiles() const;
 
 
+        //- Has face areas: S()
+        bool hasFaceAreas() const noexcept { return bool(SPtr_); }
+
+        //- Has face centres: areaCentres()
+        bool hasAreaCentres() const noexcept { return bool(faceCentresPtr_); }
+
+        //- Has edge centres: edgeCentres()
+        bool hasAdgeCentres() const noexcept { return bool(edgeCentresPtr_); }
+
+        //- Has edge length vectors: Le()
+        bool hasLe() const noexcept { return bool(LePtr_); }
+
+        //- Has edge length magnitudes: magLe()
+        bool hasMagLe() const noexcept { return bool(magLePtr_); }
+
+        //- Has face area normals: faceAreaNormals()
+        bool hasFaceAreaNormals() const noexcept
+        {
+            return bool(faceAreaNormalsPtr_);
+        }
+
+        //- Has edge area normals: edgeAreaNormals()
+        bool hasEdgeAreaNormals() const noexcept
+        {
+            return bool(edgeAreaNormalsPtr_);
+        }
+
+        //- Has point area normals: pointAreaNormals()
+        bool hasPointAreaNormals() const noexcept
+        {
+            return bool(pointAreaNormalsPtr_);
+        }
+
+        //- Has face curvatures: faceCurvatures()
+        bool hasFaceCurvatures() const noexcept
+        {
+            return bool(faceCurvaturesPtr_);
+        }
+
+
     // Mesh motion and morphing
 
         //- Is mesh moving

src/thermophysicalModels/thermophysicalProperties/thermophysicalFunctions/integratedNonUniformTable/integratedNonUniformTableThermophysicalFunction.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020 OpenFOAM Foundation
+    Copyright (C) 2020-2022 OpenFOAM Foundation
     Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
@@ -66,8 +66,8 @@ integratedNonUniformTable
 
     for (label i = 1; i < intf_.size(); ++i)
     {
-        intf_[i] = intf_[i-1] + intfdT(0, values()[i].first());
+        intf_[i] = intfdT(0, values()[i].first());
-        intfByT_[i] = intfByT_[i-1] + intfByTdT(0, values()[i].first());
+        intfByT_[i] = intfByTdT(0, values()[i].first());
     }
 
     const scalar intfStd = intfdT(Pstd, Tstd);