ENH: distributedTriSurfaceMesh: improvements

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
+        if (modes_[shellI] != DISTANCE && isA<triSurfaceMesh>(s))
-        (
-            modes_[shellI] != DISTANCE
-        &&  isA<triSurfaceMesh>(s)
-        && !isA<distributedTriSurfaceMesh>(s)
-        )
         {
             hasSurface = true;
 
@@ -196,35 +191,39 @@ void Foam::shellSurfaces::orient()
         {
             const searchableSurface& s = allGeometry_[shells_[shellI]];
 
-            if
+            if (modes_[shellI] != DISTANCE && isA<triSurfaceMesh>(s))
-            (
-                modes_[shellI] != DISTANCE
-            &&  isA<triSurfaceMesh>(s)
-            && !isA<distributedTriSurfaceMesh>(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.
+                //Pout<< "outsidePt:" << outsidePt << endl;
-                bool anyFlipped = orientedSurface::orient
+                //Pout<< "info     :" << info[0] << endl;
-                (
+                //Pout<< "normal   :" << normal[0] << endl;
-                    shell,
+                //Pout<< "region   :" << region[0] << endl;
-                    outsidePt,
+
-                    true
+                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;

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;
     }

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;

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;
+            }
+        }
+    }
+}
+
+
 // ************************************************************************* //

src/meshTools/triSurface/triSurfaceSearch/triSurfaceRegionSearch.H

@@ -128,6 +128,11 @@ public:
                 const labelList& regionIndices,
                 List<pointIndexHit>& info
             ) const;
+
+        // Edit
+
+            //- Flip orientation
+            void flip();
 };
 
 

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
 (

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
         {

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
+            //- Read my additional data from dictionary. Additional flag to
-            bool read();
+            //- 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,

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
+        // Add any missing properties (but not override existing ones)
-        dictionary dict;
+        dict_.add("bounds", meshBb);
-        dict.add("bounds", meshBb);
+        dict_.add
-        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"))
 {}
 
 

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_;
 

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}
     }
 };