diff --git a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
index a6547fa1142cc9ce4ce14ffe8a302042d99090c7..221c5c6581ad0e75806ac2b1027ffaebce962beb 100644
--- a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
+++ b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
@@ -193,6 +193,15 @@ void Foam::shellSurfaces::orient()
if (modes_[shellI] != DISTANCE && isA<triSurfaceMesh>(s))
{
+ triSurfaceMesh& shell = const_cast<triSurfaceMesh&>
+ (
+ refCast<const triSurfaceMesh>(s)
+ );
+
+ // Make sure that normals are consistent. Does not change
+ // anything if surface is consistent.
+ orientedSurface::orientConsistent(shell);
+
List<pointIndexHit> info;
vectorField normal;
labelList region;
@@ -213,15 +222,10 @@ void Foam::shellSurfaces::orient()
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_)
{
Info<< "shellSurfaces : Flipped orientation of surface "
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index f207200fabca45b1836e841a6a99cf70c771e626..eb64fa51a0d738ee53a4d409cfd56cd6161dee3d 100644
--- a/src/meshTools/Make/files
+++ b/src/meshTools/Make/files
@@ -112,6 +112,7 @@ searchableSurfaces/searchableSurfacesQueries/searchableSurfacesQueries.C
searchableSurfaces/searchableSurfaceWithGaps/searchableSurfaceWithGaps.C
searchableSurfaces/subTriSurfaceMesh/subTriSurfaceMesh.C
searchableSurfaces/triSurfaceMesh/triSurfaceMesh.C
+searchableSurfaces/triSurfaceMesh/smoothTriSurfaceMesh.C
coordSet/coordSet.C
diff --git a/src/meshTools/searchableSurfaces/triSurfaceMesh/smoothTriSurfaceMesh.C b/src/meshTools/searchableSurfaces/triSurfaceMesh/smoothTriSurfaceMesh.C
new file mode 100644
index 0000000000000000000000000000000000000000..0c6dbec6522a2a4734298224815009974e64930c
--- /dev/null
+++ b/src/meshTools/searchableSurfaces/triSurfaceMesh/smoothTriSurfaceMesh.C
@@ -0,0 +1,318 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2016-2018 OpenFOAM Foundation
+ Copyright (C) 2025 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "smoothTriSurfaceMesh.H"
+#include "addToRunTimeSelectionTable.H"
+#include "unitConversion.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(smoothTriSurfaceMesh, 0);
+addToRunTimeSelectionTable(searchableSurface, smoothTriSurfaceMesh, dict);
+
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+void Foam::smoothTriSurfaceMesh::calcFeatureEdges(const scalar featureAngle)
+{
+ if (featureAngle <= -180)
+ {
+ return;
+ }
+
+ const triSurface& s = *this;
+
+ scalar cosAngle = Foam::cos(degToRad(featureAngle));
+
+ const labelListList& edgeFaces = s.edgeFaces();
+ const edgeList& edges = s.edges();
+
+ forAll(edgeFaces, edgei)
+ {
+ const labelList& eFaces = edgeFaces[edgei];
+
+ if (eFaces.size() > 2)
+ {
+ isBorderEdge_[edgei] = true;
+ }
+ else if (eFaces.size() == 2)
+ {
+ const vector n0(s[eFaces[0]].unitNormal(points()));
+ const vector n1(s[eFaces[1]].unitNormal(points()));
+ if ((n0&n1) < cosAngle)
+ {
+ isBorderEdge_[edgei] = true;
+ }
+ }
+ }
+
+ forAll(isBorderEdge_, edgei)
+ {
+ if (isBorderEdge_[edgei])
+ {
+ const edge& e = edges[edgei];
+ isPointOnBorderEdge_[e[0]] = true;
+ isPointOnBorderEdge_[e[1]] = true;
+ }
+ }
+}
+
+
+Foam::vector Foam::smoothTriSurfaceMesh::pointNormal
+(
+ const label startFacei,
+ const label localPointi
+) const
+{
+ const triSurface& s = *this;
+
+ if (!isPointOnBorderEdge_[localPointi])
+ {
+ const labelList& pFaces = pointFaces()[localPointi];
+ vector pn(Zero);
+ forAll(pFaces, i)
+ {
+ label facei = pFaces[i];
+ pn += s[facei].unitNormal(points());
+ }
+ return pn.normalise();
+ }
+
+
+ // Calculate the local point normal on the face. This routine
+ // only gets called if the point is on a border edge so we can
+ // walk and always hit a border edge.
+
+ const edgeList& edges = triSurface::edges();
+ const labelList& fEdges = faceEdges()[startFacei];
+
+ // Get the two edges on the point
+ label e0 = -1;
+ label e1 = -1;
+ forAll(fEdges, i)
+ {
+ const edge& e = edges[fEdges[i]];
+ if (e.otherVertex(localPointi) == -1)
+ {
+ e0 = fEdges[fEdges.rcIndex(i)];
+ e1 = fEdges[fEdges.fcIndex(i)];
+ break;
+ }
+ }
+
+ label facei = startFacei;
+ label edgei = e0;
+
+ // Initialise normal
+ vector n(s[facei].unitNormal(points()));
+
+ while (!isBorderEdge_[edgei])
+ {
+ // Cross edge to next face
+ const labelList& eFaces = edgeFaces()[edgei];
+ if (eFaces.size() != 2)
+ {
+ break;
+ }
+
+ forAll(eFaces, i)
+ {
+ if (eFaces[i] != facei)
+ {
+ facei = eFaces[i];
+ break;
+ }
+ }
+
+ if (facei == startFacei)
+ {
+ break;
+ }
+
+ n += s[facei].unitNormal(points());
+
+ // Cross face to next edge
+ const labelList& fEdges = faceEdges()[facei];
+
+ forAll(fEdges, i)
+ {
+ label ei = fEdges[i];
+ if (ei != edgei && edges[ei].otherVertex(localPointi) != -1)
+ {
+ edgei = ei;
+ break;
+ }
+ }
+ }
+
+
+ // Walk in other direction
+ {
+
+ facei = startFacei;
+ edgei = e1;
+
+ while (!isBorderEdge_[edgei])
+ {
+ // Cross edge to next face
+ const labelList& eFaces = edgeFaces()[edgei];
+ if (eFaces.size() != 2)
+ {
+ break;
+ }
+
+ forAll(eFaces, i)
+ {
+ if (eFaces[i] != facei)
+ {
+ facei = eFaces[i];
+ break;
+ }
+ }
+
+ if (facei == startFacei)
+ {
+ break;
+ }
+
+ n += s[facei].unitNormal(points());
+
+ // Cross face to next edge
+ const labelList& fEdges = faceEdges()[facei];
+
+ forAll(fEdges, i)
+ {
+ label ei = fEdges[i];
+ if (ei != edgei && edges[ei].otherVertex(localPointi) != -1)
+ {
+ edgei = ei;
+ break;
+ }
+ }
+ }
+ }
+
+ return n.normalise();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::smoothTriSurfaceMesh::smoothTriSurfaceMesh
+(
+ const IOobject& io,
+ const scalar featureAngle
+)
+:
+ triSurfaceMesh(io),
+ isBorderEdge_(nEdges()),
+ isPointOnBorderEdge_(nPoints())
+{
+ calcFeatureEdges(featureAngle);
+}
+
+
+Foam::smoothTriSurfaceMesh::smoothTriSurfaceMesh
+(
+ const IOobject& io,
+ const dictionary& dict
+)
+:
+ triSurfaceMesh(io, dict),
+ isBorderEdge_(nEdges()),
+ isPointOnBorderEdge_(nPoints())
+{
+ if (dict.found("featureAngle"))
+ {
+ calcFeatureEdges(readScalar(dict.lookup("featureAngle")));
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::smoothTriSurfaceMesh::~smoothTriSurfaceMesh()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::smoothTriSurfaceMesh::getNormal
+(
+ const List<pointIndexHit>& info,
+ vectorField& normal
+) const
+{
+ const triSurface& s = static_cast<const triSurface&>(*this);
+ const pointField& pts = s.points();
+
+ DebugPout<< "smoothTriSurfaceMesh::getNormal :"
+ << " surface:" << this->searchableSurface::name()
+ << " tris:" << s.size()
+ << " feat edges:" << isBorderEdge_.count()
+ << endl;
+
+ normal.setSize(info.size());
+
+ forAll(info, i)
+ {
+ if (info[i].hit())
+ {
+ label facei = info[i].index();
+
+ // Get local coordinates in triangle
+ barycentric2D coordinates
+ (
+ s[facei].tri(pts).pointToBarycentric(info[i].hitPoint())
+ );
+
+ // Average point normals
+ const triFace& localTri = s.localFaces()[facei];
+
+ vector n(Zero);
+ forAll(localTri, fp)
+ {
+ n += coordinates[fp]*pointNormal(facei, localTri[fp]);
+ }
+ normal[i] = n.normalise();
+ }
+ else
+ {
+ // Set to what?
+ normal[i] = Zero;
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/searchableSurfaces/triSurfaceMesh/smoothTriSurfaceMesh.H b/src/meshTools/searchableSurfaces/triSurfaceMesh/smoothTriSurfaceMesh.H
new file mode 100644
index 0000000000000000000000000000000000000000..772704b3c58c2509b7592b11685e0b8163d579e3
--- /dev/null
+++ b/src/meshTools/searchableSurfaces/triSurfaceMesh/smoothTriSurfaceMesh.H
@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2016-2018 OpenFOAM Foundation
+ Copyright (C) 2025 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software: you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
+
+Class
+ Foam::smoothTriSurfaceMesh
+
+Description
+ Variant of triSurfaceMesh that interpolates face-normals to obtain
+ point normals.
+
+ Note: when constructing from dictionary has optional parameters:
+ - scale : scaling factor.
+ - tolerance : relative tolerance for doing intersections
+ (see triangle::intersection)
+ - minQuality: discard triangles with low quality when getting normal
+ - featureAngle: limit interpolation to current 'region' only.
+ Region limited by feature edges.
+
+
+SourceFiles
+ smoothTriSurfaceMesh.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef smoothTriSurfaceMesh_H
+#define smoothTriSurfaceMesh_H
+
+#include "triSurfaceMesh.H"
+#include "PackedBoolList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class smoothTriSurfaceMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class smoothTriSurfaceMesh
+:
+ public triSurfaceMesh
+{
+private:
+
+ //- Edges with discontinuity
+ bitSet isBorderEdge_;
+
+ //- Points on border edges
+ bitSet isPointOnBorderEdge_;
+
+ //- Mark feature edges
+ void calcFeatureEdges(const scalar featureAngle);
+
+ //- Calculate local point normal
+ vector pointNormal(const label, const label) const;
+
+
+ //- Disallow default bitwise copy construct
+ smoothTriSurfaceMesh(const smoothTriSurfaceMesh&) = delete;
+
+ //- Disallow default bitwise assignment
+ void operator=(const smoothTriSurfaceMesh&) = delete;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("smoothTriSurfaceMesh");
+
+
+ // Constructors
+
+ //- Construct read, specify feature angle
+ smoothTriSurfaceMesh(const IOobject& io, const scalar featureAngle);
+
+ //- Construct from IO and dictionary (used by searchableSurface).
+ // Dictionary may contain a 'scale' entry (eg, 0.001: mm -> m)
+ smoothTriSurfaceMesh
+ (
+ const IOobject& io,
+ const dictionary& dict
+ );
+
+
+ //- Destructor
+ virtual ~smoothTriSurfaceMesh();
+
+
+ // Member Functions
+
+ //- From a set of points and indices get the normal
+ virtual void getNormal
+ (
+ const List<pointIndexHit>&,
+ vectorField& normal
+ ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/meshTools/triSurface/orientedSurface/orientedSurface.H b/src/meshTools/triSurface/orientedSurface/orientedSurface.H
index e3e3bf38e0fd00dd79233dacf09fff8b71117c4d..2fda96f3083577036a6b2f86067d17877cabd238 100644
--- a/src/meshTools/triSurface/orientedSurface/orientedSurface.H
+++ b/src/meshTools/triSurface/orientedSurface/orientedSurface.H
@@ -126,10 +126,6 @@ private:
// anything flipped.
static bool flipSurface(triSurface& s, const labelList& flipState);
- //- Make surface surface has consistent orientation across connected
- // triangles.
- static bool orientConsistent(triSurface& s);
-
public:
@@ -176,6 +172,10 @@ public:
const point& samplePoint,
const bool orientOutside
);
+
+ //- Make sure surface has consistent orientation across connected
+ // triangles. Does not change anything if it has.
+ static bool orientConsistent(triSurface& s);
};
diff --git a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
index 888e4dbf49c9a74f3b31694e56fba2d7639b86bd..120dd894fb3d864b26a27a404b81951c648b8e6d 100644
--- a/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
+++ b/src/parallel/distributed/distributedTriSurfaceMesh/distributedTriSurfaceMesh.C
@@ -44,6 +44,7 @@ License
#include "OBJstream.H"
#include "labelBits.H"
#include "profiling.H"
+#include "orientedSurface.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -3114,6 +3115,16 @@ Foam::distributedTriSurfaceMesh::distributedTriSurfaceMesh(const IOobject& io)
bounds().reduce();
+ if (readFromMaster)
+ {
+ // Make sure orientation is consistent in case we
+ // want to use for inside/outside testing
+ DebugInFunction
+ << "Forcing consistent normals since undecomposed" << endl;
+ const triSurface& surf = *this;
+ orientedSurface::orientConsistent(const_cast<triSurface&>(surf));
+ }
+
if (readFromMaster && !decomposeUsingBbs_)
{
// No fill-in so store normals