diff --git a/applications/test/polygonTriangulate/Make/files b/applications/test/polygonTriangulate/Make/files
new file mode 100644
index 0000000000..0698af40d0
--- /dev/null
+++ b/applications/test/polygonTriangulate/Make/files
@@ -0,0 +1,3 @@
+Test-polygonTriangulate.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-polygonTriangulate
diff --git a/applications/test/polygonTriangulate/Make/options b/applications/test/polygonTriangulate/Make/options
new file mode 100644
index 0000000000..7ce182425d
--- /dev/null
+++ b/applications/test/polygonTriangulate/Make/options
@@ -0,0 +1,5 @@
+EXE_INC = \
+ -I$(LIB_SRC)/fileFormats/lnInclude
+
+EXE_LIBS = \
+ -lfileFormats
diff --git a/applications/test/polygonTriangulate/Test-polygonTriangulate.C b/applications/test/polygonTriangulate/Test-polygonTriangulate.C
new file mode 100644
index 0000000000..4f8ad2024e
--- /dev/null
+++ b/applications/test/polygonTriangulate/Test-polygonTriangulate.C
@@ -0,0 +1,121 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration | Website: https://openfoam.org
+ \\ / A nd | Copyright (C) 2021 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+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 (1);
+
+ // Initialise the random number generator
+ label seed;
+ if (args.optionFound("seed"))
+ {
+ seed = args.optionRead("seed");
+ }
+ else
+ {
+ seed = clock::getTime();
+ Info<< "Seeding random number generator with value " << seed
+ << nl << endl;
+ }
+ Random rndGen(seed);
+
+ // Get controls
+ const bool simple = args.optionFound("simple");
+ const bool nonOptimal = args.optionFound("nonOptimal");
+
+ // Generate a random polygon
+ const List polygon =
+ polygonTriangulate::randomPolygon
+ (
+ rndGen,
+ n,
+ args.optionLookupOrDefault("error", 0)
+ );
+
+ // Write the polygon
+ {
+ OBJstream os(args[0] + "_polygon.obj");
+ os.write(face(identity(polygon.size())), polygon, false);
+ }
+
+ // Triangulate the polygon
+ const List triPoints =
+ polygonTriangulate().triangulate(polygon, simple, !nonOptimal);
+
+ // Write the triangulation
+ faceList triFacePoints(triPoints.size());
+ forAll(triPoints, trii)
+ {
+ triFacePoints[trii] = triPoints[trii];
+ }
+ {
+ OBJstream os(args[0] + "_triangulation.obj");
+ os.write(triFacePoints, pointField(polygon), false);
+ }
+
+ Info<< "End" << nl << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK/vtkTopo.C b/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK/vtkTopo.C
index e314eb8e6b..2b90ced102 100644
--- a/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK/vtkTopo.C
+++ b/applications/utilities/postProcessing/dataConversion/foamToVTK/foamToVTK/vtkTopo.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -33,6 +33,7 @@ Description
#include "cellShape.H"
#include "cellModeller.H"
#include "Swap.H"
+#include "polygonTriangulate.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -91,11 +92,14 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
{
const face& f = mesh_.faces()[cFaces[cFacei]];
- label nQuads = 0;
- label nTris = 0;
- f.nTrianglesQuads(mesh_.points(), nTris, nQuads);
-
- nAddCells += nQuads + nTris;
+ if (f.size() != 4)
+ {
+ nAddCells += f.nTriangles();
+ }
+ else
+ {
+ nAddCells ++;
+ }
}
nAddCells--;
@@ -119,6 +123,9 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
// Label of vtk type
cellTypes_.setSize(cellShapes.size() + nAddCells);
+ // Create a triangulation engine
+ polygonTriangulate triEngine;
+
// Set counters for additional points and additional cells
label addPointi = 0, addCelli = 0;
@@ -205,17 +212,24 @@ Foam::vtkTopo::vtkTopo(const polyMesh& mesh)
const face& f = mesh_.faces()[cFaces[cFacei]];
const bool isOwner = (owner[cFaces[cFacei]] == celli);
- // Number of triangles and quads in decomposition
- label nTris = 0;
- label nQuads = 0;
- f.nTrianglesQuads(mesh_.points(), nTris, nQuads);
-
- // Do actual decomposition into triFcs and quadFcs.
- faceList triFcs(nTris);
- faceList quadFcs(nQuads);
- label trii = 0;
- label quadi = 0;
- f.trianglesQuads(mesh_.points(), trii, quadi, triFcs, quadFcs);
+ // Do decomposition into triFcs and quadFcs.
+ faceList triFcs(f.size() == 4 ? 0 : f.nTriangles());
+ faceList quadFcs(f.size() == 4 ? 1 : 0);
+ if (f.size() != 4)
+ {
+ triEngine.triangulate
+ (
+ UIndirectList(mesh.points(), f)
+ );
+ forAll(triEngine.triPoints(), trii)
+ {
+ triFcs[trii] = triEngine.triPoints(trii, f);
+ }
+ }
+ else
+ {
+ quadFcs[0] = f;
+ }
forAll(quadFcs, quadI)
{
diff --git a/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamMeshVolume.C b/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamMeshVolume.C
index 895e3533c8..930364601f 100644
--- a/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamMeshVolume.C
+++ b/applications/utilities/postProcessing/graphics/PVReaders/vtkPVFoam/vtkPVFoamMeshVolume.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -31,6 +31,7 @@ License
#include "cellModeller.H"
#include "vtkOpenFOAMPoints.H"
#include "Swap.H"
+#include "polygonTriangulate.H"
// VTK includes
#include "vtkCellArray.h"
@@ -104,11 +105,14 @@ vtkUnstructuredGrid* Foam::vtkPVFoam::volumeVTKMesh
{
const face& f = mesh.faces()[cFaces[cFacei]];
- label nQuads = 0;
- label nTris = 0;
- f.nTrianglesQuads(mesh.points(), nTris, nQuads);
-
- nAddCells += nQuads + nTris;
+ if (f.size() != 4)
+ {
+ nAddCells += f.nTriangles();
+ }
+ else
+ {
+ nAddCells ++;
+ }
}
nAddCells--;
@@ -158,6 +162,9 @@ vtkUnstructuredGrid* Foam::vtkPVFoam::volumeVTKMesh
vtkmesh->Allocate(mesh.nCells() + nAddCells);
+ // Create a triangulation engine
+ polygonTriangulate triEngine;
+
// Set counters for additional points and additional cells
label addPointi = 0, addCelli = 0;
@@ -369,17 +376,24 @@ vtkUnstructuredGrid* Foam::vtkPVFoam::volumeVTKMesh
const face& f = mesh.faces()[cFaces[cFacei]];
const bool isOwner = (owner[cFaces[cFacei]] == celli);
- // Number of triangles and quads in decomposition
- label nTris = 0;
- label nQuads = 0;
- f.nTrianglesQuads(mesh.points(), nTris, nQuads);
-
- // Do actual decomposition into triFcs and quadFcs.
- faceList triFcs(nTris);
- faceList quadFcs(nQuads);
- label trii = 0;
- label quadi = 0;
- f.trianglesQuads(mesh.points(), trii, quadi, triFcs, quadFcs);
+ // Do decomposition into triFcs and quadFcs.
+ faceList triFcs(f.size() == 4 ? 0 : f.nTriangles());
+ faceList quadFcs(f.size() == 4 ? 1 : 0);
+ if (f.size() != 4)
+ {
+ triEngine.triangulate
+ (
+ UIndirectList(mesh.points(), f)
+ );
+ forAll(triEngine.triPoints(), trii)
+ {
+ triFcs[trii] = triEngine.triPoints(trii, f);
+ }
+ }
+ else
+ {
+ quadFcs[0] = f;
+ }
forAll(quadFcs, quadI)
{
diff --git a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
index 3a3ec81cf4..a1da8660b4 100644
--- a/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
+++ b/applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -52,6 +52,7 @@ Description
#include "uindirectPrimitivePatch.H"
#include "DynamicField.H"
#include "scalarListIOList.H"
+#include "polygonTriangulate.H"
using namespace Foam;
@@ -78,31 +79,30 @@ triSurface triangulate
label newPatchi = 0;
label localTriFacei = 0;
+ polygonTriangulate triEngine;
+
forAllConstIter(labelHashSet, includePatches, iter)
{
const label patchi = iter.key();
const polyPatch& patch = bMesh[patchi];
const pointField& points = patch.points();
- label nTriTotal = 0;
-
forAll(patch, patchFacei)
{
const face& f = patch[patchFacei];
- faceList triFaces(f.nTriangles(points));
+ triEngine.triangulate(UIndirectList(points, f));
- label nTri = 0;
-
- f.triangles(points, nTri, triFaces);
-
- forAll(triFaces, triFacei)
+ forAll(triEngine.triPoints(), triFacei)
{
- const face& f = triFaces[triFacei];
-
- triangles.append(labelledTri(f[0], f[1], f[2], newPatchi));
-
- nTriTotal++;
+ triangles.append
+ (
+ labelledTri
+ (
+ triEngine.triPoints(triFacei, f),
+ newPatchi
+ )
+ );
triSurfaceToAgglom[localTriFacei++] = globalNumbering.toGlobal
(
diff --git a/src/OpenFOAM/Make/files b/src/OpenFOAM/Make/files
index d077674309..38f0be1ecd 100644
--- a/src/OpenFOAM/Make/files
+++ b/src/OpenFOAM/Make/files
@@ -710,6 +710,7 @@ $(interpolationWeights)/splineInterpolationWeights/splineInterpolationWeights.C
algorithms/indexedOctree/indexedOctreeName.C
algorithms/indexedOctree/treeDataCell.C
algorithms/indexedOctree/volumeType.C
+algorithms/polygonTriangulate/polygonTriangulate.C
algorithms/dynamicIndexedOctree/dynamicIndexedOctreeName.C
diff --git a/src/OpenFOAM/algorithms/polygonTriangulate/polygonTriangulate.C b/src/OpenFOAM/algorithms/polygonTriangulate/polygonTriangulate.C
new file mode 100644
index 0000000000..c79333ac57
--- /dev/null
+++ b/src/OpenFOAM/algorithms/polygonTriangulate/polygonTriangulate.C
@@ -0,0 +1,919 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration | Website: https://openfoam.org
+ \\ / A nd | Copyright (C) 2021 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+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
+Foam::scalar Foam::polygonTriangulate::area
+(
+ const triFace& triPoints,
+ const PointField& points,
+ const vector& normal
+)
+{
+ const point& o = points[triPoints[0]];
+ const vector ao = points[triPoints[1]] - o;
+ const vector ab = points[triPoints[2]] - o;
+
+ return (ao ^ ab) & normal;
+}
+
+
+template
+Foam::scalar Foam::polygonTriangulate::quality
+(
+ const triFace& triPoints,
+ const PointField& points,
+ const vector& normal
+)
+{
+ const scalar An = area(triPoints, points, normal);
+
+ scalar PSqr = 0;
+ forAll(triPoints, triEdgei)
+ {
+ const point& p0 = points[triPoints[triEdgei]];
+ const point& p1 = points[triPoints[(triEdgei + 1) % 3]];
+ PSqr += magSqr(p0 - p1);
+ }
+
+ static const scalar equilateralAnByPSqr = sqrt(scalar(3))/12;
+
+ return PSqr != 0 ? An/PSqr/equilateralAnByPSqr : 0;
+}
+
+
+template
+bool Foam::polygonTriangulate::intersection
+(
+ const edge& edgePointsA,
+ const edge& edgePointsB,
+ const PointField& points,
+ const vector& normal
+)
+{
+ const vector tau0 = perpendicular(normal);
+ const vector tau1 = normal ^ tau0;
+
+ const point& pointA0 = points[edgePointsA[0]];
+ const point& pointA1 = points[edgePointsA[1]];
+ const point& pointB0 = points[edgePointsB[0]];
+ const point& pointB1 = points[edgePointsB[1]];
+
+ const point2D point2DA0(tau0 & pointA0, tau1 & pointA0);
+ const point2D point2DA1(tau0 & pointA1, tau1 & pointA1);
+ const point2D point2DB0(tau0 & pointB0, tau1 & pointB0);
+ const point2D point2DB1(tau0 & pointB1, tau1 & pointB1);
+
+ const tensor2D M =
+ tensor2D(point2DA0 - point2DA1, point2DB1 - point2DB0).T();
+
+ const scalar detM = det(M);
+ const tensor2D detMInvM = cof(M);
+
+ const vector2D magDetMTu = sign(detM)*detMInvM & (point2DA0 - point2DB0);
+
+ return 0 <= cmptMin(magDetMTu) && cmptMax(magDetMTu) <= mag(detM);
+}
+
+
+template
+Foam::label Foam::polygonTriangulate::nIntersections
+(
+ const edge& edgePoints,
+ const PointField& points,
+ const vector& normal
+)
+{
+ label n = 0;
+
+ forAll(points, pi)
+ {
+ const edge otherEdgePoints(pi, points.fcIndex(pi));
+
+ if (edgePoints.commonVertex(otherEdgePoints) == -1)
+ {
+ n += intersection(edgePoints, otherEdgePoints, points, normal);
+ }
+ }
+
+ return n;
+}
+
+
+template
+Foam::scalar Foam::polygonTriangulate::angle
+(
+ const label pointi,
+ const PointField& points,
+ const vector& normal
+)
+{
+ const point& o = points[pointi];
+ const vector oa = points[points.rcIndex(pointi)] - o;
+ const vector ob = points[points.fcIndex(pointi)] - o;
+
+ const vector oaNegNNOa = oa - normal*(normal & oa);
+ const vector obNegNNOb = ob - normal*(normal & ob);
+
+ return
+ atan2(normal & (oaNegNNOa ^ obNegNNOb), - oaNegNNOa & obNegNNOb)
+ + constant::mathematical::pi;
+}
+
+
+template
+bool Foam::polygonTriangulate::ear
+(
+ const label pointi,
+ const PointField& points,
+ const vector& normal
+)
+{
+ const point& o = points[pointi];
+ const vector oa = points[points.rcIndex(pointi)] - o;
+ const vector ob = points[points.fcIndex(pointi)] - o;
+
+ const tensor A = tensor(ob, oa, normal).T();
+ const tensor T(A.y() ^ A.z(), A.z() ^ A.x(), A.x() ^ A.y());
+ const scalar detA = det(A);
+
+ for
+ (
+ label pointj = points.fcIndex(points.fcIndex(pointi));
+ pointj != points.rcIndex(pointi);
+ pointj = points.fcIndex(pointj)
+ )
+ {
+ const vector detAY = (points[pointj] - o) & T;
+
+ if (detAY.x() > 0 && detAY.y() > 0 && detAY.x() + detAY.y() < detA)
+ {
+ return false;
+ }
+ }
+
+ return true;
+}
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+Foam::List Foam::polygonTriangulate::randomPolygon
+(
+ Random& rndGen,
+ const label n,
+ const scalar error
+)
+{
+ // Get random points on a unit disk
+ List points(n);
+ forAll(points, pointi)
+ {
+ const scalar theta =
+ 2*constant::mathematical::pi*rndGen.sample01();
+ const scalar r = sqrt(rndGen.sample01());
+ points[pointi] = point(r*cos(theta), r*sin(theta), 0);
+ }
+
+ // Initialise intersected polygon
+ labelList pointis(identity(n));
+
+ // Reorder until non-intersected
+ bool incomplete = true;
+ while (incomplete)
+ {
+ incomplete = false;
+
+ for (label piA0 = 0; piA0 < n; ++ piA0)
+ {
+ const label piA1 = points.fcIndex(piA0);
+
+ for (label piB0 = piA0 + 2; piB0 < (piA0 == 0 ? n - 1 : n); ++ piB0)
+ {
+ const label piB1 = points.fcIndex(piB0);
+
+ if
+ (
+ intersection
+ (
+ edge(pointis[piA0], pointis[piA1]),
+ edge(pointis[piB0], pointis[piB1]),
+ points,
+ vector(0, 0, 1)
+ )
+ )
+ {
+ SubList subOrder(pointis, piB0 + 1 - piA1, piA1);
+ inplaceReverseList(subOrder);
+ incomplete = true;
+ }
+
+ if (incomplete) break;
+ }
+
+ if (incomplete) break;
+ }
+ }
+
+ // Add noise and potential self-intersection
+ forAll(points, pointi)
+ {
+ const scalar theta =
+ 2*constant::mathematical::pi*rndGen.sample01();
+ const scalar r = sqrt(rndGen.sample01());
+ points[pointi] =
+ (1 - error)*points[pointi]
+ + error*point(r*cos(theta), r*sin(theta), rndGen.sample01());
+ }
+
+ // Reorder and return
+ return List(points, pointis);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+template
+void Foam::polygonTriangulate::optimiseTriangulation
+(
+ const label trii,
+ const PointField& points,
+ const vector& normal,
+ UList& triPoints,
+ UList>& triEdges,
+ UList& edgeTris
+)
+{
+ const triFace& t = triPoints[trii];
+ const scalar q = quality(t, points, normal);
+
+ forAll(triEdges[trii], triEdgei)
+ {
+ const label edgei = triEdges[trii][triEdgei];
+
+ const label otherTrii = edgeTris[edgei][edgeTris[edgei][0] == trii];
+
+ if (otherTrii == -1) continue;
+
+ const label otherTriEdgei = findIndex(triEdges[otherTrii], edgei);
+
+ const triFace& otherT = triPoints[otherTrii];
+ const scalar otherQ = quality(otherT, points, normal);
+
+ const label piA = triPoints[trii][(triEdgei + 1) % 3];
+ const label piB = triPoints[trii][(triEdgei + 2) % 3];
+ const label piC = triPoints[otherTrii][(otherTriEdgei + 1) % 3];
+ const label piD = triPoints[otherTrii][(otherTriEdgei + 2) % 3];
+
+ const triFace flipT0(piA, piB, piD);
+ const triFace flipT1(piC, piD, piB);
+
+ if
+ (
+ triPointsSet_.found(flipT0)
+ || triPointsSet_.found(flipT1)
+ ) continue;
+
+ const scalar flipQ0 = quality(flipT0, points, normal);
+ const scalar flipQ1 = quality(flipT1, points, normal);
+
+ if (min(q, otherQ) < min(flipQ0, flipQ1))
+ {
+ triPointsSet_.set(t);
+ triPointsSet_.set(otherT);
+
+ const label eiA = triEdges[trii][(triEdgei + 1) % 3];
+ const label eiB = triEdges[trii][(triEdgei + 2) % 3];
+ const label eiC = triEdges[otherTrii][(otherTriEdgei + 1) % 3];
+ const label eiD = triEdges[otherTrii][(otherTriEdgei + 2) % 3];
+
+ triPoints[trii] = {piA, piB, piD};
+ triEdges[trii] = {eiA, edgei, eiD};
+ edgeTris[eiD][edgeTris[eiD][0] != otherTrii] = trii;
+
+ triPoints[otherTrii] = {piC, piD, piB};
+ triEdges[otherTrii] = {eiC, edgei, eiB};
+ edgeTris[eiB][edgeTris[eiB][0] != trii] = otherTrii;
+
+ optimiseTriangulation
+ (
+ trii,
+ points,
+ normal,
+ triPoints,
+ triEdges,
+ edgeTris
+ );
+ optimiseTriangulation
+ (
+ otherTrii,
+ points,
+ normal,
+ triPoints,
+ triEdges,
+ edgeTris
+ );
+
+ break;
+ }
+ }
+}
+
+
+void Foam::polygonTriangulate::simpleTriangulate
+(
+ const UList& allPoints,
+ const vector& normal,
+ UList& triPoints,
+ UList>& triEdges,
+ UList& edgeTris,
+ const bool optimal
+)
+{
+ // Get the polygon size
+ const label n = allPoints.size();
+
+ // Clear the triangulation
+ triPoints = triFace(-1, -1, -1);
+ triEdges = FixedList({-1, -1, -1});
+ edgeTris = labelPair(-1, -1);
+
+ // Initialise workspace
+ pointis_.setSize(n);
+ edges_.setSize(n);
+ angle_.setSize(n);
+ ear_.setSize(n);
+ forAll(pointis_, i)
+ {
+ pointis_[i] = i;
+ edges_[i] = i;
+ angle_[i] = angle(i, allPoints, normal);
+ ear_[i] = ear(i, allPoints, normal);
+ }
+
+ // Create the indirect point list for the current un-triangulated polygon
+ UIndirectList points(allPoints, pointis_);
+
+ // Generate triangles
+ forAll(triPoints, trii)
+ {
+ // Find the ear with the smallest angle
+ scalar minEarAngle = vGreat;
+ label minEarAnglei = -1;
+ for (label i = 0; i < pointis_.size(); ++ i)
+ {
+ if (ear_[i] && minEarAngle > angle_[i])
+ {
+ minEarAngle = angle_[i];
+ minEarAnglei = i;
+ }
+ }
+
+ // If nothing is an ear, then this face is degenerate. Default to
+ // the smallest angle, ignoring the ear test.
+ if (minEarAnglei == -1)
+ {
+ minEarAnglei = findMin(angle_);
+ }
+
+ // Get the adjacent points
+ label minEarAnglei0 = pointis_.rcIndex(minEarAnglei);
+ label minEarAnglei1 = pointis_.fcIndex(minEarAnglei);
+
+ // Add a triangle
+ triPoints[trii] = triFace
+ (
+ pointis_[minEarAnglei0],
+ pointis_[minEarAnglei],
+ pointis_[minEarAnglei1]
+ );
+ triEdges[trii] = FixedList
+ ({
+ edges_[minEarAnglei0],
+ edges_[minEarAnglei],
+ trii == triPoints.size() - 1 ? edges_[minEarAnglei1] : n + trii
+ });
+ forAll(triEdges[trii], triEdgei)
+ {
+ const label edgei = triEdges[trii][triEdgei];
+ edgeTris[edgei][edgeTris[edgei][0] != -1] = trii;
+ }
+
+ // Remove the cut off point from the polygon
+ edges_[minEarAnglei0] = triEdges[trii][2];
+ for (label i = minEarAnglei; i < pointis_.size() - 1; ++ i)
+ {
+ pointis_[i] = pointis_[i + 1];
+ edges_[i] = edges_[i + 1];
+ angle_[i] = angle_[i + 1];
+ ear_[i] = ear_[i + 1];
+ }
+ pointis_.resize(pointis_.size() - 1);
+ edges_.resize(edges_.size() - 1);
+ angle_.resize(angle_.size() - 1);
+ ear_.resize(ear_.size() - 1);
+
+ // Update connected points
+ if (minEarAnglei0 > minEarAnglei) -- minEarAnglei0;
+ angle_[minEarAnglei0] = angle(minEarAnglei0, points, normal);
+ ear_[minEarAnglei0] = ear(minEarAnglei0, points, normal);
+ if (minEarAnglei1 > minEarAnglei) -- minEarAnglei1;
+ angle_[minEarAnglei1] = angle(minEarAnglei1, points, normal);
+ ear_[minEarAnglei1] = ear(minEarAnglei1, points, normal);
+
+ // Optimise the quality
+ if (optimal)
+ {
+ triPointsSet_.clear();
+ optimiseTriangulation
+ (
+ trii,
+ allPoints,
+ normal,
+ triPoints,
+ triEdges,
+ edgeTris
+ );
+ }
+ }
+}
+
+
+void Foam::polygonTriangulate::partitionTriangulate
+(
+ const UList& points,
+ const vector& normal,
+ const label spanPointi,
+ const label spanEdgei,
+ UList& triPoints,
+ UList>& triEdges,
+ UList& edgeTris,
+ const bool simple,
+ const bool optimal
+)
+{
+ // Get the polygon size
+ const label n = points.size();
+
+ // Get the points of the spanning triangle
+ const label piA = (spanEdgei + 1) % n;
+ const label piP = spanPointi;
+ const label piB = spanEdgei;
+
+ // Determine the sizes of the polygons connected to either side of the
+ // spanning triangle
+ const label nA = (piP - piA + n) % n + 1;
+ const label nB = (piB - piP + n) % n + 1;
+
+ // Get the edges of the spanning triangle
+ const label eiA = nA >= 3 ? n : piA;
+ const label eiB = nB >= 3 ? n + (nA >= 3) : piP;
+ const label eiE = piB;
+
+ // Add the spanning triangle
+ triPoints[0] = triFace(piA, piP, piB);
+ triEdges[0] = FixedList({eiA, eiB, eiE});
+
+ // Rotate the point list so that it can be passed to sub-triangulations
+ // without indirect addressing
+ inplaceRotateList(const_cast&>(points), - piA);
+
+ // Triangulate the sub-polygon connected to edge A
+ if (nA >= 3)
+ {
+ // Subset the polygon and triangulate
+ SubList triPointsA(triPoints, nA - 2, 1);
+ SubList> triEdgesA(triEdges, nA - 2, 1);
+ SubList edgeTrisA(edgeTris, 2*nA - 3);
+ triangulate
+ (
+ SubList(points, nA),
+ normal,
+ triPointsA,
+ triEdgesA,
+ edgeTrisA,
+ simple,
+ optimal
+ );
+
+ // Map the point labels back to the full polygon
+ forAll(triPointsA, triAi)
+ {
+ forAll(triPointsA[triAi], triPointAi)
+ {
+ label& pointi = triPointsA[triAi][triPointAi];
+ pointi = (pointi + piA) % n;
+ }
+ }
+
+ // Map the edge labels back to the full polygon
+ forAll(triEdgesA, triAi)
+ {
+ forAll(triEdgesA[triAi], triEdgeAi)
+ {
+ label& edgei = triEdgesA[triAi][triEdgeAi];
+ edgei =
+ edgei >= nA ? max(eiA, eiB) + edgei - nA + 1
+ : edgei == nA - 1 ? eiA
+ : (piA + edgei) % n;
+ }
+ }
+
+ // Map the tri labels back to the full polygon
+ forAll(edgeTrisA, edgeAi)
+ {
+ forAll(edgeTrisA[edgeAi], edgeTriAi)
+ {
+ label& trii = edgeTrisA[edgeAi][edgeTriAi];
+ trii = trii == -1 ? -1 : trii + 1;
+ }
+ }
+ }
+
+ // Triangulate the sub-polygon connected to edge B
+ if (nB >= 3)
+ {
+ // Subset the polygon and triangulate
+ SubList triPointsB(triPoints, nB - 2, nA - 1);
+ SubList> triEdgesB(triEdges, nB - 2, nA - 1);
+ SubList edgeTrisB(edgeTris, 2*nB - 3, 2*nA - 3);
+ triangulate
+ (
+ SubList(points, nB, nA - 1),
+ normal,
+ triPointsB,
+ triEdgesB,
+ edgeTrisB,
+ simple,
+ optimal
+ );
+
+ // Map the point labels back to the full polygon
+ forAll(triPointsB, triBi)
+ {
+ forAll(triPointsB[triBi], triPointBi)
+ {
+ label& pointi = triPointsB[triBi][triPointBi];
+ pointi = (pointi + piA + nA - 1) % n;
+ }
+ }
+
+ // Map the edge labels back to the full polygon
+ forAll(triEdgesB, triBi)
+ {
+ forAll(triEdgesB[triBi], triEdgeBi)
+ {
+ label& edgei = triEdgesB[triBi][triEdgeBi];
+ edgei =
+ edgei >= nB ? max(eiA, eiB) + nA - 3 + edgei - nB + 1
+ : edgei == nB - 1 ? eiB
+ : (piA + nA - 1 + edgei) % n;
+ }
+ }
+
+ // Map the tri labels back to the full polygon
+ forAll(edgeTrisB, edgeBi)
+ {
+ forAll(edgeTrisB[edgeBi], edgeTriBi)
+ {
+ label& trii = edgeTrisB[edgeBi][edgeTriBi];
+ trii = trii == -1 ? -1 : trii + nA - 1;
+ }
+ }
+ }
+
+ // Rotate the point list back
+ inplaceRotateList(const_cast&>(points), piA);
+
+ // Reorder the edge-tris
+ {
+ // Swap the A-internal-edges and B-outer-edges to get all outer
+ // edges and all internal edges in contiguous blocks
+ SubList l(edgeTris, nA + nB - 3, nA - 1);
+ inplaceRotateList(l, - nA + 2);
+ }
+ {
+ // Rotate the outer edges right by one so that the intersection
+ // edge (currently at the end) moves adjacent to the outer edges
+ SubList l(edgeTris, nA + nB - 3, nA + nB - 2);
+ inplaceRotateList(l, 1);
+ }
+ {
+ // Rotate the outer edges so that they are in sequence with the
+ // original point ordering
+ SubList l(edgeTris, n);
+ inplaceRotateList(l, - n + piA);
+ }
+ if (nA >= 3 && nB >= 3)
+ {
+ // Move edge-B leftwards adjacent to edge-A as this is created
+ // before any of the other internal edges
+ const label fromi = n + nA - 2, toi = n + 1;
+ const labelPair temp = edgeTris[fromi];
+ for (label i = fromi; i > toi; -- i)
+ {
+ edgeTris[i] = edgeTris[i-1];
+ }
+ edgeTris[toi] = temp;
+ }
+
+ // Add associations to the intersection triangle
+ edgeTris[eiA][edgeTris[eiA][0] != -1] = 0;
+ edgeTris[eiB][edgeTris[eiB][0] != -1] = 0;
+ edgeTris[eiE] = labelPair(0, -1);
+}
+
+
+void Foam::polygonTriangulate::complexTriangulate
+(
+ const UList& points,
+ const vector& normal,
+ UList& triPoints,
+ UList>& triEdges,
+ UList& edgeTris,
+ const bool optimal
+)
+{
+ // Get the polygon size
+ const label n = points.size();
+
+ // Detect self intersections. When one is found, remove one of the
+ // intersected edges by adding a spanning triangle and recurse. Pick the
+ // spanning triangle that results in the smallest negative area.
+ for (label piA0 = 0; piA0 < n; ++ piA0)
+ {
+ const label piA1 = points.fcIndex(piA0);
+
+ for (label piB0 = piA0 + 2; piB0 < (piA0 == 0 ? n - 1 : n); ++ piB0)
+ {
+ const label piB1 = points.fcIndex(piB0);
+
+ if
+ (
+ intersection
+ (
+ edge(piA0, piA1),
+ edge(piB0, piB1),
+ points,
+ normal
+ )
+ )
+ {
+ // Get a bunch of unique spanning triangles that remove one of
+ // the intersected edges
+ FixedList spanPointAndEdgeis
+ ({
+ labelPair(piA0, piB0),
+ labelPair(piA1, piB0),
+ labelPair(piB0, piA0),
+ labelPair(piB1, piA0),
+ labelPair(points.rcIndex(piA0), piA0),
+ labelPair(points.fcIndex(piA1), piA0),
+ labelPair(points.rcIndex(piB0), piB0),
+ labelPair(points.fcIndex(piB1), piB0)
+ });
+ forAll(spanPointAndEdgeis, spani)
+ {
+ for (label spanj = spani + 1; spanj < 8; ++ spanj)
+ {
+ if
+ (
+ spanPointAndEdgeis[spani]
+ == spanPointAndEdgeis[spanj]
+ )
+ {
+ spanPointAndEdgeis[spanj] = labelPair(-1, -1);
+ }
+ }
+ }
+
+ // Find the spanning triangle that results in the smallest
+ // negative triangulation area
+ scalar spanSumNegA = - vGreat;
+ label spanPointi = -1, spanEdgei = -1;
+ forAll(spanPointAndEdgeis, spani)
+ {
+ if (spanPointAndEdgeis[spani] == labelPair(-1, -1))
+ continue;
+
+ const label piP = spanPointAndEdgeis[spani].first();
+ const label ei = spanPointAndEdgeis[spani].second();
+ const label piA0 = points.rcIndex(ei);
+ const label piA = ei;
+ const label piB = points.fcIndex(ei);
+ const label piB1 = points.fcIndex(piB);
+
+ // Only consider this spanning triangle if it decreases
+ // the number of intersections
+ const label netNIntersections =
+ - nIntersections(edge(piA, piB), points, normal)
+ + (piP == piA0 ? -1 : +1)
+ *nIntersections(edge(piP, piA), points, normal)
+ + (piP == piB1 ? -1 : +1)
+ *nIntersections(edge(piP, piB), points, normal);
+ if (netNIntersections >= 0) continue;
+
+ // Triangulate
+ partitionTriangulate
+ (
+ points,
+ normal,
+ piP,
+ ei,
+ triPoints,
+ triEdges,
+ edgeTris,
+ false,
+ false
+ );
+
+ // Compute the negative area. If it is smaller than the
+ // previous value then update the spanning triangle. If it
+ // is zero then break and use this triangle.
+ scalar sumNegA = 0;
+ forAll(triPoints, trii)
+ {
+ const scalar a = area(triPoints[trii], points, normal);
+ sumNegA += negPart(a);
+ }
+ if (sumNegA > spanSumNegA)
+ {
+ spanSumNegA = sumNegA;
+ spanPointi = piP;
+ spanEdgei = ei;
+ }
+ if (spanSumNegA == 0)
+ {
+ break;
+ }
+ }
+
+ // If a suitable spanning triangle was found then partition the
+ // triangulation
+ if (spanPointi != -1)
+ {
+ partitionTriangulate
+ (
+ points,
+ normal,
+ spanPointi,
+ spanEdgei,
+ triPoints,
+ triEdges,
+ edgeTris,
+ false,
+ optimal
+ );
+
+ return;
+ }
+ }
+ }
+ }
+
+ // If there are no intersections then do simple polygon triangulation
+ simpleTriangulate(points, normal, triPoints, triEdges, edgeTris, optimal);
+}
+
+
+void Foam::polygonTriangulate::triangulate
+(
+ const UList& points,
+ const vector& normal,
+ UList& triPoints,
+ UList>& triEdges,
+ UList& edgeTris,
+ const bool simple,
+ const bool optimal
+)
+{
+ if (simple)
+ {
+ simpleTriangulate
+ (
+ points,
+ normal,
+ triPoints,
+ triEdges,
+ edgeTris,
+ optimal
+ );
+ }
+ else
+ {
+ complexTriangulate
+ (
+ points,
+ normal,
+ triPoints,
+ triEdges,
+ edgeTris,
+ optimal
+ );
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::polygonTriangulate::polygonTriangulate()
+:
+ pointis_(),
+ edges_(),
+ angle_(),
+ ear_(),
+
+ triPointsSet_(),
+
+ points_(),
+ triPoints_(),
+ triEdges_(),
+ edgeTris_()
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::polygonTriangulate::~polygonTriangulate()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+const Foam::UList& Foam::polygonTriangulate::triangulate
+(
+ const UList& points,
+ const vector& normal,
+ const bool simple,
+ const bool optimal
+)
+{
+ // Get the polygon size
+ const label n = points.size();
+
+ // Resize workspace
+ triPoints_.resize(n - 2);
+ triEdges_.resize(n - 2);
+ edgeTris_.resize(2*n - 3);
+
+ // Run
+ triangulate
+ (
+ points,
+ normal,
+ triPoints_,
+ triEdges_,
+ edgeTris_,
+ simple,
+ optimal
+ );
+
+ return triPoints_;
+}
+
+
+const Foam::UList& Foam::polygonTriangulate::triangulate
+(
+ const UList& points,
+ const bool simple,
+ const bool optimal
+)
+{
+ return
+ triangulate
+ (
+ points,
+ normalised(face::area(points)),
+ simple,
+ optimal
+ );
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/algorithms/polygonTriangulate/polygonTriangulate.H b/src/OpenFOAM/algorithms/polygonTriangulate/polygonTriangulate.H
new file mode 100644
index 0000000000..2003861672
--- /dev/null
+++ b/src/OpenFOAM/algorithms/polygonTriangulate/polygonTriangulate.H
@@ -0,0 +1,351 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration | Website: https://openfoam.org
+ \\ / A nd | Copyright (C) 2021 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+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 pointis_;
+
+ //- Edge indices of the un-triangulated polygon
+ DynamicList edges_;
+
+ //- The interior angles of the vertices of the un-triangulated polygon
+ DynamicList angle_;
+
+ //- Flags to denote whether points of the un-triangulated polygon are
+ // "ears"; i.e., that the triangle formed by adding a diagonal between
+ // two adjacent points contains no other point
+ DynamicList ear_;
+
+ //- Set of tri-faces used to mark previously considered configurations
+ // in optimisation and prevent unnecessary repetition
+ HashSet> triPointsSet_;
+
+ //- Point array. Used if the input points cannot be SubList-d.
+ DynamicList points_;
+
+ //- Tri-points
+ DynamicList triPoints_;
+
+ //- Tri-edge associations
+ DynamicList> triEdges_;
+
+ //- Edge-tri associations
+ DynamicList edgeTris_;
+
+
+ // Private Static Member Functions
+
+ // Renumber a polygon label to global indexing
+ static inline label renumberPolyToGlobal
+ (
+ const label triPoly,
+ const UList& polyToGlobal
+ );
+
+ // Renumber a container of polygon labels to global indexing
+ template
+ static inline Type renumberPolyToGlobal
+ (
+ const Type& triPoly,
+ const UList& polyToGlobal
+ );
+
+ //- Return the area of a triangle projected in a normal direction
+ template
+ static scalar area
+ (
+ const triFace& triPoints,
+ const PointField& points,
+ const vector& normal
+ );
+
+ //- Return the quality of a triangle projected in a normal direction
+ template
+ static scalar quality
+ (
+ const triFace& triPoints,
+ const PointField& points,
+ const vector& normal
+ );
+
+ //- Return whether or not two edges intersect in the plane defined by a
+ // normal
+ template
+ static bool intersection
+ (
+ const edge& edgePointsA,
+ const edge& edgePointsB,
+ const PointField& points,
+ const vector& normal
+ );
+
+ //- Return the number of intersections of an edge within a polygon
+ template
+ static label nIntersections
+ (
+ const edge& edgePoints,
+ const PointField& points,
+ const vector& normal
+ );
+
+ //- Return the interior angle of the point in a polygon in the plane
+ // defined by a normal
+ template
+ static scalar angle
+ (
+ const label pointi,
+ const PointField& points,
+ const vector& normal
+ );
+
+ //- Return whether a point in a polygon is an "ear"; i.e., that the
+ // triangle formed by adding a diagonal between the two adjacent
+ // points contains no other point
+ template
+ static bool ear
+ (
+ const label pointi,
+ const PointField& points,
+ const vector& normal
+ );
+
+
+ // Private Member Functions
+
+ //- Optimise the triangulation quality by flipping edges
+ template
+ void optimiseTriangulation
+ (
+ const label trii,
+ const PointField& points,
+ const vector& normal,
+ UList& triPoints,
+ UList>& triEdges,
+ UList& edgeTris
+ );
+
+ //- Perform triangulation of a polygon by ear clipping, assuming that
+ // the polygon is simple/non-intersecting
+ void simpleTriangulate
+ (
+ const UList& points,
+ const vector& normal,
+ UList& triPoints,
+ UList>& triEdges,
+ UList& edgeTris,
+ const bool optimal
+ );
+
+ //- Perform triangulation of a polygon, first by adding a spanning
+ // triangle between a point and an edge, then recursing into the
+ // sub-polygons on either side
+ void partitionTriangulate
+ (
+ const UList& points,
+ const vector& normal,
+ const label spanPointi,
+ const label spanEdgei,
+ UList& triPoints,
+ UList>& triEdges,
+ UList& edgeTris,
+ const bool simple,
+ const bool optimal
+ );
+
+ //- Perform triangulation of the given polygon, detecting self
+ // intersections and resolving them by adding spanning triangles prior
+ // to simple triangulation
+ void complexTriangulate
+ (
+ const UList& points,
+ const vector& normal,
+ UList& triPoints,
+ UList>& triEdges,
+ UList& edgeTris,
+ const bool optimal
+ );
+
+ //- Perform triangulation of the given polygon
+ void triangulate
+ (
+ const UList& points,
+ const vector& normal,
+ UList& triPoints,
+ UList>& triEdges,
+ UList& edgeTris,
+ const bool simple = false,
+ const bool optimal = true
+ );
+
+
+public:
+
+ // Static Member Functions
+
+ //- Generate a random polygon for testing purposes
+ static List randomPolygon
+ (
+ Random& rndGen,
+ const label n,
+ const scalar error
+ );
+
+
+
+ // Constructors
+
+ //- Construct null
+ polygonTriangulate();
+
+ //- Disallow default bitwise copy construction
+ polygonTriangulate(const polygonTriangulate&) = delete;
+
+
+ //- Destructor
+ ~polygonTriangulate();
+
+
+ // Member Functions
+
+ // Access
+
+ //- Get the triangles' points
+ inline const UList& triPoints() const;
+
+ //- Get the triangles' renumbered points
+ inline List triPoints
+ (
+ const UList& polyPoints
+ ) const;
+
+ //- Get a triangle's renumbered points
+ inline triFace triPoints
+ (
+ const label trii,
+ const UList& polyPoints
+ ) const;
+
+ //- Get the triangles' edges
+ inline const UList>& triEdges() const;
+
+ //- Get the triangles' renumbered edges
+ inline List> triEdges
+ (
+ const UList& polyEdges
+ ) const;
+
+ //- Get a triangle's renumbered edges
+ inline FixedList triEdges
+ (
+ const label trii,
+ const UList& polyEdges
+ ) const;
+
+
+ //- Perform triangulation of the given polygon
+ const UList& triangulate
+ (
+ const UList& points,
+ const vector& normal,
+ const bool simple = false,
+ const bool optimal = true
+ );
+
+ //- Perform triangulation of the given polygon
+ const UList& triangulate
+ (
+ const UList& points,
+ const bool simple = false,
+ const bool optimal = true
+ );
+
+ //- Perform triangulation of the given polygon
+ inline const UList& triangulate
+ (
+ const UIndirectList& points,
+ const vector& normal,
+ const bool simple = false,
+ const bool optimal = true
+ );
+
+ //- Perform triangulation of the given polygon
+ inline const UList& triangulate
+ (
+ const UIndirectList& points,
+ const bool simple = false,
+ const bool optimal = true
+ );
+
+
+ // Member Operators
+
+ //- Disallow default bitwise assignment
+ void operator=(const polygonTriangulate&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "polygonTriangulateI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/algorithms/polygonTriangulate/polygonTriangulateI.H b/src/OpenFOAM/algorithms/polygonTriangulate/polygonTriangulateI.H
new file mode 100644
index 0000000000..28d9377557
--- /dev/null
+++ b/src/OpenFOAM/algorithms/polygonTriangulate/polygonTriangulateI.H
@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration | Website: https://openfoam.org
+ \\ / A nd | Copyright (C) 2021 OpenFOAM Foundation
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+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 & polyToGlobal
+)
+{
+ return triPoly < polyToGlobal.size() ? polyToGlobal[triPoly] : -1;
+}
+
+template
+inline Type Foam::polygonTriangulate::renumberPolyToGlobal
+(
+ const Type& triPoly,
+ const UList& polyToGlobal
+)
+{
+ Type result;
+ result.resize(triPoly.size());
+ forAll(triPoly, i)
+ {
+ result[i] = renumberPolyToGlobal(triPoly[i], polyToGlobal);
+ }
+ return result;
+}
+
+
+// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
+
+inline const Foam::UList&
+Foam::polygonTriangulate::triPoints() const
+{
+ return triPoints_;
+}
+
+
+inline Foam::List Foam::polygonTriangulate::triPoints
+(
+ const UList& polyPoints
+) const
+{
+ return renumberPolyToGlobal(triPoints_, polyPoints);
+}
+
+
+inline Foam::triFace Foam::polygonTriangulate::triPoints
+(
+ const label trii,
+ const UList& polyPoints
+) const
+{
+ return renumberPolyToGlobal(triPoints_[trii], polyPoints);
+}
+
+
+inline const Foam::UList>&
+Foam::polygonTriangulate::triEdges() const
+{
+ return triEdges_;
+}
+
+
+inline Foam::List>
+Foam::polygonTriangulate::triEdges
+(
+ const UList& polyEdges
+) const
+{
+ return renumberPolyToGlobal(triEdges_, polyEdges);
+}
+
+
+inline Foam::FixedList Foam::polygonTriangulate::triEdges
+(
+ const label trii,
+ const UList& polyEdges
+) const
+{
+ return renumberPolyToGlobal(triEdges_[trii], polyEdges);
+}
+
+
+inline const Foam::UList& Foam::polygonTriangulate::triangulate
+(
+ const UIndirectList& points,
+ const vector& normal,
+ const bool simple,
+ const bool optimal
+)
+{
+ points_ = points;
+
+ return triangulate(points_, normal, simple, optimal);
+}
+
+
+inline const Foam::UList& Foam::polygonTriangulate::triangulate
+(
+ const UIndirectList& points,
+ const bool simple,
+ const bool optimal
+)
+{
+ points_ = points;
+
+ return triangulate(points_, simple, optimal);
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.C b/src/OpenFOAM/meshes/meshShapes/face/face.C
index 314421a7e3..108e1ded9e 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/face.C
+++ b/src/OpenFOAM/meshes/meshShapes/face/face.C
@@ -35,256 +35,6 @@ License
const char* const Foam::face::typeName = "face";
-// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-
-Foam::tmp
-Foam::face::calcEdges(const pointField& points) const
-{
- tmp tedges(new vectorField(size()));
- vectorField& edges = tedges.ref();
-
- forAll(*this, i)
- {
- label ni = fcIndex(i);
-
- point thisPt = points[operator[](i)];
- point nextPt = points[operator[](ni)];
-
- vector vec(nextPt - thisPt);
- vec /= Foam::mag(vec) + vSmall;
-
- edges[i] = vec;
- }
-
- return tedges;
-}
-
-
-Foam::scalar Foam::face::edgeCos
-(
- const vectorField& edges,
- const label index
-) const
-{
- label leftEdgeI = left(index);
- label rightEdgeI = right(index);
-
- // Note negate on left edge to get correct left-pointing edge.
- return -(edges[leftEdgeI] & edges[rightEdgeI]);
-}
-
-
-Foam::label Foam::face::mostConcaveAngle
-(
- const pointField& points,
- const vectorField& edges,
- scalar& maxAngle
-) const
-{
- vector a(area(points));
-
- label index = 0;
- maxAngle = -great;
-
- forAll(edges, i)
- {
- label leftEdgeI = left(i);
- label rightEdgeI = right(i);
-
- vector edgeNormal = edges[rightEdgeI] ^ edges[leftEdgeI];
-
- scalar edgeCos = edges[leftEdgeI] & edges[rightEdgeI];
- scalar edgeAngle = acos(max(-1.0, min(1.0, edgeCos)));
-
- scalar angle;
-
- if ((edgeNormal & a) > 0)
- {
- // Concave angle.
- angle = constant::mathematical::pi + edgeAngle;
- }
- else
- {
- // Convex angle. Note '-' to take into account that rightEdge
- // and leftEdge are head-to-tail connected.
- angle = constant::mathematical::pi - edgeAngle;
- }
-
- if (angle > maxAngle)
- {
- maxAngle = angle;
- index = i;
- }
- }
-
- return index;
-}
-
-
-Foam::label Foam::face::split
-(
- const face::splitMode mode,
- const pointField& points,
- label& triI,
- label& quadI,
- faceList& triFaces,
- faceList& quadFaces
-) const
-{
- label oldIndices = (triI + quadI);
-
- if (size() <= 2)
- {
- FatalErrorInFunction
- << "Serious problem: asked to split a face with < 3 vertices"
- << abort(FatalError);
- }
- if (size() == 3)
- {
- // Triangle. Just copy.
- if (mode == COUNTTRIANGLE || mode == COUNTQUAD)
- {
- triI++;
- }
- else
- {
- triFaces[triI++] = *this;
- }
- }
- else if (size() == 4)
- {
- if (mode == COUNTTRIANGLE)
- {
- triI += 2;
- }
- if (mode == COUNTQUAD)
- {
- quadI++;
- }
- else if (mode == SPLITTRIANGLE)
- {
- // Start at point with largest internal angle.
- const vectorField edges(calcEdges(points));
-
- scalar minAngle;
- label startIndex = mostConcaveAngle(points, edges, minAngle);
-
- label nextIndex = fcIndex(startIndex);
- label splitIndex = fcIndex(nextIndex);
-
- // Create triangles
- face triFace(3);
- triFace[0] = operator[](startIndex);
- triFace[1] = operator[](nextIndex);
- triFace[2] = operator[](splitIndex);
-
- triFaces[triI++] = triFace;
-
- triFace[0] = operator[](splitIndex);
- triFace[1] = operator[](fcIndex(splitIndex));
- triFace[2] = operator[](startIndex);
-
- triFaces[triI++] = triFace;
- }
- else
- {
- quadFaces[quadI++] = *this;
- }
- }
- else
- {
- // General case. Like quad: search for largest internal angle.
-
- const vectorField edges(calcEdges(points));
-
- scalar minAngle = 1;
- label startIndex = mostConcaveAngle(points, edges, minAngle);
-
- scalar bisectAngle = minAngle/2;
- vector rightEdge = edges[right(startIndex)];
-
- //
- // Look for opposite point which as close as possible bisects angle
- //
-
- // split candidate starts two points away.
- label index = fcIndex(fcIndex(startIndex));
-
- label minIndex = index;
- scalar minDiff = constant::mathematical::pi;
-
- for (label i = 0; i < size() - 3; i++)
- {
- vector splitEdge
- (
- points[operator[](index)]
- - points[operator[](startIndex)]
- );
- splitEdge /= Foam::mag(splitEdge) + vSmall;
-
- const scalar splitCos = splitEdge & rightEdge;
- const scalar splitAngle = acos(max(-1.0, min(1.0, splitCos)));
- const scalar angleDiff = fabs(splitAngle - bisectAngle);
-
- if (angleDiff < minDiff)
- {
- minDiff = angleDiff;
- minIndex = index;
- }
-
- // Go to next candidate
- index = fcIndex(index);
- }
-
-
- // Split into two subshapes.
- // face1: startIndex to minIndex
- // face2: minIndex to startIndex
-
- // Get sizes of the two subshapes
- label diff = 0;
- if (minIndex > startIndex)
- {
- diff = minIndex - startIndex;
- }
- else
- {
- // Folded around
- diff = minIndex + size() - startIndex;
- }
-
- label nPoints1 = diff + 1;
- label nPoints2 = size() - diff + 1;
-
- // Collect face1 points
- face face1(nPoints1);
-
- index = startIndex;
- for (label i = 0; i < nPoints1; i++)
- {
- face1[i] = operator[](index);
- index = fcIndex(index);
- }
-
- // Collect face2 points
- face face2(nPoints2);
-
- index = minIndex;
- for (label i = 0; i < nPoints2; i++)
- {
- face2[i] = operator[](index);
- index = fcIndex(index);
- }
-
- // Split faces
- face1.split(mode, points, triI, quadI, triFaces, quadFaces);
- face2.split(mode, points, triI, quadI, triFaces, quadFaces);
- }
-
- return (triI + quadI - oldIndices);
-}
-
-
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::face::face(const triFace& f)
@@ -714,53 +464,6 @@ int Foam::face::edgeDirection(const edge& e) const
}
-Foam::label Foam::face::nTriangles(const pointField&) const
-{
- return nTriangles();
-}
-
-
-Foam::label Foam::face::triangles
-(
- const pointField& points,
- label& triI,
- faceList& triFaces
-) const
-{
- label quadI = 0;
- faceList quadFaces;
-
- return split(SPLITTRIANGLE, points, triI, quadI, triFaces, quadFaces);
-}
-
-
-Foam::label Foam::face::nTrianglesQuads
-(
- const pointField& points,
- label& triI,
- label& quadI
-) const
-{
- faceList triFaces;
- faceList quadFaces;
-
- return split(COUNTQUAD, points, triI, quadI, triFaces, quadFaces);
-}
-
-
-Foam::label Foam::face::trianglesQuads
-(
- const pointField& points,
- label& triI,
- label& quadI,
- faceList& triFaces,
- faceList& quadFaces
-) const
-{
- return split(SPLITQUAD, points, triI, quadI, triFaces, quadFaces);
-}
-
-
Foam::label Foam::longestEdge(const face& f, const pointField& pts)
{
const edgeList& eds = f.edges();
diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.H b/src/OpenFOAM/meshes/meshShapes/face/face.H
index c3c7a4155f..2f41fecf4f 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/face.H
+++ b/src/OpenFOAM/meshes/meshShapes/face/face.H
@@ -77,58 +77,6 @@ class face
:
public labelList
{
- // Private Member Functions
-
- //- Edge to the right of face vertex i
- inline label right(const label i) const;
-
- //- Edge to the left of face vertex i
- inline label left(const label i) const;
-
- //- Construct list of edge vectors for face
- tmp calcEdges
- (
- const pointField& points
- ) const;
-
- //- Cos between neighbouring edges
- scalar edgeCos
- (
- const vectorField& edges,
- const label index
- ) const;
-
- //- Find index of largest internal angle on face
- label mostConcaveAngle
- (
- const pointField& points,
- const vectorField& edges,
- scalar& edgeCos
- ) const;
-
- //- Enumeration listing the modes for split()
- enum splitMode
- {
- COUNTTRIANGLE, // count if split into triangles
- COUNTQUAD, // count if split into triangles&quads
- SPLITTRIANGLE, // split into triangles
- SPLITQUAD // split into triangles&quads
- };
-
- //- Split face into triangles or triangles&quads.
- // Stores results quadFaces[quadI], triFaces[triI]
- // Returns number of new faces created
- label split
- (
- const splitMode mode,
- const pointField& points,
- label& triI,
- label& quadI,
- faceList& triFaces,
- faceList& quadFaces
- ) const;
-
-
public:
//- Return types for classify
@@ -330,56 +278,8 @@ public:
// - -1: reverse (clockwise) on the face
int edgeDirection(const edge&) const;
- // Face splitting utilities
-
- //- Number of triangles after splitting
- inline label nTriangles() const;
-
- //- Number of triangles after splitting
- label nTriangles(const pointField& points) const;
-
- //- Split into triangles using existing points.
- // Result in triFaces[triI..triI+nTri].
- // Splits intelligently to maximise triangle quality.
- // Returns number of faces created.
- label triangles
- (
- const pointField& points,
- label& triI,
- faceList& triFaces
- ) const;
-
- //- Split into triangles using existing points.
- // Append to DynamicList.
- // Returns number of faces created.
- template
- label triangles
- (
- const pointField& points,
- DynamicList& triFaces
- ) const;
-
- //- Number of triangles and quads after splitting
- // Returns the sum of both
- label nTrianglesQuads
- (
- const pointField& points,
- label& nTris,
- label& nQuads
- ) const;
-
- //- Split into triangles and quads.
- // Results in triFaces (starting at triI) and quadFaces
- // (starting at quadI).
- // Returns number of new faces created.
- label trianglesQuads
- (
- const pointField& points,
- label& triI,
- label& quadI,
- faceList& triFaces,
- faceList& quadFaces
- ) const;
+ //- Size of the face's triangulation
+ inline label nTriangles() const;
//- Compare faces
// 0: different
diff --git a/src/OpenFOAM/meshes/meshShapes/face/faceI.H b/src/OpenFOAM/meshes/meshShapes/face/faceI.H
index 7b8cf873a4..ce1cedc9b7 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/faceI.H
+++ b/src/OpenFOAM/meshes/meshShapes/face/faceI.H
@@ -23,20 +23,6 @@ License
\*---------------------------------------------------------------------------*/
-// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-
-inline Foam::label Foam::face::right(const label i) const
-{
- return i;
-}
-
-
-inline Foam::label Foam::face::left(const label i) const
-{
- return rcIndex(i);
-}
-
-
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
inline Foam::face::face()
diff --git a/src/OpenFOAM/meshes/meshShapes/face/faceTemplates.C b/src/OpenFOAM/meshes/meshShapes/face/faceTemplates.C
index c088675b1b..fda721585e 100644
--- a/src/OpenFOAM/meshes/meshShapes/face/faceTemplates.C
+++ b/src/OpenFOAM/meshes/meshShapes/face/faceTemplates.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -136,24 +136,6 @@ Foam::vector Foam::face::area(const PointField& ps)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-template
-Foam::label Foam::face::triangles
-(
- const pointField& points,
- DynamicList& triFaces
-) const
-{
- label triI = triFaces.size();
- label quadI = 0;
- faceList quadFaces;
-
- // adjust the addressable size (and allocate space if needed)
- triFaces.setSize(triI + nTriangles());
-
- return split(SPLITTRIANGLE, points, triI, quadI, triFaces, quadFaces);
-}
-
-
template
Type Foam::face::average
(
diff --git a/src/lagrangian/parcel/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C b/src/lagrangian/parcel/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
index 348695cc7e..32abf356d5 100644
--- a/src/lagrangian/parcel/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
+++ b/src/lagrangian/parcel/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2012-2020 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2012-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -131,7 +131,6 @@ void Foam::ParticleCollector::initPolygons
label pointOffset = 0;
points_.setSize(nPoints);
faces_.setSize(polygons.size());
- faceTris_.setSize(polygons.size());
area_.setSize(polygons.size());
forAll(faces_, facei)
{
@@ -139,11 +138,6 @@ void Foam::ParticleCollector::initPolygons
face f(identity(polyPoints.size()) + pointOffset);
UIndirectList(points_, f) = polyPoints;
area_[facei] = f.mag(points_);
-
- DynamicList tris;
- f.triangles(points_, tris);
- faceTris_[facei].transfer(tris);
-
faces_[facei].transfer(f);
pointOffset += polyPoints.size();
@@ -527,7 +521,6 @@ Foam::ParticleCollector::ParticleCollector
removeCollected_(this->coeffDict().lookup("removeCollected")),
points_(),
faces_(),
- faceTris_(),
nSector_(0),
radius_(),
coordSys_("coordSys", vector::zero, axesRotation(sphericalTensor::I)),
@@ -613,7 +606,6 @@ Foam::ParticleCollector::ParticleCollector
removeCollected_(pc.removeCollected_),
points_(pc.points_),
faces_(pc.faces_),
- faceTris_(pc.faceTris_),
nSector_(pc.nSector_),
radius_(pc.radius_),
coordSys_(pc.coordSys_),
diff --git a/src/lagrangian/parcel/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.H b/src/lagrangian/parcel/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.H
index e265c24aaa..17841d0616 100644
--- a/src/lagrangian/parcel/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.H
+++ b/src/lagrangian/parcel/submodels/CloudFunctionObjects/ParticleCollector/ParticleCollector.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2012-2020 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2012-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -144,11 +144,6 @@ private:
List faces_;
- // Polygon collector
-
- //- Triangulation of faces
- List> faceTris_;
-
// Concentric circles collector
//- Number of sectors per circle
diff --git a/src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchInjection/patchInjectionBase.C b/src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchInjection/patchInjectionBase.C
index c4131bcf4b..dff1942bb5 100644
--- a/src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchInjection/patchInjectionBase.C
+++ b/src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchInjection/patchInjectionBase.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2013-2020 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2013-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -30,6 +30,7 @@ License
#include "triPointRef.H"
#include "volFields.H"
#include "polyMeshTetDecomposition.H"
+#include "polygonTriangulate.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@@ -94,24 +95,25 @@ void Foam::patchInjectionBase::updateMesh(const polyMesh& mesh)
// Triangulate the patch faces and create addressing
DynamicList triToFace(2*patch.size());
DynamicList triMagSf(2*patch.size());
- DynamicList triFace(2*patch.size());
- DynamicList tris(5);
+ DynamicList triFace(2*patch.size());
// Set zero value at the start of the tri area list
triMagSf.append(0.0);
+ // Create a triangulation engine
+ polygonTriangulate triEngine;
+
forAll(patch, facei)
{
const face& f = patch[facei];
- tris.clear();
- f.triangles(points, tris);
+ triEngine.triangulate(UIndirectList(points, f));
- forAll(tris, i)
+ forAll(triEngine.triPoints(), i)
{
triToFace.append(facei);
- triFace.append(tris[i]);
- triMagSf.append(tris[i].mag(points));
+ triFace.append(triEngine.triPoints(i, f));
+ triMagSf.append(triEngine.triPoints(i, f).mag(points));
}
}
diff --git a/src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchInjection/patchInjectionBase.H b/src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchInjection/patchInjectionBase.H
index 3e17f46025..ace6002363 100644
--- a/src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchInjection/patchInjectionBase.H
+++ b/src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchInjection/patchInjectionBase.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2013-2020 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2013-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -45,6 +45,7 @@ SourceFiles
#include "scalarList.H"
#include "vectorList.H"
#include "faceList.H"
+#include "triFaceList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -82,7 +83,7 @@ protected:
labelList cellOwners_;
//- Decomposed patch faces as a list of triangles
- faceList triFace_;
+ triFaceList triFace_;
//- Addressing from per triangle to patch face
labelList triToFace_;
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
index de267db135..cfc2890cc8 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -147,7 +147,7 @@ Foam::tmp Foam::AMIInterpolation::patchMagSf
const pointField& patchPoints = patch.localPoints();
- faceList patchFaceTris;
+ triFaceList patchFaceTris;
forAll(result, patchFacei)
{
@@ -162,12 +162,7 @@ Foam::tmp Foam::AMIInterpolation::patchMagSf
forAll(patchFaceTris, i)
{
result[patchFacei] +=
- triPointRef
- (
- patchPoints[patchFaceTris[i][0]],
- patchPoints[patchFaceTris[i][1]],
- patchPoints[patchFaceTris[i][2]]
- ).mag();
+ patchFaceTris[i].tri(patchPoints).mag();
}
}
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/sweptFaceAreaWeightAMI/sweptFaceAreaWeightAMI.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/sweptFaceAreaWeightAMI/sweptFaceAreaWeightAMI.C
index 870c66da33..1524715140 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/sweptFaceAreaWeightAMI/sweptFaceAreaWeightAMI.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIMethod/sweptFaceAreaWeightAMI/sweptFaceAreaWeightAMI.C
@@ -395,7 +395,7 @@ Foam::scalar Foam::sweptFaceAreaWeightAMI::interArea
// Triangulate the faces
const faceAreaIntersect::triangulationMode triMode = this->triMode_;
- faceList srcFaceTris, tgtFaceTris;
+ triFaceList srcFaceTris, tgtFaceTris;
faceAreaIntersect::triangulate(srcFace, srcPoints, triMode, srcFaceTris);
faceAreaIntersect::triangulate(tgtFace, tgtPoints, triMode, tgtFaceTris);
@@ -403,7 +403,7 @@ Foam::scalar Foam::sweptFaceAreaWeightAMI::interArea
scalar areaMag = Zero;
// Loop the target triangles
- forAllConstIter(faceList, tgtFaceTris, tgtIter)
+ forAllConstIter(triFaceList, tgtFaceTris, tgtIter)
{
const FixedList
tgtTri =
@@ -414,7 +414,7 @@ Foam::scalar Foam::sweptFaceAreaWeightAMI::interArea
};
// Loop the source triangles
- forAllConstIter(faceList, srcFaceTris, srcIter)
+ forAllConstIter(triFaceList, srcFaceTris, srcIter)
{
FixedList
srcTri =
diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
index c6ef121ac7..bf2da47524 100644
--- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
+++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "faceAreaIntersect.H"
+#include "polygonTriangulate.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -330,7 +331,7 @@ void Foam::faceAreaIntersect::triangulate
const face& f,
const pointField& points,
const triangulationMode& triMode,
- faceList& faceTris
+ triFaceList& faceTris
)
{
faceTris.resize(f.nTriangles());
@@ -341,7 +342,6 @@ void Foam::faceAreaIntersect::triangulate
{
for (label i = 0; i < f.nTriangles(); ++ i)
{
- faceTris[i] = face(3);
faceTris[i][0] = f[0];
faceTris[i][1] = f[i + 1];
faceTris[i][2] = f[i + 2];
@@ -351,18 +351,9 @@ void Foam::faceAreaIntersect::triangulate
}
case tmMesh:
{
- const label nFaceTris = f.nTriangles();
-
- label nFaceTris1 = 0;
- const label nFaceTris2 = f.triangles(points, nFaceTris1, faceTris);
-
- if (nFaceTris != nFaceTris1 || nFaceTris != nFaceTris2)
- {
- FatalErrorInFunction
- << "The numbers of reported triangles in the face do not "
- << "match that generated by the triangulation"
- << exit(FatalError);
- }
+ polygonTriangulate triEngine;
+ triEngine.triangulate(UIndirectList(points, f));
+ faceTris = triEngine.triPoints(f);
}
}
}
@@ -377,7 +368,7 @@ Foam::scalar Foam::faceAreaIntersect::calc
)
{
// split faces into triangles
- faceList trisA, trisB;
+ triFaceList trisA, trisB;
triangulate(faceA, pointsA_, triMode, trisA);
triangulate(faceB, pointsB_, triMode, trisB);
diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
index e654eceadb..4f41182e96 100644
--- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
+++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -42,6 +42,7 @@ SourceFiles
#include "plane.H"
#include "face.H"
#include "NamedEnum.H"
+#include "triFaceList.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -92,7 +93,7 @@ private:
inline triPoints getTriPoints
(
const pointField& points,
- const face& f,
+ const triFace& f,
const bool reverse
) const;
@@ -162,7 +163,7 @@ public:
const face& f,
const pointField& points,
const triangulationMode& triMode,
- faceList& faceTris
+ triFaceList& faceTris
);
//- Return area of intersection of faceA with faceB
diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
index 0dcf2d0597..569fccb6c5 100644
--- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
+++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersectI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -44,7 +44,7 @@ inline void Foam::faceAreaIntersect::setTriPoints
inline Foam::faceAreaIntersect::triPoints Foam::faceAreaIntersect::getTriPoints
(
const pointField& points,
- const face& f,
+ const triFace& f,
const bool reverse
) const
{
diff --git a/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C b/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C
index 13ed53727a..51b5e68f2c 100644
--- a/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C
+++ b/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.C
@@ -26,7 +26,7 @@ License
#include "intersectedSurface.H"
#include "surfaceIntersection.H"
#include "faceList.H"
-#include "faceTriangulation.H"
+#include "polygonTriangulate.H"
#include "treeBoundBox.H"
#include "OFstream.H"
#include "error.H"
@@ -1185,6 +1185,9 @@ Foam::intersectedSurface::intersectedSurface
// Start in newTris for decomposed face.
labelList startTriI(surf.size(), 0);
+ // Create a triangulation engine
+ polygonTriangulate triEngine;
+
forAll(surf, facei)
{
startTriI[facei] = newTris.size();
@@ -1209,56 +1212,21 @@ Foam::intersectedSurface::intersectedSurface
{
const face& newF = newFaces[newFacei];
-// {
-// fileName fName
-// (
-// "face_"
-// + Foam::name(facei)
-// + "_subFace_"
-// + Foam::name(newFacei)
-// + ".obj"
-// );
-// Pout<< "Writing original face:" << facei << " subFace:"
-// << newFacei << " to " << fName << endl;
-//
-// OFstream str(fName);
-//
-// forAll(newF, fp)
-// {
-// meshTools::writeOBJ(str, eSurf.points()[newF[fp]]);
-// }
-// str << 'l';
-// forAll(newF, fp)
-// {
-// str << ' ' << fp+1;
-// }
-// str<< " 1" << nl;
-// }
-
-
const vector& n = surf.faceNormals()[facei];
const label region = surf[facei].region();
- faceTriangulation tris(eSurf.points(), newF, n);
+ triEngine.triangulate
+ (
+ UIndirectList(eSurf.points(), newF),
+ n
+ );
- forAll(tris, triI)
+ forAll(triEngine.triPoints(), trii)
{
- const triFace& t = tris[triI];
-
- forAll(t, i)
- {
- if (t[i] < 0 || t[i] >= eSurf.points().size())
- {
- FatalErrorInFunction
- << "Face triangulation of face " << facei
- << " uses points outside range 0.."
- << eSurf.points().size()-1 << endl
- << "Triangulation:"
- << tris << abort(FatalError);
- }
- }
-
- newTris.append(labelledTri(t[0], t[1], t[2], region));
+ newTris.append
+ (
+ labelledTri(triEngine.triPoints(trii, newF), region)
+ );
}
}
}
diff --git a/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.H b/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.H
index b6fd2fffac..b2cc10dc43 100644
--- a/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.H
+++ b/src/meshTools/triSurface/booleanOps/intersectedSurface/intersectedSurface.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -39,8 +39,7 @@ Description
nSurfaceEdges)
- construct face-edge addressing for above edges
- for each face do a right-handed walk to reconstruct faces (splitFace)
- - retriangulate resulting faces (might be non-convex so use
- faceTriangulation which does proper bisection)
+ - retriangulate resulting faces
The resulting surface will have the points from the surface first
in the point list (0 .. nSurfacePoints-1)
diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
index 3afaa3038f..5b18cfbd31 100644
--- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
+++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -31,6 +31,7 @@ License
#include "polyMesh.H"
#include "plane.H"
#include "geompack.H"
+#include "polygonTriangulate.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -2342,6 +2343,8 @@ Foam::triSurface Foam::triSurfaceTools::triangulate
mesh.nFaces() - mesh.nInternalFaces()
);
+ polygonTriangulate triEngine;
+
label newPatchi = 0;
forAllConstIter(labelHashSet, includePatches, iter)
@@ -2356,17 +2359,14 @@ Foam::triSurface Foam::triSurfaceTools::triangulate
{
const face& f = patch[patchFacei];
- faceList triFaces(f.nTriangles(points));
+ triEngine.triangulate(UIndirectList(points, f));
- label nTri = 0;
-
- f.triangles(points, nTri, triFaces);
-
- forAll(triFaces, triFacei)
+ forAll(triEngine.triPoints(), triFacei)
{
- const face& f = triFaces[triFacei];
-
- triangles.append(labelledTri(f[0], f[1], f[2], newPatchi));
+ triangles.append
+ (
+ labelledTri(triEngine.triPoints(triFacei, f), newPatchi)
+ );
nTriTotal++;
}
@@ -2429,6 +2429,8 @@ Foam::triSurface Foam::triSurfaceTools::triangulate
mesh.nFaces() - mesh.nInternalFaces()
);
+ polygonTriangulate triEngine;
+
label newPatchi = 0;
forAllConstIter(labelHashSet, includePatches, iter)
@@ -2445,17 +2447,14 @@ Foam::triSurface Foam::triSurfaceTools::triangulate
if (bBox.containsAny(points, f))
{
- faceList triFaces(f.nTriangles(points));
+ triEngine.triangulate(UIndirectList(points, f));
- label nTri = 0;
-
- f.triangles(points, nTri, triFaces);
-
- forAll(triFaces, triFacei)
+ forAll(triEngine.triPoints(), triFacei)
{
- const face& f = triFaces[triFacei];
-
- triangles.append(labelledTri(f[0], f[1], f[2], newPatchi));
+ triangles.append
+ (
+ labelledTri(triEngine.triPoints(triFacei, f), newPatchi)
+ );
nTriTotal++;
}
diff --git a/src/sampling/cuttingPlane/cuttingPlane.C b/src/sampling/cuttingPlane/cuttingPlane.C
index 9dd3f3791d..0d0bd028a7 100644
--- a/src/sampling/cuttingPlane/cuttingPlane.C
+++ b/src/sampling/cuttingPlane/cuttingPlane.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -27,6 +27,7 @@ License
#include "primitiveMesh.H"
#include "linePointRef.H"
#include "meshTools.H"
+#include "polygonTriangulate.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -280,6 +281,9 @@ void Foam::cuttingPlane::walkCellCuts
// scratch space for calculating the face vertices
DynamicList faceVerts(10);
+ // Create a triangulation engine
+ polygonTriangulate triEngine;
+
forAll(cutCells_, i)
{
label celli = cutCells_[i];
@@ -331,9 +335,13 @@ void Foam::cuttingPlane::walkCellCuts
// the cut faces are usually quite ugly, so optionally triangulate
if (triangulate)
{
- label nTri = f.triangles(cutPoints, dynCutFaces);
- while (nTri--)
+ triEngine.triangulate
+ (
+ UIndirectList(cutPoints, f)
+ );
+ forAll(triEngine.triPoints(), i)
{
+ dynCutFaces.append(triEngine.triPoints(i, f));
dynCutCells.append(celli);
}
}
diff --git a/src/sampling/sampledSurface/thresholdCellFaces/thresholdCellFaces.C b/src/sampling/sampledSurface/thresholdCellFaces/thresholdCellFaces.C
index c17e34c9a4..35ee56454b 100644
--- a/src/sampling/sampledSurface/thresholdCellFaces/thresholdCellFaces.C
+++ b/src/sampling/sampledSurface/thresholdCellFaces/thresholdCellFaces.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -25,6 +25,7 @@ License
#include "thresholdCellFaces.H"
#include "emptyPolyPatch.H"
+#include "polygonTriangulate.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -86,6 +87,8 @@ void Foam::thresholdCellFaces::calculate
labelList oldToNewPoints(origPoints.size(), -1);
+ // Create a triangulation engine
+ polygonTriangulate triEngine;
// internal faces only
for (label facei = 0; facei < mesh_.nInternalFaces(); ++facei)
@@ -149,9 +152,13 @@ void Foam::thresholdCellFaces::calculate
if (triangulate)
{
- label count = surfFace.triangles(origPoints, surfFaces);
- while (count-- > 0)
+ triEngine.triangulate
+ (
+ UIndirectList(origPoints, surfFace)
+ );
+ forAll(triEngine.triPoints(), i)
{
+ surfFaces.append(triEngine.triPoints(i, surfFace));
surfCells.append(cellId);
}
}
@@ -207,9 +214,13 @@ void Foam::thresholdCellFaces::calculate
if (triangulate)
{
- label count = f.triangles(origPoints, surfFaces);
- while (count-- > 0)
+ triEngine.triangulate
+ (
+ UIndirectList(origPoints, f)
+ );
+ forAll(triEngine.triPoints(), i)
{
+ surfFaces.append(triEngine.triPoints(i, f));
surfCells.append(cellId);
}
}
diff --git a/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriter.C b/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriter.C
index 6455a938be..2efc345062 100644
--- a/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriter.C
+++ b/src/sampling/sampledSurface/writers/nastran/nastranSurfaceWriter.C
@@ -24,6 +24,7 @@ License
\*---------------------------------------------------------------------------*/
#include "nastranSurfaceWriter.H"
+#include "polygonTriangulate.H"
#include "makeSurfaceWriterMethods.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -283,6 +284,8 @@ void Foam::nastranSurfaceWriter::writeGeometry
label nFace = 1;
+ polygonTriangulate triEngine;
+
forAll(faces, facei)
{
const face& f = faces[facei];
@@ -300,14 +303,12 @@ void Foam::nastranSurfaceWriter::writeGeometry
else
{
// decompose poly face into tris
- label nTri = 0;
- faceList triFaces;
- f.triangles(points, nTri, triFaces);
-
- forAll(triFaces, triI)
+ triEngine.triangulate(UIndirectList(points, f));
+ forAll(triEngine.triPoints(), triI)
{
- writeFace("CTRIA3", triFaces[triI], nFace, os);
- decomposedFaces[facei].append(triFaces[triI]);
+ const face f(triEngine.triPoints(triI, f));
+ writeFace("CTRIA3", f, nFace, os);
+ decomposedFaces[facei].append(f);
}
}
}
diff --git a/src/surfMesh/MeshedSurface/MeshedSurface.C b/src/surfMesh/MeshedSurface/MeshedSurface.C
index 23010eb468..c39be166d7 100644
--- a/src/surfMesh/MeshedSurface/MeshedSurface.C
+++ b/src/surfMesh/MeshedSurface/MeshedSurface.C
@@ -33,6 +33,7 @@ License
#include "polyMesh.H"
#include "surfMesh.H"
#include "primitivePatch.H"
+#include "polygonTriangulate.H"
#include "addToRunTimeSelectionTable.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -864,22 +865,18 @@ Foam::label Foam::MeshedSurface::triangulate
else
{
// triangulate with points
- List tmpTri(maxTri);
+ polygonTriangulate triEngine;
label newFacei = 0;
forAll(faceLst, facei)
{
- // 'face' not ''
const face& f = faceLst[facei];
- label nTmp = 0;
- f.triangles(this->points(), nTmp, tmpTri);
- for (label triI = 0; triI < nTmp; triI++)
+ triEngine.triangulate(UIndirectList(this->points(), f));
+
+ forAll(triEngine.triPoints(), triI)
{
- newFaces[newFacei] = Face
- (
- static_cast(tmpTri[triI])
- );
+ newFaces[newFacei] = triEngine.triPoints(triI, f);
faceMap[newFacei] = facei;
newFacei++;
}
diff --git a/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C b/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C
index 7240692273..707844ede3 100644
--- a/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C
+++ b/src/surfMesh/surfaceFormats/starcd/STARCDsurfaceFormat.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2019 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -25,6 +25,7 @@ License
#include "STARCDsurfaceFormat.H"
#include "ListOps.H"
+#include "polygonTriangulate.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@@ -128,6 +129,9 @@ bool Foam::fileFormats::STARCDsurfaceFormat::read
readHeader(is, "PROSTAR_CELL");
+ // Create a triangulation engine
+ polygonTriangulate triEngine;
+
DynamicList dynFaces;
DynamicList dynZones;
DynamicList dynNames;
@@ -201,24 +205,14 @@ bool Foam::fileFormats::STARCDsurfaceFormat::read
SubList vertices(vertexLabels, vertexLabels.size());
if (mustTriangulate && nLabels > 3)
{
- face f(vertices);
+ triEngine.triangulate
+ (
+ UIndirectList(this->points(), vertices)
+ );
- faceList triFaces(f.nTriangles());
- label nTri = 0;
- f.triangles(this->points(), nTri, triFaces);
-
- forAll(triFaces, facei)
+ forAll(triEngine.triPoints(), trii)
{
- // a triangular face, but not yet a triFace
- dynFaces.append
- (
- triFace
- (
- static_cast(triFaces[facei])
- )
- );
- dynZones.append(zoneI);
- dynSizes[zoneI]++;
+ dynFaces.append(triEngine.triPoints(trii, vertices));
}
}
else
diff --git a/src/triSurface/Make/files b/src/triSurface/Make/files
index d1183dd516..c296dcca02 100644
--- a/src/triSurface/Make/files
+++ b/src/triSurface/Make/files
@@ -1,7 +1,6 @@
triSurfaceTools = triSurface/triSurfaceTools
geometricSurfacePatch = triSurface/geometricSurfacePatch
-faceTriangulation/faceTriangulation.C
meshTriangulation/meshTriangulation.C
triSurface/triSurface.C
diff --git a/src/triSurface/faceTriangulation/faceTriangulation.C b/src/triSurface/faceTriangulation/faceTriangulation.C
deleted file mode 100644
index f25df15c13..0000000000
--- a/src/triSurface/faceTriangulation/faceTriangulation.C
+++ /dev/null
@@ -1,656 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
- \\/ M anipulation |
--------------------------------------------------------------------------------
-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 Foam::faceTriangulation::calcEdges
-(
- const face& f,
- const pointField& points
-)
-{
- tmp tedges(new vectorField(f.size()));
- vectorField& edges = tedges.ref();
-
- forAll(f, i)
- {
- point thisPt = points[f[i]];
- point nextPt = points[f[f.fcIndex(i)]];
-
- vector vec(nextPt - thisPt);
- vec /= mag(vec) + vSmall;
-
- edges[i] = vec;
- }
-
- return tedges;
-}
-
-
-// Calculates half angle components of angle from e0 to e1
-void Foam::faceTriangulation::calcHalfAngle
-(
- const vector& normal,
- const vector& e0,
- const vector& e1,
- scalar& cosHalfAngle,
- scalar& sinHalfAngle
-)
-{
- // truncate cos to +-1 to prevent negative numbers
- scalar cos = max(-1, min(1, e0 & e1));
-
- scalar sin = (e0 ^ e1) & normal;
-
- if (sin < -rootVSmall)
- {
- // 3rd or 4th quadrant
- cosHalfAngle = - Foam::sqrt(0.5*(1 + cos));
- sinHalfAngle = Foam::sqrt(0.5*(1 - cos));
- }
- else
- {
- // 1st or 2nd quadrant
- cosHalfAngle = Foam::sqrt(0.5*(1 + cos));
- sinHalfAngle = Foam::sqrt(0.5*(1 - cos));
- }
-}
-
-
-// Calculate intersection point between edge p1-p2 and ray (in 2D).
-// Return true and intersection point if intersection between p1 and p2.
-Foam::pointHit Foam::faceTriangulation::rayEdgeIntersect
-(
- const vector& normal,
- const point& rayOrigin,
- const vector& rayDir,
- const point& p1,
- const point& p2,
- scalar& posOnEdge
-)
-{
- // Start off from miss
- pointHit result(p1);
-
- // Construct plane normal to rayDir and intersect
- const vector y = normal ^ rayDir;
-
- posOnEdge = plane(rayOrigin, y).normalIntersect(p1, (p2-p1));
-
- // Check intersection to left of p1 or right of p2
- if ((posOnEdge < 0) || (posOnEdge > 1))
- {
- // Miss
- }
- else
- {
- // Check intersection behind rayOrigin
- point intersectPt = p1 + posOnEdge * (p2 - p1);
-
- if (((intersectPt - rayOrigin) & rayDir) < 0)
- {
- // Miss
- }
- else
- {
- // Hit
- result.setHit();
- result.setPoint(intersectPt);
- result.setDistance(mag(intersectPt - rayOrigin));
- }
- }
- return result;
-}
-
-
-// Return true if triangle given its three points (anticlockwise ordered)
-// contains point
-bool Foam::faceTriangulation::triangleContainsPoint
-(
- const vector& n,
- const point& p0,
- const point& p1,
- const point& p2,
- const point& pt
-)
-{
- scalar area01Pt = triPointRef(p0, p1, pt).area() & n;
- scalar area12Pt = triPointRef(p1, p2, pt).area() & n;
- scalar area20Pt = triPointRef(p2, p0, pt).area() & n;
-
- if ((area01Pt > 0) && (area12Pt > 0) && (area20Pt > 0))
- {
- return true;
- }
- else if ((area01Pt < 0) && (area12Pt < 0) && (area20Pt < 0))
- {
- FatalErrorInFunction << abort(FatalError);
- return false;
- }
- else
- {
- return false;
- }
-}
-
-
-// Starting from startIndex find diagonal. Return in index1, index2.
-// Index1 always startIndex except when convex polygon
-void Foam::faceTriangulation::findDiagonal
-(
- const pointField& points,
- const face& f,
- const vectorField& edges,
- const vector& normal,
- const label startIndex,
- label& index1,
- label& index2
-)
-{
- const point& startPt = points[f[startIndex]];
-
- // Calculate angle at startIndex
- const vector& rightE = edges[right(f.size(), startIndex)];
- const vector leftE = -edges[left(f.size(), startIndex)];
-
- // Construct ray which bisects angle
- scalar cosHalfAngle = great;
- scalar sinHalfAngle = great;
- calcHalfAngle(normal, rightE, leftE, cosHalfAngle, sinHalfAngle);
- vector rayDir
- (
- cosHalfAngle*rightE
- + sinHalfAngle*(normal ^ rightE)
- );
- // rayDir should be normalised already but is not due to rounding errors
- // so normalise.
- rayDir /= mag(rayDir) + vSmall;
-
-
- //
- // Check all edges (apart from rightE and leftE) for nearest intersection
- //
-
- label faceVertI = f.fcIndex(startIndex);
-
- pointHit minInter(false, Zero, great, true);
- label minIndex = -1;
- scalar minPosOnEdge = great;
-
- for (label i = 0; i < f.size() - 2; i++)
- {
- scalar posOnEdge;
- pointHit inter =
- rayEdgeIntersect
- (
- normal,
- startPt,
- rayDir,
- points[f[faceVertI]],
- points[f[f.fcIndex(faceVertI)]],
- posOnEdge
- );
-
- if (inter.hit() && inter.distance() < minInter.distance())
- {
- minInter = inter;
- minIndex = faceVertI;
- minPosOnEdge = posOnEdge;
- }
-
- faceVertI = f.fcIndex(faceVertI);
- }
-
-
- if (minIndex == -1)
- {
- // WarningInFunction
- // << "Could not find intersection starting from " << f[startIndex]
- // << " for face " << f << endl;
-
- index1 = -1;
- index2 = -1;
- return;
- }
-
- const label leftIndex = minIndex;
- const label rightIndex = f.fcIndex(minIndex);
-
- // Now ray intersects edge from leftIndex to rightIndex.
- // Check for intersection being one of the edge points. Make sure never
- // to return two consecutive points.
-
- if (mag(minPosOnEdge) < edgeRelTol && f.fcIndex(startIndex) != leftIndex)
- {
- index1 = startIndex;
- index2 = leftIndex;
-
- return;
- }
- if
- (
- mag(minPosOnEdge - 1) < edgeRelTol
- && f.fcIndex(rightIndex) != startIndex
- )
- {
- index1 = startIndex;
- index2 = rightIndex;
-
- return;
- }
-
- // Select visible vertex that minimises
- // angle to bisection. Visibility checking by checking if inside triangle
- // formed by startIndex, leftIndex, rightIndex
-
- const point& leftPt = points[f[leftIndex]];
- const point& rightPt = points[f[rightIndex]];
-
- minIndex = -1;
- scalar maxCos = -great;
-
- // all vertices except for startIndex and ones to left and right of it.
- faceVertI = f.fcIndex(f.fcIndex(startIndex));
- for (label i = 0; i < f.size() - 3; i++)
- {
- const point& pt = points[f[faceVertI]];
-
- if
- (
- (faceVertI == leftIndex)
- || (faceVertI == rightIndex)
- || (triangleContainsPoint(normal, startPt, leftPt, rightPt, pt))
- )
- {
- // pt inside triangle (so perhaps visible)
- // Select based on minimal angle (so guaranteed visible).
- vector edgePt0 = pt - startPt;
- edgePt0 /= mag(edgePt0);
-
- scalar cos = rayDir & edgePt0;
- if (cos > maxCos)
- {
- maxCos = cos;
- minIndex = faceVertI;
- }
- }
- faceVertI = f.fcIndex(faceVertI);
- }
-
- if (minIndex == -1)
- {
- // no vertex found. Return startIndex and one of the intersected edge
- // endpoints.
- index1 = startIndex;
-
- if (f.fcIndex(startIndex) != leftIndex)
- {
- index2 = leftIndex;
- }
- else
- {
- index2 = rightIndex;
- }
-
- return;
- }
-
- index1 = startIndex;
- index2 = minIndex;
-}
-
-
-// Find label of vertex to start splitting from. Is:
-// 1] flattest concave angle
-// 2] flattest convex angle if no concave angles.
-Foam::label Foam::faceTriangulation::findStart
-(
- const face& f,
- const vectorField& edges,
- const vector& normal
-)
-{
- const label size = f.size();
-
- scalar minCos = great;
- label minIndex = -1;
-
- forAll(f, fp)
- {
- const vector& rightEdge = edges[right(size, fp)];
- const vector leftEdge = -edges[left(size, fp)];
-
- if (((rightEdge ^ leftEdge) & normal) < rootVSmall)
- {
- scalar cos = rightEdge & leftEdge;
- if (cos < minCos)
- {
- minCos = cos;
- minIndex = fp;
- }
- }
- }
-
- if (minIndex == -1)
- {
- // No concave angle found. Get flattest convex angle
- minCos = great;
-
- forAll(f, fp)
- {
- const vector& rightEdge = edges[right(size, fp)];
- const vector leftEdge = -edges[left(size, fp)];
-
- scalar cos = rightEdge & leftEdge;
- if (cos < minCos)
- {
- minCos = cos;
- minIndex = fp;
- }
- }
- }
-
- return minIndex;
-}
-
-
-// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
-
-// Split face f into triangles. Handles all simple (convex & concave)
-// polygons.
-bool Foam::faceTriangulation::split
-(
- const bool fallBack,
- const pointField& points,
- const face& f,
- const vector& normal,
- label& triI
-)
-{
- const label size = f.size();
-
- if (size <= 2)
- {
- WarningInFunction
- << "Illegal face:" << f
- << " with points " << UIndirectList(points, f)()
- << endl;
-
- return false;
- }
- else if (size == 3)
- {
- // Triangle. Just copy.
- triFace& tri = operator[](triI++);
- tri[0] = f[0];
- tri[1] = f[1];
- tri[2] = f[2];
-
- return true;
- }
- else
- {
- // General case. Start splitting for -flattest concave angle
- // -or flattest convex angle if no concave angles.
-
- tmp tedges(calcEdges(f, points));
- const vectorField& edges = tedges();
-
- label startIndex = findStart(f, edges, normal);
-
- // Find diagonal to split face across
- label index1 = -1;
- label index2 = -1;
-
- forAll(f, iter)
- {
- findDiagonal
- (
- points,
- f,
- edges,
- normal,
- startIndex,
- index1,
- index2
- );
-
- if (index1 != -1 && index2 != -1)
- {
- // Found correct diagonal
- break;
- }
-
- // Try splitting from next startingIndex.
- startIndex = f.fcIndex(startIndex);
- }
-
- if (index1 == -1 || index2 == -1)
- {
- if (fallBack)
- {
- // Do naive triangulation. Find smallest angle to start
- // triangulating from.
- label maxIndex = -1;
- scalar maxCos = -great;
-
- forAll(f, fp)
- {
- const vector& rightEdge = edges[right(size, fp)];
- const vector leftEdge = -edges[left(size, fp)];
-
- scalar cos = rightEdge & leftEdge;
- if (cos > maxCos)
- {
- maxCos = cos;
- maxIndex = fp;
- }
- }
-
- WarningInFunction
- << "Cannot find valid diagonal on face " << f
- << " with points " << UIndirectList(points, f)()
- << nl
- << "Returning naive triangulation starting from "
- << f[maxIndex] << " which might not be correct for a"
- << " concave or warped face" << endl;
-
-
- label fp = f.fcIndex(maxIndex);
-
- for (label i = 0; i < size-2; i++)
- {
- label nextFp = f.fcIndex(fp);
-
- triFace& tri = operator[](triI++);
- tri[0] = f[maxIndex];
- tri[1] = f[fp];
- tri[2] = f[nextFp];
-
- fp = nextFp;
- }
-
- return true;
- }
- else
- {
- WarningInFunction
- << "Cannot find valid diagonal on face " << f
- << " with points " << UIndirectList(points, f)()
- << nl
- << "Returning empty triFaceList" << endl;
-
- return false;
- }
- }
-
-
- // Split into two subshapes.
- // face1: index1 to index2
- // face2: index2 to index1
-
- // Get sizes of the two subshapes
- label diff = 0;
- if (index2 > index1)
- {
- diff = index2 - index1;
- }
- else
- {
- // folded round
- diff = index2 + size - index1;
- }
-
- label nPoints1 = diff + 1;
- label nPoints2 = size - diff + 1;
-
- if (nPoints1 == size || nPoints2 == size)
- {
- FatalErrorInFunction
- << "Illegal split of face:" << f
- << " with points " << UIndirectList(points, f)()
- << " at indices " << index1 << " and " << index2
- << abort(FatalError);
- }
-
-
- // Collect face1 points
- face face1(nPoints1);
-
- label faceVertI = index1;
- for (int i = 0; i < nPoints1; i++)
- {
- face1[i] = f[faceVertI];
- faceVertI = f.fcIndex(faceVertI);
- }
-
- // Collect face2 points
- face face2(nPoints2);
-
- faceVertI = index2;
- for (int i = 0; i < nPoints2; i++)
- {
- face2[i] = f[faceVertI];
- faceVertI = f.fcIndex(faceVertI);
- }
-
- // Decompose the split faces
- // Pout<< "Split face:" << f << " into " << face1 << " and " << face2
- // << endl;
- // string oldPrefix(Pout.prefix());
- // Pout.prefix() = " " + oldPrefix;
-
- bool splitOk =
- split(fallBack, points, face1, normal, triI)
- && split(fallBack, points, face2, normal, triI);
-
- // Pout.prefix() = oldPrefix;
-
- return splitOk;
- }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
-
-Foam::faceTriangulation::faceTriangulation()
-:
- triFaceList()
-{}
-
-
-Foam::faceTriangulation::faceTriangulation
-(
- const pointField& points,
- const face& f,
- const bool fallBack
-)
-:
- triFaceList(f.size()-2)
-{
- const vector avgNormal = f.normal(points);
-
- label triI = 0;
-
- bool valid = split(fallBack, points, f, avgNormal, triI);
-
- if (!valid)
- {
- setSize(0);
- }
-}
-
-
-Foam::faceTriangulation::faceTriangulation
-(
- const pointField& points,
- const face& f,
- const vector& n,
- const bool fallBack
-)
-:
- triFaceList(f.size()-2)
-{
- label triI = 0;
-
- bool valid = split(fallBack, points, f, n, triI);
-
- if (!valid)
- {
- setSize(0);
- }
-}
-
-
-Foam::faceTriangulation::faceTriangulation(Istream& is)
-:
- triFaceList(is)
-{}
-
-
-// ************************************************************************* //
diff --git a/src/triSurface/faceTriangulation/faceTriangulation.H b/src/triSurface/faceTriangulation/faceTriangulation.H
deleted file mode 100644
index 09bfd1618e..0000000000
--- a/src/triSurface/faceTriangulation/faceTriangulation.H
+++ /dev/null
@@ -1,202 +0,0 @@
-/*---------------------------------------------------------------------------*\
- ========= |
- \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
- \\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
- \\/ M anipulation |
--------------------------------------------------------------------------------
-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 calcEdges(const face&, const pointField&);
-
- //- Calculates half angle components of angle from e0 to e1
- // in plane given by normal.
- static void calcHalfAngle
- (
- const vector& normal,
- const vector& e0,
- const vector& e1,
- scalar& cosHalfAngle,
- scalar& sinHalfAngle
- );
-
- //- Calculate intersection point between edge p1-p2 and ray (in 2D).
- // Return true and intersection point if intersection between p1 and p2.
- static pointHit rayEdgeIntersect
- (
- const vector& normal,
- const point& rayOrigin,
- const vector& rayDir,
- const point& p1,
- const point& p2,
- scalar& posOnEdge
- );
-
- // Return true if triangle given its three points
- // (anticlockwise ordered) contains point
- static bool triangleContainsPoint
- (
- const vector& n,
- const point& p0,
- const point& p1,
- const point& p2,
- const point& pt
- );
-
- //- Starting from startIndex find diagonal. Return in index1, index2.
- // Index1 always startIndex except when convex polygon
- static void findDiagonal
- (
- const pointField& points,
- const face& f,
- const vectorField& edges,
- const vector& normal,
- const label startIndex,
- label& index1,
- label& index2
- );
-
- //- Find label of vertex to start splitting from. This will be the
- // vertex with edge angle:
- // 1] flattest concave angle
- // 2] flattest convex angle if no concave angles.
- static label findStart
- (
- const face& f,
- const vectorField& edges,
- const vector& normal
- );
-
-
- // Private Member Functions
-
- //- Split face f into triangles. Handles all simple (convex & concave)
- // polygons. Returns false if could not produce valid split.
- bool split
- (
- const bool fallBack,
- const pointField& points,
- const face& f,
- const vector& normal,
- label& triI
- );
-
-public:
-
- // Constructors
-
- //- Construct null
- faceTriangulation();
-
- //- Construct from face and points. Decomposition based on average
- // normal. After construction *this is size 0 or holds the triangles.
- // If fallBack and triangulation fails does naive triangulation
- // and never returns 0 size.
- faceTriangulation
- (
- const pointField& points,
- const face& f,
- const bool fallBack = false
- );
-
- //- Construct from face and points and user supplied (unit) normal.
- // After construction *this is size 0 or holds the triangles.
- // If fallBack and triangulation fails does naive triangulation
- // and never returns 0 size.
- faceTriangulation
- (
- const pointField& points,
- const face& f,
- const vector& n,
- const bool fallBack = false
- );
-
- //- Construct from Istream
- faceTriangulation(Istream&);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //
diff --git a/src/triSurface/meshTriangulation/meshTriangulation.C b/src/triSurface/meshTriangulation/meshTriangulation.C
index d9a55de017..cd4c3da8d5 100644
--- a/src/triSurface/meshTriangulation/meshTriangulation.C
+++ b/src/triSurface/meshTriangulation/meshTriangulation.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -25,7 +25,7 @@ License
#include "meshTriangulation.H"
#include "polyMesh.H"
-#include "faceTriangulation.H"
+#include "polygonTriangulate.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@@ -209,7 +209,7 @@ Foam::meshTriangulation::meshTriangulation
{
if (faceIsCut[facei])
{
- nTotTri += faces[facei].nTriangles(points);
+ nTotTri += faces[facei].nTriangles();
}
}
}
@@ -352,38 +352,28 @@ Foam::meshTriangulation::meshTriangulation
// Triangulation using existing vertices
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ // Create a triangulation engine
+ polygonTriangulate triEngine;
+
// Triangulate internal faces
forAll(faceIsCut, facei)
{
if (faceIsCut[facei] && isInternalFace(mesh, includedCell, facei))
{
- // Face was internal to the mesh and will be 'internal' to
- // the surface.
+ triEngine.triangulate
+ (
+ UIndirectList(points, faces[facei])
+ );
- // Triangulate face. Fall back to naive triangulation if failed.
- faceTriangulation faceTris(points, faces[facei], true);
-
- if (faceTris.empty())
- {
- WarningInFunction
- << "Could not find triangulation for face " << facei
- << " vertices " << faces[facei] << " coords "
- << IndirectList(points, faces[facei])() << endl;
- }
- else
- {
- // Copy triangles. Make them internalFacesPatch
- insertTriangles
- (
- faceTris,
- facei,
- internalFacesPatch,
- false, // no reverse
-
- triangles,
- triI
- );
- }
+ insertTriangles
+ (
+ triEngine.triPoints(faces[facei]),
+ facei,
+ internalFacesPatch,
+ false, // no reverse
+ triangles,
+ triI
+ );
}
}
nInternalFaces_ = triI;
@@ -423,30 +413,20 @@ Foam::meshTriangulation::meshTriangulation
reverse = false;
}
- // Triangulate face
- faceTriangulation faceTris(points, faces[facei], true);
+ triEngine.triangulate
+ (
+ UIndirectList(points, faces[facei])
+ );
- if (faceTris.empty())
- {
- WarningInFunction
- << "Could not find triangulation for face " << facei
- << " vertices " << faces[facei] << " coords "
- << IndirectList(points, faces[facei])() << endl;
- }
- else
- {
- // Copy triangles. Optionally reverse them
- insertTriangles
- (
- faceTris,
- facei,
- patchi,
- reverse, // whether to reverse
-
- triangles,
- triI
- );
- }
+ insertTriangles
+ (
+ triEngine.triPoints(faces[facei]),
+ facei,
+ patchi,
+ reverse, // whether to reverse
+ triangles,
+ triI
+ );
}
}
}
diff --git a/src/triSurface/triSurface/interfaces/OFF/readOFF.C b/src/triSurface/triSurface/interfaces/OFF/readOFF.C
index 6814664047..b4cb756db4 100644
--- a/src/triSurface/triSurface/interfaces/OFF/readOFF.C
+++ b/src/triSurface/triSurface/interfaces/OFF/readOFF.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | Website: https://openfoam.org
- \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2021 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -27,6 +27,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "triSurface.H"
+#include "polygonTriangulate.H"
#include "IFstream.H"
#include "IStringStream.H"
@@ -76,6 +77,9 @@ bool Foam::triSurface::readOFF(const fileName& OFFfileName)
points[pointi] = point(x, y, z);
}
+ // Create a triangulation engine
+ polygonTriangulate triEngine;
+
// Read faces & triangulate them,
DynamicList tris(nElems);
@@ -107,17 +111,11 @@ bool Foam::triSurface::readOFF(const fileName& OFFfileName)
}
else
{
- faceList triFaces(f.nTriangles(points));
+ triEngine.triangulate(UIndirectList(points, f));
- label nTri = 0;
-
- f.triangles(points, nTri, triFaces);
-
- forAll(triFaces, triFacei)
+ forAll(triEngine.triPoints(), trii)
{
- const face& f = triFaces[triFacei];
-
- tris.append(labelledTri(f[0], f[1], f[2], 0));
+ tris.append(labelledTri(triEngine.triPoints(trii, f), 0));
}
}
}