RELEASE: Updated version to v2206

CONTRIBUTORS.md

@@ -5,8 +5,6 @@ It is likely incomplete...
 
 ## Contributors (alphabetical by surname)
 
-- Tetsuo Aoyagi
-- Akira Azami
 - William Bainbridge
 - Gabriel Barajas
 - Kutalmis Bercin
@@ -21,7 +19,6 @@ It is likely incomplete...
 - Bernhard Gschaider
 - Andrew Heather
 - David Hill
-- Yoshiaki Inoue
 - Mattijs Janssens
 - Andrew Jackson
 - Hrvoje Jasak

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/DTRMParticle/DTRMParticle.H

@@ -34,8 +34,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_DTRMParticle_H
-#define Foam_DTRMParticle_H
+#ifndef DTRMParticle_H
+#define DTRMParticle_H
 
 #include "particle.H"
 #include "IOstream.H"
@@ -50,13 +50,12 @@ SourceFiles
 namespace Foam
 {
 
-// Forward Declarations
 class DTRMParticle;
-
-Ostream& operator<<(Ostream&, const DTRMParticle&);
-
 using namespace Foam::radiation;
 
+// Forward declaration of friend functions
+Ostream& operator<<(Ostream&, const DTRMParticle&);
+
 /*---------------------------------------------------------------------------*\
                         Class DTRMParticle Declaration
 \*---------------------------------------------------------------------------*/
@@ -131,7 +130,7 @@ public:
             );
 
 
-        // Member Functions
+        // Member functions
 
             inline const interpolationCell<scalar>& aInterp() const;
             inline const interpolationCell<scalar>& eInterp() const;
@@ -233,34 +232,37 @@ public:
     // Access
 
         //- Return const access to the initial position
-        const point& p0() const noexcept { return p0_; }
+        inline const point& p0() const;
 
         //- Return const access to the target position
-        const point& p1() const noexcept { return p1_; }
+        inline const point& p1() const;
 
         //- Return const access to the initial intensity
-        scalar I0() const noexcept { return I0_; }
+        inline scalar I0() const;
 
         //- Return const access to the current intensity
-        scalar I() const noexcept { return I_; }
+        inline scalar I() const;
 
         //- Return const access dA
-        scalar dA() const noexcept { return dA_; }
+        inline scalar dA() const;
 
 
     // Edit
 
         //- Return access to the target position
-        point& p1() noexcept { return p1_; }
+        inline point& p1();
 
         //- Return access to the initial intensity
-        scalar& I0() noexcept { return I0_; }
+        inline scalar& I0();
 
         //- Return access to the current intensity
-        scalar& I() noexcept { return I_; }
+        inline scalar& I();
 
         //- Return access to dA
-        scalar& dA() noexcept { return dA_; }
+        inline scalar& dA();
+
+        //- Return access to reflectedId
+        inline label& reflectedId();
 
 
    // Tracking

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/laserDTRM/DTRMParticle/DTRMParticleI.H

@@ -107,4 +107,58 @@ inline Foam::scalar& Foam::DTRMParticle::trackingData::Q(label celli)
 }
 
 
+inline const Foam::point& Foam::DTRMParticle::p0() const
+{
+    return p0_;
+}
+
+
+inline const Foam::point& Foam::DTRMParticle::p1() const
+{
+    return p1_;
+}
+
+
+inline Foam::scalar Foam::DTRMParticle::I0() const
+{
+    return I0_;
+}
+
+
+inline Foam::scalar Foam::DTRMParticle::I() const
+{
+    return I_;
+}
+
+
+inline Foam::scalar Foam::DTRMParticle::dA() const
+{
+    return dA_;
+}
+
+
+inline Foam::scalar& Foam::DTRMParticle::dA()
+{
+    return dA_;
+}
+
+
+inline Foam::point& Foam::DTRMParticle::p1()
+{
+    return p1_;
+}
+
+
+inline Foam::scalar& Foam::DTRMParticle::I0()
+{
+    return I0_;
+}
+
+
+inline Foam::scalar& Foam::DTRMParticle::I()
+{
+    return I_;
+}
+
+
 // ************************************************************************* //

applications/test/barycentric/Make/files

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

applications/test/barycentric/Make/options

@@ -1,2 +0,0 @@
-/* EXE_INC = */
-/* EXE_LIBS = */

applications/test/barycentric/Test-barycentric.C

@@ -1,93 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2022 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-barycentric
-
-Description
-    Some simple tests for barycentric coordinates and transforms
-
-\*---------------------------------------------------------------------------*/
-
-#include "barycentricTensor.H"
-#include "tetrahedron.H"
-#include "vectorField.H"
-#include "IOstreams.H"
-
-using namespace Foam;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-//  Main program:
-
-int main(int argc, char *argv[])
-{
-    // Tets to test
-    tetPoints tetA
-    (
-        point(0, 0, 0),
-        point(1, 0, 0),
-        point(1, 1, 0),
-        point(1, 1, 1)
-    );
-
-    const barycentricTensor baryT(tetA[0], tetA[1], tetA[2], tetA[3]);
-
-    Info<< nl << "Tet: " << tetA << nl;
-    Info<< "tens:" << baryT << nl;
-
-    for
-    (
-        const barycentric& bary :
-        List<barycentric>
-        ({
-            {0.25, 0.25, 0.25, 0.25},
-            {1, 0, 0, 0},
-            {0, 1, 0, 0},
-            {0, 0, 1, 0},
-            {0, 0, 0, 1},
-            {0, 0, 0, 0}  // Not really valid
-        })
-    )
-    {
-        vector v(tetA.tet().barycentricToPoint(bary));
-        barycentric b(tetA.tet().pointToBarycentric(v));
-
-        Info<< nl
-            << "bary: " << bary << nl
-            << "vec:  " << v << nl
-            // << "Vec:  " << baryT.inner(bary) << nl
-            << "Vec:  " << (baryT & bary) << nl
-            << "bary: " << b << nl
-            // This won't work (needs a differently defined tensor)
-            // << "Bary: " << (v & baryT) << nl
-            ;
-    }
-
-    Info<< "\nEnd\n" << nl;
-
-    return 0;
-}
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/test/momentOfInertia/Test-momentOfInertia.C

@@ -38,7 +38,7 @@ Description
 #include "polyMesh.H"
 #include "ListOps.H"
 #include "face.H"
-#include "tetrahedron.H"
+#include "tetPointRef.H"
 #include "triFaceList.H"
 #include "OFstream.H"
 #include "meshTools.H"

applications/test/tetTetOverlap/Test-tetTetOverlap.C

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2017 OpenFOAM Foundation
-    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,7 +31,7 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "tetrahedron.H"
+#include "tetPointRef.H"
 #include "OFstream.H"
 #include "meshTools.H"
 #include "cut.H"
@@ -45,12 +44,12 @@ void writeOBJ
 (
     Ostream& os,
     label& vertI,
-    const tetPoints& tet
+    const FixedList<point, 4>& tet
 )
 {
-    for (const point& p : tet)
+    forAll(tet, fp)
     {
-        meshTools::writeOBJ(os, p);
+        meshTools::writeOBJ(os, tet[fp]);
     }
     os  << "l " << vertI+1 << ' ' << vertI+2 << nl
         << "l " << vertI+1 << ' ' << vertI+3 << nl
@@ -62,27 +61,33 @@ void writeOBJ
 }
 
 
+tetPointRef makeTetPointRef(const FixedList<point, 4>& p)
+{
+    return tetPointRef(p[0], p[1], p[2], p[3]);
+}
+
+
 int main(int argc, char *argv[])
 {
     // Tets to test
-    tetPoints tetA
-    (
+    FixedList<point, 4> tetA
+    ({
         point(0, 0, 0),
         point(1, 0, 0),
         point(1, 1, 0),
         point(1, 1, 1)
-    );
-    tetPoints tetB
-    (
+    });
+    FixedList<point, 4> tetB
+    ({
         point(0.1, 0.1, 0.1),
         point(1.1, 0.1, 0.1),
         point(1.1, 1.1, 0.1),
         point(1.1, 1.1, 1.1)
-    );
+    });
 
 
     // Do intersection
-    typedef DynamicList<tetPoints> tetList;
+    typedef DynamicList<FixedList<point, 4>> tetList;
     tetList tetsIn1, tetsIn2, tetsOut;
     cut::appendOp<tetList> tetOpIn1(tetsIn1);
     cut::appendOp<tetList> tetOpIn2(tetsIn2);
@@ -150,25 +155,25 @@ int main(int argc, char *argv[])
 
 
     // Check the volumes
-    Info<< "Vol A: " << tetA.tet().mag() << endl;
+    Info<< "Vol A: " << makeTetPointRef(tetA).mag() << endl;
 
     scalar volIn = 0;
-    for (const auto& t : tetsIn)
+    forAll(tetsIn, i)
     {
-        volIn += t.tet().mag();
+        volIn += makeTetPointRef(tetsIn[i]).mag();
     }
     Info<< "Vol A inside B: " << volIn << endl;
 
     scalar volOut = 0;
-    for (const auto& t : tetsOut)
+    forAll(tetsOut, i)
     {
-        volOut += t.tet().mag();
+        volOut += makeTetPointRef(tetsOut[i]).mag();
     }
     Info<< "Vol A outside B: " << volOut << endl;
 
     Info<< "Sum inside and outside: " << volIn + volOut << endl;
 
-    if (mag(volIn + volOut - tetA.tet().mag()) > SMALL)
+    if (mag(volIn + volOut - makeTetPointRef(tetA).mag()) > SMALL)
     {
         FatalErrorInFunction
             << "Tet volumes do not sum up to input tet."

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellShapeControlMesh/cellShapeControlMesh.C

@@ -31,7 +31,7 @@ License
 #include "pointIOField.H"
 #include "scalarIOField.H"
 #include "triadIOField.H"
-#include "tetrahedron.H"
+#include "tetPointRef.H"
 #include "plane.H"
 #include "transform.H"
 #include "meshTools.H"

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/fileControl/fileControl.C

@@ -27,6 +27,7 @@ License
 
 #include "fileControl.H"
 #include "addToRunTimeSelectionTable.H"
+#include "tetPointRef.H"
 #include "scalarList.H"
 #include "vectorTools.H"
 #include "pointIOField.H"
@@ -49,6 +50,9 @@ addToRunTimeSelectionTable
 }
 
 
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::fileControl::fileControl

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/cellShapeControl/cellSizeAndAlignmentControl/searchableSurfaceControl/searchableSurfaceControl.C

@@ -31,6 +31,7 @@ License
 #include "cellSizeFunction.H"
 #include "triSurfaceMesh.H"
 #include "searchableBox.H"
+#include "tetPointRef.H"
 #include "vectorTools.H"
 #include "quaternion.H"
 

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellChecks.C

@@ -26,7 +26,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "plane.H"
-#include "tetrahedron.H"
+#include "tetPointRef.H"
 #include "pointConversion.H"
 #include "CGALTriangulation3DKernel.H"
 

applications/utilities/mesh/manipulation/checkMesh/writeFields.C

@@ -4,7 +4,7 @@
 #include "polyMeshTools.H"
 #include "zeroGradientFvPatchFields.H"
 #include "syncTools.H"
-#include "tetrahedron.H"
+#include "tetPointRef.H"
 #include "regionSplit.H"
 #include "wallDist.H"
 #include "cellAspectRatio.H"

applications/utilities/mesh/manipulation/createPatch/createPatch.C

@@ -606,7 +606,7 @@ void syncPoints
                 pointField nbrPatchInfo(procPatch.nPoints());
                 {
                     // We do not know the number of points on the other side
-                    // so cannot use UIPstream::read
+                    // so cannot use Pstream::read.
                     IPstream fromNbr
                     (
                         Pstream::commsTypes::blocking,

applications/utilities/mesh/manipulation/moveMesh/moveMesh.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,30 +55,17 @@ int main(int argc, char *argv[])
     (
         "deltaT",
         "time",
-        "Override deltaT (eg, for accelerated motion)"
-    );
-
-    argList::addOption
-    (
-        "endTime",
-        "time",
-        "Override endTime (eg, for shorter tests)"
+        "Override deltaT for accelerated motion"
     );
 
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createNamedMesh.H"
 
-    scalar timeVal = 0;
-    if (args.readIfPresent("deltaT", timeVal))
+    scalar deltaT = 0;
+    if (args.readIfPresent("deltaT", deltaT))
     {
-        runTime.setDeltaT(timeVal);
-    }
-
-    if (args.readIfPresent("endTime", timeVal))
-    {
-        runTime.stopAt(Time::stopAtControls::saEndTime);
-        runTime.setEndTime(timeVal);
+        runTime.setDeltaT(deltaT);
     }
 
     autoPtr<motionSolver> motionPtr = motionSolver::New(mesh);

applications/utilities/parallelProcessing/redistributePar/passivePositionParticle.H

@@ -35,8 +35,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_passivePositionParticle_H
-#define Foam_passivePositionParticle_H
+#ifndef passivePositionParticle_H
+#define passivePositionParticle_H
 
 #include "passiveParticle.H"
 
@@ -130,9 +130,7 @@ public:
         };
 
 
-    // Member Functions
-
-        const point& cachedPosition() const noexcept
+        const point& cachedPosition() const
         {
             return cachedPosition_;
         }

etc/bashrc

@@ -55,7 +55,7 @@
 
 # [WM_PROJECT_VERSION] - A human-readable version name
 # A development version is often named 'com' - as in www.openfoam.com
-export WM_PROJECT_VERSION=com
+export WM_PROJECT_VERSION=v2206
 
 #------------------------------------------------------------------------------
 # Configuration environment variables.

etc/caseDicts/postProcessing/visualization/streamlines.cfg

@@ -13,7 +13,7 @@ executeControl  writeTime;
 writeControl    writeTime;
 
 setFormat       vtk;
-direction       forward;   // (forward | backward | bidirectional)
+trackForward    true;
 
 lifeTime        10000;
 nSubCycle       5;

etc/cshrc

@@ -55,7 +55,7 @@
 
 # [WM_PROJECT_VERSION] - A human-readable version name
 # A development version is often named 'com' - as in www.openfoam.com
-setenv WM_PROJECT_VERSION com
+setenv WM_PROJECT_VERSION v2206
 
 #------------------------------------------------------------------------------
 # Configuration environment variables.

src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.C

@@ -27,7 +27,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "dynamicIndexedOctree.H"
-#include "line.H"
+#include "linePointRef.H"
 #include "OFstream.H"
 #include "ListOps.H"
 

src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H

@@ -35,8 +35,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_dynamicIndexedOctree_H
-#define Foam_dynamicIndexedOctree_H
+#ifndef dynamicIndexedOctree_H
+#define dynamicIndexedOctree_H
 
 #include "treeBoundBox.H"
 #include "pointIndexHit.H"
@@ -54,7 +54,7 @@ namespace Foam
 
 typedef DynamicList<autoPtr<DynamicList<label>>> contentListList;
 
-// Forward Declarations
+// Forward declaration of classes
 template<class Type> class dynamicIndexedOctree;
 
 template<class Type> Ostream& operator<<

src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicTreeDataPoint.H

@@ -39,12 +39,12 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_dynamicTreeDataPoint_H
-#define Foam_dynamicTreeDataPoint_H
+#ifndef dynamicTreeDataPoint_H
+#define dynamicTreeDataPoint_H
 
 #include "pointField.H"
 #include "treeBoundBox.H"
-#include "line.H"
+#include "linePointRef.H"
 #include "volumeType.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -52,7 +52,7 @@ SourceFiles
 namespace Foam
 {
 
-// Forward Declarations
+// Forward declaration of classes
 template<class Type> class dynamicIndexedOctree;
 
 /*---------------------------------------------------------------------------*\

src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C

@@ -27,7 +27,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "indexedOctree.H"
-#include "line.H"
+#include "linePointRef.H"
 #include "OFstream.H"
 #include "ListOps.H"
 #include "memInfo.H"

src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H

@@ -34,8 +34,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_indexedOctree_H
-#define Foam_indexedOctree_H
+#ifndef indexedOctree_H
+#define indexedOctree_H
 
 #include "treeBoundBox.H"
 #include "pointIndexHit.H"
@@ -51,7 +51,7 @@ SourceFiles
 namespace Foam
 {
 
-// Forward Declarations
+// Forward declaration of classes
 template<class Type> class indexedOctree;
 template<class Type> Ostream& operator<<(Ostream&, const indexedOctree<Type>&);
 class Istream;

src/OpenFOAM/db/IOobjects/decomposedBlockData/decomposedBlockData.C

@@ -789,7 +789,7 @@ bool Foam::decomposedBlockData::writeBlocks
             for (label proci = 1; proci < nProcs; ++proci)
             {
                 elems.resize(recvSizes[proci]);
-                UIPstream::read
+                IPstream::read
                 (
                     UPstream::commsTypes::scheduled,
                     proci,

src/OpenFOAM/db/Time/Time.C

@@ -967,34 +967,6 @@ bool Foam::Time::end() const
 }
 
 
-bool Foam::Time::reverseRun() const
-{
-    // Must not solve for the 0 time step
-    bool isRunning = value() > (endTime_ + 0.5*deltaT_);
-
-    return isRunning;
-}
-
-
-bool Foam::Time::reverseLoop()
-{
-    const bool isRunning = reverseRun();
-
-    if (isRunning)
-    {
-        operator--();
-    }
-
-    return isRunning;
-}
-
-
-bool Foam::Time::reverseEnd() const
-{
-    return value() < (endTime_ - 0.5*deltaT_);
-}
-
-
 bool Foam::Time::stopAt(const stopAtControls stopCtrl) const
 {
     if (stopCtrl == stopAtControls::saUnknown)
@@ -1393,242 +1365,4 @@ Foam::Time& Foam::Time::operator++(int)
 }
 
 
-Foam::Time&
-Foam::Time::operator-=(const dimensionedScalar& deltaT)
-{
-    return operator-=(deltaT.value());
-}
-
-
-Foam::Time& Foam::Time::operator-=(const scalar deltaT)
-{
-    setDeltaT(deltaT);
-    return operator--();
-}
-
-
-Foam::Time& Foam::Time::operator--()
-{
-    deltaT0_ = deltaTSave_;
-    deltaTSave_ = deltaT_;
-
-    // Save old time value and name
-    const scalar oldTimeValue = timeToUserTime(value());
-    const word oldTimeName = dimensionedScalar::name();
-
-    // Decrease time
-    setTime(value() - deltaT_, timeIndex_ - 1);
-
-    if (!subCycling_)
-    {
-        // If the time is very close to zero reset to zero
-        if (mag(value()) < 10*SMALL*deltaT_)
-        {
-            setTime(0.0, timeIndex_);
-        }
-
-        if (sigStopAtWriteNow_.active() || sigWriteNow_.active())
-        {
-            // A signal might have been sent on one processor only
-            // Reduce so all decide the same.
-
-            label flag = 0;
-            if (sigStopAtWriteNow_.active() && stopAt_ == saWriteNow)
-            {
-                flag += 1;
-            }
-            if (sigWriteNow_.active() && writeOnce_)
-            {
-                flag += 2;
-            }
-            reduce(flag, maxOp<label>());
-
-            if (flag & 1)
-            {
-                stopAt_ = saWriteNow;
-            }
-            if (flag & 2)
-            {
-                writeOnce_ = true;
-            }
-        }
-
-        writeTime_ = false;
-
-        switch (writeControl_)
-        {
-            case wcNone:
-            case wcUnknown:
-            break;
-
-            case wcTimeStep:
-                // Avoid writing just because timIndex_ is zero
-                writeTime_ = 
-                    !(timeIndex_ % label(writeInterval_)) && timeIndex_;
-            break;
-
-            case wcRunTime:
-            case wcAdjustableRunTime:
-            {
-                const label writeIndex = label
-                (
-                    ((value() - startTime_) - 0.5*deltaT_)
-                  / writeInterval_
-                );
-
-                if (writeIndex < writeTimeIndex_)
-                {
-                    writeTime_ = true;
-                    writeTimeIndex_ = writeIndex;
-                }
-            }
-            break;
-
-            case wcCpuTime:
-            {
-                const label writeIndex = label
-                (
-                    returnReduce(elapsedCpuTime(), maxOp<double>())
-                  / writeInterval_
-                );
-                if (writeIndex > writeTimeIndex_)
-                {
-                    writeTime_ = true;
-                    writeTimeIndex_ = writeIndex;
-                }
-            }
-            break;
-
-            case wcClockTime:
-            {
-                const label writeIndex = label
-                (
-                    returnReduce(elapsedClockTime(), maxOp<double>())
-                  / writeInterval_
-                );
-                if (writeIndex > writeTimeIndex_)
-                {
-                    writeTime_ = true;
-                    writeTimeIndex_ = writeIndex;
-                }
-            }
-            break;
-        }
-
-
-        // Check if endTime needs adjustment to stop at the next
-        // reverseRun()/reverseEnd()
-        if (!reverseEnd())
-        {
-            if (stopAt_ == saNoWriteNow)
-            {
-                endTime_ = value();
-            }
-            else if (stopAt_ == saWriteNow)
-            {
-                endTime_ = value();
-                writeTime_ = true;
-            }
-            else if (stopAt_ == saNextWrite && writeTime_ == true)
-            {
-                endTime_ = value();
-            }
-        }
-
-        // Override writeTime if one-shot writing
-        if (writeOnce_)
-        {
-            writeTime_ = true;
-            writeOnce_ = false;
-        }
-
-        // Adjust the precision of the time directory name if necessary
-        if (writeTime_)
-        {
-            // Tolerance used when testing time equivalence
-            const scalar timeTol =
-                max(min(pow(10.0, -precision_), 0.1*deltaT_), SMALL);
-
-            // User-time equivalent of deltaT
-            const scalar userDeltaT = timeToUserTime(deltaT_);
-
-            // Time value obtained by reading timeName
-            scalar timeNameValue = -VGREAT;
-
-            // Check that new time representation differs from old one
-            // reinterpretation of the word
-            if
-            (
-                readScalar(dimensionedScalar::name(), timeNameValue)
-             && (mag(-timeNameValue + oldTimeValue - userDeltaT) > timeTol)
-            )
-            {
-                int oldPrecision = precision_;
-                while
-                (
-                    precision_ < maxPrecision_
-                 && readScalar(dimensionedScalar::name(), timeNameValue)
-                 && (mag(-timeNameValue + oldTimeValue - userDeltaT) > timeTol)
-                )
-                {
-                    precision_++;
-                    setTime(value(), timeIndex());
-                }
-
-                if (precision_ != oldPrecision)
-                {
-                    WarningInFunction
-                        << "Increased the timePrecision from " << oldPrecision
-                        << " to " << precision_
-                        << " to distinguish between timeNames at time "
-                        << dimensionedScalar::name()
-                        << endl;
-
-                    if (precision_ == maxPrecision_)
-                    {
-                        // Reached maxPrecision limit
-                        WarningInFunction
-                            << "Current time name " << dimensionedScalar::name()
-                            << nl
-                            << "    The maximum time precision has been reached"
-                               " which might result in overwriting previous"
-                               " results."
-                            << endl;
-                    }
-
-                    // Check if round-off error caused time-reversal
-                    scalar oldTimeNameValue = -VGREAT;
-                    if
-                    (
-                        readScalar(oldTimeName, oldTimeNameValue)
-                     && (
-                            sign(timeNameValue - oldTimeNameValue)
-                         != sign(deltaT_)
-                        )
-                    )
-                    {
-                        WarningInFunction
-                            << "Current time name " << dimensionedScalar::name()
-                            << " is set to an instance prior to the "
-                               "previous one "
-                            << oldTimeName << nl
-                            << "    This might result in temporal "
-                               "discontinuities."
-                            << endl;
-                    }
-                }
-            }
-        }
-    }
-
-    return *this;
-}
-
-
-Foam::Time& Foam::Time::operator--(int)
-{
-    return operator--();
-}
-
-
 // ************************************************************************* //

src/OpenFOAM/db/Time/Time.H

@@ -7,7 +7,6 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
     Copyright (C) 2016-2020 OpenCFD Ltd.
-    Copyright (C) 2022      PCOpt/NTUA
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -572,43 +571,6 @@ public:
             //      not yield the same result
             virtual bool end() const;
 
-            //- Return true if run should continue.
-            //  Does not invoke any functionObject methods
-            //  \note
-            //  For correct behaviour, the following style of time-loop
-            //  is recommended:
-            //  \code
-            //      while (runTime.reverseRun())
-            //      {
-            //          --runTime;
-            //          solve;
-            //          runTime.write();
-            //      }
-            //  \endcode
-            virtual bool reverseRun() const;
-
-            //- Return true if run should continue and if so decrease time
-            //  Does not invoke any functionObject methods
-            //  \note
-            //  For correct behaviour, the following style of time-loop
-            //  is recommended:
-            //  \code
-            //      while (runTime.reverseLoop())
-            //      {
-            //          solve;
-            //          runTime.write();
-            //      }
-            //  \endcode
-            virtual bool reverseLoop();
-
-            //- Return true if end of run,
-            //  does not invoke any functionObject methods
-            //  \note
-            //      The rounding heuristics near endTime mean that
-            //      \code reverseRun() \endcode and \code !reverseEnd() \endcode
-            //      may not yield the same result
-            virtual bool reverseEnd() const;
-
 
         // Edit
 
@@ -691,20 +653,6 @@ public:
 
         //- Postfix increment, this is identical to the prefix increment
         virtual Time& operator++(int);
-
-        //- Set deltaT to that specified and decrease time via operator--()
-        virtual Time& operator-=(const dimensionedScalar& deltaT);
-
-        //- Set deltaT to that specified and decrease time via operator--()
-        virtual Time& operator-=(const scalar deltaT);
-
-        //- prefix decrease,
-        //  also invokes the functionobjectlist::start() or
-        //  functionobjectlist::execute() method, depending on the time-index
-        virtual Time& operator--();
-
-        //- Postfix decrease, this is identical to the decrease increment
-        virtual Time& operator--(int);
 };
 
 

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C

@@ -545,8 +545,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 (
     const IOobject& io,
     const Mesh& mesh,
-    const dictionary& dict,
-    const bool readOldTime
+    const dictionary& dict
 )
 :
     Internal(io, mesh, dimless, false),
@@ -567,11 +566,6 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
             << exit(FatalIOError);
     }
 
-    if (readOldTime)
-    {
-        readOldTimeIfPresent();
-    }
-
     DebugInFunction
         << "Finishing dictionary-construct" << nl << this->info() << endl;
 }
@@ -629,8 +623,7 @@ template<class Type, template<class> class PatchField, class GeoMesh>
 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 (
     const IOobject& io,
-    const GeometricField<Type, PatchField, GeoMesh>& gf,
-    const bool readOldTime
+    const GeometricField<Type, PatchField, GeoMesh>& gf
 )
 :
     Internal(io, gf),
@@ -643,7 +636,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
         << "Copy construct, resetting IO params" << nl
         << this->info() << endl;
 
-    if (!readIfPresent() && gf.field0Ptr_ && readOldTime)
+    if (!readIfPresent() && gf.field0Ptr_)
     {
         field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
         (

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H

@@ -296,8 +296,7 @@ public:
         (
             const IOobject& io,
             const Mesh& mesh,
-            const dictionary& dict,
-            const bool readOldTime = false
+            const dictionary& dict
         );
 
         //- Copy construct
@@ -316,8 +315,7 @@ public:
         GeometricField
         (
             const IOobject& io,
-            const GeometricField<Type, PatchField, GeoMesh>& gf,
-            const bool readOldTime = true
+            const GeometricField<Type, PatchField, GeoMesh>& gf
         );
 
         //- Construct from tmp\<GeometricField\> resetting IO parameters

src/OpenFOAM/fields/pointPatchFields/constraint/processorCyclic/processorCyclicPointPatchField.C

@@ -119,7 +119,7 @@ void Foam::processorCyclicPointPatchField<Type>::initSwapAddSeparated
         if (commsType == Pstream::commsTypes::nonBlocking)
         {
             receiveBuf_.setSize(pf.size());
-            UIPstream::read
+            IPstream::read
             (
                 commsType,
                 procPatch_.neighbProcNo(),
@@ -129,7 +129,7 @@ void Foam::processorCyclicPointPatchField<Type>::initSwapAddSeparated
                 procPatch_.comm()
             );
         }
-        UOPstream::write
+        OPstream::write
         (
             commsType,
             procPatch_.neighbProcNo(),
@@ -155,7 +155,7 @@ void Foam::processorCyclicPointPatchField<Type>::swapAddSeparated
         if (commsType != Pstream::commsTypes::nonBlocking)
         {
             receiveBuf_.setSize(this->size());
-            UIPstream::read
+            IPstream::read
             (
                 commsType,
                 procPatch_.neighbProcNo(),

src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C

@@ -86,7 +86,7 @@ Foam::LUscalarMatrix::LUscalarMatrix
 
         if (Pstream::master(comm_))
         {
-            for (const int proci : Pstream::subProcs(comm_))
+            for (const int slave : Pstream::subProcs(comm_))
             {
                 lduMatrices.set
                 (
@@ -96,7 +96,7 @@ Foam::LUscalarMatrix::LUscalarMatrix
                         IPstream
                         (
                             Pstream::commsTypes::scheduled,
-                            proci,
+                            slave,
                             0,          // bufSize
                             Pstream::msgType(),
                             comm_

src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrixTemplates.C

@@ -54,17 +54,17 @@ void Foam::LUscalarMatrix::solve
 
             SubList<Type>(X, x.size()) = x;
 
-            for (const int proci : Pstream::subProcs(comm_))
+            for (const int slave : Pstream::subProcs(comm_))
             {
-                UIPstream::read
+                IPstream::read
                 (
                     Pstream::commsTypes::scheduled,
-                    proci,
+                    slave,
                     reinterpret_cast<char*>
                     (
-                        &(X[procOffsets_[proci]])
+                        &(X[procOffsets_[slave]])
                     ),
-                    (procOffsets_[proci+1]-procOffsets_[proci])*sizeof(Type),
+                    (procOffsets_[slave+1]-procOffsets_[slave])*sizeof(Type),
                     Pstream::msgType(),
                     comm_
                 );
@@ -72,7 +72,7 @@ void Foam::LUscalarMatrix::solve
         }
         else
         {
-            UOPstream::write
+            OPstream::write
             (
                 Pstream::commsTypes::scheduled,
                 Pstream::masterNo(),
@@ -89,17 +89,17 @@ void Foam::LUscalarMatrix::solve
 
             x = SubList<Type>(X, x.size());
 
-            for (const int proci : Pstream::subProcs(comm_))
+            for (const int slave : Pstream::subProcs(comm_))
             {
-                UOPstream::write
+                OPstream::write
                 (
                     Pstream::commsTypes::scheduled,
-                    proci,
+                    slave,
                     reinterpret_cast<const char*>
                     (
-                        &(X[procOffsets_[proci]])
+                        &(X[procOffsets_[slave]])
                     ),
-                    (procOffsets_[proci+1]-procOffsets_[proci])*sizeof(Type),
+                    (procOffsets_[slave+1]-procOffsets_[slave])*sizeof(Type),
                     Pstream::msgType(),
                     comm_
                 );
@@ -107,7 +107,7 @@ void Foam::LUscalarMatrix::solve
         }
         else
         {
-            UIPstream::read
+            IPstream::read
             (
                 Pstream::commsTypes::scheduled,
                 Pstream::masterNo(),

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrix.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,7 +54,6 @@ SourceFiles
 #define Foam_LduMatrix_H
 
 #include "lduMesh.H"
-#include "lduMatrix.H"
 #include "Field.H"
 #include "FieldField.H"
 #include "LduInterfaceFieldPtrsList.H"
@@ -70,7 +69,8 @@ namespace Foam
 
 // Forward Declarations
 
-template<class Type, class DType, class LUType> class LduMatrix;
+template<class Type, class DType, class LUType>
+class LduMatrix;
 
 template<class Type, class DType, class LUType>
 Ostream& operator<<
@@ -119,13 +119,16 @@ public:
 
         // Protected Data
 
+            //- Default maximum number of iterations in the solver
+            static const label defaultMaxIter_ = 1000;
+
             word fieldName_;
             const LduMatrix<Type, DType, LUType>& matrix_;
 
-            //- Dictionary of solution controls
+            //- Dictionary of controls
             dictionary controlDict_;
 
-            //- Verbosity level for solver output statements
+            //- Level of verbosity in the solver output statements
             int log_;
 
             //- Minimum number of iterations in the solver
@@ -134,9 +137,6 @@ public:
             //- Maximum number of iterations in the solver
             label maxIter_;
 
-            //- The matrix normalisation type
-            lduMatrix::normTypes normType_;
-
             //- Final convergence tolerance
             Type tolerance_;
 
@@ -146,7 +146,7 @@ public:
 
         // Protected Member Functions
 
-            //- Read the control parameters from controlDict_
+            //- Read the control parameters from the controlDict_
             virtual void readControls();
 
 
@@ -206,7 +206,6 @@ public:
 
         // Constructors
 
-            //- Construct for given field name, matrix and controls
             solver
             (
                 const word& fieldName,
@@ -226,8 +225,9 @@ public:
             );
 
 
-        //- Destructor
-        virtual ~solver() = default;
+        // Destructor
+
+            virtual ~solver() = default;
 
 
         // Member Functions
@@ -244,33 +244,21 @@ public:
 
 
             //- Read and reset the solver parameters from the given dictionary
-            virtual void read(const dictionary&);
+            virtual void read(const dictionary& solverDict);
 
             virtual SolverPerformance<Type> solve
             (
                 Field<Type>& psi
             ) const = 0;
 
-            //- Return the matrix norm using the specified norm method
-            Type normFactor
-            (
-                const Field<Type>& psi,
-                const Field<Type>& Apsi,
-                Field<Type>& tmpField,
-                const lduMatrix::normTypes normType
-            ) const;
-
             //- Return the matrix norm used to normalise the residual for the
-            //- stopping criterion
+            //  stopping criterion
             Type normFactor
             (
                 const Field<Type>& psi,
                 const Field<Type>& Apsi,
                 Field<Type>& tmpField
-            ) const
-            {
-                return this->normFactor(psi, Apsi, tmpField, normType_);
-            }
+            ) const;
     };
 
 
@@ -326,7 +314,6 @@ public:
 
         // Constructors
 
-            //- Construct for given field name and matrix
             smoother
             (
                 const word& fieldName,
@@ -345,8 +332,9 @@ public:
             );
 
 
-        //- Destructor
-        virtual ~smoother() = default;
+        // Destructor
+
+            virtual ~smoother() = default;
 
 
         // Member Functions
@@ -417,8 +405,10 @@ public:
 
         // Constructors
 
-            //- Construct for given solver
-            preconditioner(const solver& sol)
+            preconditioner
+            (
+                const solver& sol
+            )
             :
                 solver_(sol)
             {}
@@ -434,15 +424,16 @@ public:
             );
 
 
-        //- Destructor
-        virtual ~preconditioner() = default;
+        // Destructor
+
+            virtual ~preconditioner() = default;
 
 
         // Member functions
 
             //- Read and reset the preconditioner parameters
-            //- from the given dictionary
-            virtual void read(const dictionary&)
+            //  from the given dictionary
+            virtual void read(const dictionary& preconditionerDict)
             {}
 
             //- Return wA the preconditioned form of residual rA
@@ -453,7 +444,7 @@ public:
             ) const = 0;
 
             //- Return wT the transpose-matrix preconditioned form of
-            //- residual rT.
+            //  residual rT.
             //  This is only required for preconditioning asymmetric matrices.
             virtual void preconditionT
             (
@@ -489,16 +480,17 @@ public:
         LduMatrix(const lduMesh&, Istream&);
 
 
-    //- Destructor
-    ~LduMatrix();
+    // Destructor
+
+        ~LduMatrix();
 
 
-    // Member Functions
+    // Member functions
 
         // Access to addressing
 
             //- Return the LDU mesh from which the addressing is obtained
-            const lduMesh& mesh() const noexcept
+            const lduMesh& mesh() const
             {
                 return lduMesh_;
             }
@@ -562,43 +554,43 @@ public:
             }
 
 
-            bool hasDiag() const noexcept
+            bool hasDiag() const
             {
                 return (diagPtr_);
             }
 
-            bool hasUpper() const noexcept
+            bool hasUpper() const
             {
                 return (upperPtr_);
             }
 
-            bool hasLower() const noexcept
+            bool hasLower() const
             {
                 return (lowerPtr_);
             }
 
-            bool hasSource() const noexcept
+            bool hasSource() const
             {
                 return (sourcePtr_);
             }
 
-            bool diagonal() const noexcept
+            bool diagonal() const
             {
                 return (diagPtr_ && !lowerPtr_ && !upperPtr_);
             }
 
-            bool symmetric() const noexcept
+            bool symmetric() const
             {
                 return (diagPtr_ && (!lowerPtr_ && upperPtr_));
             }
 
-            bool asymmetric() const noexcept
+            bool asymmetric() const
             {
                 return (diagPtr_ && lowerPtr_ && upperPtr_);
             }
 
 
-        // Operations
+        // operations
 
             void sumDiag();
             void negSumDiag();
@@ -648,7 +640,7 @@ public:
             tmp<Field<Type>> faceH(const tmp<Field<Type>>&) const;
 
 
-    // Member Operators
+    // Member operators
 
         void operator=(const LduMatrix<Type, DType, LUType>&);
 
@@ -661,7 +653,7 @@ public:
         void operator*=(scalar);
 
 
-    // Ostream Operator
+    // Ostream operator
 
         friend Ostream& operator<< <Type, DType, LUType>
         (

src/OpenFOAM/matrices/LduMatrix/LduMatrix/LduMatrixSolver.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -131,9 +131,8 @@ Foam::LduMatrix<Type, DType, LUType>::solver::solver
 
     log_(1),
     minIter_(0),
-    maxIter_(lduMatrix::defaultMaxIter),
-    normType_(lduMatrix::normTypes::DEFAULT_NORM),
-    tolerance_(lduMatrix::defaultTolerance*pTraits<Type>::one),
+    maxIter_(defaultMaxIter_),
+    tolerance_(1e-6*pTraits<Type>::one),
     relTol_(Zero)
 {
     readControls();
@@ -146,8 +145,6 @@ template<class Type, class DType, class LUType>
 void Foam::LduMatrix<Type, DType, LUType>::solver::readControls()
 {
     controlDict_.readIfPresent("log", log_);
-    normType_ = lduMatrix::normTypes::DEFAULT_NORM;
-    lduMatrix::normTypesNames_.readIfPresent("norm", controlDict_, normType_);
     controlDict_.readIfPresent("minIter", minIter_);
     controlDict_.readIfPresent("maxIter", maxIter_);
     controlDict_.readIfPresent("tolerance", tolerance_);
@@ -171,45 +168,21 @@ Type Foam::LduMatrix<Type, DType, LUType>::solver::normFactor
 (
     const Field<Type>& psi,
     const Field<Type>& Apsi,
-    Field<Type>& tmpField,
-    const lduMatrix::normTypes normType
+    Field<Type>& tmpField
 ) const
 {
-    switch (normType)
-    {
-        case lduMatrix::normTypes::NO_NORM :
-        {
-            break;
-        }
+    // --- Calculate A dot reference value of psi
+    matrix_.sumA(tmpField);
+    cmptMultiply(tmpField, tmpField, gAverage(psi));
 
-        case lduMatrix::normTypes::DEFAULT_NORM :
-        case lduMatrix::normTypes::L1_SCALED_NORM :
-        {
-            // --- Calculate A dot reference value of psi
-            matrix_.sumA(tmpField);
-            cmptMultiply(tmpField, tmpField, gAverage(psi));
+    return stabilise
+    (
+        gSum(cmptMag(Apsi - tmpField) + cmptMag(matrix_.source() - tmpField)),
+        SolverPerformance<Type>::small_
+    );
 
-            return stabilise
-            (
-                gSum
-                (
-                    cmptMag(Apsi - tmpField)
-                  + cmptMag(matrix_.source() - tmpField)
-                ),
-                SolverPerformance<Type>::small_
-            );
-
-            // Equivalent at convergence:
-            // return stabilise
-            // (
-            //     2*gSumCmptMag(matrix_.source()), matrix_.small_
-            // );
-            break;
-        }
-    }
-
-    // Fall-through: no norm
-    return pTraits<Type>::one;
+    // At convergence this simpler method is equivalent to the above
+    // return stabilise(2*gSumCmptMag(matrix_.source()), matrix_.small_);
 }
 
 

src/OpenFOAM/matrices/LduMatrix/Solvers/DiagonalSolver/DiagonalSolver.H

@@ -34,8 +34,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_DiagonalSolver_H
-#define Foam_DiagonalSolver_H
+#ifndef DiagonalSolver_H
+#define DiagonalSolver_H
 
 #include "LduMatrix.H"
 
@@ -53,9 +53,7 @@ class DiagonalSolver
 :
     public LduMatrix<Type, DType, LUType>::solver
 {
-public:
-
-    // Generated Methods
+    // Private Member Functions
 
         //- No copy construct
         DiagonalSolver(const DiagonalSolver&) = delete;
@@ -64,6 +62,8 @@ public:
         void operator=(const DiagonalSolver&) = delete;
 
 
+public:
+
     //- Runtime type information
     TypeName("diagonal");
 

src/OpenFOAM/matrices/LduMatrix/Solvers/SmoothSolver/SmoothSolver.H

@@ -37,8 +37,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_SmoothSolver_H
-#define Foam_SmoothSolver_H
+#ifndef SmoothSolver_H
+#define SmoothSolver_H
 
 #include "lduMatrix.H"
 
@@ -56,9 +56,10 @@ class SmoothSolver
 :
     public LduMatrix<Type, DType, LUType>::solver
 {
+
 protected:
 
-    // Protected Data
+    // Protected data
 
         //- Number of sweeps before the evaluation of residual
         label nSweeps_;

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduCalculatedProcessorField/lduCalculatedProcessorField.C

@@ -32,11 +32,13 @@ License
 template<class Type>
 Foam::lduCalculatedProcessorField<Type>::lduCalculatedProcessorField
 (
-    const lduInterface& interface
+    const lduInterface& interface,
+    const Field<Type>& iF
 )
 :
     LduInterfaceField<Type>(interface),
     procInterface_(refCast<const lduPrimitiveProcessorInterface>(interface)),
+    field_(iF),
     sendBuf_(procInterface_.faceCells().size()),
     receiveBuf_(procInterface_.faceCells().size()),
     scalarSendBuf_(procInterface_.faceCells().size()),
@@ -54,6 +56,7 @@ Foam::lduCalculatedProcessorField<Type>::lduCalculatedProcessorField
 :
     LduInterfaceField<Type>(refCast<const lduInterface>(ptf)),
     procInterface_(ptf.procInterface_),
+    field_(ptf.field_),
     sendBuf_(procInterface_.faceCells().size()),
     receiveBuf_(procInterface_.faceCells().size()),
     scalarSendBuf_(procInterface_.faceCells().size()),

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/lduCalculatedProcessorField/lduCalculatedProcessorField.H

@@ -30,7 +30,8 @@ Group
     grpGenericBoundaryConditions
 
 Description
-    A lduProcessorField type bypassing coupledFvPatchField
+    A lduProcessorField type bypassing coupledFvPatchField and holding
+    a reference to the Field<Type>.
 
     Used to add updateInterfaceMatrix capabilities to a lduMatrix
     which is fully uncoupled from the fvMesh.
@@ -71,6 +72,8 @@ protected:
         //- Local reference cast into the interface
          const lduPrimitiveProcessorInterface& procInterface_;
 
+        //- Local Field
+         const Field<Type>& field_;
 
         // Sending and receiving
 
@@ -115,8 +118,8 @@ public:
         //- Construct from patch and internal field
         lduCalculatedProcessorField
         (
-            const lduInterface& interface
-            //const Field<Type>&
+            const lduInterface& interface,
+            const Field<Type>&
         );
 
         //- Construct as copy

src/OpenFOAM/matrices/lduMatrix/lduAddressing/lduInterface/processorLduInterfaceTemplates.C

@@ -47,7 +47,7 @@ void Foam::processorLduInterface::send
      || commsType == Pstream::commsTypes::scheduled
     )
     {
-        UOPstream::write
+        OPstream::write
         (
             commsType,
             neighbProcNo(),
@@ -61,7 +61,7 @@ void Foam::processorLduInterface::send
     {
         resizeBuf(receiveBuf_, nBytes);
 
-        UIPstream::read
+        IPstream::read
         (
             commsType,
             neighbProcNo(),
@@ -77,7 +77,7 @@ void Foam::processorLduInterface::send
             static_cast<void*>(sendBuf_.data()), f.cdata(), nBytes
         );
 
-        UOPstream::write
+        OPstream::write
         (
             commsType,
             neighbProcNo(),
@@ -109,7 +109,7 @@ void Foam::processorLduInterface::receive
      || commsType == Pstream::commsTypes::scheduled
     )
     {
-        UIPstream::read
+        IPstream::read
         (
             commsType,
             neighbProcNo(),
@@ -181,7 +181,7 @@ void Foam::processorLduInterface::compressedSend
          || commsType == Pstream::commsTypes::scheduled
         )
         {
-            UOPstream::write
+            OPstream::write
             (
                 commsType,
                 neighbProcNo(),
@@ -195,7 +195,7 @@ void Foam::processorLduInterface::compressedSend
         {
             resizeBuf(receiveBuf_, nBytes);
 
-            UIPstream::read
+            IPstream::read
             (
                 commsType,
                 neighbProcNo(),
@@ -205,7 +205,7 @@ void Foam::processorLduInterface::compressedSend
                 comm()
             );
 
-            UOPstream::write
+            OPstream::write
             (
                 commsType,
                 neighbProcNo(),
@@ -252,7 +252,7 @@ void Foam::processorLduInterface::compressedReceive
         {
             resizeBuf(receiveBuf_, nBytes);
 
-            UIPstream::read
+            IPstream::read
             (
                 commsType,
                 neighbProcNo(),

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.C

@@ -41,18 +41,7 @@ namespace Foam
 }
 
 
-const Foam::scalar Foam::lduMatrix::defaultTolerance = 1e-6;
-
-const Foam::Enum
-<
-    Foam::lduMatrix::normTypes
->
-Foam::lduMatrix::normTypesNames_
-({
-    { normTypes::NO_NORM, "none" },
-    { normTypes::DEFAULT_NORM, "default" },
-    { normTypes::L1_SCALED_NORM, "L1_scaled" },
-});
+const Foam::label Foam::lduMatrix::solver::defaultMaxIter_ = 1000;
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrix.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,8 +50,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_lduMatrix_H
-#define Foam_lduMatrix_H
+#ifndef lduMatrix_H
+#define lduMatrix_H
 
 #include "lduMesh.H"
 #include "primitiveFieldsFwd.H"
@@ -62,7 +62,6 @@ SourceFiles
 #include "runTimeSelectionTables.H"
 #include "solverPerformance.H"
 #include "InfoProxy.H"
-#include "Enum.H"
 #include "profilingTrigger.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -71,6 +70,7 @@ namespace Foam
 {
 
 // Forward Declarations
+
 class lduMatrix;
 
 Ostream& operator<<(Ostream&, const lduMatrix&);
@@ -95,26 +95,6 @@ class lduMatrix
 
 public:
 
-    // Public Types
-
-        //- Enumerated matrix normalisation types
-        enum class normTypes : char
-        {
-            NO_NORM,            //!< "none" norm (returns 1)
-            DEFAULT_NORM,       //!< "default" norm (== L1_scaled)
-            L1_SCALED_NORM,     //!< "L1_scaled" norm
-        };
-
-        //- Names for the normTypes
-        static const Enum<normTypes> normTypesNames_;
-
-        //- Default maximum number of iterations for solvers (1000)
-        static constexpr const label defaultMaxIter = 1000;
-
-        //- Default (absolute) tolerance (1e-6)
-        static const scalar defaultTolerance;
-
-
     //- Abstract base-class for lduMatrix solvers
     class solver
     {
@@ -122,16 +102,19 @@ public:
 
         // Protected Data
 
+            //- Default maximum number of iterations in the solver
+            static const label defaultMaxIter_;
+
             word fieldName_;
             const lduMatrix& matrix_;
             const FieldField<Field, scalar>& interfaceBouCoeffs_;
             const FieldField<Field, scalar>& interfaceIntCoeffs_;
             lduInterfaceFieldPtrsList interfaces_;
 
-            //- Dictionary of solution controls
+            //- Dictionary of controls
             dictionary controlDict_;
 
-            //- Verbosity level for solver output statements
+            //- Level of verbosity in the solver output statements
             int log_;
 
             //- Minimum number of iterations in the solver
@@ -140,22 +123,18 @@ public:
             //- Maximum number of iterations in the solver
             label maxIter_;
 
-            //- The normalisation type
-            lduMatrix::normTypes normType_;
-
             //- Final convergence tolerance
             scalar tolerance_;
 
             //- Convergence tolerance relative to the initial
             scalar relTol_;
 
-            //- Profiling instrumentation
             profilingTrigger profiling_;
 
 
         // Protected Member Functions
 
-            //- Read the control parameters from controlDict_
+            //- Read the control parameters from the controlDict_
             virtual void readControls();
 
 
@@ -216,7 +195,6 @@ public:
 
         // Constructors
 
-            //- Construct solver for given field name, matrix etc
             solver
             (
                 const word& fieldName,
@@ -294,16 +272,6 @@ public:
                 const direction cmpt=0
             ) const;
 
-            //- Return the matrix norm using the specified norm method
-            solveScalarField::cmptType normFactor
-            (
-                const solveScalarField& psi,
-                const solveScalarField& source,
-                const solveScalarField& Apsi,
-                solveScalarField& tmpField,
-                const lduMatrix::normTypes normType
-            ) const;
-
             //- Return the matrix norm used to normalise the residual for the
             //- stopping criterion
             solveScalarField::cmptType normFactor
@@ -312,10 +280,7 @@ public:
                 const solveScalarField& source,
                 const solveScalarField& Apsi,
                 solveScalarField& tmpField
-            ) const
-            {
-                return this->normFactor(psi, source, Apsi, tmpField, normType_);
-            }
+            ) const;
     };
 
 
@@ -389,7 +354,6 @@ public:
 
         // Constructors
 
-            //- Construct for given field name, matrix etc
             smoother
             (
                 const word& fieldName,
@@ -515,8 +479,10 @@ public:
 
         // Constructors
 
-            //- Construct for given solver
-            explicit preconditioner(const solver& sol)
+            preconditioner
+            (
+                const solver& sol
+            )
             :
                 solver_(sol)
             {}
@@ -598,7 +564,7 @@ public:
         // Access to addressing
 
             //- Return the LDU mesh from which the addressing is obtained
-            const lduMesh& mesh() const noexcept
+            const lduMesh& mesh() const
             {
                 return lduMesh_;
             }
@@ -640,38 +606,38 @@ public:
             const scalarField& diag() const;
             const scalarField& upper() const;
 
-            bool hasDiag() const noexcept
+            bool hasDiag() const
             {
                 return (diagPtr_);
             }
 
-            bool hasUpper() const noexcept
+            bool hasUpper() const
             {
                 return (upperPtr_);
             }
 
-            bool hasLower() const noexcept
+            bool hasLower() const
             {
                 return (lowerPtr_);
             }
 
-            bool diagonal() const noexcept
+            bool diagonal() const
             {
                 return (diagPtr_ && !lowerPtr_ && !upperPtr_);
             }
 
-            bool symmetric() const noexcept
+            bool symmetric() const
             {
                 return (diagPtr_ && (!lowerPtr_ && upperPtr_));
             }
 
-            bool asymmetric() const noexcept
+            bool asymmetric() const
             {
                 return (diagPtr_ && lowerPtr_ && upperPtr_);
             }
 
 
-        // Operations
+        // operations
 
             void sumDiag();
             void negSumDiag();

src/OpenFOAM/matrices/lduMatrix/lduMatrix/lduMatrixSolver.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2022 OpenCFD Ltd.
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -152,14 +152,6 @@ Foam::lduMatrix::solver::solver
     interfaceIntCoeffs_(interfaceIntCoeffs),
     interfaces_(interfaces),
     controlDict_(solverControls),
-
-    log_(1),
-    minIter_(0),
-    maxIter_(lduMatrix::defaultMaxIter),
-    normType_(lduMatrix::normTypes::DEFAULT_NORM),
-    tolerance_(lduMatrix::defaultTolerance),
-    relTol_(Zero),
-
     profiling_("lduMatrix::solver." + fieldName)
 {
     readControls();
@@ -170,19 +162,11 @@ Foam::lduMatrix::solver::solver
 
 void Foam::lduMatrix::solver::readControls()
 {
-    log_ = 1;
-    minIter_ = 0;
-    maxIter_ = lduMatrix::defaultMaxIter;
-    normType_ = lduMatrix::normTypes::DEFAULT_NORM;
-    tolerance_ = lduMatrix::defaultTolerance;
-    relTol_ = 0;
-
-    controlDict_.readIfPresent("log", log_);
-    lduMatrix::normTypesNames_.readIfPresent("norm", controlDict_, normType_);
-    controlDict_.readIfPresent("minIter", minIter_);
-    controlDict_.readIfPresent("maxIter", maxIter_);
-    controlDict_.readIfPresent("tolerance", tolerance_);
-    controlDict_.readIfPresent("relTol", relTol_);
+    log_ = controlDict_.getOrDefault<int>("log", 1);
+    minIter_ = controlDict_.getOrDefault<label>("minIter", 0);
+    maxIter_ = controlDict_.getOrDefault<label>("maxIter", defaultMaxIter_);
+    tolerance_ = controlDict_.getOrDefault<scalar>("tolerance", 1e-6);
+    relTol_ = controlDict_.getOrDefault<scalar>("relTol", 0);
 }
 
 
@@ -215,40 +199,24 @@ Foam::solveScalarField::cmptType Foam::lduMatrix::solver::normFactor
     const solveScalarField& psi,
     const solveScalarField& source,
     const solveScalarField& Apsi,
-    solveScalarField& tmpField,
-    const lduMatrix::normTypes normType
+    solveScalarField& tmpField
 ) const
 {
-    switch (normType)
-    {
-        case lduMatrix::normTypes::NO_NORM :
-        {
-            break;
-        }
+    // --- Calculate A dot reference value of psi
+    matrix_.sumA(tmpField, interfaceBouCoeffs_, interfaces_);
 
-        case lduMatrix::normTypes::DEFAULT_NORM :
-        case lduMatrix::normTypes::L1_SCALED_NORM :
-        {
-            // --- Calculate A dot reference value of psi
-            matrix_.sumA(tmpField, interfaceBouCoeffs_, interfaces_);
+    tmpField *= gAverage(psi, matrix_.mesh().comm());
 
-            tmpField *= gAverage(psi, matrix_.mesh().comm());
+    return
+        gSum
+        (
+            (mag(Apsi - tmpField) + mag(source - tmpField))(),
+            matrix_.mesh().comm()
+        )
+      + solverPerformance::small_;
 
-            return
-                gSum
-                (
-                    (mag(Apsi - tmpField) + mag(source - tmpField))(),
-                    matrix_.mesh().comm()
-                ) + solverPerformance::small_;
-
-            // Equivalent at convergence:
-            // return 2*gSumMag(source) + solverPerformance::small_;
-            break;
-        }
-    }
-
-    // Fall-through: no norm
-    return solveScalarField::cmptType(1);
+    // At convergence this simpler method is equivalent to the above
+    // return 2*gSumMag(source) + solverPerformance::small_;
 }
 
 

src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.C

@@ -67,6 +67,12 @@ Foam::GAMGPreconditioner::GAMGPreconditioner
 }
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::GAMGPreconditioner::~GAMGPreconditioner()
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 void Foam::GAMGPreconditioner::readControls()
@@ -83,7 +89,7 @@ void Foam::GAMGPreconditioner::precondition
     const direction cmpt
 ) const
 {
-    wA = Zero;
+    wA = 0.0;
     solveScalarField AwA(wA.size());
     solveScalarField finestCorrection(wA.size());
     solveScalarField finestResidual(rA_ss);

src/OpenFOAM/matrices/lduMatrix/preconditioners/GAMGPreconditioner/GAMGPreconditioner.H

@@ -61,8 +61,7 @@ class GAMGPreconditioner
     public lduMatrix::preconditioner
 {
 protected:
-
-    // Protected Data
+    // Protected data
 
         //- Number of V-cycles to perform
         label nVcycles_;
@@ -87,7 +86,7 @@ public:
 
 
     //- Destructor
-    virtual ~GAMGPreconditioner() = default;
+    virtual ~GAMGPreconditioner();
 
 
     // Member Functions

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGProcAgglomerations/procFacesGAMGProcAgglomeration/procFacesGAMGProcAgglomeration.H

@@ -51,7 +51,6 @@ SourceFiles
 namespace Foam
 {
 
-// Forward Declarations
 class GAMGAgglomeration;
 class lduMesh;
 class lduPrimitiveMesh;
@@ -64,7 +63,7 @@ class procFacesGAMGProcAgglomeration
 :
     public GAMGProcAgglomeration
 {
-    // Private Data
+    // Private data
 
         //- When to processor agglomerate
         const label nAgglomeratingCells_;
@@ -122,8 +121,9 @@ public:
 
     // Member Functions
 
-        //- Modify agglomeration. Return true if modified
+       //- Modify agglomeration. Return true if modified
         virtual bool agglomerate();
+
 };
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2021-2022 OpenCFD Ltd.
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -69,6 +69,7 @@ Foam::GAMGSolver::GAMGSolver
 
     // Default values for all controls
     // which may be overridden by those in controlDict
+    cacheAgglomeration_(true),
     nPreSweeps_(0),
     preSweepsLevelMultiplier_(1),
     maxPreSweeps_(4),
@@ -76,12 +77,9 @@ Foam::GAMGSolver::GAMGSolver
     postSweepsLevelMultiplier_(1),
     maxPostSweeps_(4),
     nFinestSweeps_(2),
-
-    cacheAgglomeration_(true),
     interpolateCorrection_(false),
     scaleCorrection_(matrix.symmetric()),
     directSolveCoarsest_(false),
-
     agglomeration_(GAMGAgglomeration::New(matrix_, controlDict_)),
 
     matrixLevels_(agglomeration_.size()),
@@ -391,7 +389,14 @@ void Foam::GAMGSolver::readControls()
 
 const Foam::lduMatrix& Foam::GAMGSolver::matrixLevel(const label i) const
 {
-    return i ? matrixLevels_[i-1] : matrix_;
+    if (i == 0)
+    {
+        return matrix_;
+    }
+    else
+    {
+        return matrixLevels_[i - 1];
+    }
 }
 
 
@@ -400,7 +405,14 @@ const Foam::lduInterfaceFieldPtrsList& Foam::GAMGSolver::interfaceLevel
     const label i
 ) const
 {
-    return i ? interfaceLevels_[i-1] : interfaces_;
+    if (i == 0)
+    {
+        return interfaces_;
+    }
+    else
+    {
+        return interfaceLevels_[i - 1];
+    }
 }
 
 
@@ -410,7 +422,14 @@ Foam::GAMGSolver::interfaceBouCoeffsLevel
     const label i
 ) const
 {
-    return i ? interfaceLevelsBouCoeffs_[i-1] : interfaceBouCoeffs_;
+    if (i == 0)
+    {
+        return interfaceBouCoeffs_;
+    }
+    else
+    {
+        return interfaceLevelsBouCoeffs_[i - 1];
+    }
 }
 
 
@@ -420,7 +439,14 @@ Foam::GAMGSolver::interfaceIntCoeffsLevel
     const label i
 ) const
 {
-    return i ? interfaceLevelsIntCoeffs_[i-1] : interfaceIntCoeffs_;
+    if (i == 0)
+    {
+        return interfaceIntCoeffs_;
+    }
+    else
+    {
+        return interfaceLevelsIntCoeffs_[i - 1];
+    }
 }
 
 

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolver.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -56,11 +56,12 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_GAMGSolver_H
-#define Foam_GAMGSolver_H
+#ifndef GAMGSolver_H
+#define GAMGSolver_H
 
 #include "GAMGAgglomeration.H"
 #include "lduMatrix.H"
+#include "labelField.H"
 #include "primitiveFields.H"
 #include "LUscalarMatrix.H"
 
@@ -77,7 +78,9 @@ class GAMGSolver
 :
     public lduMatrix::solver
 {
-    // Private Data
+    // Private data
+
+        bool cacheAgglomeration_;
 
         //- Number of pre-smoothing sweeps
         label nPreSweeps_;
@@ -100,9 +103,6 @@ class GAMGSolver
         //- Number of smoothing sweeps on finest mesh
         label nFinestSweeps_;
 
-        //- Cache the agglomeration (default: true)
-        bool cacheAgglomeration_;
-
         //- Choose if the corrections should be interpolated after injection.
         //  By default corrections are not interpolated.
         bool interpolateCorrection_;

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/GAMGSolverAgglomerateMatrix.C

@@ -318,7 +318,7 @@ void Foam::GAMGSolver::gatherMatrices
         {
             label otherI = proci-1;
 
-            IPstream fromProc
+            IPstream fromSlave
             (
                 Pstream::commsTypes::scheduled,
                 procIDs[proci],
@@ -327,14 +327,14 @@ void Foam::GAMGSolver::gatherMatrices
                 meshComm
             );
 
-            otherMats.set(otherI, new lduMatrix(dummyMesh, fromProc));
+            otherMats.set(otherI, new lduMatrix(dummyMesh, fromSlave));
 
             // Receive number of/valid interfaces
             boolList& procTransforms = otherTransforms[otherI];
             List<label>& procRanks = otherRanks[otherI];
 
-            fromProc >> procTransforms;
-            fromProc >> procRanks;
+            fromSlave >> procTransforms;
+            fromSlave >> procRanks;
 
             // Size coefficients
             otherBouCoeffs.set
@@ -354,12 +354,12 @@ void Foam::GAMGSolver::gatherMatrices
                     otherBouCoeffs[otherI].set
                     (
                         intI,
-                        new scalarField(fromProc)
+                        new scalarField(fromSlave)
                     );
                     otherIntCoeffs[otherI].set
                     (
                         intI,
-                        new scalarField(fromProc)
+                        new scalarField(fromSlave)
                     );
                 }
             }

src/OpenFOAM/matrices/lduMatrix/solvers/GAMG/interfaceFields/processorGAMGInterfaceField/processorGAMGInterfaceField.C

@@ -110,7 +110,7 @@ void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate
         // Fast path.
         scalarReceiveBuf_.setSize(scalarSendBuf_.size());
         outstandingRecvRequest_ = UPstream::nRequests();
-        UIPstream::read
+        IPstream::read
         (
             Pstream::commsTypes::nonBlocking,
             procInterface_.neighbProcNo(),
@@ -121,7 +121,7 @@ void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate
         );
 
         outstandingSendRequest_ = UPstream::nRequests();
-        UOPstream::write
+        OPstream::write
         (
             Pstream::commsTypes::nonBlocking,
             procInterface_.neighbProcNo(),

src/OpenFOAM/matrices/lduMatrix/solvers/PBiCG/PBiCG.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -224,13 +224,11 @@ Foam::solverPerformance Foam::PBiCG::solve
     }
 
     // Recommend PBiCGStab if PBiCG fails to converge
-    const label upperMaxIters = max(maxIter_, lduMatrix::defaultMaxIter);
-
-    if (solverPerf.nIterations() > upperMaxIters)
+    if (solverPerf.nIterations() > max(defaultMaxIter_, maxIter_))
     {
         FatalErrorInFunction
-            << "PBiCG has failed to converge within the maximum iterations: "
-            << upperMaxIters << nl
+            << "PBiCG has failed to converge within the maximum number"
+               " of iterations " << max(defaultMaxIter_, maxIter_) << nl
             << "    Please try the more robust PBiCGStab solver."
             << exit(FatalError);
     }

src/OpenFOAM/matrices/lduMatrix/solvers/diagonalSolver/diagonalSolver.H

@@ -37,8 +37,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_diagonalSolver_H
-#define Foam_diagonalSolver_H
+#ifndef diagonalSolver_H
+#define diagonalSolver_H
 
 #include "lduMatrix.H"
 
@@ -55,9 +55,7 @@ class diagonalSolver
 :
     public lduMatrix::solver
 {
-public:
-
-    // Generated Methods
+    // Private Member Functions
 
         //- No copy construct
         diagonalSolver(const diagonalSolver&) = delete;
@@ -66,6 +64,8 @@ public:
         void operator=(const diagonalSolver&) = delete;
 
 
+public:
+
     //- Runtime type information
     TypeName("diagonal");
 

src/OpenFOAM/matrices/lduMatrix/solvers/smoothSolver/smoothSolver.H

@@ -43,8 +43,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_smoothSolver_H
-#define Foam_smoothSolver_H
+#ifndef smoothSolver_H
+#define smoothSolver_H
 
 #include "lduMatrix.H"
 
@@ -63,7 +63,7 @@ class smoothSolver
 {
 protected:
 
-    // Protected Data
+    // Protected data
 
         //- Number of sweeps before the evaluation of residual
         label nSweeps_;
@@ -95,7 +95,6 @@ public:
     //- Destructor
     virtual ~smoothSolver() = default;
 
-
     // Member Functions
 
         //- Solve the matrix with this solver

src/OpenFOAM/meshes/lduMesh/lduMesh.H

@@ -33,8 +33,8 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_lduMesh_H
-#define Foam_lduMesh_H
+#ifndef lduMesh_H
+#define lduMesh_H
 
 #include "lduAddressing.H"
 #include "lduInterfacePtrsList.H"
@@ -46,8 +46,11 @@ Description
 namespace Foam
 {
 
-// Forward Declarations
 class objectRegistry;
+
+
+// Forward declaration of friend functions and operators
+
 class lduMesh;
 
 Ostream& operator<<(Ostream&, const InfoProxy<lduMesh>&);
@@ -59,12 +62,15 @@ Ostream& operator<<(Ostream&, const InfoProxy<lduMesh>&);
 
 class lduMesh
 {
+
 public:
 
     //- Runtime type information
     TypeName("lduMesh");
 
 
+    // Constructors
+
     //- Destructor
     virtual ~lduMesh() = default;
 
@@ -108,7 +114,7 @@ public:
             }
 
 
-    // Ostream Operator
+    // Ostream operator
 
         friend Ostream& operator<<(Ostream&, const InfoProxy<lduMesh>&);
 };

src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C

@@ -37,6 +37,25 @@ License
 namespace Foam
 {
     defineTypeNameAndDebug(lduPrimitiveMesh, 0);
+
+    //- Less operator for pairs of \<processor\>\<index\>
+    class procLess
+    {
+        const labelPairList& list_;
+
+    public:
+
+        procLess(const labelPairList& list)
+        :
+            list_(list)
+        {}
+
+        bool operator()(const label a, const label b)
+        {
+            return list_[a].first() < list_[b].first();
+        }
+    };
+
 }
 
 
@@ -1073,15 +1092,15 @@ void Foam::lduPrimitiveMesh::gather
 
     if (Pstream::myProcNo(comm) == procIDs[0])
     {
-        // Master.
         otherMeshes.setSize(procIDs.size()-1);
 
+        // Slave meshes
         for (label i = 1; i < procIDs.size(); ++i)
         {
             //Pout<< "on master :"
-            //    << " receiving from proc " << procIDs[i] << endl;
+            //    << " receiving from slave " << procIDs[i] << endl;
 
-            IPstream fromProc
+            IPstream fromSlave
             (
                 Pstream::commsTypes::scheduled,
                 procIDs[i],
@@ -1090,10 +1109,10 @@ void Foam::lduPrimitiveMesh::gather
                 comm
             );
 
-            label nCells = readLabel(fromProc);
-            labelList lowerAddr(fromProc);
-            labelList upperAddr(fromProc);
-            boolList validInterface(fromProc);
+            label nCells = readLabel(fromSlave);
+            labelList lowerAddr(fromSlave);
+            labelList upperAddr(fromSlave);
+            boolList validInterface(fromSlave);
 
 
             // Construct mesh without interfaces
@@ -1116,7 +1135,7 @@ void Foam::lduPrimitiveMesh::gather
             {
                 if (validInterface[intI])
                 {
-                    word coupleType(fromProc);
+                    word coupleType(fromSlave);
 
                     newInterfaces.set
                     (
@@ -1126,7 +1145,7 @@ void Foam::lduPrimitiveMesh::gather
                             coupleType,
                             intI,
                             otherMeshes[i-1].rawInterfaces(),
-                            fromProc
+                            fromSlave
                         ).ptr()
                     );
                 }

src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H

@@ -35,8 +35,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_lduPrimitiveMesh_H
-#define Foam_lduPrimitiveMesh_H
+#ifndef lduPrimitiveMesh_H
+#define lduPrimitiveMesh_H
 
 #include "lduMesh.H"
 #include "labelList.H"
@@ -55,7 +55,7 @@ class lduPrimitiveMesh
     public lduMesh,
     public lduAddressing
 {
-    // Private Data
+    // Private data
 
         //- Lower addressing
         labelList lowerAddr_;

src/OpenFOAM/meshes/meshShapes/cell/cell.C

@@ -27,7 +27,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "cell.H"
-#include "pyramid.H"
+#include "pyramidPointFaceRef.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/meshShapes/cellModel/cellModel.C

@@ -27,7 +27,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "cellModel.H"
-#include "pyramid.H"
+#include "pyramidPointFaceRef.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/meshShapes/edge/edge.H

@@ -45,11 +45,11 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_edge_H
-#define Foam_edge_H
+#ifndef edge_H
+#define edge_H
 
 #include "labelPair.H"
-#include "line.H"
+#include "linePointRef.H"
 #include "pointField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/meshes/meshShapes/face/face.C

@@ -28,7 +28,7 @@ License
 
 #include "face.H"
 #include "triFace.H"
-#include "triangle.H"
+#include "triPointRef.H"
 #include "mathematicalConstants.H"
 #include "Circulator.H"
 #include <algorithm>

src/OpenFOAM/meshes/meshShapes/face/face.H

@@ -43,8 +43,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_face_H
-#define Foam_face_H
+#ifndef face_H
+#define face_H
 
 #include "pointField.H"
 #include "labelList.H"
@@ -435,12 +435,6 @@ public:
         static bool sameVertices(const face& a, const face& b);
 
 
-    // Member Operators
-
-        //- Increment (offset) vertices by given amount
-        inline void operator+=(const label vertexOffset);
-
-
     // Hashing
 
         //- The symmetric hash value for face.
@@ -500,11 +494,19 @@ template<> struct Hash<face> : face::hasher {};
 template<>
 struct offsetOp<face>
 {
-    face operator()(const face& x, const label offset) const
+    inline face operator()
+    (
+        const face& x,
+        const label offset
+    ) const
     {
-        face f(x);
-        f += offset;
-        return f;
+        face result(x.size());
+
+        forAll(x, i)
+        {
+            result[i] = x[i] + offset;
+        }
+        return result;
     }
 };
 

src/OpenFOAM/meshes/meshShapes/face/faceI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -190,20 +190,6 @@ inline Foam::label Foam::face::nTriangles() const
 }
 
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-inline void Foam::face::operator+=(const label vertexOffset)
-{
-    if (vertexOffset)
-    {
-        for (label& vrt : static_cast<labelList&>(*this))
-        {
-            vrt += vertexOffset;
-        }
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Global Operators   * * * * * * * * * * * * * //
 
 inline bool Foam::operator==(const face& a, const face& b)

src/OpenFOAM/meshes/meshShapes/face/faceIntersection.C

@@ -27,7 +27,7 @@ License
 
 #include "face.H"
 #include "pointHit.H"
-#include "triangle.H"
+#include "triPointRef.H"
 #include "line.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/meshes/meshShapes/face/oppositeFace.H

@@ -30,10 +30,13 @@ Description
     Class containing opposite face for a prismatic cell with addressing
     and a possibility of failure.
 
+SourceFiles
+    oppositeFace.C
+
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_oppositeFace_H
-#define Foam_oppositeFace_H
+#ifndef oppositeFace_H
+#define oppositeFace_H
 
 #include "face.H"
 
@@ -50,7 +53,7 @@ class oppositeFace
 :
     public face
 {
-    // Private Data
+    // Private data
 
         //- Master face index
         const label masterIndex_;
@@ -81,19 +84,19 @@ public:
     // Member Functions
 
         //- Master face index
-        label masterIndex() const noexcept
+        inline label masterIndex() const
         {
             return masterIndex_;
         }
 
-        //- Opposite face index
-        label oppositeIndex() const noexcept
+        //- Slave face index
+        inline label oppositeIndex() const
         {
             return oppositeIndex_;
         }
 
         //- Does the opposite face exist?
-        bool found() const noexcept
+        inline bool found() const
         {
             return oppositeIndex_ >= 0;
         }

src/OpenFOAM/meshes/meshShapes/labelledTri/labelledTri.H

@@ -35,8 +35,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_labelledTri_H
-#define Foam_labelledTri_H
+#ifndef labelledTri_H
+#define labelledTri_H
 
 #include "triFace.H"
 #include "ListListOps.H"
@@ -91,9 +91,9 @@ public:
         //- and optional region index (0 if unspecified)
         inline labelledTri
         (
-            const label p0,
-            const label p1,
-            const label p2,
+            const label a,
+            const label b,
+            const label c,
             const label region = 0
         );
 
@@ -179,9 +179,13 @@ struct offsetOp<labelledTri>
 {
     labelledTri operator()(const labelledTri& x, const label offset) const
     {
-        labelledTri f(x);
-        f += offset;
-        return f;
+        labelledTri result(x);
+
+        forAll(x, xi)
+        {
+            result[xi] = x[xi] + offset;
+        }
+        return result;
     }
 };
 

src/OpenFOAM/meshes/meshShapes/labelledTri/labelledTriI.H

@@ -77,13 +77,13 @@ inline Foam::labelledTri::labelledTri
 
 inline Foam::labelledTri::labelledTri
 (
-    const label p0,
-    const label p1,
-    const label p2,
+    const label a,
+    const label b,
+    const label c,
     const label region
 )
 :
-    triFace(p0, p1, p2),
+    triFace(a, b, c),
     index_(region)
 {}
 

src/OpenFOAM/meshes/meshShapes/tetCell/tetCell.H

@@ -40,15 +40,15 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_tetCell_H
-#define Foam_tetCell_H
+#ifndef tetCell_H
+#define tetCell_H
 
 #include "FixedList.H"
 #include "triFace.H"
 #include "faceList.H"
 #include "edgeList.H"
 #include "pointField.H"
-#include "tetrahedron.H"
+#include "tetPointRef.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -85,10 +85,10 @@ public:
         //- Construct from four point labels
         inline tetCell
         (
-            const label p0,
-            const label p1,
-            const label p2,
-            const label p3
+            const label a,
+            const label b,
+            const label c,
+            const label d
         );
 
         //- Construct from an initializer list of four point labels

src/OpenFOAM/meshes/meshShapes/tetCell/tetCellI.H

@@ -36,16 +36,16 @@ inline Foam::tetCell::tetCell()
 
 inline Foam::tetCell::tetCell
 (
-    const label p0,
-    const label p1,
-    const label p2,
-    const label p3
+    const label a,
+    const label b,
+    const label c,
+    const label d
 )
 {
-    operator[](0) = p0;
-    operator[](1) = p1;
-    operator[](2) = p2;
-    operator[](3) = p3;
+    operator[](0) = a;
+    operator[](1) = b;
+    operator[](2) = c;
+    operator[](3) = d;
 }
 
 

src/OpenFOAM/meshes/meshShapes/triFace/triFace.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2022 OpenCFD Ltd.
+    Copyright (C) 2017-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,7 +48,7 @@ SourceFiles
 #include "pointHit.H"
 #include "intersection.H"
 #include "pointField.H"
-#include "triangle.H"
+#include "triPointRef.H"
 #include "ListListOps.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -79,7 +79,12 @@ public:
         inline triFace();
 
         //- Construct from three point labels
-        inline triFace(const label p0, const label p1, const label p2);
+        inline triFace
+        (
+            const label a,
+            const label b,
+            const label c
+        );
 
         //- Construct from an initializer list of three point labels
         inline explicit triFace(std::initializer_list<label> list);
@@ -316,12 +321,6 @@ public:
         static inline int compare(const triFace& a, const triFace& b);
 
 
-    // Member Operators
-
-        //- Increment (offset) vertices by given amount
-        inline void operator+=(const label vertexOffset);
-
-
     // Hashing
 
         //- The (commutative) hash value for triFace
@@ -378,11 +377,19 @@ template<> struct Hash<triFace> : triFace::hasher {};
 template<>
 struct offsetOp<triFace>
 {
-    triFace operator()(const triFace& x, const label offset) const
+    inline triFace operator()
+    (
+        const triFace& x,
+        const label offset
+    ) const
     {
-        triFace f(x);
-        f += offset;
-        return f;
+        triFace result;
+
+        forAll(x, i)
+        {
+            result[i] = x[i] + offset;
+        }
+        return result;
     }
 };
 

src/OpenFOAM/meshes/meshShapes/triFace/triFaceI.H

@@ -28,6 +28,7 @@ License
 
 #include "IOstreams.H"
 #include "face.H"
+#include "triPointRef.H"
 #include <algorithm>  // For std::swap
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
@@ -67,11 +68,16 @@ inline Foam::triFace::triFace()
 {}
 
 
-inline Foam::triFace::triFace(const label p0, const label p1, const label p2)
+inline Foam::triFace::triFace
+(
+    const label a,
+    const label b,
+    const label c
+)
 {
-    operator[](0) = p0;
-    operator[](1) = p1;
-    operator[](2) = p2;
+    operator[](0) = a;
+    operator[](1) = b;
+    operator[](2) = c;
 }
 
 
@@ -470,19 +476,6 @@ inline int Foam::triFace::edgeDirection(const Foam::edge& e) const
 }
 
 
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-inline void Foam::triFace::operator+=(const label vertexOffset)
-{
-    if (vertexOffset)
-    {
-        (*this)[0] += vertexOffset;
-        (*this)[1] += vertexOffset;
-        (*this)[2] += vertexOffset;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::operator==(const triFace& a, const triFace& b)

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeBaseSubset.C

@@ -158,7 +158,7 @@ void Foam::mapDistributeBase::exchangeMasks
         {
             if (proci != myRank && recvMasks[proci].size())
             {
-                UIPstream::read
+                IPstream::read
                 (
                     UPstream::commsTypes::nonBlocking,
                     proci,
@@ -174,7 +174,7 @@ void Foam::mapDistributeBase::exchangeMasks
         {
             if (proci != myRank && sendMasks[proci].size())
             {
-                UOPstream::write
+                OPstream::write
                 (
                     UPstream::commsTypes::nonBlocking,
                     proci,

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeBaseTemplates.C

@@ -520,7 +520,7 @@ void Foam::mapDistributeBase::distribute
                     sendFields[domain] =
                         accessAndFlip(field, map, subHasFlip, negOp);
 
-                    UOPstream::write
+                    OPstream::write
                     (
                         Pstream::commsTypes::nonBlocking,
                         domain,
@@ -543,7 +543,8 @@ void Foam::mapDistributeBase::distribute
                 if (domain != myRank && map.size())
                 {
                     recvFields[domain].resize(map.size());
-                    UIPstream::read
+
+                    IPstream::read
                     (
                         Pstream::commsTypes::nonBlocking,
                         domain,
@@ -555,6 +556,7 @@ void Foam::mapDistributeBase::distribute
                 }
             }
 
+
             // Set up 'send' to myself
             {
                 sendFields[myRank] =
@@ -980,7 +982,7 @@ void Foam::mapDistributeBase::distribute
                     sendFields[domain] =
                         accessAndFlip(field, map, subHasFlip, negOp);
 
-                    UOPstream::write
+                    OPstream::write
                     (
                         Pstream::commsTypes::nonBlocking,
                         domain,
@@ -1002,7 +1004,7 @@ void Foam::mapDistributeBase::distribute
 
                 if (domain != myRank && map.size())
                 {
-                    recvFields[domain].resize(map.size());
+                    recvFields[domain].setSize(map.size());
                     UIPstream::read
                     (
                         Pstream::commsTypes::nonBlocking,
@@ -1016,6 +1018,7 @@ void Foam::mapDistributeBase::distribute
             }
 
             // Set up 'send' to myself
+
             {
                 sendFields[myRank] =
                     accessAndFlip(field, subMap[myRank], subHasFlip, negOp);

src/OpenFOAM/meshes/polyMesh/polyMeshCheck/polyMeshTools.C

@@ -28,7 +28,7 @@ License
 
 #include "polyMeshTools.H"
 #include "syncTools.H"
-#include "pyramid.H"
+#include "pyramidPointFaceRef.H"
 #include "primitiveMeshTools.H"
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/polyMeshTetDecomposition.H

@@ -39,12 +39,13 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_polyMeshTetDecomposition_H
-#define Foam_polyMeshTetDecomposition_H
+#ifndef polyMeshTetDecomposition_H
+#define polyMeshTetDecomposition_H
 
 #include "polyMesh.H"
 #include "coupledPolyPatch.H"
 #include "syncTools.H"
+#include "tetPointRef.H"
 #include "tetIndices.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -53,7 +54,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                  Class polyMeshTetDecomposition Declaration
+                Class polyMeshTetDecomposition Declaration
 \*---------------------------------------------------------------------------*/
 
 class polyMeshTetDecomposition

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndices.H

@@ -59,8 +59,8 @@ SourceFiles
 #define Foam_tetIndices_H
 
 #include "label.H"
-#include "tetrahedron.H"
-#include "triangle.H"
+#include "tetPointRef.H"
+#include "triPointRef.H"
 #include "polyMesh.H"
 #include "triFace.H"
 #include "face.H"
@@ -76,6 +76,7 @@ class tetIndices;
 Istream& operator>>(Istream&, tetIndices&);
 Ostream& operator<<(Ostream&, const tetIndices&);
 
+
 /*---------------------------------------------------------------------------*\
                          Class tetIndices Declaration
 \*---------------------------------------------------------------------------*/
@@ -148,49 +149,35 @@ public:
         label& tetPt() noexcept { return tetPti_; }
 
 
-    // Mesh-related
+    // Searching
 
-        //- The indices corresponding to the tri on the face for this tet.
-        //- The normal of the tri points out of the cell.
+        //- Return the indices corresponding to the tri on the face for
+        //  this tet. The normal of the tri points out of the cell
         inline triFace faceTriIs
         (
             const polyMesh& mesh,
             const bool warn = true
         ) const;
 
-        //- Local indices corresponding to the tri on the face for this tet.
-        //- The normal of the tri points out of the cell.
+        //- Return the local indices corresponding to the tri on the face
+        //- for this tet. The normal of the tri points out of the cell
         inline triFace triIs
         (
             const polyMesh& mesh,
             const bool warn = true
         ) const;
 
-
-        //- The tet geometry for this tet,
-        //- where point0 is the cell centre.
+        //- Return the geometry corresponding to this tet
         inline tetPointRef tet(const polyMesh& mesh) const;
 
-        //- The tet geometry for this tet (using old positions),
-        //- where point0 is the cell centre.
-        inline tetPointRef oldTet(const polyMesh& mesh) const;
-
-        //- The triangle geometry for the face for this tet.
-        //- The normal of the tri points out of the cell
+        //- Return the geometry corresponding to the tri on the face for
+        //- this tet. The normal of the tri points out of the cell
         inline triPointRef faceTri(const polyMesh& mesh) const;
 
-        //- The triangle geometry for the face for this tet
-        //- (using old positions)
+        //- Return the geometry corresponding to the tri on the face for
+        //- this tet using the old positions
         inline triPointRef oldFaceTri(const polyMesh& mesh) const;
 
-        //- The x/y/z position for given barycentric coordinates
-        //- (where point0 is the cell centre).
-        inline point barycentricToPoint
-        (
-            const polyMesh& mesh,
-            const barycentric& bary
-        ) const;
-
 
     // Other
 

src/OpenFOAM/meshes/polyMesh/polyMeshTetDecomposition/tetIndicesI.H

@@ -57,9 +57,9 @@ inline Foam::triFace Foam::tetIndices::faceTriIs
     const bool warn
 ) const
 {
-    const Foam::face& f = mesh.faces()[facei_];
+    const Foam::face& f = mesh.faces()[face()];
 
-    label faceBasePtI = mesh.tetBasePtIs()[facei_];
+    label faceBasePtI = mesh.tetBasePtIs()[face()];
 
     if (faceBasePtI < 0)
     {
@@ -67,25 +67,24 @@ inline Foam::triFace Foam::tetIndices::faceTriIs
 
         if (warn && nWarnings_ < maxNWarnings)
         {
+            ++nWarnings_;
             WarningInFunction
-                << "No base point for face " << facei_ << ", " << f
+                << "No base point for face " << face() << ", " << f
                 << ", produces a valid tet decomposition." << endl;
-
-            if (++nWarnings_ == maxNWarnings)
+            if (nWarnings_ == maxNWarnings)
             {
                 Warning
-                    << "Suppressing further warnings." << endl;
+                    << "Suppressing any further warnings." << endl;
             }
         }
     }
 
-    const label facePtI = (tetPti_ + faceBasePtI) % f.size();
-    const label faceOtherPtI = f.fcIndex(facePtI);
+    label facePtI = (tetPti_ + faceBasePtI) % f.size();
+    label faceOtherPtI = f.fcIndex(facePtI);
 
-    if (mesh.faceOwner()[facei_] != celli_)
+    if (mesh.faceOwner()[face()] != cell())
     {
-        // Neighbour: return flipped face
-        return triFace(f[faceBasePtI], f[faceOtherPtI], f[facePtI]);
+        std::swap(facePtI, faceOtherPtI);
     }
 
     return triFace(f[faceBasePtI], f[facePtI], f[faceOtherPtI]);
@@ -98,9 +97,9 @@ inline Foam::triFace Foam::tetIndices::triIs
     const bool warn
 ) const
 {
-    const Foam::face& f = mesh.faces()[facei_];
+    const Foam::face& f = mesh.faces()[face()];
 
-    label faceBasePtI = mesh.tetBasePtIs()[facei_];
+    label faceBasePtI = mesh.tetBasePtIs()[face()];
 
     if (faceBasePtI < 0)
     {
@@ -108,25 +107,24 @@ inline Foam::triFace Foam::tetIndices::triIs
 
         if (warn && nWarnings_ < maxNWarnings)
         {
+            ++nWarnings_;
             WarningInFunction
-                << "No base point for face " << facei_ << ", " << f
+                << "No base point for face " << face() << ", " << f
                 << ", produces a valid tet decomposition." << endl;
-
-            if (++nWarnings_ == maxNWarnings)
+            if (nWarnings_ == maxNWarnings)
             {
                 Warning
-                    << "Suppressing further warnings." << endl;
+                    << "Suppressing any further warnings." << endl;
             }
         }
     }
 
-    const label facePtI = (tetPti_ + faceBasePtI) % f.size();
-    const label faceOtherPtI = f.fcIndex(facePtI);
+    label facePtI = (tetPti_ + faceBasePtI) % f.size();
+    label faceOtherPtI = f.fcIndex(facePtI);
 
-    if (mesh.faceOwner()[facei_] != celli_)
+    if (mesh.faceOwner()[face()] != cell())
     {
-        // Neighbour: return flipped face
-        return triFace(faceBasePtI, faceOtherPtI, facePtI);
+        std::swap(facePtI, faceOtherPtI);
     }
 
     return triFace(faceBasePtI, facePtI, faceOtherPtI);
@@ -135,60 +133,30 @@ inline Foam::triFace Foam::tetIndices::triIs
 
 inline Foam::tetPointRef Foam::tetIndices::tet(const polyMesh& mesh) const
 {
-    const pointField& pts = mesh.points();
+    const pointField& meshPoints = mesh.points();
     const triFace tri = faceTriIs(mesh);
 
     return tetPointRef
     (
-        mesh.cellCentres()[celli_],
-        pts[tri[0]],
-        pts[tri[1]],
-        pts[tri[2]]
-    );
-}
-
-
-inline Foam::tetPointRef Foam::tetIndices::oldTet(const polyMesh& mesh) const
-{
-    const pointField& pts = mesh.oldPoints();
-    const triFace tri = faceTriIs(mesh);
-
-    return tetPointRef
-    (
-        mesh.oldCellCentres()[celli_],
-        pts[tri[0]],
-        pts[tri[1]],
-        pts[tri[2]]
-    );
-}
-
-
-inline Foam::point Foam::tetIndices::barycentricToPoint
-(
-    const polyMesh& mesh,
-    const barycentric& bary
-) const
-{
-    const pointField& pts = mesh.points();
-    const triFace tri = faceTriIs(mesh);
-
-    const point& a = mesh.cellCentres()[celli_];
-    const point& b = pts[tri[0]];
-    const point& c = pts[tri[1]];
-    const point& d = pts[tri[2]];
-
-    return point
-    (
-        (bary.a()*a.x() + bary.b()*b.x() + bary.c()*c.x() + bary.d()*d.x()),
-        (bary.a()*a.y() + bary.b()*b.y() + bary.c()*c.y() + bary.d()*d.y()),
-        (bary.a()*a.z() + bary.b()*b.z() + bary.c()*c.z() + bary.d()*d.z())
+        mesh.cellCentres()[cell()],
+        meshPoints[tri[0]],
+        meshPoints[tri[1]],
+        meshPoints[tri[2]]
     );
 }
 
 
 inline Foam::triPointRef Foam::tetIndices::faceTri(const polyMesh& mesh) const
 {
-    return faceTriIs(mesh).tri(mesh.points());
+    const pointField& meshPoints = mesh.points();
+    const triFace tri = faceTriIs(mesh);
+
+    return triPointRef
+    (
+        meshPoints[tri[0]],
+        meshPoints[tri[1]],
+        meshPoints[tri[2]]
+    );
 }
 
 
@@ -197,7 +165,15 @@ inline Foam::triPointRef Foam::tetIndices::oldFaceTri
     const polyMesh& mesh
 ) const
 {
-    return faceTriIs(mesh).tri(mesh.oldPoints());
+    const pointField& meshOldPoints = mesh.oldPoints();
+    const triFace tri = faceTriIs(mesh);
+
+    return triPointRef
+    (
+        meshOldPoints[tri[0]],
+        meshOldPoints[tri[1]],
+        meshOldPoints[tri[2]]
+    );
 }
 
 

src/OpenFOAM/meshes/polyMesh/syncTools/syncToolsTemplates.C

@@ -1029,7 +1029,7 @@ void Foam::syncTools::syncBoundaryFaceList
                         pp.start()-boundaryOffset
                     );
 
-                    UIPstream::read
+                    IPstream::read
                     (
                         Pstream::commsTypes::nonBlocking,
                         procPatch.neighbProcNo(),
@@ -1055,7 +1055,7 @@ void Foam::syncTools::syncBoundaryFaceList
                         pp.start()-boundaryOffset
                     );
 
-                    UOPstream::write
+                    OPstream::write
                     (
                         Pstream::commsTypes::nonBlocking,
                         procPatch.neighbProcNo(),
@@ -1261,7 +1261,7 @@ void Foam::syncTools::syncFaceList
                 recvInfos.set(patchi, new PackedList<Width>(patchSize));
                 PackedList<Width>& recvInfo = recvInfos[patchi];
 
-                UIPstream::read
+                IPstream::read
                 (
                     Pstream::commsTypes::nonBlocking,
                     procPatch.neighbProcNo(),
@@ -1296,7 +1296,7 @@ void Foam::syncTools::syncFaceList
                 );
                 PackedList<Width>& sendInfo = sendInfos[patchi];
 
-                UOPstream::write
+                OPstream::write
                 (
                     Pstream::commsTypes::nonBlocking,
                     procPatch.neighbProcNo(),

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C

@@ -27,7 +27,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "primitiveMesh.H"
-#include "pyramid.H"
+#include "pyramidPointFaceRef.H"
 #include "ListOps.H"
 #include "unitConversion.H"
 #include "SortableList.H"

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshTools.C

@@ -29,7 +29,7 @@ License
 #include "primitiveMeshTools.H"
 #include "primitiveMesh.H"
 #include "syncTools.H"
-#include "pyramid.H"
+#include "pyramidPointFaceRef.H"
 #include "PrecisionAdaptor.H"
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //

src/OpenFOAM/meshes/primitiveShapes/cut/cut.H

@@ -33,13 +33,13 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_cut_H
-#define Foam_cut_H
+#ifndef cut_H
+#define cut_H
 
 #include "plane.H"
 #include "FixedList.H"
-#include "tetrahedron.H"
-#include "triangle.H"
+#include "tetPointRef.H"
+#include "triPointRef.H"
 #include "zero.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/OpenFOAM/meshes/primitiveShapes/line/line.H

@@ -35,13 +35,12 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_line_H
-#define Foam_line_H
+#ifndef line_H
+#define line_H
 
-#include "point.H"
-#include "point2D.H"
 #include "vector.H"
 #include "PointHit.H"
+#include "point2D.H"
 #include "FixedList.H"
 #include "UList.H"
 
@@ -61,12 +60,6 @@ template<class Point, class PointRef>
 inline Ostream& operator<<(Ostream& os, const line<Point, PointRef>& l);
 
 
-// Common Typedefs
-
-//- A line using referred points
-typedef line<point, const point&> linePointRef;
-
-
 /*---------------------------------------------------------------------------*\
                            Class line Declaration
 \*---------------------------------------------------------------------------*/
@@ -107,19 +100,19 @@ public:
     // Access
 
         //- Return first point
-        inline PointRef first() const noexcept { return a_; }
+        inline PointRef first() const noexcept;
 
         //- Return second (last) point
-        inline PointRef second() const noexcept { return b_; }
+        inline PointRef second() const noexcept;
 
         //- Return last (second) point
-        inline PointRef last() const noexcept { return b_; }
+        inline PointRef last() const noexcept;
 
-        //- Return start (first) point
-        inline PointRef start() const noexcept { return a_; }
+        //- Return first point
+        inline PointRef start() const noexcept;
 
-        //- Return end (second) point
-        inline PointRef end() const noexcept { return b_; }
+        //- Return second (last) point
+        inline PointRef end() const noexcept;
 
 
     // Properties

src/OpenFOAM/meshes/primitiveShapes/line/lineI.H

@@ -60,6 +60,40 @@ inline Foam::line<Point, PointRef>::line(Istream& is)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Point, class PointRef>
+inline PointRef Foam::line<Point, PointRef>::first() const noexcept
+{
+    return a_;
+}
+
+
+template<class Point, class PointRef>
+inline PointRef Foam::line<Point, PointRef>::second() const noexcept
+{
+    return b_;
+}
+
+
+template<class Point, class PointRef>
+inline PointRef Foam::line<Point, PointRef>::last() const noexcept
+{
+    return b_;
+}
+
+
+template<class Point, class PointRef>
+inline PointRef Foam::line<Point, PointRef>::start() const noexcept
+{
+    return first();
+}
+
+template<class Point, class PointRef>
+inline PointRef Foam::line<Point, PointRef>::end() const noexcept
+{
+    return second();
+}
+
+
 template<class Point, class PointRef>
 inline Point Foam::line<Point, PointRef>::centre() const
 {

src/OpenFOAM/meshes/primitiveShapes/line/linePoint2DRef.H

@@ -5,8 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2022 PCOpt/NTUA
-    Copyright (C) 2022      FOSS GP
+    Copyright (C) 2011 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,22 +23,33 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
+Typedef
+    Foam::linePoint2DRef
+
+Description
+    A line using referred 2D points
+
 \*---------------------------------------------------------------------------*/
 
-#include "shortGeometricField.H"
-#include "HashTable.H"
+#ifndef linePoint2DRef_H
+#define linePoint2DRef_H
+
+#include "point2D.H"
+#include "line.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-makeCompressedGeometricField(shortGeometricField)
+typedef line<point2D, const point2D&> linePoint2DRef;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
 
+#endif
+
 // ************************************************************************* //

src/OpenFOAM/meshes/primitiveShapes/line/linePointRef.H

@@ -1,11 +1,55 @@
-// Compatibility include.
-// linePointRef typedef included in line.H (JUL-2022)
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
 
-#ifndef FoamCompat_linePointRef_H
-#define FoamCompat_linePointRef_H
+    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/>.
+
+Typedef
+    Foam::linePointRef
+
+Description
+    A line using referred points
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef linePointRef_H
+#define linePointRef_H
+
+#include "point.H"
 #include "line.H"
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+typedef line<point, const point&> linePointRef;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
 #endif
 
 // ************************************************************************* //

src/OpenFOAM/meshes/primitiveShapes/point/point.H

@@ -31,8 +31,8 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_point_H
-#define Foam_point_H
+#ifndef point_H
+#define point_H
 
 #include "vector.H"
 

src/OpenFOAM/meshes/primitiveShapes/pyramid/pyramid.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2020-2022 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,11 +37,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_pyramid_H
-#define Foam_pyramid_H
-
-#include "point.H"
-#include "face.H"
+#ifndef pyramid_H
+#define pyramid_H
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,39 +47,34 @@ namespace Foam
 
 // Forward Declarations
 
-template<class Point, class PointRef, class PolygonRef> class pyramid;
+template<class Point, class PointRef, class polygonRef>
+class pyramid;
 
-template<class Point, class PointRef, class PolygonRef>
+template<class Point, class PointRef, class polygonRef>
 inline Istream& operator>>
 (
     Istream& is,
-    pyramid<Point, PointRef, PolygonRef>& p
+    pyramid<Point, PointRef, polygonRef>& p
 );
 
-template<class Point, class PointRef, class PolygonRef>
+template<class Point, class PointRef, class polygonRef>
 inline Ostream& operator<<
 (
     Ostream& os,
-    const pyramid<Point, PointRef, PolygonRef>& p
+    const pyramid<Point, PointRef, polygonRef>& p
 );
 
 
-// Common Typedefs
-
-//- A pyramid using referred point and face
-typedef pyramid<point, const point&, const face&> pyramidPointFaceRef;
-
-
 /*---------------------------------------------------------------------------*\
                            Class pyramid Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class Point, class PointRef, class PolygonRef>
+template<class Point, class PointRef, class polygonRef>
 class pyramid
 {
     // Private Data
 
-        PolygonRef base_;
+        polygonRef base_;
 
         PointRef apex_;
 
@@ -98,7 +90,7 @@ public:
     // Constructors
 
         //- Construct from base polygon and apex point
-        inline pyramid(PolygonRef base, const Point& apex);
+        inline pyramid(polygonRef base, const Point& apex);
 
         //- Construct from Istream
         inline explicit pyramid(Istream& is);
@@ -106,36 +98,36 @@ public:
 
     // Member Functions
 
-    // Access
+        // Access
 
-        //- The apex point
-        const Point& apex() const { return apex_; }
+            //- Return apex point
+            inline const Point& apex() const;
 
-        //- The base polygon
-        PolygonRef base() const { return base_; }
+            //- Return base polygon
+            inline polygonRef base() const;
 
 
-    // Properties
+        // Properties
 
-        //- Return centre (centroid)
-        inline Point centre(const UList<point>& points) const;
+            //- Return centre (centroid)
+            inline Point centre(const UList<point>& points) const;
 
-        //- Return height vector
-        inline vector height(const UList<point>& points) const;
+            //- Return height vector
+            inline vector height(const UList<point>& points) const;
 
-        //- Return scalar magnitude - returns volume of pyramid
-        inline scalar mag(const UList<point>& points) const;
+            //- Return scalar magnitude - returns volume of pyramid
+            inline scalar mag(const UList<point>& points) const;
 
 
     // IOstream Operators
 
-        friend Istream& operator>> <Point, PointRef, PolygonRef>
+        friend Istream& operator>> <Point, PointRef, polygonRef>
         (
             Istream& is,
             pyramid& p
         );
 
-        friend Ostream& operator<< <Point, PointRef, PolygonRef>
+        friend Ostream& operator<< <Point, PointRef, polygonRef>
         (
             Ostream& os,
             const pyramid& p

src/OpenFOAM/meshes/primitiveShapes/pyramid/pyramidI.H

@@ -29,10 +29,10 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class Point, class PointRef, class PolygonRef>
-inline Foam::pyramid<Point, PointRef, PolygonRef>::pyramid
+template<class Point, class PointRef, class polygonRef>
+inline Foam::pyramid<Point, PointRef, polygonRef>::pyramid
 (
-    PolygonRef base,
+    polygonRef base,
     const Point& apex
 )
 :
@@ -41,17 +41,31 @@ inline Foam::pyramid<Point, PointRef, PolygonRef>::pyramid
 {}
 
 
-template<class Point, class PointRef, class PolygonRef>
-inline Foam::pyramid<Point, PointRef, PolygonRef>::pyramid(Istream& is)
+template<class Point, class PointRef, class polygonRef>
+inline Foam::pyramid<Point, PointRef, polygonRef>::pyramid(Istream& is)
 {
-    is >> *this;
+    is >> base_ >> apex_;
+    is.check(FUNCTION_NAME);
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Point, class PointRef, class PolygonRef>
-inline Point Foam::pyramid<Point, PointRef, PolygonRef>::centre
+template<class Point, class PointRef, class polygonRef>
+inline const Point& Foam::pyramid<Point, PointRef, polygonRef>::apex() const
+{
+    return apex_;
+}
+
+template<class Point, class PointRef, class polygonRef>
+inline polygonRef Foam::pyramid<Point, PointRef, polygonRef>::base() const
+{
+    return base_;
+}
+
+
+template<class Point, class PointRef, class polygonRef>
+inline Point Foam::pyramid<Point, PointRef, polygonRef>::centre
 (
     const UList<point>& points
 ) const
@@ -60,8 +74,8 @@ inline Point Foam::pyramid<Point, PointRef, PolygonRef>::centre
 }
 
 
-template<class Point, class PointRef, class PolygonRef>
-inline Foam::vector Foam::pyramid<Point, PointRef, PolygonRef>::height
+template<class Point, class PointRef, class polygonRef>
+inline Foam::vector Foam::pyramid<Point, PointRef, polygonRef>::height
 (
     const UList<point>& points
 ) const
@@ -71,8 +85,8 @@ inline Foam::vector Foam::pyramid<Point, PointRef, PolygonRef>::height
 }
 
 
-template<class Point, class PointRef, class PolygonRef>
-inline Foam::scalar Foam::pyramid<Point, PointRef, PolygonRef>::mag
+template<class Point, class PointRef, class polygonRef>
+inline Foam::scalar Foam::pyramid<Point, PointRef, polygonRef>::mag
 (
     const UList<point>& points
 ) const
@@ -83,11 +97,11 @@ inline Foam::scalar Foam::pyramid<Point, PointRef, PolygonRef>::mag
 
 // * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * * //
 
-template<class Point, class PointRef, class PolygonRef>
+template<class Point, class PointRef, class polygonRef>
 inline Foam::Istream& Foam::operator>>
 (
     Istream& is,
-    pyramid<Point, PointRef, PolygonRef>& p
+    pyramid<Point, PointRef, polygonRef>& p
 )
 {
     is  >> p.base_ >> p.apex_;
@@ -96,11 +110,11 @@ inline Foam::Istream& Foam::operator>>
 }
 
 
-template<class Point, class PointRef, class PolygonRef>
+template<class Point, class PointRef, class polygonRef>
 inline Foam::Ostream& Foam::operator<<
 (
     Ostream& os,
-    const pyramid<Point, PointRef, PolygonRef>& p
+    const pyramid<Point, PointRef, polygonRef>& p
 )
 {
     os  << p.base_ << tab << p.apex_ << nl;

src/OpenFOAM/meshes/primitiveShapes/pyramid/pyramidPointFaceRef.H

@@ -1,11 +1,56 @@
-// Compatibility include
-// pyramidPointFaceRef typedef included in tetrahedron.H (SEP-2017)
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
 
-#ifndef FoamCompat_pyramidPointFaceRef_H
-#define FoamCompat_pyramidPointFaceRef_H
+    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/>.
+
+Typedef
+    Foam::pyramidPointFaceRef
+
+Description
+    A pyramid using referred points and faces
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pyramidPointFaceRef_H
+#define pyramidPointFaceRef_H
+
+#include "point.H"
+#include "face.H"
 #include "pyramid.H"
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+typedef pyramid<point, const point&, const face&> pyramidPointFaceRef;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
 #endif
 
 // ************************************************************************* //

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPointRef.H

@@ -1,11 +1,55 @@
-// Compatibility include.
-// tetPointRef typedef included in tetrahedron.H (SEP-2017)
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
 
-#ifndef FoamCompat_tetPointRef_H
-#define FoamCompat_tetPointRef_H
+    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/>.
+
+Typedef
+    Foam::tetPointRef
+
+Description
+    A tetrahedron using referred points
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef tetPointRef_H
+#define tetPointRef_H
+
+#include "point.H"
 #include "tetrahedron.H"
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+typedef tetrahedron<point, const point&> tetPointRef;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
 #endif
 
 // ************************************************************************* //

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetPoints.H

@@ -1,10 +1,126 @@
-// Compatibility include.
-// tetPoints defined in tetrahedron.H (JUL-2022)
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2012 OpenFOAM Foundation
+    Copyright (C) 2017-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
 
-#ifndef FoamCompat_tetPoints_H
-#define FoamCompat_tetPoints_H
+    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/>.
+
+Class
+    Foam::tetPoints
+
+Description
+    Tet storage. Null constructable (unfortunately tetrahedron<point, point>
+    is not)
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef tetPoints_H
+#define tetPoints_H
 
 #include "tetrahedron.H"
+#include "FixedList.H"
+#include "treeBoundBox.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class tetPoints Declaration
+\*---------------------------------------------------------------------------*/
+
+class tetPoints
+:
+    public FixedList<point, 4>
+{
+public:
+
+    // Constructors
+
+        //- Default construct
+        inline tetPoints()
+        {}
+
+        //- Construct from four points
+        inline tetPoints
+        (
+            const point& a,
+            const point& b,
+            const point& c,
+            const point& d
+        )
+        {
+            operator[](0) = a;
+            operator[](1) = b;
+            operator[](2) = c;
+            operator[](3) = d;
+        }
+
+        //- Copy construct from subset of points
+        inline tetPoints
+        (
+            const UList<point>& points,
+            const FixedList<label, 4>& indices
+        )
+        :
+            FixedList<point, 4>(points, indices)
+        {}
+
+
+    // Member Functions
+
+        //- Return the tetrahedron
+        inline tetPointRef tet() const
+        {
+            return tetPointRef
+            (
+                operator[](0),
+                operator[](1),
+                operator[](2),
+                operator[](3)
+            );
+        }
+
+        //- Calculate the bounding box
+        inline treeBoundBox bounds() const
+        {
+            treeBoundBox bb(operator[](0));
+            for (label i = 1; i < size(); ++i)
+            {
+                bb.add(operator[](i));
+            }
+            return bb;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
 

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.C

@@ -26,6 +26,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "tetrahedron.H"
+#include "triPointRef.H"
 #include "scalarField.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedron.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,17 +39,17 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_tetrahedron_H
-#define Foam_tetrahedron_H
+#ifndef tetrahedron_H
+#define tetrahedron_H
 
 #include "point.H"
-#include "pointHit.H"
 #include "primitiveFieldsFwd.H"
+#include "pointHit.H"
 #include "Random.H"
 #include "FixedList.H"
 #include "UList.H"
-#include "triangle.H"
-#include "treeBoundBox.H"
+#include "triPointRef.H"
+#include "boundBox.H"
 #include "barycentric.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -59,103 +59,28 @@ namespace Foam
 
 // Forward Declarations
 class plane;
+class tetPoints;
 
 template<class Point, class PointRef> class tetrahedron;
 
 template<class Point, class PointRef>
-inline Istream& operator>>(Istream&, tetrahedron<Point, PointRef>&);
+inline Istream& operator>>
+(
+    Istream&,
+    tetrahedron<Point, PointRef>&
+);
 
 template<class Point, class PointRef>
-inline Ostream& operator<<(Ostream&, const tetrahedron<Point, PointRef>&);
+inline Ostream& operator<<
+(
+    Ostream&,
+    const tetrahedron<Point, PointRef>&
+);
 
-
-// Common Typedefs
-
-//- A tetrahedron using referred points
 typedef tetrahedron<point, const point&> tetPointRef;
 
-
 /*---------------------------------------------------------------------------*\
-                         Class tetPoints Declaration
-\*---------------------------------------------------------------------------*/
-
-//- Tet point storage. Default constructable (tetrahedron is not)
-class tetPoints
-:
-    public FixedList<point, 4>
-{
-public:
-
-    // Generated Methods
-
-        //- Default construct
-        tetPoints() = default;
-
-
-    // Constructors
-
-        //- Construct from four points
-        inline tetPoints
-        (
-            const point& p0,
-            const point& p1,
-            const point& p2,
-            const point& p3
-        );
-
-        //- Construct from point references
-        inline explicit tetPoints(const tetPointRef& pts);
-
-        //- Construct from four points
-        inline tetPoints(const FixedList<point, 4>& pts);
-
-        //- Copy construct from subset of points
-        inline tetPoints
-        (
-            const UList<point>& points,
-            const FixedList<label, 4>& indices
-        );
-
-
-    // Member Functions
-
-        //- The first vertex
-        const point& a() const { return operator[](0); }
-
-        //- The second vertex
-        const point& b() const { return operator[](1); }
-
-        //- The third vertex
-        const point& c() const { return operator[](2); }
-
-        //- The fourth vertex
-        const point& d() const { return operator[](3); }
-
-        //- The first vertex
-        point& a() { return operator[](0); }
-
-        //- The second vertex
-        point& b() { return operator[](1); }
-
-        //- The third vertex
-        point& c() { return operator[](2); }
-
-        //- The fourth vertex
-        point& d() { return operator[](3); }
-
-        //- Return as tetrahedron reference
-        inline tetPointRef tet() const;
-
-        //- Invert tetrahedron by swapping third and fourth vertices
-        inline void flip();
-
-        //- The bounding box for the tetrahedron
-        inline treeBoundBox bounds() const;
-};
-
-
-/*---------------------------------------------------------------------------*\
-                         Class tetrahedron Declaration
+                           class tetrahedron Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class Point, class PointRef>
@@ -163,10 +88,7 @@ class tetrahedron
 {
 public:
 
-    // Public Typedefs
-
-        //- The point type
-        typedef Point point_type;
+    // Public typedefs
 
         //- Storage type for tets originating from intersecting tets.
         //  (can possibly be smaller than 200)
@@ -208,13 +130,10 @@ public:
 
 private:
 
-    // Private Data
+    // Private data
 
         PointRef a_, b_, c_, d_;
 
-
-    // Private Member Functions
-
         inline static point planeIntersection
         (
             const FixedList<scalar, 4>&,
@@ -242,7 +161,7 @@ private:
 
 public:
 
-    // Member Constants
+    // Member constants
 
         enum
         {
@@ -253,18 +172,15 @@ public:
 
     // Constructors
 
-        //- Construct from four points
+        //- Construct from points
         inline tetrahedron
         (
-            const Point& p0,
-            const Point& p1,
-            const Point& p2,
-            const Point& p3
+            const Point& a,
+            const Point& b,
+            const Point& c,
+            const Point& d
         );
 
-        //- Construct from four points
-        inline tetrahedron(const FixedList<Point, 4>& pts);
-
         //- Construct from four points in the list of points
         inline tetrahedron
         (
@@ -273,27 +189,24 @@ public:
         );
 
         //- Construct from Istream
-        inline explicit tetrahedron(Istream&);
+        inline tetrahedron(Istream&);
 
 
     // Member Functions
 
-    // Access
+        // Access
 
-        //- Return vertex a
-        const Point& a() const noexcept { return a_; }
+            //- Return vertices
+            inline const Point& a() const;
 
-        //- Return vertex b
-        const Point& b() const noexcept { return b_; }
+            inline const Point& b() const;
 
-        //- Return vertex c
-        const Point& c() const noexcept { return c_; }
+            inline const Point& c() const;
 
-        //- Return vertex d
-        const Point& d() const noexcept { return d_; }
+            inline const Point& d() const;
 
-        //- Return i-th face
-        inline triPointRef tri(const label facei) const;
+            //- Return i-th face
+            inline triPointRef tri(const label facei) const;
 
 
         // Properties
@@ -323,12 +236,12 @@ public:
             inline scalar circumRadius() const;
 
             //- Return quality: Ratio of tetrahedron and circum-sphere
-            //- volume, scaled so that a regular tetrahedron has a
-            //- quality of 1
+            //  volume, scaled so that a regular tetrahedron has a
+            //  quality of 1
             inline scalar quality() const;
 
             //- Return a random point in the tetrahedron from a
-            //- uniform distribution
+            //  uniform distribution
             inline Point randomPoint(Random& rndGen) const;
 
             //- Calculate the point from the given barycentric coordinates.

src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,76 +28,24 @@ License
 
 #include "triangle.H"
 #include "IOstreams.H"
+#include "tetPoints.H"
 #include "plane.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-inline Foam::tetPoints::tetPoints
-(
-    const point& p0,
-    const point& p1,
-    const point& p2,
-    const point& p3
-)
-{
-    operator[](0) = p0;
-    operator[](1) = p1;
-    operator[](2) = p2;
-    operator[](3) = p3;
-}
-
-
-inline Foam::tetPoints::tetPoints(const tetPointRef& pts)
-{
-    operator[](0) = pts.a();
-    operator[](1) = pts.b();
-    operator[](2) = pts.c();
-    operator[](3) = pts.d();
-}
-
-
-inline Foam::tetPoints::tetPoints(const FixedList<point, 4>& pts)
-:
-    FixedList<point, 4>(pts)
-{}
-
-
-inline Foam::tetPoints::tetPoints
-(
-    const UList<point>& points,
-    const FixedList<label, 4>& indices
-)
-:
-    FixedList<point, 4>(points, indices)
-{}
-
-
 template<class Point, class PointRef>
 inline Foam::tetrahedron<Point, PointRef>::tetrahedron
 (
-    const Point& p0,
-    const Point& p1,
-    const Point& p2,
-    const Point& p3
+    const Point& a,
+    const Point& b,
+    const Point& c,
+    const Point& d
 )
 :
-    a_(p0),
-    b_(p1),
-    c_(p2),
-    d_(p3)
-{}
-
-
-template<class Point, class PointRef>
-inline Foam::tetrahedron<Point, PointRef>::tetrahedron
-(
-    const FixedList<Point, 4>& pts
-)
-:
-    a_(pts[0]),
-    b_(pts[1]),
-    c_(pts[2]),
-    d_(pts[3])
+    a_(a),
+    b_(b),
+    c_(c),
+    d_(d)
 {}
 
 
@@ -124,31 +72,33 @@ inline Foam::tetrahedron<Point, PointRef>::tetrahedron(Istream& is)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::tetPointRef Foam::tetPoints::tet() const
+template<class Point, class PointRef>
+inline const Point& Foam::tetrahedron<Point, PointRef>::a() const
 {
-    return tetPointRef(a(), b(), c(), d());
+    return a_;
 }
 
 
-inline void Foam::tetPoints::flip()
+template<class Point, class PointRef>
+inline const Point& Foam::tetrahedron<Point, PointRef>::b() const
 {
-    // swap pt2 <-> pt3
-    point t(c());
-    c() = d();
-    d() = t;
+    return b_;
 }
 
 
-inline Foam::treeBoundBox Foam::tetPoints::bounds() const
+template<class Point, class PointRef>
+inline const Point& Foam::tetrahedron<Point, PointRef>::c() const
 {
-    treeBoundBox bb;
-    bb.add(static_cast<const FixedList<point, 4>&>(*this));
-    return bb;
+    return c_;
 }
 
 
-// Warning. Ordering of faces needs to be the same for a tetrahedron class,
-// tetrahedron cell shape model and a tetCell
+template<class Point, class PointRef>
+inline const Point& Foam::tetrahedron<Point, PointRef>::d() const
+{
+    return d_;
+}
+
 
 template<class Point, class PointRef>
 inline Foam::triPointRef Foam::tetrahedron<Point, PointRef>::tri
@@ -307,7 +257,7 @@ inline Point Foam::tetrahedron<Point, PointRef>::barycentricToPoint
     const barycentric& bary
 ) const
 {
-    return Point(bary.a()*a_ + bary.b()*b_ + bary.c()*c_ + bary.d()*d_);
+    return bary[0]*a_ + bary[1]*b_ + bary[2]*c_ + bary[3]*d_;
 }
 
 
@@ -360,7 +310,7 @@ inline Foam::scalar Foam::tetrahedron<Point, PointRef>::pointToBarycentric
     bary[0] = res.x();
     bary[1] = res.y();
     bary[2] = res.z();
-    bary[3] = 1 - bary[0] - bary[1] - bary[2];
+    bary[3] = 1 - cmptSum(res);
 
     return detT;
 }
@@ -384,83 +334,63 @@ inline Foam::pointHit Foam::tetrahedron<Point, PointRef>::nearestPoint
 
     bool inside = true;
 
-    // Side a
+    if (((p - b_) & Sa()) >= 0)
     {
-        const triangle<Point, PointRef> tria(b_, c_, d_);
+        // p is outside halfspace plane of tri
+        pointHit info = triangle<Point, PointRef>(b_, c_, d_).nearestPoint(p);
 
-        if (((p - b_) & tria.areaNormal()) >= 0)
+        inside = false;
+
+        if (info.distance() < minOutsideDistance)
         {
-            // p is outside halfspace plane of tri
-            pointHit info = tria.nearestPoint(p);
+            closestPt = info.rawPoint();
 
-            inside = false;
-
-            if (info.distance() < minOutsideDistance)
-            {
-                closestPt = info.rawPoint();
-
-                minOutsideDistance = info.distance();
-            }
+            minOutsideDistance = info.distance();
         }
     }
 
-    // Side b
+    if (((p - a_) & Sb()) >= 0)
     {
-        const triangle<Point, PointRef> tria(a_, d_, c_);
+        // p is outside halfspace plane of tri
+        pointHit info = triangle<Point, PointRef>(a_, d_, c_).nearestPoint(p);
 
-        if (((p - a_) & tria.areaNormal()) >= 0)
+        inside = false;
+
+        if (info.distance() < minOutsideDistance)
         {
-            // p is outside halfspace plane of tri
-            pointHit info = tria.nearestPoint(p);
+            closestPt = info.rawPoint();
 
-            inside = false;
-
-            if (info.distance() < minOutsideDistance)
-            {
-                closestPt = info.rawPoint();
-
-                minOutsideDistance = info.distance();
-            }
+            minOutsideDistance = info.distance();
         }
     }
 
-    // Side c
+    if (((p - a_) & Sc()) >= 0)
     {
-        const triangle<Point, PointRef> tria(a_, b_, d_);
+        // p is outside halfspace plane of tri
+        pointHit info = triangle<Point, PointRef>(a_, b_, d_).nearestPoint(p);
 
-        if (((p - a_) & tria.areaNormal()) >= 0)
+        inside = false;
+
+        if (info.distance() < minOutsideDistance)
         {
-            // p is outside halfspace plane of tri
-            pointHit info = tria.nearestPoint(p);
+            closestPt = info.rawPoint();
 
-            inside = false;
-
-            if (info.distance() < minOutsideDistance)
-            {
-                closestPt = info.rawPoint();
-
-                minOutsideDistance = info.distance();
-            }
+            minOutsideDistance = info.distance();
         }
     }
 
-    // Side c
+    if (((p - a_) & Sd()) >= 0)
     {
-        const triangle<Point, PointRef> tria(a_, c_, b_);
+        // p is outside halfspace plane of tri
+        pointHit info = triangle<Point, PointRef>(a_, c_, b_).nearestPoint(p);
 
-        if (((p - a_) & tria.areaNormal()) >= 0)
+        inside = false;
+
+        if (info.distance() < minOutsideDistance)
         {
-            // p is outside halfspace plane of tri
-            pointHit info = tria.nearestPoint(p);
+            closestPt = info.rawPoint();
 
-            inside = false;
-
-            if (info.distance() < minOutsideDistance)
-            {
-                closestPt = info.rawPoint();
-
-                minOutsideDistance = info.distance();
-            }
+            minOutsideDistance = info.distance();
         }
     }
 
@@ -482,7 +412,7 @@ inline Foam::pointHit Foam::tetrahedron<Point, PointRef>::nearestPoint
 
 
 template<class Point, class PointRef>
-bool Foam::tetrahedron<Point, PointRef>::inside(const point& p) const
+bool Foam::tetrahedron<Point, PointRef>::inside(const point& pt) const
 {
     // For robustness, assuming that the point is in the tet unless
     // "definitively" shown otherwise by obtaining a positive dot
@@ -499,41 +429,55 @@ bool Foam::tetrahedron<Point, PointRef>::inside(const point& p) const
     // planeBase[2] = tetBasePt = b_
     // planeBase[3] = tetBasePt = b_
 
-    // Side a
-    {
-        const triangle<Point, PointRef> tria(b_, c_, d_);
+    vector n = Zero;
 
-        if (((p - b_) & tria.unitNormal()) > SMALL)
+    {
+        // 0, a
+        const point& basePt = b_;
+
+        n = Sa();
+        n /= (Foam::mag(n) + VSMALL);
+
+        if (((pt - basePt) & n) > SMALL)
         {
             return false;
         }
     }
 
-    // Side b
     {
-        const triangle<Point, PointRef> tria(a_, d_, c_);
+        // 1, b
+        const point& basePt = c_;
 
-        if (((p - a_) & tria.unitNormal()) > SMALL)
+        n = Sb();
+        n /= (Foam::mag(n) + VSMALL);
+
+        if (((pt - basePt) & n) > SMALL)
         {
             return false;
         }
     }
 
-    // Side c
     {
-        const triangle<Point, PointRef> tria(a_, b_, d_);
+        // 2, c
+        const point& basePt = b_;
 
-        if (((p - a_) & tria.unitNormal()) > SMALL)
+        n = Sc();
+        n /= (Foam::mag(n) + VSMALL);
+
+        if (((pt - basePt) & n) > SMALL)
         {
             return false;
         }
     }
 
-    // Side d
     {
-        const triangle<Point, PointRef> tria(a_, c_, b_);
+        // 3, d
+        const point& basePt = b_;
 
-        if (((p - a_) & tria.unitNormal()) > SMALL)
+        n = Sd();
+        n /= (Foam::mag(n) + VSMALL);
+
+        if (((pt - basePt) & n) > SMALL)
         {
             return false;
         }
@@ -1111,7 +1055,6 @@ inline Foam::Ostream& Foam::operator<<
     const tetrahedron<Point, PointRef>& t
 )
 {
-    // Format like FixedList
     os  << nl
         << token::BEGIN_LIST
         << t.a_ << token::SPACE

src/OpenFOAM/meshes/primitiveShapes/triangle/triPointRef.H

@@ -1,11 +1,55 @@
-// Compatibility include.
-// triPointRef typedef included in triangle.H (NOV-2015)
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
 
-#ifndef FoamCompat_triPointRef_H
-#define FoamCompat_triPointRef_H
+    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/>.
+
+Typedef
+    Foam::triPointRef
+
+Description
+    A triangle using referred points
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef triPointRef_H
+#define triPointRef_H
+
+#include "point.H"
 #include "triangle.H"
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+typedef triangle<point, const point&> triPointRef;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
 #endif
 
 // ************************************************************************* //

src/OpenFOAM/meshes/primitiveShapes/triangle/triPoints.H

@@ -1,10 +1,111 @@
-// Compatibility include.
-// triPoints defined in triangle.H (JUL-2022)
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2015 OpenFOAM Foundation
+    Copyright (C) 2017-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
 
-#ifndef FoamCompat_triPoints_H
-#define FoamCompat_triPoints_H
+    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.
 
-#include "triangle.H"
+    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/>.
+
+Class
+    Foam::triPoints
+
+Description
+    Triangle storage. Null constructable (unfortunately triangle<point, point>
+    is not)
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef triPoints_H
+#define triPoints_H
+
+#include "FixedList.H"
+#include "treeBoundBox.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class triPoints Declaration
+\*---------------------------------------------------------------------------*/
+
+class triPoints
+:
+    public FixedList<point, 3>
+{
+public:
+
+    // Constructors
+
+        //- Default construct
+        inline triPoints()
+        {}
+
+        //- Construct from points
+        inline triPoints
+        (
+            const point& a,
+            const point& b,
+            const point& c
+        )
+        {
+            operator[](0) = a;
+            operator[](1) = b;
+            operator[](2) = c;
+        }
+
+        //- Copy construct from subset of points
+        inline triPoints
+        (
+            const UList<point>& points,
+            const FixedList<label, 3>& indices
+        )
+        :
+            FixedList<point, 3>(points, indices)
+        {}
+
+
+    // Member Functions
+
+        //- Calculate the bounding box
+        inline treeBoundBox bounds() const
+        {
+            treeBoundBox bb(operator[](0));
+            for (label i = 1; i < size(); ++i)
+            {
+                bb.add(operator[](i));
+            }
+            return bb;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
 

src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C

@@ -26,6 +26,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "triangle.H"
+#include "triPoints.H"
 #include "plane.H"
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,20 +35,18 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_triangle_H
-#define Foam_triangle_H
+#ifndef triangle_H
+#define triangle_H
 
 #include "intersection.H"
-#include "point.H"
 #include "vector.H"
 #include "tensor.H"
 #include "pointHit.H"
 #include "Random.H"
 #include "FixedList.H"
 #include "UList.H"
-#include "line.H"
+#include "linePointRef.H"
 #include "barycentric2D.H"
-#include "treeBoundBox.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -57,6 +55,7 @@ namespace Foam
 
 // Forward Declarations
 class plane;
+class triPoints;
 
 template<class Point, class PointRef> class triangle;
 
@@ -73,73 +72,6 @@ inline Ostream& operator<<(Ostream&, const triangle<Point, PointRef>&);
 typedef triangle<point, const point&> triPointRef;
 
 
-/*---------------------------------------------------------------------------*\
-                          Class triPoints Declaration
-\*---------------------------------------------------------------------------*/
-
-//- Triangle point storage. Default constructable (triangle is not)
-class triPoints
-:
-    public FixedList<point, 3>
-{
-public:
-
-    // Generated Methods
-
-        //- Default construct
-        triPoints() = default;
-
-
-    // Constructors
-
-        //- Construct from three points
-        inline triPoints(const point& p0, const point& p1, const point& p2);
-
-        //- Construct from point references
-        inline explicit triPoints(const triPointRef& pts);
-
-        //- Construct from three points
-        inline triPoints(const FixedList<point, 3>& pts);
-
-        //- Copy construct from subset of points
-        inline triPoints
-        (
-            const UList<point>& points,
-            const FixedList<label, 3>& indices
-        );
-
-
-    // Member Functions
-
-        //- The first vertex
-        const point& a() const { return operator[](0); }
-
-        //- The second vertex
-        const point& b() const { return operator[](1); }
-
-        //- The third vertex
-        const point& c() const { return operator[](2); }
-
-        //- The first vertex
-        point& a() { return operator[](0); }
-
-        //- The second vertex
-        point& b() { return operator[](1); }
-
-        //- The third vertex
-        point& c() { return operator[](2); }
-
-        //- Return as triangle reference
-        inline triPointRef tri() const;
-
-        //- Flip triangle orientation by swapping second and third vertices
-        inline void flip();
-
-        //- The bounding box for the tetrahedron
-        inline treeBoundBox bounds() const;
-};
-
-
 /*---------------------------------------------------------------------------*\
                           Class triangle Declaration
 \*---------------------------------------------------------------------------*/
@@ -158,7 +90,7 @@ public:
         //- with another triangle
         typedef FixedList<triPoints, 27> triIntersectionList;
 
-        //- Proximity classifications
+        //- The proximity classifications
         enum proxType
         {
             NONE = 0,   //!< Unknown proximity
@@ -167,10 +99,10 @@ public:
         };
 
 
-    // Public Classes
+    // Public classes
 
-        // Classes for use in sliceWithPlane.
-        // What to do with decomposition of triangle.
+        // Classes for use in sliceWithPlane. What to do with decomposition
+        // of triangle.
 
             //- Dummy
             class dummyOp
@@ -237,10 +169,10 @@ public:
     // Constructors
 
         //- Construct from three points
-        inline triangle(const Point& p0, const Point& p1, const Point& p2);
+        inline triangle(const Point& a, const Point& b, const Point& c);
 
         //- Construct from three points
-        inline triangle(const FixedList<Point, 3>& pts);
+        inline triangle(const FixedList<Point, 3>& tri);
 
         //- Construct from three points in the list of points
         //  The indices could be from triFace etc.
@@ -256,16 +188,16 @@ public:
 
     // Member Functions
 
-    // Access
+        // Access
 
-        //- The first vertex
-        const Point& a() const noexcept { return a_; }
+            //- Return first vertex
+            inline const Point& a() const;
 
-        //- The second vertex
-        const Point& b() const noexcept { return b_; }
+            //- Return second vertex
+            inline const Point& b() const;
 
-        //- The third vertex
-        const Point& c() const noexcept { return c_; }
+            //- Return third vertex
+            inline const Point& c() const;
 
 
         // Properties
@@ -313,7 +245,7 @@ public:
             ) const;
 
             //- Return a random point on the triangle from a uniform
-            //- distribution
+            //  distribution
             inline Point randomPoint(Random& rndGen) const;
 
             //- Calculate the point from the given barycentric coordinates.

src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,70 +28,34 @@ License
 
 #include "IOstreams.H"
 #include "pointHit.H"
+#include "triPoints.H"
 #include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
 
-inline Foam::triPoints::triPoints
+template<class Point, class PointRef>
+inline Foam::triangle<Point, PointRef>::triangle
 (
-    const point& p0,
-    const point& p1,
-    const point& p2
-)
-{
-    operator[](0) = p0;
-    operator[](1) = p1;
-    operator[](2) = p2;
-}
-
-
-inline Foam::triPoints::triPoints(const triPointRef& pts)
-{
-    operator[](0) = pts.a();
-    operator[](1) = pts.b();
-    operator[](2) = pts.c();
-}
-
-
-inline Foam::triPoints::triPoints(const FixedList<point, 3>& pts)
-:
-    FixedList<point, 3>(pts)
-{}
-
-
-inline Foam::triPoints::triPoints
-(
-    const UList<point>& points,
-    const FixedList<label, 3>& indices
+    const Point& a,
+    const Point& b,
+    const Point& c
 )
 :
-    FixedList<point, 3>(points, indices)
+    a_(a),
+    b_(b),
+    c_(c)
 {}
 
 
 template<class Point, class PointRef>
 inline Foam::triangle<Point, PointRef>::triangle
 (
-    const Point& p0,
-    const Point& p1,
-    const Point& p2
+    const FixedList<Point, 3>& tri
 )
 :
-    a_(p0),
-    b_(p1),
-    c_(p2)
-{}
-
-
-template<class Point, class PointRef>
-inline Foam::triangle<Point, PointRef>::triangle
-(
-    const FixedList<Point, 3>& pts
-)
-:
-    a_(pts[0]),
-    b_(pts[1]),
-    c_(pts[2])
+    a_(tri[0]),
+    b_(tri[1]),
+    c_(tri[2])
 {}
 
 
@@ -108,6 +72,7 @@ inline Foam::triangle<Point, PointRef>::triangle
 {}
 
 
+
 template<class Point, class PointRef>
 inline Foam::triangle<Point, PointRef>::triangle(Istream& is)
 {
@@ -117,26 +82,22 @@ inline Foam::triangle<Point, PointRef>::triangle(Istream& is)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::triPointRef Foam::triPoints::tri() const
+template<class Point, class PointRef>
+inline const Point& Foam::triangle<Point, PointRef>::a() const
 {
-    return triPointRef(a(), b(), c());
+    return a_;
 }
 
-
-inline void Foam::triPoints::flip()
+template<class Point, class PointRef>
+inline const Point& Foam::triangle<Point, PointRef>::b() const
 {
-    // swap pt1 <-> pt2
-    point t(b());
-    b() = c();
-    c() = t;
+    return b_;
 }
 
-
-inline Foam::treeBoundBox Foam::triPoints::bounds() const
+template<class Point, class PointRef>
+inline const Point& Foam::triangle<Point, PointRef>::c() const
 {
-    treeBoundBox bb;
-    bb.add(static_cast<const FixedList<point, 3>&>(*this));
-    return bb;
+    return c_;
 }
 
 
@@ -955,7 +916,6 @@ inline Foam::Ostream& Foam::operator<<
     const triangle<Point, PointRef>& t
 )
 {
-    // Format like FixedList
     os  << nl
         << token::BEGIN_LIST
         << t.a_ << token::SPACE

src/OpenFOAM/parallel/globalIndex/globalIndexTemplates.C

@@ -100,7 +100,7 @@ void Foam::globalIndex::gatherValues
         {
             if (is_contiguous<Type>::value)
             {
-                UIPstream::read
+                IPstream::read
                 (
                     commsType,
                     procIDs[i],
@@ -123,7 +123,7 @@ void Foam::globalIndex::gatherValues
 
         if (is_contiguous<Type>::value)
         {
-            UOPstream::write
+            OPstream::write
             (
                 commsType,
                 procIDs[0],
@@ -198,7 +198,7 @@ void Foam::globalIndex::gather
             }
             else if (is_contiguous<Type>::value)
             {
-                UIPstream::read
+                IPstream::read
                 (
                     commsType,
                     procIDs[i],
@@ -223,7 +223,7 @@ void Foam::globalIndex::gather
         }
         else if (is_contiguous<Type>::value)
         {
-            UOPstream::write
+            OPstream::write
             (
                 commsType,
                 procIDs[0],
@@ -889,7 +889,7 @@ void Foam::globalIndex::scatter
             }
             else if (is_contiguous<Type>::value)
             {
-                UOPstream::write
+                OPstream::write
                 (
                     commsType,
                     procIDs[i],
@@ -926,7 +926,7 @@ void Foam::globalIndex::scatter
         }
         else if (is_contiguous<Type>::value)
         {
-            UIPstream::read
+            IPstream::read
             (
                 commsType,
                 procIDs[0],

src/OpenFOAM/primitives/Barycentric/Barycentric.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 OpenFOAM Foundation
-    Copyright (C) 2019-2022 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,16 +28,16 @@ Class
     Foam::Barycentric
 
 Description
-    Templated 3D Barycentric derived from VectorSpace.
-    Has 4 components, one of which is redundant.
+    Templated 3D Barycentric derived from VectorSpace. Has 4 components, one of
+    which is redundant.
 
 SourceFiles
     BarycentricI.H
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Foam_Barycentric_H
-#define Foam_Barycentric_H
+#ifndef Barycentric_H
+#define Barycentric_H
 
 #include "contiguous.H"
 #include "VectorSpace.H"
@@ -74,7 +74,7 @@ public:
     enum components { A, B, C, D };
 
 
-    // Generated Methods: copy construct/assignment
+    // Generated Methods
 
         //- Default construct
         Barycentric() = default;
@@ -94,29 +94,20 @@ public:
             const Cmpt& vd
         );
 
-        //- Construct from three components, calculate fourth component
-        inline Barycentric(const Cmpt& va, const Cmpt& vb, const Cmpt& vc);
-
 
     // Member Functions
 
-    // Component access
+        // Access
 
-        inline const Cmpt& a() const;
-        inline const Cmpt& b() const;
-        inline const Cmpt& c() const;
-        inline const Cmpt& d() const;
+            inline const Cmpt& a() const;
+            inline const Cmpt& b() const;
+            inline const Cmpt& c() const;
+            inline const Cmpt& d() const;
 
-        inline Cmpt& a();
-        inline Cmpt& b();
-        inline Cmpt& c();
-        inline Cmpt& d();
-
-
-    // Operations
-
-        //- Scalar-product of \c this with another Barycentric.
-        inline Cmpt inner(const Barycentric<Cmpt>& b2) const;
+            inline Cmpt& a();
+            inline Cmpt& b();
+            inline Cmpt& c();
+            inline Cmpt& d();
 };
 
 

src/OpenFOAM/primitives/Barycentric/BarycentricI.H

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016-2017 OpenFOAM Foundation
-    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -51,18 +50,6 @@ inline Foam::Barycentric<Cmpt>::Barycentric
 }
 
 
-template<class Cmpt>
-inline Foam::Barycentric<Cmpt>::Barycentric
-(
-    const Cmpt& va,
-    const Cmpt& vb,
-    const Cmpt& vc
-)
-:
-    Barycentric<Cmpt>(va, vb, vc, (pTraits<Cmpt>::one - va - vb - vc))
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class Cmpt>
@@ -121,17 +108,6 @@ inline Cmpt& Foam::Barycentric<Cmpt>::d()
 }
 
 
-// * * * * * * * * * * * * * * * Member Operations * * * * * * * * * * * * * //
-
-template<class Cmpt>
-inline Cmpt Foam::Barycentric<Cmpt>::inner(const Barycentric<Cmpt>& b2) const
-{
-    const Barycentric<Cmpt>& b1 = *this;
-
-    return (b1.a()*b2.a() + b1.b()*b2.b() + b1.c()*b2.c() + b1.d()*b2.d());
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -142,7 +118,7 @@ namespace Foam
 template<class Cmpt>
 inline Cmpt operator&(const Barycentric<Cmpt>& b1, const Barycentric<Cmpt>& b2)
 {
-    return b1.inner(b2);
+    return b1.a()*b2.a() + b1.b()*b2.b() + b1.c()*b2.c() + b1.d()*b2.d();
 }