From c9c85d9afebd385dd1b5dd81b499f1c7953dd6d9 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 31 Mar 2021 12:11:20 +0100
Subject: [PATCH] BUG: primitiveMesh: incorrect uncached edgeFaces. Fixes
#2047.
It was only looking for faces that were used in both
endpoints but not actually checking whether they were indeed
an edge (== consecutive vertex) in all faces. So if one
face had an additional crossing edge and another didn't it
would find more edgeFaces than the proper
'primitiveMesh::edgeFaces()' routine.
This occasionally happened inside snappyHexMesh
(e.g. motorBike tutorial)
---
.../mesh/advanced/modifyMesh/modifyMesh.C | 56 ++++++++++++++++++-
etc/caseDicts/annotated/modifyMeshDict | 10 ++++
.../primitiveMesh/primitiveMeshEdgeFaces.C | 35 +++++++++---
3 files changed, 92 insertions(+), 9 deletions(-)
diff --git a/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C b/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C
index c3b0c7a320c..048359375f2 100644
--- a/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C
+++ b/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C
@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
- Copyright (C) 2016-2017 OpenCFD Ltd.
+ Copyright (C) 2016-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -365,12 +365,16 @@ int main(int argc, char *argv[])
(
dict.lookup("facesToTriangulate")
);
+ // Optional
+ List<Pair<point>> facesToSplit;
+ dict.readIfPresent("facesToSplit", facesToSplit);
bool cutBoundary =
(
pointsToMove.size()
|| edgesToSplit.size()
|| facesToTriangulate.size()
+ || facesToSplit.size()
);
List<Pair<point>> edgesToCollapse(dict.lookup("edgesToCollapse"));
@@ -388,6 +392,7 @@ int main(int argc, char *argv[])
<< " Boundary cutting module:" << nl
<< " points to move :" << pointsToMove.size() << nl
<< " edges to split :" << edgesToSplit.size() << nl
+ << " faces to split :" << facesToSplit.size() << nl
<< " faces to triangulate:" << facesToTriangulate.size() << nl
<< " Cell splitting module:" << nl
<< " cells to split :" << cellsToPyramidise.size() << nl
@@ -481,6 +486,53 @@ int main(int argc, char *argv[])
}
+ Info<< nl << "Looking up faces to split ..." << nl << endl;
+ Map<labelPair> faceToSplit(facesToSplit.size());
+ for (const Pair<point>& pts : facesToSplit)
+ {
+ label facei = findFace(mesh, allBoundary, pts.first());
+ if (facei == -1)
+ {
+ Info<< "Could not insert mesh face " << facei
+ << " for input point " << pts.first() << nl
+ << "Perhaps the face is already marked for splitting?" << endl;
+
+ validInputs = false;
+ }
+ else
+ {
+ // Find nearest points on face
+ const primitivePatch pp
+ (
+ SubList<face>(mesh.faces(), 1, facei),
+ mesh.points()
+ );
+
+ const label p0 = findPoint(pp, pts.first());
+ const label p1 = findPoint(pp, pts.second());
+
+ const face& f = mesh.faces()[facei];
+
+ if
+ (
+ p0 != -1
+ && p1 != -1
+ && faceToSplit.insert(facei, labelPair(f.find(p0), f.find(p1)))
+ )
+ {}
+ else
+ {
+ Info<< "Could not insert mesh face " << facei
+ << " for input coordinates " << pts
+ << " with vertices " << p0 << " and " << p1 << nl
+ << "Perhaps the face is already marked for splitting?"
+ << endl;
+
+ }
+ }
+ }
+DebugVar(faceToSplit);
+
Info<< nl << "Looking up cells to convert to pyramids around"
<< " cell centre ..." << nl << endl;
@@ -652,7 +704,7 @@ int main(int argc, char *argv[])
(
pointToPos,
edgeToCuts,
- Map<labelPair>(0), // Faces to split diagonally
+ faceToSplit, // Faces to split diagonally
faceToDecompose, // Faces to triangulate
meshMod
);
diff --git a/etc/caseDicts/annotated/modifyMeshDict b/etc/caseDicts/annotated/modifyMeshDict
index 4080792b681..00b6660413e 100644
--- a/etc/caseDicts/annotated/modifyMeshDict
+++ b/etc/caseDicts/annotated/modifyMeshDict
@@ -30,6 +30,16 @@ edgesToSplit
(( -0.17692 -0.45312 0.74516)( -0.18 -0.45 0.742))
);
+// Split a boundary face:
+// The two coordinates indicate an edge on a face:
+// - find the face where the first coordinate is on (so the point should be
+// on the plane of the face
+// - find nearest vertex on the face for both coordinates
+facesToSplit
+(
+ ((0.99 0 0.01)(0.99 0 0.99)) // edge (1 0 0)(1 0 1) on bottom face
+);
+
// Triangulate a boundary face:
// First coord is a point on the face to triangulate. It will introduce a
// point on the face, triangulate and move the point to the second coordinate.
diff --git a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeFaces.C b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeFaces.C
index c7e7cd62ddc..a7e84341cf5 100644
--- a/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeFaces.C
+++ b/src/OpenFOAM/meshes/primitiveMesh/primitiveMeshEdgeFaces.C
@@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2015 OpenFOAM Foundation
+ Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -68,7 +69,9 @@ const Foam::labelList& Foam::primitiveMesh::edgeFaces
}
else
{
- // Use the fact that pointEdges are sorted in incrementing edge order
+ // Use the fact that pointFaces are sorted in incrementing edge order
+ // (since they get constructed by inverting the faces which walks
+ // in increasing face order)
const edge& e = edges()[edgeI];
const labelList& pFaces0 = pointFaces()[e[0]];
const labelList& pFaces1 = pointFaces()[e[1]];
@@ -80,20 +83,38 @@ const Foam::labelList& Foam::primitiveMesh::edgeFaces
while (i0 < pFaces0.size() && i1 < pFaces1.size())
{
- if (pFaces0[i0] < pFaces1[i1])
+ const label f0 = pFaces0[i0];
+ const label f1 = pFaces1[i1];
+
+ if (f0 < f1)
{
++i0;
}
- else if (pFaces0[i0] > pFaces1[i1])
+ else if (f0 > f1)
{
++i1;
}
else
{
- // Equal. Append.
- storage.append(pFaces0[i0]);
- ++i0;
- ++i1;
+ // Equal face. Check if indeed on consecutive vertices on both
+ // faces since it could be that there is an 'ear' where one
+ // side is a triangle and the other side is part of a bigger
+ // face (e.g. quad). Now all vertex-vertex pairs on the
+ // triangle are edges but there is no cross connection on the
+ // bigger face.
+ const face& f = faces()[f0];
+ const label fp0 = f.find(e[0]);
+
+ if (f[f.fcIndex(fp0)] == e[1] || f[f.rcIndex(fp0)] == e[1])
+ {
+ storage.append(f0);
+ ++i0;
+ ++i1;
+ }
+ else
+ {
+ ++i1;
+ }
}
}
--
GitLab