diff --git a/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C b/applications/utilities/mesh/advanced/modifyMesh/modifyMesh.C index c3b0c7a320..048359375f 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> facesToSplit; + dict.readIfPresent("facesToSplit", facesToSplit); bool cutBoundary = ( pointsToMove.size() || edgesToSplit.size() || facesToTriangulate.size() + || facesToSplit.size() ); List> 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 faceToSplit(facesToSplit.size()); + for (const Pair& 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(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(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 4080792b68..00b6660413 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 c7e7cd62dd..a7e84341cf 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; + } } }