From cd133a28f589c5eb48e53737e9a71f530b1192c0 Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Wed, 20 Dec 2023 15:54:25 +0000
Subject: [PATCH 1/4] fvMeshTopoChangers::refiner: Fix hang in parallel

This has required switching the order of dependence of the fvMesh
geometry (i.e., Sf, Cf, and C), and the interpolation weights. The
weights are now considered dependent on the geometry fields, and the
geometry fields must therefore construct without reference to the
weights.
---
 .../surfaceInterpolation/surfaceInterpolation.C      | 12 ++++--------
 .../refiner/fvMeshTopoChangersRefiner.C              |  5 +----
 2 files changed, 5 insertions(+), 12 deletions(-)

diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
index 146dc770a8..44d46ada8b 100644
--- a/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
+++ b/src/finiteVolume/interpolation/surfaceInterpolation/surfaceInterpolation/surfaceInterpolation.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -158,15 +158,11 @@ void Foam::surfaceInterpolation::makeWeights() const
     surfaceScalarField& weights = *weights_;
 
     // Set local references to mesh data
-    // (note that we should not use fvMesh sliced fields at this point yet
-    //  since this causes a loop when generating weighting factors in
-    //  coupledFvPatchField evaluation phase)
     const labelUList& owner = mesh_.owner();
     const labelUList& neighbour = mesh_.neighbour();
-
-    const vectorField& Cf = mesh_.faceCentres();
-    const vectorField& C = mesh_.cellCentres();
-    const vectorField& Sf = mesh_.faceAreas();
+    const surfaceVectorField& Sf = mesh_.Sf();
+    const surfaceVectorField& Cf = mesh_.Cf();
+    const volVectorField& C = mesh_.C();
 
     // ... and reference to the internal field of the weighting factors
     scalarField& w = weights.primitiveFieldRef();
diff --git a/src/fvMeshTopoChangers/refiner/fvMeshTopoChangersRefiner.C b/src/fvMeshTopoChangers/refiner/fvMeshTopoChangersRefiner.C
index 8a8b3bc2e5..87091f9ae2 100644
--- a/src/fvMeshTopoChangers/refiner/fvMeshTopoChangersRefiner.C
+++ b/src/fvMeshTopoChangers/refiner/fvMeshTopoChangersRefiner.C
@@ -575,10 +575,7 @@ void Foam::fvMeshTopoChangers::refiner::refineUfs
         {
             const surfaceVectorField UfU
             (
-                fvc::interpolate
-                (
-                    mesh().lookupObject<volVectorField>(Uname)
-                )
+                fvc::interpolate(mesh().lookupObject<volVectorField>(Uname))
             );
 
             // Recalculate new internal faces.

From f9b4027b05a0f77e59ad3e1d7ec5346bd386e50b Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Thu, 21 Dec 2023 09:10:07 +0000
Subject: [PATCH 2/4] vtkWriteOps: Protect against float overflow

---
 src/fileFormats/vtk/vtkWriteOps.C             | 26 +++++-----
 src/fileFormats/vtk/vtkWriteOps.H             | 19 +++++--
 src/fileFormats/vtk/vtkWriteOpsI.H            | 49 +++++++++++++++++++
 .../vtk/vtkWritePolyDataTemplates.C           | 13 +++--
 4 files changed, 86 insertions(+), 21 deletions(-)
 create mode 100644 src/fileFormats/vtk/vtkWriteOpsI.H

diff --git a/src/fileFormats/vtk/vtkWriteOps.C b/src/fileFormats/vtk/vtkWriteOps.C
index 67c2216d3d..a3dffc236c 100644
--- a/src/fileFormats/vtk/vtkWriteOps.C
+++ b/src/fileFormats/vtk/vtkWriteOps.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -90,7 +90,7 @@ void Foam::vtkWriteOps::write
         os.write
         (
             reinterpret_cast<char*>(fField.begin()),
-            fField.size()*sizeof(float)
+            fField.size()*sizeof(floatScalar)
         );
 
         os  << std::endl;
@@ -228,7 +228,7 @@ void Foam::vtkWriteOps::writePointDataHeader
 
 void Foam::vtkWriteOps::insert(const scalar src, DynamicList<floatScalar>& dest)
 {
-    dest.append(float(src));
+    dest.append(cast(src));
 }
 
 
@@ -240,7 +240,7 @@ void Foam::vtkWriteOps::insert
 {
     for (direction cmpt = 0; cmpt < vector::nComponents; ++cmpt)
     {
-        dest.append(float(src[cmpt]));
+        dest.append(cast(src[cmpt]));
     }
 }
 
@@ -253,7 +253,7 @@ void Foam::vtkWriteOps::insert
 {
     for (direction cmpt = 0; cmpt < sphericalTensor::nComponents; ++cmpt)
     {
-        dest.append(float(src[cmpt]));
+        dest.append(cast(src[cmpt]));
     }
 }
 
@@ -264,12 +264,12 @@ void Foam::vtkWriteOps::insert
     DynamicList<floatScalar>& dest
 )
 {
-    dest.append(float(src.xx()));
-    dest.append(float(src.yy()));
-    dest.append(float(src.zz()));
-    dest.append(float(src.xy()));
-    dest.append(float(src.yz()));
-    dest.append(float(src.xz()));
+    dest.append(cast(src.xx()));
+    dest.append(cast(src.yy()));
+    dest.append(cast(src.zz()));
+    dest.append(cast(src.xy()));
+    dest.append(cast(src.yz()));
+    dest.append(cast(src.xz()));
 }
 
 
@@ -281,7 +281,7 @@ void Foam::vtkWriteOps::insert
 {
     for (direction cmpt = 0; cmpt < tensor::nComponents; ++cmpt)
     {
-        dest.append(float(src[cmpt]));
+        dest.append(cast(src[cmpt]));
     }
 }
 
@@ -301,7 +301,7 @@ void Foam::vtkWriteOps::insert
 {
     forAll(map, i)
     {
-        dest.append(float(source[map[i]]));
+        dest.append(cast(source[map[i]]));
     }
 }
 
diff --git a/src/fileFormats/vtk/vtkWriteOps.H b/src/fileFormats/vtk/vtkWriteOps.H
index 494c54f8fc..5643d3b8ff 100644
--- a/src/fileFormats/vtk/vtkWriteOps.H
+++ b/src/fileFormats/vtk/vtkWriteOps.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,6 +58,17 @@ namespace vtkWriteOps
     //- Swap halves of word
     void swapWords(const label nWords, label* words32);
 
+    //- Cast an integer to the corresponding VTK write type. Does nothing.
+    inline label cast(const label& x);
+
+    //- Cast a float to the corresponding VTK write type. Does nothing.
+    inline floatScalar cast(const floatScalar& x);
+
+    //- Cast a scalar to the corresponding VTK write type. Clips to
+    //  pTraits<floatScalar>::max.
+    template<class Scalar>
+    inline floatScalar cast(const Scalar& x);
+
     //- Write header
     void writeHeader
     (
@@ -66,6 +77,7 @@ namespace vtkWriteOps
         const std::string& title
     );
 
+    //- Write cell data header
     void writeCellDataHeader
     (
         std::ostream&,
@@ -73,6 +85,7 @@ namespace vtkWriteOps
         const label nFields
     );
 
+    //- Write point data header
     void writePointDataHeader
     (
         std::ostream&,
@@ -80,7 +93,6 @@ namespace vtkWriteOps
         const label nFields
     );
 
-
     //- Write floats ascii or binary.
     //  If binary optionally in-place swaps argument
     void write(std::ostream& os, const bool binary, List<floatScalar>& fField);
@@ -97,7 +109,6 @@ namespace vtkWriteOps
     //  If binary optionally in-place swaps argument
     void write(std::ostream&, const bool, DynamicList<label>&);
 
-
     //- Append scalar to given DynamicList
     void insert(const scalar, DynamicList<floatScalar>&);
 
@@ -148,6 +159,8 @@ namespace vtkWriteOps
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "vtkWriteOpsI.H"
+
 #ifdef NoRepository
     #include "vtkWriteOpsTemplates.C"
 #endif
diff --git a/src/fileFormats/vtk/vtkWriteOpsI.H b/src/fileFormats/vtk/vtkWriteOpsI.H
new file mode 100644
index 0000000000..c48b14a8a6
--- /dev/null
+++ b/src/fileFormats/vtk/vtkWriteOpsI.H
@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "vtkWriteOps.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline Foam::label Foam::vtkWriteOps::cast(const label& x)
+{
+    return x;
+}
+
+
+inline Foam::floatScalar Foam::vtkWriteOps::cast(const floatScalar& x)
+{
+    return x;
+}
+
+
+template<class Scalar>
+inline Foam::floatScalar Foam::vtkWriteOps::cast(const Scalar& x)
+{
+    return sign(x)*floatScalar(min(mag(x), floatScalarVGreat));
+}
+
+
+// ************************************************************************* //
diff --git a/src/fileFormats/vtk/vtkWritePolyDataTemplates.C b/src/fileFormats/vtk/vtkWritePolyDataTemplates.C
index 94596eebeb..f3ba744ba9 100644
--- a/src/fileFormats/vtk/vtkWritePolyDataTemplates.C
+++ b/src/fileFormats/vtk/vtkWritePolyDataTemplates.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2021-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -60,10 +60,13 @@ void Foam::vtkWritePolyData::writeFieldTypeValues
                 for (direction cmpt = 0; cmpt < nCmpt; ++ cmpt)
                 {
                     data[i ++] =
-                        component
+                        vtkWriteOps::cast
                         (
-                            fieldTypeValues[fieldi][fieldValuei],
-                            cmpt
+                            component
+                            (
+                                fieldTypeValues[fieldi][fieldValuei],
+                                cmpt
+                            )
                         );
                 }
             }
@@ -109,7 +112,7 @@ void Foam::vtkWritePolyData::write
             const point& p = points[pointi];
             forAll(p, i)
             {
-                coordinates[3*pointi + i] = float(p[i]);
+                coordinates[3*pointi + i] = vtkWriteOps::cast(p[i]);
             }
         }
         vtkWriteOps::write(os, binary, coordinates);

From 705a4999bd9a2d01a0802081836a120b8f235cee Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Thu, 21 Dec 2023 10:37:03 +0000
Subject: [PATCH 3/4] triIntersect: Stabilise source-point-to-target-triangle
 projections

---
 src/meshTools/triIntersect/triIntersect.C | 24 ++++++++++++++++++++++-
 1 file changed, 23 insertions(+), 1 deletion(-)

diff --git a/src/meshTools/triIntersect/triIntersect.C b/src/meshTools/triIntersect/triIntersect.C
index c280069f52..6f3103b9a6 100644
--- a/src/meshTools/triIntersect/triIntersect.C
+++ b/src/meshTools/triIntersect/triIntersect.C
@@ -1251,11 +1251,33 @@ Foam::barycentric2D Foam::triIntersect::srcPointTgtTriIntersection
 
     const scalar maxMagDetAY = mag(detAY[findMax(cmptMag(detAY))]);
 
-    const vector y =
+    vector y =
         maxMagDetAY/vGreat < mag(detA)
       ? detAY/detA
       : detAY/maxMagDetAY*vGreat;
 
+    // Project into the source triangle
+    if (cmptMin(y) < 0)
+    {
+        const direction iMin = findMin(y);
+        const direction iMax = findMax(y);
+        const direction iMid = 3 - iMin - iMax;
+
+        if (y[iMid] < 0)
+        {
+            y[iMin] = 0;
+            y[iMax] = 1;
+            y[iMid] = 0;
+        }
+        else
+        {
+            const scalar t = y[iMax] + y[iMid];
+            y[iMin] = 0;
+            y[iMax] /= t;
+            y[iMid] /= t;
+        }
+    }
+
     return barycentric2D(y.x(), y.y(), y.z());
 }
 

From d600bd11e42cdd29aa00995ad4da5032c3027755 Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Thu, 21 Dec 2023 12:57:55 +0000
Subject: [PATCH 4/4] triIntersect: Move stabilisation in to geometry
 construction

Source-point-to-target-triangle projections might be done in order to
determine the location outside of the target triangle, as well as the
location inside. Forcing the source point inside the target triangle is
therefore only appropriate when constructing the intersection geometry.
---
 src/meshTools/triIntersect/triIntersect.C | 46 +++++++++++------------
 1 file changed, 23 insertions(+), 23 deletions(-)

diff --git a/src/meshTools/triIntersect/triIntersect.C b/src/meshTools/triIntersect/triIntersect.C
index 6f3103b9a6..6ab0e490b1 100644
--- a/src/meshTools/triIntersect/triIntersect.C
+++ b/src/meshTools/triIntersect/triIntersect.C
@@ -641,9 +641,31 @@ void generateGeometryForLocations
         {
             const point& srcP = srcPs[l.srcPointi()];
             const vector& srcN = srcNs[l.srcPointi()];
-            const barycentric2D tgtTs =
+            barycentric2D tgtTs =
                 srcPointTgtTriIntersection(srcP, srcN, tgtPs);
 
+            // Force inside the target triangle
+            if (cmptMin(tgtTs) < 0)
+            {
+                const direction iMin = findMin(tgtTs);
+                const direction iMax = findMax(tgtTs);
+                const direction iMid = 3 - iMin - iMax;
+
+                if (tgtTs[iMid] < 0)
+                {
+                    tgtTs[iMin] = 0;
+                    tgtTs[iMax] = 1;
+                    tgtTs[iMid] = 0;
+                }
+                else
+                {
+                    const scalar t = tgtTs[iMax] + tgtTs[iMid];
+                    tgtTs[iMin] = 0;
+                    tgtTs[iMax] /= t;
+                    tgtTs[iMid] /= t;
+                }
+            }
+
             srcPoints[pointi] = srcP;
             srcPointNormals[pointi] = srcN;
             tgtPoints[pointi] = tgtTriInterpolate(tgtTs, tgtPs);
@@ -1256,28 +1278,6 @@ Foam::barycentric2D Foam::triIntersect::srcPointTgtTriIntersection
       ? detAY/detA
       : detAY/maxMagDetAY*vGreat;
 
-    // Project into the source triangle
-    if (cmptMin(y) < 0)
-    {
-        const direction iMin = findMin(y);
-        const direction iMax = findMax(y);
-        const direction iMid = 3 - iMin - iMax;
-
-        if (y[iMid] < 0)
-        {
-            y[iMin] = 0;
-            y[iMax] = 1;
-            y[iMid] = 0;
-        }
-        else
-        {
-            const scalar t = y[iMax] + y[iMid];
-            y[iMin] = 0;
-            y[iMax] /= t;
-            y[iMid] /= t;
-        }
-    }
-
     return barycentric2D(y.x(), y.y(), y.z());
 }