diff --git a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
index 198be01bcf..a6547fa114 100644
--- a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
+++ b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -168,12 +168,7 @@ void Foam::shellSurfaces::orient()
     {
         const searchableSurface& s = allGeometry_[shells_[shellI]];
 
-        if
-        (
-            modes_[shellI] != DISTANCE
-        &&  isA<triSurfaceMesh>(s)
-        && !isA<distributedTriSurfaceMesh>(s)
-        )
+        if (modes_[shellI] != DISTANCE && isA<triSurfaceMesh>(s))
         {
             hasSurface = true;
 
@@ -196,35 +191,39 @@ void Foam::shellSurfaces::orient()
         {
             const searchableSurface& s = allGeometry_[shells_[shellI]];
 
-            if
-            (
-                modes_[shellI] != DISTANCE
-            &&  isA<triSurfaceMesh>(s)
-            && !isA<distributedTriSurfaceMesh>(s)
-            )
+            if (modes_[shellI] != DISTANCE && isA<triSurfaceMesh>(s))
             {
-                triSurfaceMesh& shell = const_cast<triSurfaceMesh&>
+                List<pointIndexHit> info;
+                vectorField normal;
+                labelList region;
+                s.findNearest
                 (
-                    refCast<const triSurfaceMesh>(s)
+                    pointField(1, outsidePt),
+                    scalarField(1, GREAT),
+                    info,
+                    normal,
+                    region
                 );
 
-                // Flip surface so outsidePt is outside.
-                bool anyFlipped = orientedSurface::orient
-                (
-                    shell,
-                    outsidePt,
-                    true
-                );
+                //Pout<< "outsidePt:" << outsidePt << endl;
+                //Pout<< "info     :" << info[0] << endl;
+                //Pout<< "normal   :" << normal[0] << endl;
+                //Pout<< "region   :" << region[0] << endl;
+
+                bool anyFlipped = false;
+                if ((normal[0] & (info[0].point()-outsidePt)) > 0)
+                {
+                    triSurfaceMesh& shell = const_cast<triSurfaceMesh&>
+                    (
+                        refCast<const triSurfaceMesh>(s)
+                    );
+                    shell.flip();
+                    anyFlipped = true;
+                }
+
 
                 if (anyFlipped && !dryRun_)
                 {
-                    // orientedSurface will have done a clearOut of the surface.
-                    // we could do a clearout of the triSurfaceMeshes::trees()
-                    // but these aren't affected by orientation
-                    // (except for cached
-                    // sideness which should not be set at this point.
-                    // !!Should check!)
-
                     Info<< "shellSurfaces : Flipped orientation of surface "
                         << s.name()
                         << " so point " << outsidePt << " is outside." << endl;
diff --git a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
index 8976887905..4d13ca4deb 100644
--- a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
+++ b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2020,2022 OpenCFD Ltd.
+    Copyright (C) 2015-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -101,6 +101,7 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const
     if (debug)
     {
         Pout<< "triSurfaceMesh::isSurfaceClosed:"
+            << " surface:" << searchableSurface::name()
             << " determining closedness for surface with "
             << triSurface::size() << " triangles" << endl;
     }
@@ -223,6 +224,7 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const
                 if (debug)
                 {
                     Pout<< "triSurfaceMesh::isSurfaceClosed :"
+                        << " surface:" << searchableSurface::name()
                         << " surface is non-manifold" << endl;
                 }
                 return false;
@@ -240,6 +242,7 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const
             if (debug)
             {
                 Pout<< "triSurfaceMesh::isSurfaceClosed :"
+                    << " surface:" << searchableSurface::name()
                     << " surface is open" << endl;
             }
             return false;
@@ -274,6 +277,7 @@ bool Foam::triSurfaceMesh::isSurfaceClosed() const
     if (debug)
     {
         Pout<< "triSurfaceMesh::isSurfaceClosed :"
+            << " surface:" << searchableSurface::name()
             << " surface is closed" << endl;
     }
     return true;
@@ -738,6 +742,7 @@ void Foam::triSurfaceMesh::movePoints(const pointField& newPoints)
     if (debug)
     {
         Pout<< "triSurfaceMesh::movePoints :"
+            << " surface:" << searchableSurface::name()
             << " moving at time " << objectRegistry::time().timeName()
             << endl;
     }
@@ -773,6 +778,7 @@ Foam::triSurfaceMesh::edgeTree() const
         if (debug)
         {
             Pout<< "triSurfaceMesh::edgeTree :"
+                << " surface:" << searchableSurface::name()
                 << " constructing tree for " << nEdges() - nInternalEdges()
                 << " boundary edges" << endl;
         }
@@ -882,6 +888,7 @@ Foam::volumeType Foam::triSurfaceMesh::outsideVolumeType() const
         if (debug)
         {
             Pout<< "triSurfaceMesh::outsideVolumeType :"
+                << " surface:" << searchableSurface::name()
                 << " triggering outsidePoint:" << outsidePt
                 << " orientation" << endl;
         }
@@ -906,6 +913,49 @@ Foam::volumeType Foam::triSurfaceMesh::outsideVolumeType() const
 }
 
 
+void Foam::triSurfaceMesh::flip()
+{
+    if (debug)
+    {
+        Pout<< "triSurfaceMesh::flip :"
+            << " surface:" << searchableSurface::name()
+            << " with current orientation "
+            << volumeType::names[outsideVolType_]
+            << endl;
+    }
+
+    // Don't bother getting nearest etc. Just flip the triangles.
+
+    // triSurface
+    {
+        triSurface& s = *this;
+        s.clearOut();
+        for (auto& tri : s)
+        {
+            tri.flip();
+        }
+    }
+
+    // triSurfaceRegionSearch (if cached volume type)
+    triSurfaceRegionSearch::flip();
+
+    // edge tree not relevant
+
+    if (hasVolumeType())
+    {
+        // outsideVolType_
+        if (outsideVolType_ == volumeType::INSIDE)
+        {
+            outsideVolType_ = volumeType::OUTSIDE;
+        }
+        else if (outsideVolType_ == volumeType::OUTSIDE)
+        {
+            outsideVolType_ = volumeType::INSIDE;
+        }
+    }
+}
+
+
 void Foam::triSurfaceMesh::findNearest
 (
     const pointField& samples,
@@ -916,6 +966,7 @@ void Foam::triSurfaceMesh::findNearest
     if (debug)
     {
         Pout<< "triSurfaceMesh::findNearest :"
+            << " surface:" << searchableSurface::name()
             << " trying to find nearest for " << samples.size()
             << " samples with max sphere "
             << (samples.size() ? Foam::sqrt(max(nearestDistSqr)) : Zero)
@@ -942,6 +993,7 @@ void Foam::triSurfaceMesh::findNearest
     if (debug)
     {
         Pout<< "triSurfaceMesh::findNearest :"
+            << " surface:" << searchableSurface::name()
             << " trying to find nearest and region for " << samples.size()
             << " samples with max sphere "
             << (samples.size() ? Foam::sqrt(max(nearestDistSqr)) : Zero)
@@ -973,6 +1025,7 @@ void Foam::triSurfaceMesh::findLine
     if (debug)
     {
         Pout<< "triSurfaceMesh::findLine :"
+            << " surface:" << searchableSurface::name()
             << " intersecting with "
             << start.size() << " rays" << endl;
     }
@@ -996,6 +1049,7 @@ void Foam::triSurfaceMesh::findLineAny
     if (debug)
     {
         Pout<< "triSurfaceMesh::findLineAny :"
+            << " surface:" << searchableSurface::name()
             << " intersecting with "
             << start.size() << " rays" << endl;
     }
@@ -1019,6 +1073,7 @@ void Foam::triSurfaceMesh::findLineAll
     if (debug)
     {
         Pout<< "triSurfaceMesh::findLineAll :"
+            << " surface:" << searchableSurface::name()
             << " intersecting with "
             << start.size() << " rays" << endl;
     }
@@ -1041,6 +1096,7 @@ void Foam::triSurfaceMesh::getRegion
     if (debug)
     {
         Pout<< "triSurfaceMesh::getRegion :"
+            << " surface:" << searchableSurface::name()
             << " getting region for "
             << info.size() << " triangles" << endl;
     }
@@ -1074,6 +1130,7 @@ void Foam::triSurfaceMesh::getNormal
     if (debug)
     {
         Pout<< "triSurfaceMesh::getNormal :"
+            << " surface:" << searchableSurface::name()
             << " getting normal for "
             << info.size() << " triangles" << endl;
     }
@@ -1182,6 +1239,7 @@ void Foam::triSurfaceMesh::setField(const labelList& values)
     if (debug)
     {
         Pout<< "triSurfaceMesh::setField :"
+            << " surface:" << searchableSurface::name()
             << " finished setting field for "
             << values.size() << " triangles" << endl;
     }
@@ -1213,6 +1271,7 @@ void Foam::triSurfaceMesh::getField
     if (debug)
     {
         Pout<< "triSurfaceMesh::setField :"
+            << " surface:" << searchableSurface::name()
             << " finished getting field for "
             << info.size() << " triangles" << endl;
     }
@@ -1231,6 +1290,7 @@ void Foam::triSurfaceMesh::getVolumeType
     if (debug)
     {
         Pout<< "triSurfaceMesh::getVolumeType :"
+            << " surface:" << searchableSurface::name()
             << " finding orientation for " << points.size()
             << " samples" << endl;
     }
diff --git a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H
index e5a18bdeb4..c833e24d2b 100644
--- a/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H
+++ b/src/meshTools/searchableSurfaces/triSurfaceMesh/triSurfaceMesh.H
@@ -216,6 +216,9 @@ public:
                 return triSurface::size();
             }
 
+            //- Flip triangles, outsideVolumeType and all cached inside/outside.
+            virtual void flip();
+
             //- Get representative set of element coordinates
             //  Usually the element centres (should be of length size()).
             virtual tmp<pointField> coordinates() const;
diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C
index cb8a67db5c..eae21fb06e 100644
--- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C
+++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.C
@@ -251,4 +251,26 @@ void Foam::triSurfaceRegionSearch::findNearest
 }
 
 
+void Foam::triSurfaceRegionSearch::flip()
+{
+    triSurfaceSearch::flip();
+
+    for (auto& tree : treeByRegion_)
+    {
+        PackedList<2>& nodeTypes = tree.nodeTypes();
+        forAll(nodeTypes, i)
+        {
+            if (nodeTypes[i] == volumeType::INSIDE)
+            {
+                nodeTypes[i] = volumeType::OUTSIDE;
+            }
+            else if (nodeTypes[i] == volumeType::OUTSIDE)
+            {
+                nodeTypes[i] = volumeType::INSIDE;
+            }
+        }
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.H b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.H
index 724617d132..54143930eb 100644
--- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.H
+++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.H
@@ -128,6 +128,11 @@ public:
                 const labelList& regionIndices,
                 List<pointIndexHit>& info
             ) const;
+
+        // Edit
+
+            //- Flip orientation
+            void flip();
 };
 
 
diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
index 33fa142a61..b80944a491 100644
--- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
+++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.C
@@ -255,6 +255,26 @@ Foam::triSurfaceSearch::tree() const
 }
 
 
+void Foam::triSurfaceSearch::flip()
+{
+    if (treePtr_)
+    {
+        PackedList<2>& nodeTypes = treePtr_->nodeTypes();
+        forAll(nodeTypes, i)
+        {
+            if (nodeTypes[i] == volumeType::INSIDE)
+            {
+                nodeTypes[i] = volumeType::OUTSIDE;
+            }
+            else if (nodeTypes[i] == volumeType::OUTSIDE)
+            {
+                nodeTypes[i] = volumeType::INSIDE;
+            }
+        }
+    }
+}
+
+
 // Determine inside/outside for samples
 Foam::boolList Foam::triSurfaceSearch::calcInside
 (
diff --git a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H
index a4ce7c8c98..2deb9abf9e 100644
--- a/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H
+++ b/src/meshTools/triSurface/triSurfaceSearch/triSurfaceSearch.H
@@ -126,6 +126,9 @@ public:
         //- Demand driven construction of the octree
         const indexedOctree<treeDataTriSurface>& tree() const;
 
+        //- Flip orientation (if cached on octree)
+        void flip();
+
         //- Return reference to the surface.
         const triSurface& surface() const
         {
diff --git a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
index 8c32aface9..bdb63f152e 100644
--- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
+++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
@@ -110,6 +110,7 @@ namespace Foam
     };
 
 
+    ////- Same as above but with some more debugging
     //typedef Tuple2<Pair<point>, volumeType> NearType;
     //
     ////- Combine operator for volume types
@@ -337,73 +338,145 @@ Foam::word Foam::distributedTriSurfaceMesh::findLocalInstance
 }
 
 
-// Read my additional data from the dictionary
-bool Foam::distributedTriSurfaceMesh::read()
+bool Foam::distributedTriSurfaceMesh::readSettings(const bool undecomposed)
 {
     // Get bb of all domains.
     procBb_.resize_nocopy(Pstream::nProcs());
 
-    if (dict_.empty())
+    if (!dict_.readIfPresent("bounds", procBb_[Pstream::myProcNo()]))
     {
-        // Did not find the distributed version; assume master has loaded the
-        // triSurfaceMesh version. Make up some settings.
-
+        // Not found. Use local triangles' bounding box
         const boundBox& localBb = triSurfaceMesh::bounds();
 
         procBb_[Pstream::myProcNo()] =
             treeBoundBoxList(1, treeBoundBox(localBb));
 
         dict_.add("bounds", procBb_[Pstream::myProcNo()]);
+    }
 
-        // Wanted distribution type
+
+    decomposeUsingBbs_ = dict_.getOrDefault<bool>("decomposeUsingBbs", true);
+
+    // Wanted distribution type
+    if
+    (
+       !distributionTypeNames_.readIfPresent
+        (
+            "distributionType",
+            dict_,
+            distType_
+        )
+    )
+    {
+        // Not found. Make up most sensible data
         distType_ = DISTRIBUTED;    //INDEPENDENT;
         dict_.add("distributionType", distributionTypeNames_[distType_]);
+    }
 
-        // Merge distance
+    // Merge distance
+    if (!dict_.readIfPresent("mergeDistance", mergeDist_))
+    {
         mergeDist_ = SMALL;
         dict_.add("mergeDistance", mergeDist_);
-
-        // Force underlying triSurfaceMesh to calculate volume type
-        // (is topological walk; does not construct tree)
-        surfaceClosed_ = triSurfaceMesh::hasVolumeType();
-        Pstream::broadcast(surfaceClosed_);
-        dict_.add("closed", surfaceClosed_);
-
-        // Delay calculating outside vol type since constructs tree. Is ok
-        // after distributing since then local surfaces much smaller
-        //outsideVolType_ = volumeType::UNKNOWN;
-        //if (surfaceClosed_)
-        //{
-        //    point outsidePt(localBb.max()+localBb.centre());
-        //    List<volumeType> outsideVolTypes;
-        //    triSurfaceMesh::getVolumeType
-        //    (
-        //        pointField(1, outsidePt),
-        //        outsideVolTypes
-        //    );
-        //    outsideVolType_ = outsideVolTypes[0];
-        //}
-        //dict_.add("outsideVolumeType", volumeType::names[outsideVolType_]);
     }
-    else
+
+    // Is closed?
+    if (!dict_.readIfPresent<label>("closed", surfaceClosed_))
     {
-        dict_.readEntry("bounds", procBb_[Pstream::myProcNo()]);
+        if (undecomposed)
+        {
+            surfaceClosed_ = triSurfaceMesh::hasVolumeType();
+            Pstream::broadcast(surfaceClosed_);
+        }
+        else
+        {
+            surfaceClosed_ = false;
+        }
+        dict_.add("closed", surfaceClosed_);
+    }
 
-        // Wanted distribution type
-        distType_ = distributionTypeNames_.get("distributionType", dict_);
-
-        // Merge distance
-        dict_.readEntry("mergeDistance", mergeDist_);
-
-        // Distribution type
-        surfaceClosed_ = dict_.getOrDefault<bool>("closed", false);
-
-        outsideVolType_ = volumeType::names.getOrDefault
+    // Do we know volume inside/outside?
+    volumeType::type vt;
+    if
+    (
+        volumeType::names.readIfPresent
         (
             "outsideVolumeType",
             dict_,
-            volumeType::UNKNOWN
-        );
+            vt
+        )
+    )
+    {
+        outsideVolType_ = vt;
+    }
+    else
+    {
+        if (undecomposed && surfaceClosed_)
+        {
+            if (UPstream::master())
+            {
+                // Determine nearest to do inside/outside testing.
+                // Do linear search to avoid building tree on this large
+                // undecomposed surface. See e.g. triSurfaceTools::surfaceSide
+
+                const triSurface& surf = *this;
+                const pointField& points = surf.points();
+                const boundBox& localBb = triSurfaceMesh::bounds();
+                const point outsidePt(localBb.max()+localBb.centre());
+
+                label nearTrii = -1;
+                pointHit nearInfo;  // distance=GREAT
+                forAll(surf, trii)
+                {
+                    pointHit info = surf[trii].nearestPoint(outsidePt, points);
+
+                    if (info.distance() < nearInfo.distance())
+                    {
+                        nearTrii = trii;
+                        nearInfo = info;
+                    }
+                }
+
+                if (nearTrii != -1)
+                {
+                    const triSurfaceTools::sideType t =
+                        triSurfaceTools::surfaceSide
+                        (
+                            surf,
+                            outsidePt,
+                            nearTrii
+                        );
+
+                    if (t == triSurfaceTools::UNKNOWN)
+                    {
+                        outsideVolType_ = volumeType::UNKNOWN;
+                    }
+                    else if (t == triSurfaceTools::INSIDE)
+                    {
+                        outsideVolType_ = volumeType::INSIDE;
+                    }
+                    else if (t == triSurfaceTools::OUTSIDE)
+                    {
+                        outsideVolType_ = volumeType::OUTSIDE;
+                    }
+                    else
+                    {
+                        FatalErrorInFunction << "problem" << abort(FatalError);
+                    }
+                }
+            }
+            // Scatter master
+            label volType = outsideVolType_;
+            Pstream::broadcast(volType);
+            outsideVolType_ = volumeType(volType);
+        }
+        else
+        {
+            // Mark as to-be-done
+            outsideVolType_ = volumeType::UNKNOWN;
+        }
+
+        dict_.add("outsideVolumeType", volumeType::names[outsideVolType_]);
     }
 
     Pstream::allGatherList(procBb_);
@@ -412,6 +485,41 @@ bool Foam::distributedTriSurfaceMesh::read()
 }
 
 
+void Foam::distributedTriSurfaceMesh::calcVertexNormals
+(
+    const triSurface& surf,
+    List<FixedList<vector, 3>>& vn
+) const
+{
+    // For now don't do features so just use the point normals
+    vectorField pointNormals(surf.points().size(), vector::zero);
+    forAll(surf, facei)
+    {
+        const triSurface::face_type& f = surf[facei];
+        const vector n = f.unitNormal(surf.points());
+        forAll(f, fp)
+        {
+            const label pointi = f[fp];
+            pointNormals[pointi] += n;
+        }
+    }
+    pointNormals.normalise();
+
+
+    vn.resize_nocopy(surf.size());
+    vn = Zero;
+
+    forAll(surf, facei)
+    {
+        const triSurface::face_type& f = surf[facei];
+        forAll(f, fp)
+        {
+            vn[facei][fp] = pointNormals[f[fp]];
+        }
+    }
+}
+
+
 // Is segment fully local?
 bool Foam::distributedTriSurfaceMesh::isLocal
 (
@@ -504,28 +612,31 @@ void Foam::distributedTriSurfaceMesh::distributeSegment
     List<DynamicList<label>>& sendMap
 ) const
 {
-    // 1. Fully local already handled outside. Note: retest is cheap.
-    if (isLocal(procBb_[Pstream::myProcNo()], start, end))
+    if (decomposeUsingBbs_)
     {
-        return;
-    }
-
-
-    // 2. If fully inside one other processor, then only need to send
-    // to that one processor even if it intersects another. Rare occurrence
-    // but cheap to test.
-    forAll(procBb_, proci)
-    {
-        if (proci != Pstream::myProcNo())
+        // 1. Fully local already handled outside. Note: retest is cheap.
+        if (isLocal(procBb_[Pstream::myProcNo()], start, end))
         {
-            const List<treeBoundBox>& bbs = procBb_[proci];
+            return;
+        }
 
-            if (isLocal(bbs, start, end))
+
+        // 2. If fully inside one other processor, then only need to send
+        // to that one processor even if it intersects another. Rare occurrence
+        // but cheap to test.
+        forAll(procBb_, proci)
+        {
+            if (proci != Pstream::myProcNo())
             {
-                sendMap[proci].append(allSegments.size());
-                allSegmentMap.append(segmenti);
-                allSegments.append(segment(start, end));
-                return;
+                const List<treeBoundBox>& bbs = procBb_[proci];
+
+                if (isLocal(bbs, start, end))
+                {
+                    sendMap[proci].append(allSegments.size());
+                    allSegmentMap.append(segmenti);
+                    allSegments.append(segment(start, end));
+                    return;
+                }
             }
         }
     }
@@ -644,7 +755,7 @@ void Foam::distributedTriSurfaceMesh::findLine
 (
     const bool nearestIntersection,
     const pointField& start,
-    const pointField& end,
+    const pointField& initialEnd,
     List<pointIndexHit>& info
 ) const
 {
@@ -659,6 +770,9 @@ void Foam::distributedTriSurfaceMesh::findLine
 
     const indexedOctree<treeDataTriSurface>& octree = tree();
 
+    // end-point - updated with local hits
+    pointField end(initialEnd);
+
     // Initialise
     info.setSize(start.size());
     forAll(info, i)
@@ -691,6 +805,8 @@ void Foam::distributedTriSurfaceMesh::findLine
             if (info[i].hit())
             {
                 info[i].setIndex(triIndexer.toGlobal(info[i].index()));
+                // Update endpoint
+                end[i] = info[i].hitPoint();
             }
             nLocal++;
         }
@@ -699,113 +815,121 @@ void Foam::distributedTriSurfaceMesh::findLine
 
     if
     (
-        returnReduce(nLocal, sumOp<label>())
-      < returnReduce(start.size(), sumOp<label>())
+        decomposeUsingBbs_
+     && (
+            returnReduce(nLocal, sumOp<label>())
+         == returnReduce(start.size(), sumOp<label>())
+        )
     )
     {
-        // Not all can be resolved locally. Build segments and map,
-        // send over segments, do intersections, send back and merge.
+        // In decomposeUsingBbs_ mode we're guaranteed to have all the geometry
+        // inside the local bb. Since we've done all tests we're done.
+        return;
+    }
 
 
-        // Construct queries (segments)
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Not all can be resolved locally. Build segments and map,
+    // send over segments, do intersections, send back and merge.
 
-        // Segments to test
-        List<segment> allSegments(start.size());
-        // Original index of segment
-        labelList allSegmentMap(start.size());
 
-        const autoPtr<mapDistribute> mapPtr
+    // Construct queries (segments)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Segments to test
+    List<segment> allSegments(start.size());
+    // Original index of segment
+    labelList allSegmentMap(start.size());
+
+    const autoPtr<mapDistribute> mapPtr
+    (
+        distributeSegments
         (
-            distributeSegments
-            (
-                start,
-                end,
-                allSegments,
-                allSegmentMap
-            )
-        );
-        const mapDistribute& map = mapPtr();
+            start,
+            end,
+            allSegments,
+            allSegmentMap
+        )
+    );
+    const mapDistribute& map = mapPtr();
 
-        label nOldAllSegments = allSegments.size();
+    label nOldAllSegments = allSegments.size();
 
 
-        // Exchange the segments
-        // ~~~~~~~~~~~~~~~~~~~~~
+    // Exchange the segments
+    // ~~~~~~~~~~~~~~~~~~~~~
 
-        map.distribute(allSegments);
+    map.distribute(allSegments);
 
 
-        // Do tests i need to do
-        // ~~~~~~~~~~~~~~~~~~~~~
+    // Do tests i need to do
+    // ~~~~~~~~~~~~~~~~~~~~~
 
-        // Intersections
-        List<pointIndexHit> intersections(allSegments.size());
+    // Intersections
+    List<pointIndexHit> intersections(allSegments.size());
 
-        forAll(allSegments, i)
+    forAll(allSegments, i)
+    {
+        if (nearestIntersection)
         {
-            if (nearestIntersection)
-            {
-                intersections[i] = octree.findLine
-                (
-                    allSegments[i].first(),
-                    allSegments[i].second()
-                );
-            }
-            else
-            {
-                intersections[i] = octree.findLineAny
-                (
-                    allSegments[i].first(),
-                    allSegments[i].second()
-                );
-            }
-
-            // Convert triangle index to global numbering
-            if (intersections[i].hit())
-            {
-                intersections[i].setIndex
-                (
-                    triIndexer.toGlobal(intersections[i].index())
-                );
-            }
+            intersections[i] = octree.findLine
+            (
+                allSegments[i].first(),
+                allSegments[i].second()
+            );
+        }
+        else
+        {
+            intersections[i] = octree.findLineAny
+            (
+                allSegments[i].first(),
+                allSegments[i].second()
+            );
         }
 
-
-        // Exchange the intersections (opposite to segments)
-        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
-        map.reverseDistribute(nOldAllSegments, intersections);
-
-
-        // Extract the hits
-        // ~~~~~~~~~~~~~~~~
-
-        forAll(intersections, i)
+        // Convert triangle index to global numbering
+        if (intersections[i].hit())
         {
-            const pointIndexHit& allInfo = intersections[i];
-            label segmenti = allSegmentMap[i];
-            pointIndexHit& hitInfo = info[segmenti];
+            intersections[i].setIndex
+            (
+                triIndexer.toGlobal(intersections[i].index())
+            );
+        }
+    }
 
-            if (allInfo.hit())
+
+    // Exchange the intersections (opposite to segments)
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    map.reverseDistribute(nOldAllSegments, intersections);
+
+
+    // Extract the hits
+    // ~~~~~~~~~~~~~~~~
+
+    forAll(intersections, i)
+    {
+        const pointIndexHit& allInfo = intersections[i];
+        label segmenti = allSegmentMap[i];
+        pointIndexHit& hitInfo = info[segmenti];
+
+        if (allInfo.hit())
+        {
+            if (!hitInfo.hit())
             {
-                if (!hitInfo.hit())
+                // No intersection yet so take this one
+                hitInfo = allInfo;
+            }
+            else if (nearestIntersection)
+            {
+                // Nearest intersection
+                if
+                (
+                    start[segmenti].distSqr(allInfo.point())
+                  < start[segmenti].distSqr(hitInfo.point())
+                )
                 {
-                    // No intersection yet so take this one
                     hitInfo = allInfo;
                 }
-                else if (nearestIntersection)
-                {
-                    // Nearest intersection
-                    if
-                    (
-                        start[segmenti].distSqr(allInfo.point())
-                      < start[segmenti].distSqr(hitInfo.point())
-                    )
-                    {
-                        hitInfo = allInfo;
-                    }
-                }
             }
         }
     }
@@ -1154,7 +1278,7 @@ Foam::volumeType Foam::distributedTriSurfaceMesh::edgeSide
     const label face1
 ) const
 {
-    const triSurface& surf = static_cast<const triSurface&>(*this);
+    const triSurface& surf = *this;
     const pointField& points = surf.points();
 
     // Compare to bisector. This is actually correct since edge is
@@ -1182,7 +1306,7 @@ Foam::label Foam::distributedTriSurfaceMesh::findOtherFace
     const label nearLabel
 ) const
 {
-    const triSurface& surf = static_cast<const triSurface&>(*this);
+    const triSurface& surf = *this;
     const triSurface::face_type& nearF = surf[nearFacei];
 
     const edge e(nearF[nearLabel], nearF[nearF.fcIndex(nearLabel)]);
@@ -1255,7 +1379,8 @@ void Foam::distributedTriSurfaceMesh::surfaceSide
         Pout<< "distributedTriSurfaceMesh::surfaceSide :"
             << " on surface " << searchableSurface::name()
             << " finding surface side given points on surface for "
-            << samples.size() << " samples" << endl;
+            << samples.size() << " samples"
+            << " cached vertexnormals " << vertexNormals_.valid() << endl;
     }
 
     // Use global index to send local tri and nearest back to originating
@@ -1276,6 +1401,10 @@ void Foam::distributedTriSurfaceMesh::surfaceSide
     pointField localSamples(samples);
     map.distribute(localSamples);
 
+    // Send over nearest point since we have it anyway
+    List<pointIndexHit> localNearestInfo(nearestInfo);
+    map.distribute(localNearestInfo);
+
 
     // Do my tests
     // ~~~~~~~~~~~
@@ -1283,8 +1412,43 @@ void Foam::distributedTriSurfaceMesh::surfaceSide
     volType.setSize(triangleIndex.size());
     volType = volumeType::UNKNOWN;
 
-    const triSurface& surf = static_cast<const triSurface&>(*this);
+    const triSurface& surf = *this;
     const pointField& points = surf.points();
+
+
+    if (vertexNormals_)
+    {
+        // Use smooth interpolation for the normal
+
+        forAll(triangleIndex, i)
+        {
+            const label facei = triangleIndex[i];
+            const triSurface::face_type& f = surf[facei];
+            const auto& vn = vertexNormals_()[facei];
+
+            const point& sample = localSamples[i];
+
+            const point& nearest = localNearestInfo[i].point();
+
+            const vector sampleNearestVec(sample - nearest);
+
+            const barycentric2D w(f.tri(points).pointToBarycentric(nearest));
+
+            const vector n = w[0]*vn[0]+w[1]*vn[1]+w[2]*vn[2];
+
+            const scalar c = (sampleNearestVec & n);
+
+            if (c > 0)
+            {
+                volType[i] = volumeType::OUTSIDE;
+            }
+            else
+            {
+                volType[i] = volumeType::INSIDE;
+            }
+        }
+    }
+    else
     {
         //const labelListList& pointFaces = surf.pointFaces();
         // Construct pointFaces. Let's hope surface has compact point
@@ -1313,7 +1477,7 @@ void Foam::distributedTriSurfaceMesh::surfaceSide
 
                 // Nearest to face interior. Use faceNormal to determine side
                 //scalar c = sampleNearestVec & surf.faceNormals()[facei];
-                scalar c = sampleNearestVec & surf[facei].areaNormal(points);
+                const scalar c = (sampleNearestVec & f.areaNormal(points));
 
                 if (c > 0)
                 {
@@ -1451,13 +1615,45 @@ void Foam::distributedTriSurfaceMesh::surfaceSide
     // send the volume-query to a single processor)
 
 
-    //forAll(localSamples, i)
     //{
-    //    Pout<< "surfaceSide : for localSample:" << localSamples[i]
-    //        << " found volType:" << volumeType::names[volType[i]]
-    //        << endl;
+    //- Alternative which also sends over source sample so we get better
+    //- debug messages.
+    //    // Pack sample, nearest and volType so we get nicer error messages
+    //    typedef Tuple2<Pair<point>, volumeType> NearType;
+    //    List<NearType> nearTypes(volType.size());
+    //    forAll(localSamples, i)
+    //    {
+    //        const point& sample = localSamples[i];
+    //        const point& near = localNearestInfo[i].point();
+    //        nearTypes[i] = NearType(Pair<point>(sample, near), volType[i]);
+    //    }
+    //
+    //
+    //    const NearType zero(Pair<point>(Zero, Zero), volumeType::UNKNOWN);
+    //    mapDistributeBase::distribute
+    //    (
+    //        Pstream::commsTypes::nonBlocking,
+    //        List<labelPair>::null(),
+    //        nearestInfo.size(),
+    //        map.constructMap(),
+    //        map.constructHasFlip(),
+    //        map.subMap(),
+    //        map.subHasFlip(),
+    //        nearTypes,
+    //        zero,
+    //        NearTypeCombineOp(),
+    //        noOp(),           // no flipping
+    //        UPstream::msgType(),
+    //        map.comm()
+    //    );
+    //    volType.setSize(nearTypes.size());
+    //    forAll(nearTypes, i)
+    //    {
+    //        volType[i] = nearTypes[i].second();
+    //    }
     //}
 
+
     const volumeType zero(volumeType::UNKNOWN);
     mapDistributeBase::distribute
     (
@@ -1476,44 +1672,10 @@ void Foam::distributedTriSurfaceMesh::surfaceSide
         map.comm()
     );
 
-    //// Pack sample, nearest and volType so we get nicer error messages
-    ////typedef Tuple2<Pair<point>, volumeType>> NearType;
-    //List<NearType> nearTypes(volType.size());
-    //forAll(localSamples, i)
-    //{
-    //    const point& sample = localSamples[i];
-    //    const point& near = nearest[i];
-    //    nearTypes[i] = NearType(Pair<point>(sample, near), volType[i]);
-    //}
-    //
-    //
-    //const NearType zero(Pair<point>(Zero, Zero), volumeType::UNKNOWN);
-    //mapDistributeBase::distribute
-    //(
-    //    Pstream::commsTypes::nonBlocking,
-    //    List<labelPair>::null(),
-    //    nearestInfo.size(),
-    //    map.constructMap(),
-    //    map.constructHasFlip(),
-    //    map.subMap(),
-    //    map.subHasFlip(),
-    //    nearTypes,
-    //    zero,
-    //    NearTypeCombineOp(),
-    //    noOp(),           // no flipping
-    //    UPstream::msgType(),
-    //    map.comm()
-    //);
-    //volType.setSize(nearTypes.size());
-    //forAll(nearTypes, i)
-    //{
-    //    volType[i] = nearTypes[i].second();
-    //}
-
     if (debug)
     {
         Pout<< "distributedTriSurfaceMesh::surfaceSide :"
-            << " finished finding surface" << searchableSurface::name()
+            << " finished finding surface " << searchableSurface::name()
             << " given points on surface "
             << searchableSurface::name() << " for "
             << samples.size() << " samples" << endl;
@@ -1521,6 +1683,196 @@ void Foam::distributedTriSurfaceMesh::surfaceSide
 }
 
 
+Foam::volumeType Foam::distributedTriSurfaceMesh::markMixed
+(
+    const indexedOctree<treeDataTriSurface>& tree,
+    const label nodei,
+    const triPointRef& tri,
+    PackedList<2>& nodeTypes
+) const
+{
+    const auto& nod = tree.nodes()[nodei];
+    //const auto& contents = tree.contents();
+
+    //Pout<< indent
+    //    << "Entering node:" << nodei
+    //    << " with bb:" << nod.bb_ << endl;
+    //Pout<< incrIndent;
+
+    volumeType myType = volumeType::UNKNOWN;
+
+    for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
+    {
+        const labelBits index = nod.subNodes_[octant];
+        const treeBoundBox subBb(nod.bb_.subBbox(octant));
+
+        if (subBb.intersects(tri))
+        {
+            volumeType subType;
+
+            if (indexedOctree<treeDataTriSurface>::isNode(index))
+            {
+                //Pout<< indent
+                //    << "Testing node octant:" << octant << " with node:"
+                //    << indexedOctree<treeDataTriSurface>::getNode(index)
+                //    << " subBb:" << subBb << endl;
+
+                // tree node. Recurse.
+                subType = markMixed
+                (
+                    tree,
+                    indexedOctree<treeDataTriSurface>::getNode(index),
+                    tri,
+                    nodeTypes
+                );
+            }
+            else if (indexedOctree<treeDataTriSurface>::isContent(index))
+            {
+                // Content inside so we cannot pre-determine whether
+                // whole node box is inside or outside
+
+                //const label contenti =
+                //    indexedOctree<treeDataTriSurface>::getContent(index);
+                //Pout<< indent
+                //    << "Testing content octant:" << octant
+                //    << " with contenti:" << contenti << " subBb:" << subBb
+                //    << " holding shapes:" << flatOutput(contents[contenti])
+                //    << endl;
+
+                // Store octant type
+                subType = volumeType::MIXED;
+            }
+            else
+            {
+                // Nothing in it but intersects triangle. Can't cache
+
+                //Pout<< indent
+                //    << "Testing empty octant:" << octant
+                //    << " subBb:" << subBb << endl;
+
+                // Store octant type
+                subType = volumeType::MIXED;
+            }
+
+            // Store octant type
+            nodeTypes.set(labelBits::pack(nodei, octant), subType);
+
+            // Combine sub node types into type for treeNode. Result is
+            // 'mixed' if types differ among subnodes.
+            if (myType == volumeType::UNKNOWN)
+            {
+                myType = subType;
+            }
+            else if (subType != myType)
+            {
+                //Pout<< indent
+                //    << "for octant:" << octant
+                //    << " siwtched from myType:" << myType
+                //    << " and subType:" << subType
+                //    << " switched to MIXED" << endl;
+                myType = volumeType::MIXED;
+            }
+        }
+    }
+
+    //Pout<< indent
+    //    << "Result for node:" << nodei << " myType:" << myType << endl;
+    //Pout<< decrIndent;
+
+    return myType;
+}
+
+
+void Foam::distributedTriSurfaceMesh::markMixedOverlap
+(
+    const indexedOctree<treeDataTriSurface>& tree,
+    PackedList<2>& nodeTypes
+) const
+{
+    // At this point each triangle is on a unique processor so we can just
+    // go and
+    // check what additional triangles would be sent over based on the
+    // bounding boxes. Note that the surface at this point is already
+    // uniquely distributed (every triangle is on a single processor
+    // only) and the tree is only for those.
+
+    PstreamBuffers pBufs;
+
+    // Test & send
+    forAll(procBb_, proci)
+    {
+        if (proci != UPstream::myProcNo())
+        {
+            labelList pointMap;
+            labelList faceMap;
+            overlappingSurface
+            (
+                *this,
+                procBb_[proci],
+                pointMap,    // not used
+                faceMap
+            );
+            if (!faceMap.empty())
+            {
+                const triSurface subSurface
+                (
+                    subsetMesh
+                    (
+                        *this,
+                        faceMap,
+                        pointMap
+                    )
+                );
+                UOPstream os(proci, pBufs);
+                os << subSurface;
+            }
+        }
+    }
+
+    pBufs.finishedSends();
+
+    // Receive
+    for (const int proci : pBufs.allProcs())
+    {
+        if (proci != Pstream::myProcNo() && pBufs.recvDataCount(proci))
+        {
+            UIPstream is(proci, pBufs);
+
+            // Receive
+            const triSurface subSurface(is);
+
+            //Pout<< "** from proc:" << proci
+            //    << " received:" << subSurface.size()
+            //    << " compared to local:" << this->size() << endl;
+            //const fileName fName
+            //(
+            //    searchableSurface::time().path()
+            //   /searchableSurface::name()
+            // + "_from_proc" + Foam::name(proci) + ".obj"
+            //);
+            //Pout<< "fName:" << fName << endl;
+            //subSurface.write(fName);
+
+
+            // Set all tree nodes overlapping buffer triangles to
+            // mixed. This state will be preserved when caching the
+            // volume type (see below). This uses the actual local
+            // triangles.
+            for (const auto& tri : subSurface)
+            {
+                markMixed
+                (
+                    tree,
+                    0,          //nodei,
+                    tri.tri(subSurface.points()),
+                    nodeTypes
+                );
+            }
+        }
+    }
+}
+
+
 void Foam::distributedTriSurfaceMesh::collectLeafMids
 (
     const label nodeI,
@@ -1623,6 +1975,96 @@ Foam::volumeType Foam::distributedTriSurfaceMesh::calcVolumeType
 }
 
 
+void Foam::distributedTriSurfaceMesh::cacheVolumeType(PackedList<2>& nt) const
+{
+    // Get local tree
+    const indexedOctree<treeDataTriSurface>& t = tree();
+    const auto& nodes = t.nodes();
+
+    // Collect midpoints
+    DynamicField<point> midPoints(label(0.5*nodes.size()));
+    if (nodes.size())
+    {
+        collectLeafMids(0, midPoints);
+    }
+
+    if (debug)
+    {
+        Pout<< "distributedTriSurfaceMesh::cacheVolumeType :"
+            << " surface " << searchableSurface::name()
+            << " triggering orientation caching for "
+            << midPoints.size() << " leaf mids" << endl;
+    }
+
+    // Get volume type of mid points
+    List<volumeType> midVolTypes;
+
+    // Note: could recurse here (since now outsideVolType_ is set)
+    // but this would use the cached mid point data which we're trying
+    // to calculate. Instead bypass caching and do a full search
+    {
+        List<pointIndexHit> nearestInfo;
+        findNearest
+        (
+            midPoints,
+            scalarField(midPoints.size(), Foam::sqr(GREAT)),
+            nearestInfo
+        );
+        surfaceSide(midPoints, nearestInfo, midVolTypes);
+
+        if (debug & 2)
+        {
+            OBJstream insideStr
+            (
+                searchableSurface::time().path()
+               /searchableSurface::name() + "_cacheVolumeType_inside.obj"
+            );
+            OBJstream outsideStr
+            (
+                searchableSurface::time().path()
+               /searchableSurface::name() + "_cacheVolumeType_outside.obj"
+            );
+            forAll(midVolTypes, i)
+            {
+                const linePointRef ln(midPoints[i], nearestInfo[i].hitPoint());
+                if (midVolTypes[i] == volumeType::INSIDE)
+                {
+                    insideStr.write(ln);
+                }
+                else if (midVolTypes[i] == volumeType::OUTSIDE)
+                {
+                    outsideStr.write(ln);
+                }
+            }
+            Pout<< "Whilst caching " << searchableSurface::name()
+                << " have inside:" << insideStr.nVertices()
+                << " have outside:" << outsideStr.nVertices()
+                << endl;
+        }
+    }
+
+    // Cache on local tree
+    if (nodes.size())
+    {
+        label index = 0;
+        calcVolumeType
+        (
+            midVolTypes,
+            index,
+            nt,
+            0               // nodeI
+        );
+    }
+    if (debug)
+    {
+        Pout<< "distributedTriSurfaceMesh::cacheVolumeType :"
+            << " surface " << searchableSurface::name()
+            << " done orientation caching for "
+            << midPoints.size() << " leaf mids" << endl;
+    }
+}
+
+
 Foam::volumeType Foam::distributedTriSurfaceMesh::cachedVolumeType
 (
     const label nodeI,
@@ -1698,17 +2140,95 @@ Foam::volumeType Foam::distributedTriSurfaceMesh::cachedVolumeType
 }
 
 
+const Foam::decompositionMethod&
+Foam::distributedTriSurfaceMesh::decomposer() const
+{
+    if (!decomposer_)
+    {
+        // Check if local dictionary contains decomposition entries
+        if (dict_.found("method"))
+        {
+            if (debug)
+            {
+                Pout<< "distributedTriSurfaceMesh::decomposer() :"
+                    << " constructing decomposer from provided dictionary "
+                    << dict_ << endl;
+            }
+
+            decomposer_ = decompositionMethod::New(dict_);
+        }
+        else
+        {
+            // Use singleton decomposeParDict. Cannot use decompositionModel
+            // here since we've only got Time and not a mesh.
+
+            const auto* dictPtr =
+                searchableSurface::time().findObject<IOdictionary>
+                (
+                    // == decompositionModel::canonicalName
+                    "decomposeParDict"
+                );
+
+            if (dictPtr)
+            {
+                if (debug)
+                {
+                    Pout<< "distributedTriSurfaceMesh::decomposer() :"
+                        << " constructing decomposer from registered"
+                        << " dictionary " << *dictPtr << endl;
+                }
+
+                decomposer_ = decompositionMethod::New(*dictPtr);
+            }
+            else
+            {
+                if (!decomposeParDict_)
+                {
+                    decomposeParDict_.reset
+                    (
+                        new IOdictionary
+                        (
+                            IOobject
+                            (
+                                // == decompositionModel::canonicalName
+                                "decomposeParDict",
+                                searchableSurface::time().system(),
+                                searchableSurface::time(),
+                                IOobject::MUST_READ,
+                                IOobject::NO_WRITE
+                            )
+                        )
+                    );
+
+                    if (debug)
+                    {
+                        Pout<< "distributedTriSurfaceMesh::decomposer() :"
+                            << " reading "
+                            << decomposeParDict_().objectRelPath()
+                            << endl;
+                    }
+                }
+                decomposer_ = decompositionMethod::New(*decomposeParDict_);
+            }
+        }
+    }
+
+    return decomposer_();
+}
+
+
 // Find bounding boxes that guarantee a more or less uniform distribution
 // of triangles. Decomposition in here is only used to get the bounding
 // boxes, actual decomposition is done later on.
 // Returns a per processor a list of bounding boxes that most accurately
 // describe the shape. For now just a single bounding box per processor but
 // optimisation might be to determine a better fitting shape.
-Foam::List<Foam::List<Foam::treeBoundBox>>
-Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
+void Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
 (
-    const triSurface& s
-)
+    const triSurface& s,
+    labelList& distribution,    // per local triangle the destination processor
+    List<List<treeBoundBox>>& bbs
+) const
 {
     addProfiling
     (
@@ -1716,53 +2236,9 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
         "distributedTriSurfaceMesh::independentlyDistributedBbs"
     );
 
-    if (!decomposer_)
-    {
-        // Use singleton decomposeParDict. Cannot use decompositionModel
-        // here since we've only got Time and not a mesh.
-
-        const auto* dictPtr =
-            searchableSurface::time().findObject<IOdictionary>
-            (
-                // == decompositionModel::canonicalName
-                "decomposeParDict"
-            );
-
-        if (dictPtr)
-        {
-            decomposer_ = decompositionMethod::New(*dictPtr);
-        }
-        else
-        {
-            if (!decomposeParDict_)
-            {
-                decomposeParDict_.reset
-                (
-                    new IOdictionary
-                    (
-                        IOobject
-                        (
-                            // == decompositionModel::canonicalName
-                            "decomposeParDict",
-                            searchableSurface::time().system(),
-                            searchableSurface::time(),
-                            IOobject::MUST_READ,
-                            IOobject::NO_WRITE
-                        )
-                    )
-                );
-            }
-            decomposer_ = decompositionMethod::New(*decomposeParDict_);
-        }
-    }
-
-
     // Initialise to inverted box
-    List<List<treeBoundBox>> bbs
-    (
-        Pstream::nProcs(),
-        List<treeBoundBox>(1, treeBoundBox::null())
-    );
+    bbs.resize_nocopy(Pstream::nProcs());
+    bbs = List<treeBoundBox>(1, treeBoundBox::null());
 
     const globalIndex& triIndexer = globalTris();
 
@@ -1795,8 +2271,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
             triCentres[trii] = s[trii].centre(s.points());
         }
 
-        labelList distribution;
-        if (!isA<geomDecomp>(decomposer_()))
+        if (!isA<geomDecomp>(decomposer()))
         {
             // Use connectivity
             labelListList pointFaces;
@@ -1806,19 +2281,19 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
 
             // Do the actual decomposition
             const bool oldParRun = UPstream::parRun(false);
-            distribution = decomposer_().decompose(faceFaces, triCentres);
+            distribution = decomposer().decompose(faceFaces, triCentres);
             UPstream::parRun(oldParRun);  // Restore parallel state
         }
         else
         {
             // Do the actual decomposition
-            distribution = decomposer_().decompose(triCentres);
+            distribution = decomposer().decompose(triCentres);
         }
 
         if (debug)
         {
             Pout<< "distributedTriSurfaceMesh::independentlyDistributedBbs :"
-                << " determining processor bounding boxes for surface"
+                << " determining processor bounding boxes for surface "
                 << searchableSurface::name() << endl;
         }
 
@@ -1845,10 +2320,10 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
     {
         // Fully distributed decomposition. Does not filter duplicate
         // triangles.
-        if (!decomposer_().parallelAware())
+        if (!decomposer().parallelAware())
         {
             FatalErrorInFunction
-                << "The decomposition method " << decomposer_().typeName
+                << "The decomposition method " << decomposer().typeName
                 << " does not decompose in parallel."
                 << " Please choose one that does." << exit(FatalError);
         }
@@ -1868,7 +2343,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
             triCentres[trii] = s[trii].centre(s.points());
         }
 
-        labelList distribution = decomposer_().decompose(triCentres);
+        distribution = decomposer().decompose(triCentres);
 
         if (debug)
         {
@@ -1901,10 +2376,10 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
 //    {
 //        // Fully distributed decomposition. Does not filter duplicate
 //        // triangles.
-//        if (!decomposer_().parallelAware())
+//        if (!decomposer().parallelAware())
 //        {
 //            FatalErrorInFunction
-//                << "The decomposition method " << decomposer_().typeName
+//                << "The decomposition method " << decomposer().typeName
 //                << " does not decompose in parallel."
 //                << " Please choose one that does." << exit(FatalError);
 //        }
@@ -2060,7 +2535,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
 //        pointField masterCentres(allCentres, triMap);
 //
 //        // Do the actual decomposition
-//        labelList masterDistribution(decomposer_().decompose(masterCentres));
+//        labelList masterDistribution(decomposer().decompose(masterCentres));
 //
 //        // Make map to get the decomposition from the master of each border
 //        labelList borderGlobalMaster(borderTris.size());
@@ -2120,7 +2595,6 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
         }
 
         // Determine local decomposition
-        labelList distribution(s.size());
         {
             labelList allDistribution;
             if (Pstream::master())
@@ -2149,7 +2623,7 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
 
                 // Decompose merged centres
                 const bool oldParRun = UPstream::parRun(false);
-                labelList mergedDist(decomposer_().decompose(mergedPoints));
+                labelList mergedDist(decomposer().decompose(mergedPoints));
                 UPstream::parRun(oldParRun);  // Restore parallel state
 
                 // Convert to all
@@ -2202,7 +2676,6 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
         }
         Pstream::broadcast(bbs);
     }
-    return bbs;
 }
 
 
@@ -2576,8 +3049,10 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
     ),
     currentDistType_(FROZEN)    // only used to trigger re-distribution
 {
-    // Read from the provided dictionary
-    read();
+    // Construct from components. If empty dictionary get all the info from the
+    // supplied (undecomposed) triSurface.
+
+    readSettings(dict_.empty());
 
     bounds().reduce();
 
@@ -2639,21 +3114,34 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
     ),
     currentDistType_(FROZEN)    // only used to trigger re-distribution
 {
-    // Read from the local, decomposed dictionary
-    read();
-
-    bounds().reduce();
-
     // Try and find out where the triSurface was loaded from. On slave this
     // might give empty filename.
     const fileName actualFile(triSurfaceMesh::findFile(io, true));
 
+    // Was it read undecomposed?
+    const bool readFromMaster =
+    (
+        actualFile.empty()
+     || actualFile != io.localFilePath(triSurfaceMesh::typeName)
+    );
+
+    readSettings(readFromMaster);
+
+    bounds().reduce();
+
+    if (readFromMaster && !decomposeUsingBbs_)
+    {
+        // No fill-in so store normals
+        DebugInFunction
+            << "Determining vertex based normals since undecomposed" << endl;
+        const triSurface& surf = *this;
+        vertexNormals_ = autoPtr<List<FixedList<vector, 3>>>::New();
+        calcVertexNormals(surf, vertexNormals_());
+    }
+
     if
     (
-        (
-            actualFile.empty()
-         || actualFile != io.localFilePath(triSurfaceMesh::typeName)
-        )
+        readFromMaster
      && (distType_ == INDEPENDENT || distType_ == DISTRIBUTED)
     )
     {
@@ -2748,50 +3236,37 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh
     ),
     currentDistType_(FROZEN)    // only used to trigger re-distribution
 {
-    // Read from the local, decomposed dictionary
-    read();
-
-    // Optionally override settings from provided dictionary
-    {
-        // Wanted distribution type
-        distributionTypeNames_.readIfPresent
-        (
-            "distributionType",
-            dict,
-            distType_
-        );
-
-        // Merge distance
-        dict.readIfPresent("mergeDistance", mergeDist_);
-
-        // Distribution type
-        bool closed;
-        if (dict.readIfPresent<bool>("closed", closed))
-        {
-            surfaceClosed_ = closed;
-        }
-
-        outsideVolType_ = volumeType::names.getOrDefault
-        (
-            "outsideVolumeType",
-            dict,
-            outsideVolType_
-        );
-    }
-
-
-    bounds().reduce();
-
     // Try and find out where the triSurface was loaded from. On slave this
     // might give empty filename.
     const fileName actualFile(triSurfaceMesh::findFile(io, dict, true));
 
+    // Was it read undecomposed?
+    const bool readFromMaster =
+    (
+        actualFile.empty()
+     || actualFile != io.localFilePath(triSurfaceMesh::typeName)
+    );
+
+    // Add supplied dictionary entries (override existing)
+    dict_ <<= dict;
+
+    readSettings(readFromMaster);
+
+    bounds().reduce();
+
+    if (readFromMaster && !decomposeUsingBbs_)
+    {
+        // No fill-in so store enough normals to calculate
+        DebugInFunction
+            << "Determining vertex based normals since undecomposed" << endl;
+        const triSurface& surf = *this;
+        vertexNormals_ = autoPtr<List<FixedList<vector, 3>>>::New();
+        calcVertexNormals(surf, vertexNormals_());
+    }
+
     if
     (
-        (
-            actualFile.empty()
-         || actualFile != io.localFilePath(triSurfaceMesh::typeName)
-        )
+        readFromMaster
      && (distType_ == INDEPENDENT || distType_ == DISTRIBUTED)
     )
     {
@@ -2878,6 +3353,35 @@ const Foam::globalIndex& Foam::distributedTriSurfaceMesh::globalTris() const
 }
 
 
+void Foam::distributedTriSurfaceMesh::flip()
+{
+    if (debug)
+    {
+        Pout<< "distributedTriSurfaceMesh::flip :"
+            << " surface " << searchableSurface::name()
+            << " with current orientation "
+            << volumeType::names[outsideVolType_]
+            << " vertexNormals_ " << vertexNormals_.good()
+            << endl;
+    }
+
+    triSurfaceMesh::flip();
+
+    if (vertexNormals_)
+    {
+        List<FixedList<vector, 3>>& vn = vertexNormals_();
+
+        for (auto& normals : vn)
+        {
+            for (auto& n : normals)
+            {
+                n = -n;
+            }
+        }
+    }
+}
+
+
 //void Foam::distributedTriSurfaceMesh::findNearest
 //(
 //    const pointField& samples,
@@ -3902,7 +4406,7 @@ void Foam::distributedTriSurfaceMesh::getRegion
     // Do my tests
     // ~~~~~~~~~~~
 
-    const triSurface& s = static_cast<const triSurface&>(*this);
+    const triSurface& s = *this;
 
     region.setSize(triangleIndex.size());
 
@@ -3969,14 +4473,36 @@ void Foam::distributedTriSurfaceMesh::getNormal
     // Do my tests
     // ~~~~~~~~~~~
 
-    const triSurface& s = static_cast<const triSurface&>(*this);
+    const triSurface& s = *this;
 
     normal.setSize(triangleIndex.size());
 
-    forAll(triangleIndex, i)
+    if (vertexNormals_)
     {
-        label trii = triangleIndex[i];
-        normal[i] = s[trii].unitNormal(s.points());
+        // Use smooth interpolation for the normal
+
+        // Send over nearest point since we have it anyway
+        List<pointIndexHit> localNearestInfo(info);
+        map.distribute(localNearestInfo);
+
+        forAll(triangleIndex, i)
+        {
+            const point& nearest = localNearestInfo[i].point();
+
+            const label trii = triangleIndex[i];
+            const triPointRef tri(s[trii].tri(s.points()));
+            const auto& vn = vertexNormals_()[trii];
+            const barycentric2D w(tri.pointToBarycentric(nearest));
+            normal[i] = w[0]*vn[0]+w[1]*vn[1]+w[2]*vn[2];
+        }
+    }
+    else
+    {
+        forAll(triangleIndex, i)
+        {
+            const label trii = triangleIndex[i];
+            normal[i] = s[trii].unitNormal(s.points());
+        }
     }
 
 
@@ -4138,64 +4664,18 @@ void Foam::distributedTriSurfaceMesh::getVolumeType
          || outsideVolType_ == volumeType::OUTSIDE
         )
         {
-            // Get local tree
             const indexedOctree<treeDataTriSurface>& t = tree();
             const auto& nodes = t.nodes();
             PackedList<2>& nt = t.nodeTypes();
             nt.resize(nodes.size());
             nt = volumeType::UNKNOWN;
 
-            // Collect midpoints
-            DynamicField<point> midPoints(label(0.5*nodes.size()));
-            if (nodes.size())
+            if (!decomposeUsingBbs_)
             {
-                collectLeafMids(0, midPoints);
+                markMixedOverlap(t, nt);
             }
 
-            if (debug)
-            {
-                Pout<< "distributedTriSurfaceMesh::getVolumeType :"
-                    << " surface " << searchableSurface::name()
-                    << " triggering orientation caching for "
-                    << midPoints.size() << " leaf mids" << endl;
-            }
-
-            // Get volume type of mid points
-            List<volumeType> midVolTypes;
-
-            // Note: could recurse here (since now outsideVolType_ is set)
-            // but this would use the cached mid point data which we're trying
-            // to calculate. Instead bypass caching and do a full search
-            {
-                List<pointIndexHit> nearestInfo;
-                findNearest
-                (
-                    midPoints,
-                    scalarField(midPoints.size(), Foam::sqr(GREAT)),
-                    nearestInfo
-                );
-                surfaceSide(midPoints, nearestInfo, midVolTypes);
-            }
-
-            // Cache on local tree
-            if (nodes.size())
-            {
-                label index = 0;
-                calcVolumeType
-                (
-                    midVolTypes,
-                    index,
-                    nt,
-                    0               // nodeI
-                );
-            }
-            if (debug)
-            {
-                Pout<< "distributedTriSurfaceMesh::getVolumeType :"
-                    << " surface " << searchableSurface::name()
-                    << " done orientation caching for "
-                    << midPoints.size() << " leaf mids" << endl;
-            }
+            cacheVolumeType(nt);
         }
     }
 
@@ -4282,9 +4762,18 @@ void Foam::distributedTriSurfaceMesh::getVolumeType
     // those I need to search the nearest for
     DynamicField<point> fullSearchPoints(localPoints.size());
     DynamicList<label> fullSearchMap(localPoints.size());
+
     forAll(localPoints, i)
     {
-        if (tree().nodes().size())
+        // If no-fill in:
+        // - tree nodes on the outside will not include fill-in triangles
+        // - so might decide the wrong state (might be 'mixed' with fill-in
+        //   triangles)
+        // - for now do not cache
+        // - TBD. Note that the cached volume types already include remote
+        //   fill-in triangles. There is still a problem there though. TBD.
+
+        if (decomposeUsingBbs_ && tree().nodes().size())
         {
             volType[i] = cachedVolumeType(0, localPoints[i]);
         }
@@ -4322,9 +4811,68 @@ void Foam::distributedTriSurfaceMesh::getVolumeType
     // Combine
     forAll(fullSearchMap, i)
     {
-        volType[fullSearchMap[i]] = fullSearchType[i];
+        volumeType& vt = volType[fullSearchMap[i]];
+        if (vt == volumeType::UNKNOWN)
+        {
+            vt = fullSearchType[i];
+        }
+        else if (vt != fullSearchType[i])
+        {
+            Pout<< "At point:" << fullSearchPoints[i]
+                << " or point:" << localPoints[fullSearchMap[i]]
+                << " have cached status:" << vt
+                << " have full-search status:" << fullSearchType[i]
+                << endl;
+        }
+        //volType[fullSearchMap[i]] = fullSearchType[i];
     }
 
+    //{
+    //    //- Variant of below that produces better error messages
+    //    typedef Tuple2<Pair<point>, volumeType> NearType;
+    //    const NearType zero(Pair<point>(Zero, Zero), volumeType::UNKNOWN);
+    //
+    //    // Combine nearest and volType so we can give nicer error messages
+    //    List<NearType> nearTypes(volType.size(), zero);
+    //
+    //    NearTypeCombineOp cop;
+    //
+    //    // For all volType we do have the originating point but not the
+    //    // nearest
+    //    forAll(localPoints, i)
+    //    {
+    //        nearTypes[i] =
+    //            NearType(Pair<point>(localPoints[i], Zero), volType[i]);
+    //    }
+    //    // But for all full-search we also have the nearest
+    //    forAll(fullSearchMap, i)
+    //    {
+    //        const NearType nt
+    //        (
+    //            Pair<point>(fullSearchPoints[i], nearestInfo[i].hitPoint()),
+    //            fullSearchType[i]
+    //        );
+    //
+    //        cop(nearTypes[fullSearchMap[i]], nt);
+    //    }
+    //    mapDistributeBase::distribute
+    //    (
+    //        Pstream::commsTypes::nonBlocking,
+    //        List<labelPair>::null(),
+    //        samples.size(),
+    //        map.constructMap(),
+    //        map.constructHasFlip(),
+    //        map.subMap(),
+    //        map.subHasFlip(),
+    //        nearTypes,
+    //        zero,
+    //        cop,                //NearTypeCombineOp(),
+    //        noOp(),             // no flipping
+    //        UPstream::msgType(),
+    //        map.comm()
+    //    );
+    //}
+
     // Send back to originator. In case of multiple answers choose inside or
     // outside
     const volumeType zero(volumeType::UNKNOWN);
@@ -4345,47 +4893,6 @@ void Foam::distributedTriSurfaceMesh::getVolumeType
         map.comm()
     );
 
-
-    //// Combine nearest and volType so we can give nicer error messages
-    //List<NearType> nearTypes(volType.size());
-    //// For all volType we do have the originating point but not the nearest
-    //forAll(localPoints, i)
-    //{
-    //    nearTypes[i] =
-    //        NearType(Pair<point>(localPoints[i], Zero), volType[i]);
-    //}
-    //// But for all full-search we also have the nearest
-    //forAll(fullSearchMap, i)
-    //{
-    //    nearTypes[fullSearchMap[i]] = NearType
-    //    (
-    //        Pair<point>(fullSearchPoints[i], nearestInfo[i].hitPoint()),
-    //        fullSearchType[i]
-    //    );
-    //}
-    //const NearType zero(Pair<point>(Zero, Zero), volumeType::UNKNOWN);
-    //mapDistributeBase::distribute
-    //(
-    //    Pstream::commsTypes::nonBlocking,
-    //    List<labelPair>::null(),
-    //    samples.size(),
-    //    map.constructMap(),
-    //    map.constructHasFlip(),
-    //    map.subMap(),
-    //    map.subHasFlip(),
-    //    nearTypes,
-    //    zero,
-    //    NearTypeCombineOp(),
-    //    noOp(),           // no flipping
-    //    UPstream::msgType(),
-    //    map.comm()
-    //);
-    //volType.setSize(nearTypes.size());
-    //forAll(nearTypes, i)
-    //{
-    //    volType[i] = nearTypes[i].second();
-    //}
-
     // Add the points outside the bounding box
     for (label samplei : outsideSamples)
     {
@@ -4613,6 +5120,7 @@ void Foam::distributedTriSurfaceMesh::distribute
             << " surface " << searchableSurface::name()
             << " distributing surface according to method:"
             << distributionTypeNames_[distType_]
+            << " fill-in:" << decomposeUsingBbs_
             << " follow bbs:" << flatOutput(bbs) << endl;
     }
 
@@ -4622,7 +5130,10 @@ void Foam::distributedTriSurfaceMesh::distribute
     // Get bbs of all domains
     // ~~~~~~~~~~~~~~~~~~~~~~
 
+    // Per triangle the destination processor
+    labelList distribution;
     {
+        // Per processor the bounding box of all (destination) triangles
         List<List<treeBoundBox>> newProcBb(Pstream::nProcs());
 
         switch (distType_)
@@ -4638,7 +5149,7 @@ void Foam::distributedTriSurfaceMesh::distribute
                 {
                     return;
                 }
-                newProcBb = independentlyDistributedBbs(*this);
+                independentlyDistributedBbs(*this, distribution, newProcBb);
             break;
 
             case FROZEN:
@@ -4691,15 +5202,37 @@ void Foam::distributedTriSurfaceMesh::distribute
     labelListList faceSendMap(Pstream::nProcs());
     labelListList pointSendMap(Pstream::nProcs());
 
-    forAll(procBb_, proci)
+    if (decomposeUsingBbs_)
     {
-        overlappingSurface
-        (
-            *this,
-            procBb_[proci],
-            pointSendMap[proci],
-            faceSendMap[proci]
-        );
+        forAll(procBb_, proci)
+        {
+            overlappingSurface
+            (
+                *this,
+                procBb_[proci],
+                pointSendMap[proci],
+                faceSendMap[proci]
+            );
+        }
+    }
+    else
+    {
+        const triSurface& s = *this;
+        forAll(procBb_, proci)
+        {
+            boolList includedFace(s.size(), false);
+            forAll(s, trii)
+            {
+                includedFace[trii] = (distribution[trii] == proci);
+            }
+            subsetMesh
+            (
+                s,
+                includedFace,
+                pointSendMap[proci],
+                faceSendMap[proci]
+            );
+        }
     }
 
     if (keepNonLocal)
@@ -4707,7 +5240,7 @@ void Foam::distributedTriSurfaceMesh::distribute
         // Include in faceSendMap/pointSendMap the triangles that are
         // not mapped to any processor so they stay local.
 
-        const triSurface& s = static_cast<const triSurface&>(*this);
+        const triSurface& s = *this;
 
         boolList includedFace(s.size(), true);
 
@@ -4819,7 +5352,7 @@ void Foam::distributedTriSurfaceMesh::distribute
             UIPstream is(proci, pBufs);
 
             // Receive
-            triSurface subSurface(is);
+            const triSurface subSurface(is);
 
             // Merge into allSurf
             merge
@@ -4876,6 +5409,63 @@ void Foam::distributedTriSurfaceMesh::distribute
     distributeFields<sphericalTensor>(faceMap());
     distributeFields<symmTensor>(faceMap());
     distributeFields<tensor>(faceMap());
+    if (vertexNormals_)
+    {
+        List<FixedList<vector, 3>>& vn = vertexNormals_();
+        faceMap().distribute(vn);
+
+        if (debug & 2)
+        {
+            OBJstream str
+            (
+                searchableSurface::time().path()
+               /searchableSurface::name() + "_vertex_normals.obj"
+            );
+            const triSurface& s = *this;
+            forAll(s, trii)
+            {
+                const auto& f = s[trii];
+                const scalar l = Foam::sqrt(f.mag(s.points()));
+                forAll(f, fp)
+                {
+                    const point& pt = s.points()[f[fp]];
+                    const vector& n = vn[trii][fp];
+                    str.write(linePointRef(pt, pt+l*n));
+                }
+            }
+            Pout<< "Dumped local vertex normals to " << str.name() << endl;
+        }
+    }
+
+
+    if
+    (
+        outsideVolType_ == volumeType::INSIDE
+     || outsideVolType_ == volumeType::OUTSIDE
+    )
+    {
+        // Re-build tree & inside/outside for closed surfaces
+
+        const indexedOctree<treeDataTriSurface>& t = tree();
+        const auto& nodes = t.nodes();
+        PackedList<2>& nt = t.nodeTypes();
+        nt.resize(nodes.size());
+        nt = volumeType::UNKNOWN;
+
+        // Send over any overlapping but not included triangles
+        if (!decomposeUsingBbs_)
+        {
+            // We do not backfill the bounding box so the local tree might
+            // not have all the triangles to cache inside/outside correctly.
+            // Do a pre-pass to mark 'mixed' the tree leaves that intersect
+            // the triangles of the (would be) fill-in triangles.
+            markMixedOverlap(t, nt);
+        }
+
+        cacheVolumeType(nt);
+    }
+
+
 
     if (debug)
     {
@@ -4995,16 +5585,17 @@ void Foam::distributedTriSurfaceMesh::writeStats(Ostream& os) const
 {
     boundBox bb;
     label nPoints;
-    PatchTools::calcBounds(static_cast<const triSurface&>(*this), bb, nPoints);
+    PatchTools::calcBounds(*this, bb, nPoints);
     bb.reduce();
 
-    os  << "Triangles    : " << returnReduce(triSurface::size(), sumOp<label>())
-        << endl
-        << "Vertices     : " << returnReduce(nPoints, sumOp<label>()) << endl
-        << "Bounding Box : " << bb << endl
-        << "Closed       : " << surfaceClosed_ << endl
-        << "Outside point: " << volumeType::names[outsideVolType_] << endl
-        << "Distribution : " << distributionTypeNames_[distType_] << endl;
+    os  << "Triangles      : "
+        << returnReduce(triSurface::size(), sumOp<label>()) << endl
+        << "Vertices       : " << returnReduce(nPoints, sumOp<label>()) << endl
+        << "Vertex normals : " << vertexNormals_.valid() << endl
+        << "Bounding Box   : " << bb << endl
+        << "Closed         : " << surfaceClosed_ << endl
+        << "Outside type   : " << volumeType::names[outsideVolType_] << endl
+        << "Distribution   : " << distributionTypeNames_[distType_] << endl;
 }
 
 
diff --git a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.H b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.H
index 05adb82d5c..5c4498994b 100644
--- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.H
+++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.H
@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2015-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -109,20 +109,28 @@ private:
         //- Merging distance
         scalar mergeDist_;
 
-        autoPtr<IOdictionary> decomposeParDict_;
+        mutable autoPtr<IOdictionary> decomposeParDict_;
 
         //- Decomposition used when independently decomposing surface.
-        autoPtr<decompositionMethod> decomposer_;
+        mutable autoPtr<decompositionMethod> decomposer_;
 
         //- Bounding box settings
         localIOdictionary dict_;
 
+        //- Use bounding boxes (default) or unique decomposition of triangles
+        //- (i.e. do not duplicate triangles)
+        bool decomposeUsingBbs_;
+
         //- Bounding boxes of all processors
         List<List<treeBoundBox>> procBb_;
 
         //- Global triangle numbering
         mutable autoPtr<globalIndex> globalTris_;
 
+        //- Optional per-vertex normals. TBD: move to triSurface? or have
+        //- per-triangle 3 normals so we can interpolate and have features
+        mutable autoPtr<List<FixedList<vector, 3>>> vertexNormals_;
+
         //- The (wanted) distribution type.
         distributionType distType_;
 
@@ -139,8 +147,17 @@ private:
             //  directory or in parent directory
             static word findLocalInstance(const IOobject& io);
 
-            //- Read my additional data
-            bool read();
+            //- Read my additional data from dictionary. Additional flag to
+            //- say whether we can use master-only geometric tests.
+            bool readSettings(const bool isUndecomposed);
+
+            //- Construction helper: generate vertex normals upon reading
+            //- undecomposed surface
+            void calcVertexNormals
+            (
+                const triSurface& surf,
+                List<FixedList<vector, 3>>& vn
+            ) const;
 
 
         // Line intersection
@@ -287,6 +304,23 @@ private:
 
             // Caching of volume type (based on indexedOctree)
 
+            //- Set node type on any node containing the triangle
+            volumeType markMixed
+            (
+                const indexedOctree<treeDataTriSurface>& tree,
+                const label nodei,
+                const triPointRef& tri,
+                PackedList<2>& nodeTypes
+            ) const;
+
+            //- Set node type on any node overlapping any remote triangles.
+            //- Only valid if using unique decomposition.
+            void markMixedOverlap
+            (
+                const indexedOctree<treeDataTriSurface>& tree,
+                PackedList<2>& nodeTypes
+            ) const;
+
             //- Collect mid points of tree boxes
             void collectLeafMids
             (
@@ -303,6 +337,9 @@ private:
                 const label nodeI
             ) const;
 
+            //- Calculate inside/outside of midpoint of tree nodes
+            void cacheVolumeType(PackedList<2>& nt) const;
+
             //- Look up any cached data. Return unknown if cannot be determined.
             volumeType cachedVolumeType
             (
@@ -321,11 +358,16 @@ private:
                 labelListList& faceFaces
             );
 
+            //- Helper: get decompositionMethod
+            const decompositionMethod& decomposer() const;
+
             //- Finds new bounds based on an independent decomposition.
-            List<List<treeBoundBox>> independentlyDistributedBbs
+            void independentlyDistributedBbs
             (
-                const triSurface&
-            );
+                const triSurface& s,
+                labelList& distribution,
+                List<List<treeBoundBox>>& bbs
+            ) const;
 
             //- Does any part of triangle overlap bb.
             static bool overlaps
@@ -413,7 +455,8 @@ public:
 
     // Constructors
 
-        //- Construct from triSurface
+        //- Construct from components. Assumes triSurface is already decomposed
+        //- and dictionary contains corresponding information
         distributedTriSurfaceMesh
         (
             const IOobject&,
@@ -421,11 +464,17 @@ public:
             const dictionary& dict
         );
 
-        //- Construct read. Does findInstance to find io.local().
+        //- Construct read. Does findInstance to find io.local()
+        //-   - if found local : assume distributed
+        //-   - if found in parent : assume undistributed. Can e.g. check for
+        //-     closedness.
         distributedTriSurfaceMesh(const IOobject& io);
 
         //- Construct from dictionary (used by searchableSurface).
         //  Does read. Does findInstance to find io.local().
+        //-   - if found local : assume distributed
+        //-   - if found in parent : assume undistributed. Can e.g. check for
+        //-     closedness.
         distributedTriSurfaceMesh
         (
             const IOobject& io,
@@ -454,6 +503,9 @@ public:
                 return globalTris().totalSize();
             }
 
+            //- Flip triangles, outsideVolumeType and all cached inside/outside.
+            virtual void flip();
+
             virtual void findNearest
             (
                 const pointField& sample,
diff --git a/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.C b/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.C
index 73fe2e15fe..2fce36be06 100644
--- a/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.C
+++ b/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2022 OpenCFD Ltd.
+    Copyright (C) 2018-2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -61,10 +61,9 @@ Foam::patchDistMethods::exact::patchSurface() const
             localBb.extend(rndGen, 1E-3)
         );
 
-        // Dummy bounds dictionary
-        dictionary dict;
-        dict.add("bounds", meshBb);
-        dict.add
+        // Add any missing properties (but not override existing ones)
+        dict_.add("bounds", meshBb);
+        dict_.add
         (
             "distributionType",
             distributedTriSurfaceMesh::distributionTypeNames_
@@ -74,8 +73,7 @@ Foam::patchDistMethods::exact::patchSurface() const
                 distributedTriSurfaceMesh::DISTRIBUTED      // parallel decomp
             ]
         );
-        dict.add("mergeDistance", 1e-6*localBb.mag());
-
+        dict_.add("mergeDistance", 1e-6*localBb.mag());
 
         Info<< "Triangulating local patch faces" << nl << endl;
 
@@ -100,7 +98,7 @@ Foam::patchDistMethods::exact::patchSurface() const
                     patchIDs_,
                     mapTriToGlobal
                 ),
-                dict
+                dict_
             )
         );
 
@@ -131,7 +129,8 @@ Foam::patchDistMethods::exact::exact
     const labelHashSet& patchIDs
 )
 :
-    patchDistMethod(mesh, patchIDs)
+    patchDistMethod(mesh, patchIDs),
+    dict_(dict.subOrEmptyDict(typeName + "Coeffs"))
 {}
 
 
diff --git a/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.H b/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.H
index 4fd1f82fb1..3ac84561ce 100644
--- a/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.H
+++ b/src/parallel/distributed/patchDistMethods/exact/exactPatchDistMethod.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2018-2019 OpenCFD Ltd.
+    Copyright (C) 2018-2019,2024 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,9 +30,32 @@ Description
     Calculation of exact distance to nearest patch for all cells and
     boundary by constructing a search tree for all patch faces.
 
+Usage
+
+    \begin{verbatim}
+    wallDist
+    {
+        method          exactDistance;
+
+        // Optional entries (currently for distributedTriSurfaceMesh)
+        exactDistanceCoeffs
+        {
+            // Optional decomposition method. If not supplied will use the
+            // system/decomposeParDict.
+            method              hierarchical;
+            numberOfSubdomains  8;
+            n                   (2 2 2);
+
+            // (not)add fill-in triangles
+            decomposeUsingBbs   false;
+        }
+    }
+    \end{verbatim}
+
 See also
     Foam::patchDistMethod::meshWave
     Foam::wallDist
+    Foam::distributedTriSurfaceMesh
 
 SourceFiles
     exactPatchDistMethod.C
@@ -55,7 +78,7 @@ namespace patchDistMethods
 {
 
 /*---------------------------------------------------------------------------*\
-                          Class exact Declaration
+                            Class exact Declaration
 \*---------------------------------------------------------------------------*/
 
 class exact
@@ -64,6 +87,9 @@ class exact
 {
     // Private Member Data
 
+        //- Dictionary contents from the dictionary constructor
+        mutable dictionary dict_;
+
         //- Cache surface+searching of patch
         mutable autoPtr<distributedTriSurfaceMesh> patchSurfPtr_;
 
diff --git a/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/snappyHexMeshDict b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/snappyHexMeshDict
index 9fb6072860..bb4f915e9e 100644
--- a/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/snappyHexMeshDict
+++ b/tutorials/mesh/snappyHexMesh/distributedTriSurfaceMesh/system/snappyHexMeshDict
@@ -20,6 +20,18 @@ snap            true;
 addLayers       false;
 
 
+// Define common distributedTriSurface properties
+_decomp
+{
+    numberOfSubdomains  8;
+    method              hierarchical;
+    n                   (2 2 2);
+
+    // unique triangle decomposition; no fill-in
+    decomposeUsingBbs   false;
+}
+
+
 // Geometry. Definition of all surfaces. All surfaces are of class
 // searchableSurface.
 // Surfaces are used
@@ -32,16 +44,19 @@ geometry
     {
         type distributedTriSurfaceMesh;
         file "box.obj"; //"box_12_2.obj";
+        ${_decomp}
     }
     box_trans
     {
-        file "box_trans.obj";
         type distributedTriSurfaceMesh;
+        file "box_trans.obj";
+        ${_decomp}
     }
     shell
     {
-        file "box_scaled.obj";
         type distributedTriSurfaceMesh;
+        file "box_scaled.obj";
+        ${_decomp}
     }
 };