From b79cbfe49fc07ca0d367887178bb602de57c9a1b Mon Sep 17 00:00:00 2001
From: Mattijs Janssens <ext-mjanssens@esi-group.com>
Date: Thu, 15 May 2025 17:21:53 +0100
Subject: [PATCH] ENH: snappyHexMesh: work better on non-manifold. Fixes ##3361
---
.../shellSurfaces/shellSurfaces.C | 14 +-
src/meshTools/Make/files | 1 +
.../triSurfaceMesh/smoothTriSurfaceMesh.C | 318 ++++++++++++++++++
.../triSurfaceMesh/smoothTriSurfaceMesh.H | 132 ++++++++
.../orientedSurface/orientedSurface.H | 8 +-
.../distributedTriSurfaceMesh.C | 11 +
6 files changed, 475 insertions(+), 9 deletions(-)
create mode 100644 src/meshTools/searchableSurfaces/triSurfaceMesh/smoothTriSurfaceMesh.C
create mode 100644 src/meshTools/searchableSurfaces/triSurfaceMesh/smoothTriSurfaceMesh.H
diff --git a/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C b/src/mesh/snappyHexMesh/shellSurfaces/shellSurfaces.C
index a6547fa1142..221c5c6581a 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 f207200fabc..eb64fa51a0d 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 00000000000..0c6dbec6522
--- /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 00000000000..772704b3c58
--- /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 e3e3bf38e0f..2fda96f3083 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 888e4dbf49c..120dd894fb3 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
--
GitLab