TUT: added two unsteady adjoint optimisation tutorials

showcasing the use of compressedFullStorage and binomialCheckPointing

CONTRIBUTORS.md

@@ -5,6 +5,8 @@ It is likely incomplete...
 
 ## Contributors (alphabetical by surname)
 
+- Tetsuo Aoyagi
+- Akira Azami
 - William Bainbridge
 - Gabriel Barajas
 - Kutalmis Bercin
@@ -19,6 +21,7 @@ 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 DTRMParticle_H
-#define DTRMParticle_H
+#ifndef Foam_DTRMParticle_H
+#define Foam_DTRMParticle_H
 
 #include "particle.H"
 #include "IOstream.H"
@@ -50,12 +50,13 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class DTRMParticle;
-using namespace Foam::radiation;
 
-// Forward declaration of friend functions
 Ostream& operator<<(Ostream&, const DTRMParticle&);
 
+using namespace Foam::radiation;
+
 /*---------------------------------------------------------------------------*\
                         Class DTRMParticle Declaration
 \*---------------------------------------------------------------------------*/
@@ -130,7 +131,7 @@ public:
             );
 
 
-        // Member functions
+        // Member Functions
 
             inline const interpolationCell<scalar>& aInterp() const;
             inline const interpolationCell<scalar>& eInterp() const;
@@ -232,37 +233,34 @@ public:
     // Access
 
         //- Return const access to the initial position
-        inline const point& p0() const;
+        const point& p0() const noexcept { return p0_; }
 
         //- Return const access to the target position
-        inline const point& p1() const;
+        const point& p1() const noexcept { return p1_; }
 
         //- Return const access to the initial intensity
-        inline scalar I0() const;
+        scalar I0() const noexcept { return I0_; }
 
         //- Return const access to the current intensity
-        inline scalar I() const;
+        scalar I() const noexcept { return I_; }
 
         //- Return const access dA
-        inline scalar dA() const;
+        scalar dA() const noexcept { return dA_; }
 
 
     // Edit
 
         //- Return access to the target position
-        inline point& p1();
+        point& p1() noexcept { return p1_; }
 
         //- Return access to the initial intensity
-        inline scalar& I0();
+        scalar& I0() noexcept { return I0_; }
 
         //- Return access to the current intensity
-        inline scalar& I();
+        scalar& I() noexcept { return I_; }
 
         //- Return access to dA
-        inline scalar& dA();
-
-        //- Return access to reflectedId
-        inline label& reflectedId();
+        scalar& dA() noexcept { return dA_; }
 
 
    // Tracking

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

@@ -107,58 +107,4 @@ 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

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

applications/test/barycentric/Make/options

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

applications/test/barycentric/Test-barycentric.C

@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "tetPointRef.H"
+#include "tetrahedron.H"
 #include "triFaceList.H"
 #include "OFstream.H"
 #include "meshTools.H"

applications/test/tetTetOverlap/Test-tetTetOverlap.C

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

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

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

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

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

applications/utilities/parallelProcessing/redistributePar/passivePositionParticle.H

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

etc/caseDicts/postProcessing/visualization/streamlines.cfg

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

src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.C

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

src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H

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

src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicTreeDataPoint.H

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

src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C

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

src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H

@@ -34,8 +34,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef indexedOctree_H
-#define indexedOctree_H
+#ifndef Foam_indexedOctree_H
+#define Foam_indexedOctree_H
 
 #include "treeBoundBox.H"
 #include "pointIndexHit.H"
@@ -51,7 +51,7 @@ SourceFiles
 namespace Foam
 {
 
-// Forward declaration of classes
+// Forward Declarations
 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]);
-                IPstream::read
+                UIPstream::read
                 (
                     UPstream::commsTypes::scheduled,
                     proci,

src/OpenFOAM/db/Time/Time.C

@@ -967,6 +967,34 @@ 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)
@@ -1365,4 +1393,242 @@ 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,6 +7,7 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
     Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2022      PCOpt/NTUA
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -571,6 +572,43 @@ 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
 
@@ -653,6 +691,20 @@ 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,7 +545,8 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 (
     const IOobject& io,
     const Mesh& mesh,
-    const dictionary& dict
+    const dictionary& dict,
+    const bool readOldTime
 )
 :
     Internal(io, mesh, dimless, false),
@@ -566,6 +567,11 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
             << exit(FatalIOError);
     }
 
+    if (readOldTime)
+    {
+        readOldTimeIfPresent();
+    }
+
     DebugInFunction
         << "Finishing dictionary-construct" << nl << this->info() << endl;
 }
@@ -623,7 +629,8 @@ 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 GeometricField<Type, PatchField, GeoMesh>& gf,
+    const bool readOldTime
 )
 :
     Internal(io, gf),
@@ -636,7 +643,7 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
         << "Copy construct, resetting IO params" << nl
         << this->info() << endl;
 
-    if (!readIfPresent() && gf.field0Ptr_)
+    if (!readIfPresent() && gf.field0Ptr_ && readOldTime)
     {
         field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
         (

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

@@ -296,7 +296,8 @@ public:
         (
             const IOobject& io,
             const Mesh& mesh,
-            const dictionary& dict
+            const dictionary& dict,
+            const bool readOldTime = false
         );
 
         //- Copy construct
@@ -315,7 +316,8 @@ public:
         GeometricField
         (
             const IOobject& io,
-            const GeometricField<Type, PatchField, GeoMesh>& gf
+            const GeometricField<Type, PatchField, GeoMesh>& gf,
+            const bool readOldTime = true
         );
 
         //- 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());
-            IPstream::read
+            UIPstream::read
             (
                 commsType,
                 procPatch_.neighbProcNo(),
@@ -129,7 +129,7 @@ void Foam::processorCyclicPointPatchField<Type>::initSwapAddSeparated
                 procPatch_.comm()
             );
         }
-        OPstream::write
+        UOPstream::write
         (
             commsType,
             procPatch_.neighbProcNo(),
@@ -155,7 +155,7 @@ void Foam::processorCyclicPointPatchField<Type>::swapAddSeparated
         if (commsType != Pstream::commsTypes::nonBlocking)
         {
             receiveBuf_.setSize(this->size());
-            IPstream::read
+            UIPstream::read
             (
                 commsType,
                 procPatch_.neighbProcNo(),

src/OpenFOAM/matrices/LUscalarMatrix/LUscalarMatrix.C

@@ -86,7 +86,7 @@ Foam::LUscalarMatrix::LUscalarMatrix
 
         if (Pstream::master(comm_))
         {
-            for (const int slave : Pstream::subProcs(comm_))
+            for (const int proci : Pstream::subProcs(comm_))
             {
                 lduMatrices.set
                 (
@@ -96,7 +96,7 @@ Foam::LUscalarMatrix::LUscalarMatrix
                         IPstream
                         (
                             Pstream::commsTypes::scheduled,
-                            slave,
+                            proci,
                             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 slave : Pstream::subProcs(comm_))
+            for (const int proci : Pstream::subProcs(comm_))
             {
-                IPstream::read
+                UIPstream::read
                 (
                     Pstream::commsTypes::scheduled,
-                    slave,
+                    proci,
                     reinterpret_cast<char*>
                     (
-                        &(X[procOffsets_[slave]])
+                        &(X[procOffsets_[proci]])
                     ),
-                    (procOffsets_[slave+1]-procOffsets_[slave])*sizeof(Type),
+                    (procOffsets_[proci+1]-procOffsets_[proci])*sizeof(Type),
                     Pstream::msgType(),
                     comm_
                 );
@@ -72,7 +72,7 @@ void Foam::LUscalarMatrix::solve
         }
         else
         {
-            OPstream::write
+            UOPstream::write
             (
                 Pstream::commsTypes::scheduled,
                 Pstream::masterNo(),
@@ -89,17 +89,17 @@ void Foam::LUscalarMatrix::solve
 
             x = SubList<Type>(X, x.size());
 
-            for (const int slave : Pstream::subProcs(comm_))
+            for (const int proci : Pstream::subProcs(comm_))
             {
-                OPstream::write
+                UOPstream::write
                 (
                     Pstream::commsTypes::scheduled,
-                    slave,
+                    proci,
                     reinterpret_cast<const char*>
                     (
-                        &(X[procOffsets_[slave]])
+                        &(X[procOffsets_[proci]])
                     ),
-                    (procOffsets_[slave+1]-procOffsets_[slave])*sizeof(Type),
+                    (procOffsets_[proci+1]-procOffsets_[proci])*sizeof(Type),
                     Pstream::msgType(),
                     comm_
                 );
@@ -107,7 +107,7 @@ void Foam::LUscalarMatrix::solve
         }
         else
         {
-            IPstream::read
+            UIPstream::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-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -54,6 +54,7 @@ SourceFiles
 #define Foam_LduMatrix_H
 
 #include "lduMesh.H"
+#include "lduMatrix.H"
 #include "Field.H"
 #include "FieldField.H"
 #include "LduInterfaceFieldPtrsList.H"
@@ -69,8 +70,7 @@ 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,16 +119,13 @@ 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 controls
+            //- Dictionary of solution controls
             dictionary controlDict_;
 
-            //- Level of verbosity in the solver output statements
+            //- Verbosity level for solver output statements
             int log_;
 
             //- Minimum number of iterations in the solver
@@ -137,6 +134,9 @@ 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 the controlDict_
+            //- Read the control parameters from controlDict_
             virtual void readControls();
 
 
@@ -206,6 +206,7 @@ public:
 
         // Constructors
 
+            //- Construct for given field name, matrix and controls
             solver
             (
                 const word& fieldName,
@@ -225,9 +226,8 @@ public:
             );
 
 
-        // Destructor
-
-            virtual ~solver() = default;
+        //- Destructor
+        virtual ~solver() = default;
 
 
         // Member Functions
@@ -244,21 +244,33 @@ public:
 
 
             //- Read and reset the solver parameters from the given dictionary
-            virtual void read(const dictionary& solverDict);
+            virtual void read(const dictionary&);
 
             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;
+            ) const
+            {
+                return this->normFactor(psi, Apsi, tmpField, normType_);
+            }
     };
 
 
@@ -314,6 +326,7 @@ public:
 
         // Constructors
 
+            //- Construct for given field name and matrix
             smoother
             (
                 const word& fieldName,
@@ -332,9 +345,8 @@ public:
             );
 
 
-        // Destructor
-
-            virtual ~smoother() = default;
+        //- Destructor
+        virtual ~smoother() = default;
 
 
         // Member Functions
@@ -405,10 +417,8 @@ public:
 
         // Constructors
 
-            preconditioner
-            (
-                const solver& sol
-            )
+            //- Construct for given solver
+            preconditioner(const solver& sol)
             :
                 solver_(sol)
             {}
@@ -424,16 +434,15 @@ 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& preconditionerDict)
+            //- from the given dictionary
+            virtual void read(const dictionary&)
             {}
 
             //- Return wA the preconditioned form of residual rA
@@ -444,7 +453,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
             (
@@ -480,17 +489,16 @@ 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
+            const lduMesh& mesh() const noexcept
             {
                 return lduMesh_;
             }
@@ -554,43 +562,43 @@ public:
             }
 
 
-            bool hasDiag() const
+            bool hasDiag() const noexcept
             {
                 return (diagPtr_);
             }
 
-            bool hasUpper() const
+            bool hasUpper() const noexcept
             {
                 return (upperPtr_);
             }
 
-            bool hasLower() const
+            bool hasLower() const noexcept
             {
                 return (lowerPtr_);
             }
 
-            bool hasSource() const
+            bool hasSource() const noexcept
             {
                 return (sourcePtr_);
             }
 
-            bool diagonal() const
+            bool diagonal() const noexcept
             {
                 return (diagPtr_ && !lowerPtr_ && !upperPtr_);
             }
 
-            bool symmetric() const
+            bool symmetric() const noexcept
             {
                 return (diagPtr_ && (!lowerPtr_ && upperPtr_));
             }
 
-            bool asymmetric() const
+            bool asymmetric() const noexcept
             {
                 return (diagPtr_ && lowerPtr_ && upperPtr_);
             }
 
 
-        // operations
+        // Operations
 
             void sumDiag();
             void negSumDiag();
@@ -640,7 +648,7 @@ public:
             tmp<Field<Type>> faceH(const tmp<Field<Type>>&) const;
 
 
-    // Member operators
+    // Member Operators
 
         void operator=(const LduMatrix<Type, DType, LUType>&);
 
@@ -653,7 +661,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-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -131,8 +131,9 @@ Foam::LduMatrix<Type, DType, LUType>::solver::solver
 
     log_(1),
     minIter_(0),
-    maxIter_(defaultMaxIter_),
-    tolerance_(1e-6*pTraits<Type>::one),
+    maxIter_(lduMatrix::defaultMaxIter),
+    normType_(lduMatrix::normTypes::DEFAULT_NORM),
+    tolerance_(lduMatrix::defaultTolerance*pTraits<Type>::one),
     relTol_(Zero)
 {
     readControls();
@@ -145,6 +146,8 @@ 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_);
@@ -168,21 +171,45 @@ Type Foam::LduMatrix<Type, DType, LUType>::solver::normFactor
 (
     const Field<Type>& psi,
     const Field<Type>& Apsi,
-    Field<Type>& tmpField
+    Field<Type>& tmpField,
+    const lduMatrix::normTypes normType
 ) const
 {
-    // --- Calculate A dot reference value of psi
-    matrix_.sumA(tmpField);
-    cmptMultiply(tmpField, tmpField, gAverage(psi));
+    switch (normType)
+    {
+        case lduMatrix::normTypes::NO_NORM :
+        {
+            break;
+        }
 
-    return stabilise
-    (
-        gSum(cmptMag(Apsi - tmpField) + cmptMag(matrix_.source() - tmpField)),
-        SolverPerformance<Type>::small_
-    );
+        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));
 
-    // At convergence this simpler method is equivalent to the above
-    // return stabilise(2*gSumCmptMag(matrix_.source()), matrix_.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;
 }
 
 

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

@@ -34,8 +34,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef DiagonalSolver_H
-#define DiagonalSolver_H
+#ifndef Foam_DiagonalSolver_H
+#define Foam_DiagonalSolver_H
 
 #include "LduMatrix.H"
 
@@ -53,7 +53,9 @@ class DiagonalSolver
 :
     public LduMatrix<Type, DType, LUType>::solver
 {
-    // Private Member Functions
+public:
+
+    // Generated Methods
 
         //- No copy construct
         DiagonalSolver(const DiagonalSolver&) = delete;
@@ -62,8 +64,6 @@ class DiagonalSolver
         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 SmoothSolver_H
-#define SmoothSolver_H
+#ifndef Foam_SmoothSolver_H
+#define Foam_SmoothSolver_H
 
 #include "lduMatrix.H"
 
@@ -56,10 +56,9 @@ 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,13 +32,11 @@ License
 template<class Type>
 Foam::lduCalculatedProcessorField<Type>::lduCalculatedProcessorField
 (
-    const lduInterface& interface,
-    const Field<Type>& iF
+    const lduInterface& interface
 )
 :
     LduInterfaceField<Type>(interface),
     procInterface_(refCast<const lduPrimitiveProcessorInterface>(interface)),
-    field_(iF),
     sendBuf_(procInterface_.faceCells().size()),
     receiveBuf_(procInterface_.faceCells().size()),
     scalarSendBuf_(procInterface_.faceCells().size()),
@@ -56,7 +54,6 @@ 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,8 +30,7 @@ Group
     grpGenericBoundaryConditions
 
 Description
-    A lduProcessorField type bypassing coupledFvPatchField and holding
-    a reference to the Field<Type>.
+    A lduProcessorField type bypassing coupledFvPatchField
 
     Used to add updateInterfaceMatrix capabilities to a lduMatrix
     which is fully uncoupled from the fvMesh.
@@ -72,8 +71,6 @@ protected:
         //- Local reference cast into the interface
          const lduPrimitiveProcessorInterface& procInterface_;
 
-        //- Local Field
-         const Field<Type>& field_;
 
         // Sending and receiving
 
@@ -118,8 +115,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
     )
     {
-        OPstream::write
+        UOPstream::write
         (
             commsType,
             neighbProcNo(),
@@ -61,7 +61,7 @@ void Foam::processorLduInterface::send
     {
         resizeBuf(receiveBuf_, nBytes);
 
-        IPstream::read
+        UIPstream::read
         (
             commsType,
             neighbProcNo(),
@@ -77,7 +77,7 @@ void Foam::processorLduInterface::send
             static_cast<void*>(sendBuf_.data()), f.cdata(), nBytes
         );
 
-        OPstream::write
+        UOPstream::write
         (
             commsType,
             neighbProcNo(),
@@ -109,7 +109,7 @@ void Foam::processorLduInterface::receive
      || commsType == Pstream::commsTypes::scheduled
     )
     {
-        IPstream::read
+        UIPstream::read
         (
             commsType,
             neighbProcNo(),
@@ -181,7 +181,7 @@ void Foam::processorLduInterface::compressedSend
          || commsType == Pstream::commsTypes::scheduled
         )
         {
-            OPstream::write
+            UOPstream::write
             (
                 commsType,
                 neighbProcNo(),
@@ -195,7 +195,7 @@ void Foam::processorLduInterface::compressedSend
         {
             resizeBuf(receiveBuf_, nBytes);
 
-            IPstream::read
+            UIPstream::read
             (
                 commsType,
                 neighbProcNo(),
@@ -205,7 +205,7 @@ void Foam::processorLduInterface::compressedSend
                 comm()
             );
 
-            OPstream::write
+            UOPstream::write
             (
                 commsType,
                 neighbProcNo(),
@@ -252,7 +252,7 @@ void Foam::processorLduInterface::compressedReceive
         {
             resizeBuf(receiveBuf_, nBytes);
 
-            IPstream::read
+            UIPstream::read
             (
                 commsType,
                 neighbProcNo(),

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

@@ -41,7 +41,18 @@ namespace Foam
 }
 
 
-const Foam::label Foam::lduMatrix::solver::defaultMaxIter_ = 1000;
+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" },
+});
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -50,8 +50,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef lduMatrix_H
-#define lduMatrix_H
+#ifndef Foam_lduMatrix_H
+#define Foam_lduMatrix_H
 
 #include "lduMesh.H"
 #include "primitiveFieldsFwd.H"
@@ -62,6 +62,7 @@ SourceFiles
 #include "runTimeSelectionTables.H"
 #include "solverPerformance.H"
 #include "InfoProxy.H"
+#include "Enum.H"
 #include "profilingTrigger.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -70,7 +71,6 @@ namespace Foam
 {
 
 // Forward Declarations
-
 class lduMatrix;
 
 Ostream& operator<<(Ostream&, const lduMatrix&);
@@ -95,6 +95,26 @@ 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
     {
@@ -102,19 +122,16 @@ 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 controls
+            //- Dictionary of solution controls
             dictionary controlDict_;
 
-            //- Level of verbosity in the solver output statements
+            //- Verbosity level for solver output statements
             int log_;
 
             //- Minimum number of iterations in the solver
@@ -123,18 +140,22 @@ 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 the controlDict_
+            //- Read the control parameters from controlDict_
             virtual void readControls();
 
 
@@ -195,6 +216,7 @@ public:
 
         // Constructors
 
+            //- Construct solver for given field name, matrix etc
             solver
             (
                 const word& fieldName,
@@ -272,6 +294,16 @@ 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
@@ -280,7 +312,10 @@ public:
                 const solveScalarField& source,
                 const solveScalarField& Apsi,
                 solveScalarField& tmpField
-            ) const;
+            ) const
+            {
+                return this->normFactor(psi, source, Apsi, tmpField, normType_);
+            }
     };
 
 
@@ -354,6 +389,7 @@ public:
 
         // Constructors
 
+            //- Construct for given field name, matrix etc
             smoother
             (
                 const word& fieldName,
@@ -479,10 +515,8 @@ public:
 
         // Constructors
 
-            preconditioner
-            (
-                const solver& sol
-            )
+            //- Construct for given solver
+            explicit preconditioner(const solver& sol)
             :
                 solver_(sol)
             {}
@@ -564,7 +598,7 @@ public:
         // Access to addressing
 
             //- Return the LDU mesh from which the addressing is obtained
-            const lduMesh& mesh() const
+            const lduMesh& mesh() const noexcept
             {
                 return lduMesh_;
             }
@@ -606,38 +640,38 @@ public:
             const scalarField& diag() const;
             const scalarField& upper() const;
 
-            bool hasDiag() const
+            bool hasDiag() const noexcept
             {
                 return (diagPtr_);
             }
 
-            bool hasUpper() const
+            bool hasUpper() const noexcept
             {
                 return (upperPtr_);
             }
 
-            bool hasLower() const
+            bool hasLower() const noexcept
             {
                 return (lowerPtr_);
             }
 
-            bool diagonal() const
+            bool diagonal() const noexcept
             {
                 return (diagPtr_ && !lowerPtr_ && !upperPtr_);
             }
 
-            bool symmetric() const
+            bool symmetric() const noexcept
             {
                 return (diagPtr_ && (!lowerPtr_ && upperPtr_));
             }
 
-            bool asymmetric() const
+            bool asymmetric() const noexcept
             {
                 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-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -152,6 +152,14 @@ 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();
@@ -162,11 +170,19 @@ Foam::lduMatrix::solver::solver
 
 void Foam::lduMatrix::solver::readControls()
 {
-    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);
+    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_);
 }
 
 
@@ -199,24 +215,40 @@ Foam::solveScalarField::cmptType Foam::lduMatrix::solver::normFactor
     const solveScalarField& psi,
     const solveScalarField& source,
     const solveScalarField& Apsi,
-    solveScalarField& tmpField
+    solveScalarField& tmpField,
+    const lduMatrix::normTypes normType
 ) const
 {
-    // --- Calculate A dot reference value of psi
-    matrix_.sumA(tmpField, interfaceBouCoeffs_, interfaces_);
+    switch (normType)
+    {
+        case lduMatrix::normTypes::NO_NORM :
+        {
+            break;
+        }
 
-    tmpField *= gAverage(psi, matrix_.mesh().comm());
+        case lduMatrix::normTypes::DEFAULT_NORM :
+        case lduMatrix::normTypes::L1_SCALED_NORM :
+        {
+            // --- Calculate A dot reference value of psi
+            matrix_.sumA(tmpField, interfaceBouCoeffs_, interfaces_);
 
-    return
-        gSum
-        (
-            (mag(Apsi - tmpField) + mag(source - tmpField))(),
-            matrix_.mesh().comm()
-        )
-      + solverPerformance::small_;
+            tmpField *= gAverage(psi, matrix_.mesh().comm());
 
-    // At convergence this simpler method is equivalent to the above
-    // return 2*gSumMag(source) + 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);
 }
 
 

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

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

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

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

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

@@ -51,6 +51,7 @@ SourceFiles
 namespace Foam
 {
 
+// Forward Declarations
 class GAMGAgglomeration;
 class lduMesh;
 class lduPrimitiveMesh;
@@ -63,7 +64,7 @@ class procFacesGAMGProcAgglomeration
 :
     public GAMGProcAgglomeration
 {
-    // Private data
+    // Private Data
 
         //- When to processor agglomerate
         const label nAgglomeratingCells_;
@@ -121,9 +122,8 @@ 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 OpenCFD Ltd.
+    Copyright (C) 2021-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -69,7 +69,6 @@ Foam::GAMGSolver::GAMGSolver
 
     // Default values for all controls
     // which may be overridden by those in controlDict
-    cacheAgglomeration_(true),
     nPreSweeps_(0),
     preSweepsLevelMultiplier_(1),
     maxPreSweeps_(4),
@@ -77,9 +76,12 @@ 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()),
@@ -389,14 +391,7 @@ void Foam::GAMGSolver::readControls()
 
 const Foam::lduMatrix& Foam::GAMGSolver::matrixLevel(const label i) const
 {
-    if (i == 0)
-    {
-        return matrix_;
-    }
-    else
-    {
-        return matrixLevels_[i - 1];
-    }
+    return i ? matrixLevels_[i-1] : matrix_;
 }
 
 
@@ -405,14 +400,7 @@ const Foam::lduInterfaceFieldPtrsList& Foam::GAMGSolver::interfaceLevel
     const label i
 ) const
 {
-    if (i == 0)
-    {
-        return interfaces_;
-    }
-    else
-    {
-        return interfaceLevels_[i - 1];
-    }
+    return i ? interfaceLevels_[i-1] : interfaces_;
 }
 
 
@@ -422,14 +410,7 @@ Foam::GAMGSolver::interfaceBouCoeffsLevel
     const label i
 ) const
 {
-    if (i == 0)
-    {
-        return interfaceBouCoeffs_;
-    }
-    else
-    {
-        return interfaceLevelsBouCoeffs_[i - 1];
-    }
+    return i ? interfaceLevelsBouCoeffs_[i-1] : interfaceBouCoeffs_;
 }
 
 
@@ -439,14 +420,7 @@ Foam::GAMGSolver::interfaceIntCoeffsLevel
     const label i
 ) const
 {
-    if (i == 0)
-    {
-        return interfaceIntCoeffs_;
-    }
-    else
-    {
-        return interfaceLevelsIntCoeffs_[i - 1];
-    }
+    return i ? interfaceLevelsIntCoeffs_[i-1] : interfaceIntCoeffs_;
 }
 
 

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

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -56,12 +56,11 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef GAMGSolver_H
-#define GAMGSolver_H
+#ifndef Foam_GAMGSolver_H
+#define Foam_GAMGSolver_H
 
 #include "GAMGAgglomeration.H"
 #include "lduMatrix.H"
-#include "labelField.H"
 #include "primitiveFields.H"
 #include "LUscalarMatrix.H"
 
@@ -78,9 +77,7 @@ class GAMGSolver
 :
     public lduMatrix::solver
 {
-    // Private data
-
-        bool cacheAgglomeration_;
+    // Private Data
 
         //- Number of pre-smoothing sweeps
         label nPreSweeps_;
@@ -103,6 +100,9 @@ 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 fromSlave
+            IPstream fromProc
             (
                 Pstream::commsTypes::scheduled,
                 procIDs[proci],
@@ -327,14 +327,14 @@ void Foam::GAMGSolver::gatherMatrices
                 meshComm
             );
 
-            otherMats.set(otherI, new lduMatrix(dummyMesh, fromSlave));
+            otherMats.set(otherI, new lduMatrix(dummyMesh, fromProc));
 
             // Receive number of/valid interfaces
             boolList& procTransforms = otherTransforms[otherI];
             List<label>& procRanks = otherRanks[otherI];
 
-            fromSlave >> procTransforms;
-            fromSlave >> procRanks;
+            fromProc >> procTransforms;
+            fromProc >> procRanks;
 
             // Size coefficients
             otherBouCoeffs.set
@@ -354,12 +354,12 @@ void Foam::GAMGSolver::gatherMatrices
                     otherBouCoeffs[otherI].set
                     (
                         intI,
-                        new scalarField(fromSlave)
+                        new scalarField(fromProc)
                     );
                     otherIntCoeffs[otherI].set
                     (
                         intI,
-                        new scalarField(fromSlave)
+                        new scalarField(fromProc)
                     );
                 }
             }

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();
-        IPstream::read
+        UIPstream::read
         (
             Pstream::commsTypes::nonBlocking,
             procInterface_.neighbProcNo(),
@@ -121,7 +121,7 @@ void Foam::processorGAMGInterfaceField::initInterfaceMatrixUpdate
         );
 
         outstandingSendRequest_ = UPstream::nRequests();
-        OPstream::write
+        UOPstream::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-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -224,11 +224,13 @@ Foam::solverPerformance Foam::PBiCG::solve
     }
 
     // Recommend PBiCGStab if PBiCG fails to converge
-    if (solverPerf.nIterations() > max(defaultMaxIter_, maxIter_))
+    const label upperMaxIters = max(maxIter_, lduMatrix::defaultMaxIter);
+
+    if (solverPerf.nIterations() > upperMaxIters)
     {
         FatalErrorInFunction
-            << "PBiCG has failed to converge within the maximum number"
-               " of iterations " << max(defaultMaxIter_, maxIter_) << nl
+            << "PBiCG has failed to converge within the maximum iterations: "
+            << upperMaxIters << nl
             << "    Please try the more robust PBiCGStab solver."
             << exit(FatalError);
     }

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

@@ -37,8 +37,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef diagonalSolver_H
-#define diagonalSolver_H
+#ifndef Foam_diagonalSolver_H
+#define Foam_diagonalSolver_H
 
 #include "lduMatrix.H"
 
@@ -55,7 +55,9 @@ class diagonalSolver
 :
     public lduMatrix::solver
 {
-    // Private Member Functions
+public:
+
+    // Generated Methods
 
         //- No copy construct
         diagonalSolver(const diagonalSolver&) = delete;
@@ -64,8 +66,6 @@ class diagonalSolver
         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 smoothSolver_H
-#define smoothSolver_H
+#ifndef Foam_smoothSolver_H
+#define Foam_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,6 +95,7 @@ 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 lduMesh_H
-#define lduMesh_H
+#ifndef Foam_lduMesh_H
+#define Foam_lduMesh_H
 
 #include "lduAddressing.H"
 #include "lduInterfacePtrsList.H"
@@ -46,11 +46,8 @@ Description
 namespace Foam
 {
 
+// Forward Declarations
 class objectRegistry;
-
-
-// Forward declaration of friend functions and operators
-
 class lduMesh;
 
 Ostream& operator<<(Ostream&, const InfoProxy<lduMesh>&);
@@ -62,15 +59,12 @@ Ostream& operator<<(Ostream&, const InfoProxy<lduMesh>&);
 
 class lduMesh
 {
-
 public:
 
     //- Runtime type information
     TypeName("lduMesh");
 
 
-    // Constructors
-
     //- Destructor
     virtual ~lduMesh() = default;
 
@@ -114,7 +108,7 @@ public:
             }
 
 
-    // Ostream operator
+    // Ostream Operator
 
         friend Ostream& operator<<(Ostream&, const InfoProxy<lduMesh>&);
 };

src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.C

@@ -37,25 +37,6 @@ 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();
-        }
-    };
-
 }
 
 
@@ -1092,15 +1073,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 slave " << procIDs[i] << endl;
+            //    << " receiving from proc " << procIDs[i] << endl;
 
-            IPstream fromSlave
+            IPstream fromProc
             (
                 Pstream::commsTypes::scheduled,
                 procIDs[i],
@@ -1109,10 +1090,10 @@ void Foam::lduPrimitiveMesh::gather
                 comm
             );
 
-            label nCells = readLabel(fromSlave);
-            labelList lowerAddr(fromSlave);
-            labelList upperAddr(fromSlave);
-            boolList validInterface(fromSlave);
+            label nCells = readLabel(fromProc);
+            labelList lowerAddr(fromProc);
+            labelList upperAddr(fromProc);
+            boolList validInterface(fromProc);
 
 
             // Construct mesh without interfaces
@@ -1135,7 +1116,7 @@ void Foam::lduPrimitiveMesh::gather
             {
                 if (validInterface[intI])
                 {
-                    word coupleType(fromSlave);
+                    word coupleType(fromProc);
 
                     newInterfaces.set
                     (
@@ -1145,7 +1126,7 @@ void Foam::lduPrimitiveMesh::gather
                             coupleType,
                             intI,
                             otherMeshes[i-1].rawInterfaces(),
-                            fromSlave
+                            fromProc
                         ).ptr()
                     );
                 }

src/OpenFOAM/meshes/lduMesh/lduPrimitiveMesh.H

@@ -35,8 +35,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef lduPrimitiveMesh_H
-#define lduPrimitiveMesh_H
+#ifndef Foam_lduPrimitiveMesh_H
+#define Foam_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 "pyramidPointFaceRef.H"
+#include "pyramid.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 

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

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

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

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

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

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

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

@@ -43,8 +43,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef face_H
-#define face_H
+#ifndef Foam_face_H
+#define Foam_face_H
 
 #include "pointField.H"
 #include "labelList.H"
@@ -435,6 +435,12 @@ 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.
@@ -494,19 +500,11 @@ template<> struct Hash<face> : face::hasher {};
 template<>
 struct offsetOp<face>
 {
-    inline face operator()
-    (
-        const face& x,
-        const label offset
-    ) const
+    face operator()(const face& x, const label offset) const
     {
-        face result(x.size());
-
-        forAll(x, i)
-        {
-            result[i] = x[i] + offset;
-        }
-        return result;
+        face f(x);
+        f += offset;
+        return f;
     }
 };
 

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

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -190,6 +190,20 @@ 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 "triPointRef.H"
+#include "triangle.H"
 #include "line.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

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

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

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

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

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

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

@@ -40,15 +40,15 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef tetCell_H
-#define tetCell_H
+#ifndef Foam_tetCell_H
+#define Foam_tetCell_H
 
 #include "FixedList.H"
 #include "triFace.H"
 #include "faceList.H"
 #include "edgeList.H"
 #include "pointField.H"
-#include "tetPointRef.H"
+#include "tetrahedron.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -85,10 +85,10 @@ public:
         //- Construct from four point labels
         inline tetCell
         (
-            const label a,
-            const label b,
-            const label c,
-            const label d
+            const label p0,
+            const label p1,
+            const label p2,
+            const label p3
         );
 
         //- 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 a,
-    const label b,
-    const label c,
-    const label d
+    const label p0,
+    const label p1,
+    const label p2,
+    const label p3
 )
 {
-    operator[](0) = a;
-    operator[](1) = b;
-    operator[](2) = c;
-    operator[](3) = d;
+    operator[](0) = p0;
+    operator[](1) = p1;
+    operator[](2) = p2;
+    operator[](3) = p3;
 }
 
 

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

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

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

@@ -28,7 +28,6 @@ License
 
 #include "IOstreams.H"
 #include "face.H"
-#include "triPointRef.H"
 #include <algorithm>  // For std::swap
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
@@ -68,16 +67,11 @@ inline Foam::triFace::triFace()
 {}
 
 
-inline Foam::triFace::triFace
-(
-    const label a,
-    const label b,
-    const label c
-)
+inline Foam::triFace::triFace(const label p0, const label p1, const label p2)
 {
-    operator[](0) = a;
-    operator[](1) = b;
-    operator[](2) = c;
+    operator[](0) = p0;
+    operator[](1) = p1;
+    operator[](2) = p2;
 }
 
 
@@ -476,6 +470,19 @@ 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())
             {
-                IPstream::read
+                UIPstream::read
                 (
                     UPstream::commsTypes::nonBlocking,
                     proci,
@@ -174,7 +174,7 @@ void Foam::mapDistributeBase::exchangeMasks
         {
             if (proci != myRank && sendMasks[proci].size())
             {
-                OPstream::write
+                UOPstream::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);
 
-                    OPstream::write
+                    UOPstream::write
                     (
                         Pstream::commsTypes::nonBlocking,
                         domain,
@@ -543,8 +543,7 @@ void Foam::mapDistributeBase::distribute
                 if (domain != myRank && map.size())
                 {
                     recvFields[domain].resize(map.size());
-
-                    IPstream::read
+                    UIPstream::read
                     (
                         Pstream::commsTypes::nonBlocking,
                         domain,
@@ -556,7 +555,6 @@ void Foam::mapDistributeBase::distribute
                 }
             }
 
-
             // Set up 'send' to myself
             {
                 sendFields[myRank] =
@@ -982,7 +980,7 @@ void Foam::mapDistributeBase::distribute
                     sendFields[domain] =
                         accessAndFlip(field, map, subHasFlip, negOp);
 
-                    OPstream::write
+                    UOPstream::write
                     (
                         Pstream::commsTypes::nonBlocking,
                         domain,
@@ -1004,7 +1002,7 @@ void Foam::mapDistributeBase::distribute
 
                 if (domain != myRank && map.size())
                 {
-                    recvFields[domain].setSize(map.size());
+                    recvFields[domain].resize(map.size());
                     UIPstream::read
                     (
                         Pstream::commsTypes::nonBlocking,
@@ -1018,7 +1016,6 @@ 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 "pyramidPointFaceRef.H"
+#include "pyramid.H"
 #include "primitiveMeshTools.H"
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //

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

@@ -39,13 +39,12 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef polyMeshTetDecomposition_H
-#define polyMeshTetDecomposition_H
+#ifndef Foam_polyMeshTetDecomposition_H
+#define Foam_polyMeshTetDecomposition_H
 
 #include "polyMesh.H"
 #include "coupledPolyPatch.H"
 #include "syncTools.H"
-#include "tetPointRef.H"
 #include "tetIndices.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -54,7 +53,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 "tetPointRef.H"
-#include "triPointRef.H"
+#include "tetrahedron.H"
+#include "triangle.H"
 #include "polyMesh.H"
 #include "triFace.H"
 #include "face.H"
@@ -76,7 +76,6 @@ class tetIndices;
 Istream& operator>>(Istream&, tetIndices&);
 Ostream& operator<<(Ostream&, const tetIndices&);
 
-
 /*---------------------------------------------------------------------------*\
                          Class tetIndices Declaration
 \*---------------------------------------------------------------------------*/
@@ -149,35 +148,49 @@ public:
         label& tetPt() noexcept { return tetPti_; }
 
 
-    // Searching
+    // Mesh-related
 
-        //- Return the indices corresponding to the tri on the face for
-        //  this tet. The normal of the tri points out of the cell
+        //- 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;
 
-        //- Return the local indices corresponding to the tri on the face
-        //- for this tet. The normal of the tri points out of the cell
+        //- 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;
 
-        //- Return the geometry corresponding to this tet
+
+        //- The tet geometry for this tet,
+        //- where point0 is the cell centre.
         inline tetPointRef tet(const polyMesh& mesh) const;
 
-        //- Return the geometry corresponding to the tri on the face for
-        //- this tet. The normal of the tri points out of the cell
+        //- 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
         inline triPointRef faceTri(const polyMesh& mesh) const;
 
-        //- Return the geometry corresponding to the tri on the face for
-        //- this tet using the old positions
+        //- The triangle geometry for the face for this tet
+        //- (using 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()[face()];
+    const Foam::face& f = mesh.faces()[facei_];
 
-    label faceBasePtI = mesh.tetBasePtIs()[face()];
+    label faceBasePtI = mesh.tetBasePtIs()[facei_];
 
     if (faceBasePtI < 0)
     {
@@ -67,24 +67,25 @@ inline Foam::triFace Foam::tetIndices::faceTriIs
 
         if (warn && nWarnings_ < maxNWarnings)
         {
-            ++nWarnings_;
             WarningInFunction
-                << "No base point for face " << face() << ", " << f
+                << "No base point for face " << facei_ << ", " << f
                 << ", produces a valid tet decomposition." << endl;
-            if (nWarnings_ == maxNWarnings)
+
+            if (++nWarnings_ == maxNWarnings)
             {
                 Warning
-                    << "Suppressing any further warnings." << endl;
+                    << "Suppressing further warnings." << endl;
             }
         }
     }
 
-    label facePtI = (tetPti_ + faceBasePtI) % f.size();
-    label faceOtherPtI = f.fcIndex(facePtI);
+    const label facePtI = (tetPti_ + faceBasePtI) % f.size();
+    const label faceOtherPtI = f.fcIndex(facePtI);
 
-    if (mesh.faceOwner()[face()] != cell())
+    if (mesh.faceOwner()[facei_] != celli_)
     {
-        std::swap(facePtI, faceOtherPtI);
+        // Neighbour: return flipped face
+        return triFace(f[faceBasePtI], f[faceOtherPtI], f[facePtI]);
     }
 
     return triFace(f[faceBasePtI], f[facePtI], f[faceOtherPtI]);
@@ -97,9 +98,9 @@ inline Foam::triFace Foam::tetIndices::triIs
     const bool warn
 ) const
 {
-    const Foam::face& f = mesh.faces()[face()];
+    const Foam::face& f = mesh.faces()[facei_];
 
-    label faceBasePtI = mesh.tetBasePtIs()[face()];
+    label faceBasePtI = mesh.tetBasePtIs()[facei_];
 
     if (faceBasePtI < 0)
     {
@@ -107,24 +108,25 @@ inline Foam::triFace Foam::tetIndices::triIs
 
         if (warn && nWarnings_ < maxNWarnings)
         {
-            ++nWarnings_;
             WarningInFunction
-                << "No base point for face " << face() << ", " << f
+                << "No base point for face " << facei_ << ", " << f
                 << ", produces a valid tet decomposition." << endl;
-            if (nWarnings_ == maxNWarnings)
+
+            if (++nWarnings_ == maxNWarnings)
             {
                 Warning
-                    << "Suppressing any further warnings." << endl;
+                    << "Suppressing further warnings." << endl;
             }
         }
     }
 
-    label facePtI = (tetPti_ + faceBasePtI) % f.size();
-    label faceOtherPtI = f.fcIndex(facePtI);
+    const label facePtI = (tetPti_ + faceBasePtI) % f.size();
+    const label faceOtherPtI = f.fcIndex(facePtI);
 
-    if (mesh.faceOwner()[face()] != cell())
+    if (mesh.faceOwner()[facei_] != celli_)
     {
-        std::swap(facePtI, faceOtherPtI);
+        // Neighbour: return flipped face
+        return triFace(faceBasePtI, faceOtherPtI, facePtI);
     }
 
     return triFace(faceBasePtI, facePtI, faceOtherPtI);
@@ -133,30 +135,60 @@ inline Foam::triFace Foam::tetIndices::triIs
 
 inline Foam::tetPointRef Foam::tetIndices::tet(const polyMesh& mesh) const
 {
-    const pointField& meshPoints = mesh.points();
+    const pointField& pts = mesh.points();
     const triFace tri = faceTriIs(mesh);
 
     return tetPointRef
     (
-        mesh.cellCentres()[cell()],
-        meshPoints[tri[0]],
-        meshPoints[tri[1]],
-        meshPoints[tri[2]]
+        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())
     );
 }
 
 
 inline Foam::triPointRef Foam::tetIndices::faceTri(const polyMesh& mesh) const
 {
-    const pointField& meshPoints = mesh.points();
-    const triFace tri = faceTriIs(mesh);
-
-    return triPointRef
-    (
-        meshPoints[tri[0]],
-        meshPoints[tri[1]],
-        meshPoints[tri[2]]
-    );
+    return faceTriIs(mesh).tri(mesh.points());
 }
 
 
@@ -165,15 +197,7 @@ inline Foam::triPointRef Foam::tetIndices::oldFaceTri
     const polyMesh& mesh
 ) const
 {
-    const pointField& meshOldPoints = mesh.oldPoints();
-    const triFace tri = faceTriIs(mesh);
-
-    return triPointRef
-    (
-        meshOldPoints[tri[0]],
-        meshOldPoints[tri[1]],
-        meshOldPoints[tri[2]]
-    );
+    return faceTriIs(mesh).tri(mesh.oldPoints());
 }
 
 

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

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

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

@@ -27,7 +27,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "primitiveMesh.H"
-#include "pyramidPointFaceRef.H"
+#include "pyramid.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 "pyramidPointFaceRef.H"
+#include "pyramid.H"
 #include "PrecisionAdaptor.H"
 
 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //

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

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

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

@@ -35,12 +35,13 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef line_H
-#define line_H
+#ifndef Foam_line_H
+#define Foam_line_H
 
+#include "point.H"
+#include "point2D.H"
 #include "vector.H"
 #include "PointHit.H"
-#include "point2D.H"
 #include "FixedList.H"
 #include "UList.H"
 
@@ -60,6 +61,12 @@ 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
 \*---------------------------------------------------------------------------*/
@@ -100,19 +107,19 @@ public:
     // Access
 
         //- Return first point
-        inline PointRef first() const noexcept;
+        inline PointRef first() const noexcept { return a_; }
 
         //- Return second (last) point
-        inline PointRef second() const noexcept;
+        inline PointRef second() const noexcept { return b_; }
 
         //- Return last (second) point
-        inline PointRef last() const noexcept;
+        inline PointRef last() const noexcept { return b_; }
 
-        //- Return first point
-        inline PointRef start() const noexcept;
+        //- Return start (first) point
+        inline PointRef start() const noexcept { return a_; }
 
-        //- Return second (last) point
-        inline PointRef end() const noexcept;
+        //- Return end (second) point
+        inline PointRef end() const noexcept { return b_; }
 
 
     // Properties

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

@@ -60,40 +60,6 @@ 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/linePointRef.H

@@ -1,55 +1,11 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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.
+// Compatibility include.
+// linePointRef typedef included in line.H (JUL-2022)
 
-    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.
+#ifndef FoamCompat_linePointRef_H
+#define FoamCompat_linePointRef_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/>.
-
-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 point_H
-#define point_H
+#ifndef Foam_point_H
+#define Foam_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 OpenCFD Ltd.
+    Copyright (C) 2020-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,8 +37,11 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef pyramid_H
-#define pyramid_H
+#ifndef Foam_pyramid_H
+#define Foam_pyramid_H
+
+#include "point.H"
+#include "face.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,34 +50,39 @@ 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_;
 
@@ -90,7 +98,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);
@@ -98,36 +106,36 @@ public:
 
     // Member Functions
 
-        // Access
+    // Access
 
-            //- Return apex point
-            inline const Point& apex() const;
+        //- The apex point
+        const Point& apex() const { return apex_; }
 
-            //- Return base polygon
-            inline polygonRef base() const;
+        //- The base polygon
+        PolygonRef base() const { return base_; }
 
 
-        // 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,31 +41,17 @@ 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 >> base_ >> apex_;
-    is.check(FUNCTION_NAME);
+    is >> *this;
 }
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-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
+template<class Point, class PointRef, class PolygonRef>
+inline Point Foam::pyramid<Point, PointRef, PolygonRef>::centre
 (
     const UList<point>& points
 ) const
@@ -74,8 +60,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
@@ -85,8 +71,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
@@ -97,11 +83,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_;
@@ -110,11 +96,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,56 +1,11 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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.
+// Compatibility include
+// pyramidPointFaceRef typedef included in tetrahedron.H (SEP-2017)
 
-    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.
+#ifndef FoamCompat_pyramidPointFaceRef_H
+#define FoamCompat_pyramidPointFaceRef_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/>.
-
-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,55 +1,11 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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.
+// Compatibility include.
+// tetPointRef typedef included in tetrahedron.H (SEP-2017)
 
-    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.
+#ifndef FoamCompat_tetPointRef_H
+#define FoamCompat_tetPointRef_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/>.
-
-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,126 +1,10 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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.
+// Compatibility include.
+// tetPoints defined in tetrahedron.H (JUL-2022)
 
-    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
+#ifndef FoamCompat_tetPoints_H
+#define FoamCompat_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,7 +26,6 @@ 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 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -39,17 +39,17 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef tetrahedron_H
-#define tetrahedron_H
+#ifndef Foam_tetrahedron_H
+#define Foam_tetrahedron_H
 
 #include "point.H"
-#include "primitiveFieldsFwd.H"
 #include "pointHit.H"
+#include "primitiveFieldsFwd.H"
 #include "Random.H"
 #include "FixedList.H"
 #include "UList.H"
-#include "triPointRef.H"
-#include "boundBox.H"
+#include "triangle.H"
+#include "treeBoundBox.H"
 #include "barycentric.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -59,28 +59,103 @@ 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 tetrahedron Declaration
+                         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
 \*---------------------------------------------------------------------------*/
 
 template<class Point, class PointRef>
@@ -88,7 +163,10 @@ class tetrahedron
 {
 public:
 
-    // Public typedefs
+    // Public Typedefs
+
+        //- The point type
+        typedef Point point_type;
 
         //- Storage type for tets originating from intersecting tets.
         //  (can possibly be smaller than 200)
@@ -130,10 +208,13 @@ public:
 
 private:
 
-    // Private data
+    // Private Data
 
         PointRef a_, b_, c_, d_;
 
+
+    // Private Member Functions
+
         inline static point planeIntersection
         (
             const FixedList<scalar, 4>&,
@@ -161,7 +242,7 @@ private:
 
 public:
 
-    // Member constants
+    // Member Constants
 
         enum
         {
@@ -172,15 +253,18 @@ public:
 
     // Constructors
 
-        //- Construct from points
+        //- Construct from four points
         inline tetrahedron
         (
-            const Point& a,
-            const Point& b,
-            const Point& c,
-            const Point& d
+            const Point& p0,
+            const Point& p1,
+            const Point& p2,
+            const Point& p3
         );
 
+        //- Construct from four points
+        inline tetrahedron(const FixedList<Point, 4>& pts);
+
         //- Construct from four points in the list of points
         inline tetrahedron
         (
@@ -189,24 +273,27 @@ public:
         );
 
         //- Construct from Istream
-        inline tetrahedron(Istream&);
+        inline explicit tetrahedron(Istream&);
 
 
     // Member Functions
 
-        // Access
+    // Access
 
-            //- Return vertices
-            inline const Point& a() const;
+        //- Return vertex a
+        const Point& a() const noexcept { return a_; }
 
-            inline const Point& b() const;
+        //- Return vertex b
+        const Point& b() const noexcept { return b_; }
 
-            inline const Point& c() const;
+        //- Return vertex c
+        const Point& c() const noexcept { return c_; }
 
-            inline const Point& d() const;
+        //- Return vertex d
+        const Point& d() const noexcept { return d_; }
 
-            //- Return i-th face
-            inline triPointRef tri(const label facei) const;
+        //- Return i-th face
+        inline triPointRef tri(const label facei) const;
 
 
         // Properties
@@ -236,12 +323,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 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,24 +28,76 @@ 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& a,
-    const Point& b,
-    const Point& c,
-    const Point& d
+    const Point& p0,
+    const Point& p1,
+    const Point& p2,
+    const Point& p3
 )
 :
-    a_(a),
-    b_(b),
-    c_(c),
-    d_(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])
 {}
 
 
@@ -72,33 +124,31 @@ inline Foam::tetrahedron<Point, PointRef>::tetrahedron(Istream& is)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Point, class PointRef>
-inline const Point& Foam::tetrahedron<Point, PointRef>::a() const
+inline Foam::tetPointRef Foam::tetPoints::tet() const
 {
-    return a_;
+    return tetPointRef(a(), b(), c(), d());
 }
 
 
-template<class Point, class PointRef>
-inline const Point& Foam::tetrahedron<Point, PointRef>::b() const
+inline void Foam::tetPoints::flip()
 {
-    return b_;
+    // swap pt2 <-> pt3
+    point t(c());
+    c() = d();
+    d() = t;
 }
 
 
-template<class Point, class PointRef>
-inline const Point& Foam::tetrahedron<Point, PointRef>::c() const
+inline Foam::treeBoundBox Foam::tetPoints::bounds() const
 {
-    return c_;
+    treeBoundBox bb;
+    bb.add(static_cast<const FixedList<point, 4>&>(*this));
+    return bb;
 }
 
 
-template<class Point, class PointRef>
-inline const Point& Foam::tetrahedron<Point, PointRef>::d() const
-{
-    return d_;
-}
-
+// 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 Foam::triPointRef Foam::tetrahedron<Point, PointRef>::tri
@@ -257,7 +307,7 @@ inline Point Foam::tetrahedron<Point, PointRef>::barycentricToPoint
     const barycentric& bary
 ) const
 {
-    return bary[0]*a_ + bary[1]*b_ + bary[2]*c_ + bary[3]*d_;
+    return Point(bary.a()*a_ + bary.b()*b_ + bary.c()*c_ + bary.d()*d_);
 }
 
 
@@ -310,7 +360,7 @@ inline Foam::scalar Foam::tetrahedron<Point, PointRef>::pointToBarycentric
     bary[0] = res.x();
     bary[1] = res.y();
     bary[2] = res.z();
-    bary[3] = 1 - cmptSum(res);
+    bary[3] = 1 - bary[0] - bary[1] - bary[2];
 
     return detT;
 }
@@ -334,63 +384,83 @@ inline Foam::pointHit Foam::tetrahedron<Point, PointRef>::nearestPoint
 
     bool inside = true;
 
-    if (((p - b_) & Sa()) >= 0)
+    // Side a
     {
-        // p is outside halfspace plane of tri
-        pointHit info = triangle<Point, PointRef>(b_, c_, d_).nearestPoint(p);
+        const triangle<Point, PointRef> tria(b_, c_, d_);
 
-        inside = false;
-
-        if (info.distance() < minOutsideDistance)
+        if (((p - b_) & tria.areaNormal()) >= 0)
         {
-            closestPt = info.rawPoint();
+            // p is outside halfspace plane of tri
+            pointHit info = tria.nearestPoint(p);
 
-            minOutsideDistance = info.distance();
+            inside = false;
+
+            if (info.distance() < minOutsideDistance)
+            {
+                closestPt = info.rawPoint();
+
+                minOutsideDistance = info.distance();
+            }
         }
     }
 
-    if (((p - a_) & Sb()) >= 0)
+    // Side b
     {
-        // p is outside halfspace plane of tri
-        pointHit info = triangle<Point, PointRef>(a_, d_, c_).nearestPoint(p);
+        const triangle<Point, PointRef> tria(a_, d_, c_);
 
-        inside = false;
-
-        if (info.distance() < minOutsideDistance)
+        if (((p - a_) & tria.areaNormal()) >= 0)
         {
-            closestPt = info.rawPoint();
+            // p is outside halfspace plane of tri
+            pointHit info = tria.nearestPoint(p);
 
-            minOutsideDistance = info.distance();
+            inside = false;
+
+            if (info.distance() < minOutsideDistance)
+            {
+                closestPt = info.rawPoint();
+
+                minOutsideDistance = info.distance();
+            }
         }
     }
 
-    if (((p - a_) & Sc()) >= 0)
+    // Side c
     {
-        // p is outside halfspace plane of tri
-        pointHit info = triangle<Point, PointRef>(a_, b_, d_).nearestPoint(p);
+        const triangle<Point, PointRef> tria(a_, b_, d_);
 
-        inside = false;
-
-        if (info.distance() < minOutsideDistance)
+        if (((p - a_) & tria.areaNormal()) >= 0)
         {
-            closestPt = info.rawPoint();
+            // p is outside halfspace plane of tri
+            pointHit info = tria.nearestPoint(p);
 
-            minOutsideDistance = info.distance();
+            inside = false;
+
+            if (info.distance() < minOutsideDistance)
+            {
+                closestPt = info.rawPoint();
+
+                minOutsideDistance = info.distance();
+            }
         }
     }
 
-    if (((p - a_) & Sd()) >= 0)
+    // Side c
     {
-        // p is outside halfspace plane of tri
-        pointHit info = triangle<Point, PointRef>(a_, c_, b_).nearestPoint(p);
+        const triangle<Point, PointRef> tria(a_, c_, b_);
 
-        inside = false;
-
-        if (info.distance() < minOutsideDistance)
+        if (((p - a_) & tria.areaNormal()) >= 0)
         {
-            closestPt = info.rawPoint();
+            // p is outside halfspace plane of tri
+            pointHit info = tria.nearestPoint(p);
 
-            minOutsideDistance = info.distance();
+            inside = false;
+
+            if (info.distance() < minOutsideDistance)
+            {
+                closestPt = info.rawPoint();
+
+                minOutsideDistance = info.distance();
+            }
         }
     }
 
@@ -412,7 +482,7 @@ inline Foam::pointHit Foam::tetrahedron<Point, PointRef>::nearestPoint
 
 
 template<class Point, class PointRef>
-bool Foam::tetrahedron<Point, PointRef>::inside(const point& pt) const
+bool Foam::tetrahedron<Point, PointRef>::inside(const point& p) const
 {
     // For robustness, assuming that the point is in the tet unless
     // "definitively" shown otherwise by obtaining a positive dot
@@ -429,55 +499,41 @@ bool Foam::tetrahedron<Point, PointRef>::inside(const point& pt) const
     // planeBase[2] = tetBasePt = b_
     // planeBase[3] = tetBasePt = b_
 
-    vector n = Zero;
-
+    // Side a
     {
-        // 0, a
-        const point& basePt = b_;
+        const triangle<Point, PointRef> tria(b_, c_, d_);
 
-        n = Sa();
-        n /= (Foam::mag(n) + VSMALL);
-
-        if (((pt - basePt) & n) > SMALL)
+        if (((p - b_) & tria.unitNormal()) > SMALL)
         {
             return false;
         }
     }
 
+    // Side b
     {
-        // 1, b
-        const point& basePt = c_;
+        const triangle<Point, PointRef> tria(a_, d_, c_);
 
-        n = Sb();
-        n /= (Foam::mag(n) + VSMALL);
-
-        if (((pt - basePt) & n) > SMALL)
+        if (((p - a_) & tria.unitNormal()) > SMALL)
         {
             return false;
         }
     }
 
+    // Side c
     {
-        // 2, c
-        const point& basePt = b_;
+        const triangle<Point, PointRef> tria(a_, b_, d_);
 
-        n = Sc();
-        n /= (Foam::mag(n) + VSMALL);
-
-        if (((pt - basePt) & n) > SMALL)
+        if (((p - a_) & tria.unitNormal()) > SMALL)
         {
             return false;
         }
     }
 
+    // Side d
     {
-        // 3, d
-        const point& basePt = b_;
+        const triangle<Point, PointRef> tria(a_, c_, b_);
 
-        n = Sd();
-        n /= (Foam::mag(n) + VSMALL);
-
-        if (((pt - basePt) & n) > SMALL)
+        if (((p - a_) & tria.unitNormal()) > SMALL)
         {
             return false;
         }
@@ -1055,6 +1111,7 @@ 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,55 +1,11 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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.
+// Compatibility include.
+// triPointRef typedef included in triangle.H (NOV-2015)
 
-    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.
+#ifndef FoamCompat_triPointRef_H
+#define FoamCompat_triPointRef_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/>.
-
-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,111 +1,10 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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.
+// Compatibility include.
+// triPoints defined in triangle.H (JUL-2022)
 
-    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.
+#ifndef FoamCompat_triPoints_H
+#define FoamCompat_triPoints_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
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#include "triangle.H"
 
 #endif
 

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

@@ -26,7 +26,6 @@ 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-2020 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -35,18 +35,20 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef triangle_H
-#define triangle_H
+#ifndef Foam_triangle_H
+#define Foam_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 "linePointRef.H"
+#include "line.H"
 #include "barycentric2D.H"
+#include "treeBoundBox.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -55,7 +57,6 @@ namespace Foam
 
 // Forward Declarations
 class plane;
-class triPoints;
 
 template<class Point, class PointRef> class triangle;
 
@@ -72,6 +73,73 @@ 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
 \*---------------------------------------------------------------------------*/
@@ -90,7 +158,7 @@ public:
         //- with another triangle
         typedef FixedList<triPoints, 27> triIntersectionList;
 
-        //- The proximity classifications
+        //- Proximity classifications
         enum proxType
         {
             NONE = 0,   //!< Unknown proximity
@@ -99,10 +167,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
@@ -169,10 +237,10 @@ public:
     // Constructors
 
         //- Construct from three points
-        inline triangle(const Point& a, const Point& b, const Point& c);
+        inline triangle(const Point& p0, const Point& p1, const Point& p2);
 
         //- Construct from three points
-        inline triangle(const FixedList<Point, 3>& tri);
+        inline triangle(const FixedList<Point, 3>& pts);
 
         //- Construct from three points in the list of points
         //  The indices could be from triFace etc.
@@ -188,16 +256,16 @@ public:
 
     // Member Functions
 
-        // Access
+    // Access
 
-            //- Return first vertex
-            inline const Point& a() const;
+        //- The first vertex
+        const Point& a() const noexcept { return a_; }
 
-            //- Return second vertex
-            inline const Point& b() const;
+        //- The second vertex
+        const Point& b() const noexcept { return b_; }
 
-            //- Return third vertex
-            inline const Point& c() const;
+        //- The third vertex
+        const Point& c() const noexcept { return c_; }
 
 
         // Properties
@@ -245,7 +313,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 OpenCFD Ltd.
+    Copyright (C) 2018-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,34 +28,70 @@ License
 
 #include "IOstreams.H"
 #include "pointHit.H"
-#include "triPoints.H"
 #include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
 
-template<class Point, class PointRef>
-inline Foam::triangle<Point, PointRef>::triangle
+inline Foam::triPoints::triPoints
 (
-    const Point& a,
-    const Point& b,
-    const Point& c
+    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
 )
 :
-    a_(a),
-    b_(b),
-    c_(c)
+    FixedList<point, 3>(points, indices)
 {}
 
 
 template<class Point, class PointRef>
 inline Foam::triangle<Point, PointRef>::triangle
 (
-    const FixedList<Point, 3>& tri
+    const Point& p0,
+    const Point& p1,
+    const Point& p2
 )
 :
-    a_(tri[0]),
-    b_(tri[1]),
-    c_(tri[2])
+    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])
 {}
 
 
@@ -72,7 +108,6 @@ inline Foam::triangle<Point, PointRef>::triangle
 {}
 
 
-
 template<class Point, class PointRef>
 inline Foam::triangle<Point, PointRef>::triangle(Istream& is)
 {
@@ -82,22 +117,26 @@ inline Foam::triangle<Point, PointRef>::triangle(Istream& is)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class Point, class PointRef>
-inline const Point& Foam::triangle<Point, PointRef>::a() const
+inline Foam::triPointRef Foam::triPoints::tri() const
 {
-    return a_;
+    return triPointRef(a(), b(), c());
 }
 
-template<class Point, class PointRef>
-inline const Point& Foam::triangle<Point, PointRef>::b() const
+
+inline void Foam::triPoints::flip()
 {
-    return b_;
+    // swap pt1 <-> pt2
+    point t(b());
+    b() = c();
+    c() = t;
 }
 
-template<class Point, class PointRef>
-inline const Point& Foam::triangle<Point, PointRef>::c() const
+
+inline Foam::treeBoundBox Foam::triPoints::bounds() const
 {
-    return c_;
+    treeBoundBox bb;
+    bb.add(static_cast<const FixedList<point, 3>&>(*this));
+    return bb;
 }
 
 
@@ -916,6 +955,7 @@ 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)
             {
-                IPstream::read
+                UIPstream::read
                 (
                     commsType,
                     procIDs[i],
@@ -123,7 +123,7 @@ void Foam::globalIndex::gatherValues
 
         if (is_contiguous<Type>::value)
         {
-            OPstream::write
+            UOPstream::write
             (
                 commsType,
                 procIDs[0],
@@ -198,7 +198,7 @@ void Foam::globalIndex::gather
             }
             else if (is_contiguous<Type>::value)
             {
-                IPstream::read
+                UIPstream::read
                 (
                     commsType,
                     procIDs[i],
@@ -223,7 +223,7 @@ void Foam::globalIndex::gather
         }
         else if (is_contiguous<Type>::value)
         {
-            OPstream::write
+            UOPstream::write
             (
                 commsType,
                 procIDs[0],
@@ -889,7 +889,7 @@ void Foam::globalIndex::scatter
             }
             else if (is_contiguous<Type>::value)
             {
-                OPstream::write
+                UOPstream::write
                 (
                     commsType,
                     procIDs[i],
@@ -926,7 +926,7 @@ void Foam::globalIndex::scatter
         }
         else if (is_contiguous<Type>::value)
         {
-            IPstream::read
+            UIPstream::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-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2022 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 Barycentric_H
-#define Barycentric_H
+#ifndef Foam_Barycentric_H
+#define Foam_Barycentric_H
 
 #include "contiguous.H"
 #include "VectorSpace.H"
@@ -74,7 +74,7 @@ public:
     enum components { A, B, C, D };
 
 
-    // Generated Methods
+    // Generated Methods: copy construct/assignment
 
         //- Default construct
         Barycentric() = default;
@@ -94,20 +94,29 @@ public:
             const Cmpt& vd
         );
 
+        //- Construct from three components, calculate fourth component
+        inline Barycentric(const Cmpt& va, const Cmpt& vb, const Cmpt& vc);
+
 
     // Member Functions
 
-        // Access
+    // Component 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();
+        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;
 };
 
 

src/OpenFOAM/primitives/Barycentric/BarycentricI.H

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

src/OpenFOAM/primitives/Barycentric/BarycentricTensor.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2017 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,18 +28,17 @@ Class
     Foam::BarycentricTensor
 
 Description
-    Templated 4x3 tensor derived from VectorSpace. Has 12 components. Can
-    represent a barycentric transformation as a matrix-barycentric inner-
-    product. Can alternatively represent an inverse barycentric transformation
-    as a vector-matrix inner-product.
+    Templated 4x3 tensor derived from VectorSpace. Has 12 components.
+    Can represent a barycentric transformation as a matrix-barycentric
+    inner-product.
 
 SourceFiles
     BarycentricTensorI.H
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef BarycentricTensor_H
-#define BarycentricTensor_H
+#ifndef Foam_BarycentricTensor_H
+#define Foam_BarycentricTensor_H
 
 #include "Barycentric.H"
 #include "Tensor.H"
@@ -64,7 +63,8 @@ public:
     // Typedefs
 
         //- Equivalent type of labels used for valid component indexing
-        typedef Tensor<label> labelType;
+        //- (unused)
+        typedef BarycentricTensor<label> labelType;
 
 
     // Member Constants
@@ -77,7 +77,7 @@ public:
     enum components { XA, XB, XC, XD, YA, YB, YC, YD, ZA, ZB, ZC, ZD };
 
 
-    // Generated Methods
+    // Generated Methods: copy construct/assignment
 
         //- Default construct
         BarycentricTensor() = default;
@@ -97,6 +97,7 @@ public:
         );
 
         //- Construct given four vector components (columns)
+        //  Eg, the corners of a tetrahedron
         inline BarycentricTensor
         (
             const Vector<Cmpt>& a,
@@ -108,18 +109,44 @@ public:
 
     // Member Functions
 
-        // Row-barycentric access
+    // Component access
 
-            inline Barycentric<Cmpt> x() const;
-            inline Barycentric<Cmpt> y() const;
-            inline Barycentric<Cmpt> z() const;
+        inline const Cmpt& xa() const;
+        inline const Cmpt& xb() const;
+        inline const Cmpt& xc() const;
+        inline const Cmpt& xd() const;
 
-        // Column-vector access
+        inline const Cmpt& ya() const;
+        inline const Cmpt& yb() const;
+        inline const Cmpt& yc() const;
+        inline const Cmpt& yd() const;
 
-            inline Vector<Cmpt> a() const;
-            inline Vector<Cmpt> b() const;
-            inline Vector<Cmpt> c() const;
-            inline Vector<Cmpt> d() const;
+        inline const Cmpt& za() const;
+        inline const Cmpt& zb() const;
+        inline const Cmpt& zc() const;
+        inline const Cmpt& zd() const;
+
+
+    // Row-barycentric access
+
+        inline Barycentric<Cmpt> x() const;
+        inline Barycentric<Cmpt> y() const;
+        inline Barycentric<Cmpt> z() const;
+
+
+    // Column-vector access
+
+        inline Vector<Cmpt> a() const;
+        inline Vector<Cmpt> b() const;
+        inline Vector<Cmpt> c() const;
+        inline Vector<Cmpt> d() const;
+
+
+    // Operations
+
+        //- Tensor/barycentric inner product
+        //  (transforms barycentric coordinates to vector)
+        inline Vector<Cmpt> inner(const Barycentric<Cmpt>& bry) const;
 };
 
 

src/OpenFOAM/primitives/Barycentric/BarycentricTensorI.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2017 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -87,6 +88,90 @@ inline Foam::BarycentricTensor<Cmpt>::BarycentricTensor
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class Cmpt>
+inline const Cmpt& Foam::BarycentricTensor<Cmpt>::xa() const
+{
+    return this->v_[XA];
+}
+
+
+template<class Cmpt>
+inline const Cmpt& Foam::BarycentricTensor<Cmpt>::xb() const
+{
+    return this->v_[XB];
+}
+
+
+template<class Cmpt>
+inline const Cmpt& Foam::BarycentricTensor<Cmpt>::xc() const
+{
+    return this->v_[XC];
+}
+
+
+template<class Cmpt>
+inline const Cmpt& Foam::BarycentricTensor<Cmpt>::xd() const
+{
+    return this->v_[XD];
+}
+
+
+template<class Cmpt>
+inline const Cmpt& Foam::BarycentricTensor<Cmpt>::ya() const
+{
+    return this->v_[YA];
+}
+
+
+template<class Cmpt>
+inline const Cmpt& Foam::BarycentricTensor<Cmpt>::yb() const
+{
+    return this->v_[YB];
+}
+
+
+template<class Cmpt>
+inline const Cmpt& Foam::BarycentricTensor<Cmpt>::yc() const
+{
+    return this->v_[YC];
+}
+
+
+template<class Cmpt>
+inline const Cmpt& Foam::BarycentricTensor<Cmpt>::yd() const
+{
+    return this->v_[YD];
+}
+
+
+template<class Cmpt>
+inline const Cmpt& Foam::BarycentricTensor<Cmpt>::za() const
+{
+    return this->v_[ZA];
+}
+
+
+template<class Cmpt>
+inline const Cmpt& Foam::BarycentricTensor<Cmpt>::zb() const
+{
+    return this->v_[ZB];
+}
+
+
+template<class Cmpt>
+inline const Cmpt& Foam::BarycentricTensor<Cmpt>::zc() const
+{
+    return this->v_[ZC];
+}
+
+
+template<class Cmpt>
+inline const Cmpt& Foam::BarycentricTensor<Cmpt>::zd() const
+{
+    return this->v_[ZD];
+}
+
+
 template<class Cmpt>
 inline Foam::Barycentric<Cmpt> Foam::BarycentricTensor<Cmpt>::x() const
 {
@@ -157,6 +242,26 @@ inline Foam::Vector<Cmpt> Foam::BarycentricTensor<Cmpt>::d() const
 }
 
 
+// NB: same workaround for gcc (11+) failure on (tensor dot vector)
+// - not sure it will indeed be required here as well.
+template<class Cmpt>
+#if defined(__GNUC__) && !defined(__clang__)
+__attribute__((optimize("no-tree-vectorize")))
+#endif
+inline Foam::Vector<Cmpt>
+Foam::BarycentricTensor<Cmpt>::inner(const Barycentric<Cmpt>& bry) const
+{
+    const BarycentricTensor<Cmpt>& t = *this;
+
+    return Vector<Cmpt>
+    (
+        (bry.a()*t.xa() + bry.b()*t.xb() + bry.c()*t.xc() + bry.d()*t.xd()),
+        (bry.a()*t.ya() + bry.b()*t.yb() + bry.c()*t.yc() + bry.d()*t.yd()),
+        (bry.a()*t.za() + bry.b()*t.zb() + bry.c()*t.zc() + bry.d()*t.zd())
+    );
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -164,17 +269,20 @@ namespace Foam
 
 // * * * * * * * * * * * * * * * Global Operators  * * * * * * * * * * * * * //
 
+// Transform Barycentric coordinates to Vector
 template<class Cmpt>
 inline Vector<Cmpt> operator&
 (
-    const BarycentricTensor<Cmpt>& T,
+    const BarycentricTensor<Cmpt>& t,
     const Barycentric<Cmpt>& b
 )
 {
-    return Vector<Cmpt>(T.x() & b, T.y() & b, T.z() & b);
+    return t.inner(b);
 }
 
 
+// Transform Vector to Barycentric coordinates.
+// Caution: the tensor must be inverse one (see particle.C)
 template<class Cmpt>
 inline Barycentric<Cmpt> operator&
 (

src/OpenFOAM/primitives/Barycentric2D/Barycentric2D.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2017 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,8 +36,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Barycentric2D_H
-#define Barycentric2D_H
+#ifndef Foam_Barycentric2D_H
+#define Foam_Barycentric2D_H
 
 #include "contiguous.H"
 #include "VectorSpace.H"
@@ -86,33 +86,31 @@ public:
         inline Barycentric2D(const Foam::zero);
 
         //- Construct from components
-        inline Barycentric2D
-        (
-            const Cmpt& va,
-            const Cmpt& vb,
-            const Cmpt& vc
-        );
+        inline Barycentric2D(const Cmpt& va, const Cmpt& vb, const Cmpt& vc);
+
+        //- Construct from two components, calculating third component
+        inline Barycentric2D(const Cmpt& va, const Cmpt& vb);
 
 
     // Member Functions
 
-    // Access
+    // Component access
 
         inline const Cmpt& a() const;
         inline const Cmpt& b() const;
         inline const Cmpt& c() const;
 
-
-    // Edit
-
         inline Cmpt& a();
         inline Cmpt& b();
         inline Cmpt& c();
 
 
-    // Tests
+    // Operations, Tests
 
-        //- True if any coordinates are negative
+        //- Scalar-product of \c this with another Barycentric2D.
+        inline Cmpt inner(const Barycentric2D<Cmpt>& b2) const;
+
+        //- True if any coordinate is negative
         inline bool outside() const;
 };