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);
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.
diff --git a/src/meshTools/triIntersect/triIntersect.C b/src/meshTools/triIntersect/triIntersect.C
index c280069f52..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);
@@ -1251,7 +1273,7 @@ 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;