polygonTriangulate: Added robust polygon triangulation algorithm
The new algorithm provides robust quality triangulations of non-convex polygons. It also produces a best attempt for polygons that are badly warped or self intersecting by minimising the area in which the local normal is in the opposite direction to the overal polygon normal. It is memory efficient when applied to multiple polygons as it maintains and reuses its workspace. This algorithm replaces implementations in the face and faceTriangulation classes, which have been removed. Faces can no longer be decomposed into mixtures of tris and quadrilaterals. Polygonal faces with more than 4 sides are now decomposed into triangles in foamToVTK and in paraFoam.
This commit is contained in:
Will Bainbridge
3
applications/test/polygonTriangulate/Make/files
Normal file
3
applications/test/polygonTriangulate/Make/files
Normal file
@ -0,0 +1,3 @@
|
||||
Test-polygonTriangulate.C
|
||||
|
||||
EXE = $(FOAM_USER_APPBIN)/Test-polygonTriangulate
|
||||
5
applications/test/polygonTriangulate/Make/options
Normal file
5
applications/test/polygonTriangulate/Make/options
Normal file
@ -0,0 +1,5 @@
|
||||
EXE_INC = \
|
||||
-I$(LIB_SRC)/fileFormats/lnInclude
|
||||
|
||||
EXE_LIBS = \
|
||||
-lfileFormats
|
||||
121
applications/test/polygonTriangulate/Test-polygonTriangulate.C
Normal file
121
applications/test/polygonTriangulate/Test-polygonTriangulate.C
Normal file
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "argList.H"
|
||||
#include "clock.H"
|
||||
#include "OBJstream.H"
|
||||
#include "polygonTriangulate.H"
|
||||
|
||||
using namespace Foam;
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
argList::validArgs.append("number of edges");
|
||||
argList::addOption
|
||||
(
|
||||
"error",
|
||||
"value",
|
||||
"polygon error - default is 0"
|
||||
);
|
||||
argList::addOption
|
||||
(
|
||||
"seed",
|
||||
"value",
|
||||
"random number generator seed - default is clock time"
|
||||
);
|
||||
argList::addBoolOption
|
||||
(
|
||||
"simple",
|
||||
"assume polygon as no self intersections"
|
||||
);
|
||||
argList::addBoolOption
|
||||
(
|
||||
"nonOptimal",
|
||||
"do not optimise the triangulation quality"
|
||||
);
|
||||
|
||||
Foam::argList args(argc, argv);
|
||||
Info<< nl;
|
||||
|
||||
// Get the size of the polygon
|
||||
const label n = args.argRead<label>(1);
|
||||
|
||||
// Initialise the random number generator
|
||||
label seed;
|
||||
if (args.optionFound("seed"))
|
||||
{
|
||||
seed = args.optionRead<label>("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<point> polygon =
|
||||
polygonTriangulate::randomPolygon
|
||||
(
|
||||
rndGen,
|
||||
n,
|
||||
args.optionLookupOrDefault<scalar>("error", 0)
|
||||
);
|
||||
|
||||
// Write the polygon
|
||||
{
|
||||
OBJstream os(args[0] + "_polygon.obj");
|
||||
os.write(face(identity(polygon.size())), polygon, false);
|
||||
}
|
||||
|
||||
// Triangulate the polygon
|
||||
const List<triFace> 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;
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -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<point>(mesh.points(), f)
|
||||
);
|
||||
forAll(triEngine.triPoints(), trii)
|
||||
{
|
||||
triFcs[trii] = triEngine.triPoints(trii, f);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
quadFcs[0] = f;
|
||||
}
|
||||
|
||||
forAll(quadFcs, quadI)
|
||||
{
|
||||
|
||||
@ -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<point>(mesh.points(), f)
|
||||
);
|
||||
forAll(triEngine.triPoints(), trii)
|
||||
{
|
||||
triFcs[trii] = triEngine.triPoints(trii, f);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
quadFcs[0] = f;
|
||||
}
|
||||
|
||||
forAll(quadFcs, quadI)
|
||||
{
|
||||
|
||||
@ -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<point>(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
|
||||
(
|
||||
|
||||
@ -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
|
||||
|
||||
919
src/OpenFOAM/algorithms/polygonTriangulate/polygonTriangulate.C
Normal file
919
src/OpenFOAM/algorithms/polygonTriangulate/polygonTriangulate.C
Normal file
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "polygonTriangulate.H"
|
||||
#include "tensor2D.H"
|
||||
|
||||
// * * * * * * * * * * * Private Static Member Functions * * * * * * * * * * //
|
||||
|
||||
template<class PointField>
|
||||
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<class PointField>
|
||||
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<class PointField>
|
||||
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<class PointField>
|
||||
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<class PointField>
|
||||
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<class PointField>
|
||||
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::point> Foam::polygonTriangulate::randomPolygon
|
||||
(
|
||||
Random& rndGen,
|
||||
const label n,
|
||||
const scalar error
|
||||
)
|
||||
{
|
||||
// Get random points on a unit disk
|
||||
List<point> points(n);
|
||||
forAll(points, pointi)
|
||||
{
|
||||
const scalar theta =
|
||||
2*constant::mathematical::pi*rndGen.sample01<scalar>();
|
||||
const scalar r = sqrt(rndGen.sample01<scalar>());
|
||||
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<label> 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<scalar>();
|
||||
const scalar r = sqrt(rndGen.sample01<scalar>());
|
||||
points[pointi] =
|
||||
(1 - error)*points[pointi]
|
||||
+ error*point(r*cos(theta), r*sin(theta), rndGen.sample01<scalar>());
|
||||
}
|
||||
|
||||
// Reorder and return
|
||||
return List<point>(points, pointis);
|
||||
}
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
template<class PointField>
|
||||
void Foam::polygonTriangulate::optimiseTriangulation
|
||||
(
|
||||
const label trii,
|
||||
const PointField& points,
|
||||
const vector& normal,
|
||||
UList<triFace>& triPoints,
|
||||
UList<FixedList<label, 3>>& triEdges,
|
||||
UList<labelPair>& 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<point>& allPoints,
|
||||
const vector& normal,
|
||||
UList<triFace>& triPoints,
|
||||
UList<FixedList<label, 3>>& triEdges,
|
||||
UList<labelPair>& edgeTris,
|
||||
const bool optimal
|
||||
)
|
||||
{
|
||||
// Get the polygon size
|
||||
const label n = allPoints.size();
|
||||
|
||||
// Clear the triangulation
|
||||
triPoints = triFace(-1, -1, -1);
|
||||
triEdges = FixedList<label, 3>({-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<point> 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<label, 3>
|
||||
({
|
||||
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<point>& points,
|
||||
const vector& normal,
|
||||
const label spanPointi,
|
||||
const label spanEdgei,
|
||||
UList<triFace>& triPoints,
|
||||
UList<FixedList<label, 3>>& triEdges,
|
||||
UList<labelPair>& 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<label, 3>({eiA, eiB, eiE});
|
||||
|
||||
// Rotate the point list so that it can be passed to sub-triangulations
|
||||
// without indirect addressing
|
||||
inplaceRotateList(const_cast<UList<point>&>(points), - piA);
|
||||
|
||||
// Triangulate the sub-polygon connected to edge A
|
||||
if (nA >= 3)
|
||||
{
|
||||
// Subset the polygon and triangulate
|
||||
SubList<triFace> triPointsA(triPoints, nA - 2, 1);
|
||||
SubList<FixedList<label, 3>> triEdgesA(triEdges, nA - 2, 1);
|
||||
SubList<labelPair> edgeTrisA(edgeTris, 2*nA - 3);
|
||||
triangulate
|
||||
(
|
||||
SubList<point>(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<triFace> triPointsB(triPoints, nB - 2, nA - 1);
|
||||
SubList<FixedList<label, 3>> triEdgesB(triEdges, nB - 2, nA - 1);
|
||||
SubList<labelPair> edgeTrisB(edgeTris, 2*nB - 3, 2*nA - 3);
|
||||
triangulate
|
||||
(
|
||||
SubList<point>(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<UList<point>&>(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<labelPair> 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<labelPair> 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<labelPair> 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<point>& points,
|
||||
const vector& normal,
|
||||
UList<triFace>& triPoints,
|
||||
UList<FixedList<label, 3>>& triEdges,
|
||||
UList<labelPair>& 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<labelPair, 8> 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<point>& points,
|
||||
const vector& normal,
|
||||
UList<triFace>& triPoints,
|
||||
UList<FixedList<label, 3>>& triEdges,
|
||||
UList<labelPair>& 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::triFace>& Foam::polygonTriangulate::triangulate
|
||||
(
|
||||
const UList<point>& 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::triFace>& Foam::polygonTriangulate::triangulate
|
||||
(
|
||||
const UList<point>& points,
|
||||
const bool simple,
|
||||
const bool optimal
|
||||
)
|
||||
{
|
||||
return
|
||||
triangulate
|
||||
(
|
||||
points,
|
||||
normalised(face::area(points)),
|
||||
simple,
|
||||
optimal
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
351
src/OpenFOAM/algorithms/polygonTriangulate/polygonTriangulate.H
Normal file
351
src/OpenFOAM/algorithms/polygonTriangulate/polygonTriangulate.H
Normal file
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
Class
|
||||
Foam::polygonTriangulate
|
||||
|
||||
Description
|
||||
Triangulation of three-dimensional polygons
|
||||
|
||||
SourceFiles
|
||||
polygonTriangulate.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef polygonTriangulate_H
|
||||
#define polygonTriangulate_H
|
||||
|
||||
#include "DynamicList.H"
|
||||
#include "HashSet.H"
|
||||
#include "labelPair.H"
|
||||
#include "triFace.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class polygonTriangulate Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class polygonTriangulate
|
||||
{
|
||||
// Private Data
|
||||
|
||||
//- Indirect addressing into point field to form the current as yet
|
||||
// un-triangulated polygon
|
||||
DynamicList<label> pointis_;
|
||||
|
||||
//- Edge indices of the un-triangulated polygon
|
||||
DynamicList<label> edges_;
|
||||
|
||||
//- The interior angles of the vertices of the un-triangulated polygon
|
||||
DynamicList<scalar> 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<bool> ear_;
|
||||
|
||||
//- Set of tri-faces used to mark previously considered configurations
|
||||
// in optimisation and prevent unnecessary repetition
|
||||
HashSet<triFace, Hash<triFace>> triPointsSet_;
|
||||
|
||||
//- Point array. Used if the input points cannot be SubList-d.
|
||||
DynamicList<point> points_;
|
||||
|
||||
//- Tri-points
|
||||
DynamicList<triFace> triPoints_;
|
||||
|
||||
//- Tri-edge associations
|
||||
DynamicList<FixedList<label, 3>> triEdges_;
|
||||
|
||||
//- Edge-tri associations
|
||||
DynamicList<labelPair> edgeTris_;
|
||||
|
||||
|
||||
// Private Static Member Functions
|
||||
|
||||
// Renumber a polygon label to global indexing
|
||||
static inline label renumberPolyToGlobal
|
||||
(
|
||||
const label triPoly,
|
||||
const UList<label>& polyToGlobal
|
||||
);
|
||||
|
||||
// Renumber a container of polygon labels to global indexing
|
||||
template<class Type>
|
||||
static inline Type renumberPolyToGlobal
|
||||
(
|
||||
const Type& triPoly,
|
||||
const UList<label>& polyToGlobal
|
||||
);
|
||||
|
||||
//- Return the area of a triangle projected in a normal direction
|
||||
template<class PointField>
|
||||
static scalar area
|
||||
(
|
||||
const triFace& triPoints,
|
||||
const PointField& points,
|
||||
const vector& normal
|
||||
);
|
||||
|
||||
//- Return the quality of a triangle projected in a normal direction
|
||||
template<class PointField>
|
||||
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<class PointField>
|
||||
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<class PointField>
|
||||
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<class PointField>
|
||||
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<class PointField>
|
||||
static bool ear
|
||||
(
|
||||
const label pointi,
|
||||
const PointField& points,
|
||||
const vector& normal
|
||||
);
|
||||
|
||||
|
||||
// Private Member Functions
|
||||
|
||||
//- Optimise the triangulation quality by flipping edges
|
||||
template<class PointField>
|
||||
void optimiseTriangulation
|
||||
(
|
||||
const label trii,
|
||||
const PointField& points,
|
||||
const vector& normal,
|
||||
UList<triFace>& triPoints,
|
||||
UList<FixedList<label, 3>>& triEdges,
|
||||
UList<labelPair>& edgeTris
|
||||
);
|
||||
|
||||
//- Perform triangulation of a polygon by ear clipping, assuming that
|
||||
// the polygon is simple/non-intersecting
|
||||
void simpleTriangulate
|
||||
(
|
||||
const UList<point>& points,
|
||||
const vector& normal,
|
||||
UList<triFace>& triPoints,
|
||||
UList<FixedList<label, 3>>& triEdges,
|
||||
UList<labelPair>& 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<point>& points,
|
||||
const vector& normal,
|
||||
const label spanPointi,
|
||||
const label spanEdgei,
|
||||
UList<triFace>& triPoints,
|
||||
UList<FixedList<label, 3>>& triEdges,
|
||||
UList<labelPair>& 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<point>& points,
|
||||
const vector& normal,
|
||||
UList<triFace>& triPoints,
|
||||
UList<FixedList<label, 3>>& triEdges,
|
||||
UList<labelPair>& edgeTris,
|
||||
const bool optimal
|
||||
);
|
||||
|
||||
//- Perform triangulation of the given polygon
|
||||
void triangulate
|
||||
(
|
||||
const UList<point>& points,
|
||||
const vector& normal,
|
||||
UList<triFace>& triPoints,
|
||||
UList<FixedList<label, 3>>& triEdges,
|
||||
UList<labelPair>& edgeTris,
|
||||
const bool simple = false,
|
||||
const bool optimal = true
|
||||
);
|
||||
|
||||
|
||||
public:
|
||||
|
||||
// Static Member Functions
|
||||
|
||||
//- Generate a random polygon for testing purposes
|
||||
static List<point> 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<triFace>& triPoints() const;
|
||||
|
||||
//- Get the triangles' renumbered points
|
||||
inline List<triFace> triPoints
|
||||
(
|
||||
const UList<label>& polyPoints
|
||||
) const;
|
||||
|
||||
//- Get a triangle's renumbered points
|
||||
inline triFace triPoints
|
||||
(
|
||||
const label trii,
|
||||
const UList<label>& polyPoints
|
||||
) const;
|
||||
|
||||
//- Get the triangles' edges
|
||||
inline const UList<FixedList<label, 3>>& triEdges() const;
|
||||
|
||||
//- Get the triangles' renumbered edges
|
||||
inline List<FixedList<label, 3>> triEdges
|
||||
(
|
||||
const UList<label>& polyEdges
|
||||
) const;
|
||||
|
||||
//- Get a triangle's renumbered edges
|
||||
inline FixedList<label, 3> triEdges
|
||||
(
|
||||
const label trii,
|
||||
const UList<label>& polyEdges
|
||||
) const;
|
||||
|
||||
|
||||
//- Perform triangulation of the given polygon
|
||||
const UList<triFace>& triangulate
|
||||
(
|
||||
const UList<point>& points,
|
||||
const vector& normal,
|
||||
const bool simple = false,
|
||||
const bool optimal = true
|
||||
);
|
||||
|
||||
//- Perform triangulation of the given polygon
|
||||
const UList<triFace>& triangulate
|
||||
(
|
||||
const UList<point>& points,
|
||||
const bool simple = false,
|
||||
const bool optimal = true
|
||||
);
|
||||
|
||||
//- Perform triangulation of the given polygon
|
||||
inline const UList<triFace>& triangulate
|
||||
(
|
||||
const UIndirectList<point>& points,
|
||||
const vector& normal,
|
||||
const bool simple = false,
|
||||
const bool optimal = true
|
||||
);
|
||||
|
||||
//- Perform triangulation of the given polygon
|
||||
inline const UList<triFace>& triangulate
|
||||
(
|
||||
const UIndirectList<point>& 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
|
||||
|
||||
// ************************************************************************* //
|
||||
138
src/OpenFOAM/algorithms/polygonTriangulate/polygonTriangulateI.H
Normal file
138
src/OpenFOAM/algorithms/polygonTriangulate/polygonTriangulateI.H
Normal file
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "polygonTriangulate.H"
|
||||
|
||||
// * * * * * * * * * * * Private Static Member Functions * * * * * * * * * * //
|
||||
|
||||
inline Foam::label Foam::polygonTriangulate::renumberPolyToGlobal
|
||||
(
|
||||
const label triPoly,
|
||||
const UList<label>& polyToGlobal
|
||||
)
|
||||
{
|
||||
return triPoly < polyToGlobal.size() ? polyToGlobal[triPoly] : -1;
|
||||
}
|
||||
|
||||
template<class Type>
|
||||
inline Type Foam::polygonTriangulate::renumberPolyToGlobal
|
||||
(
|
||||
const Type& triPoly,
|
||||
const UList<label>& 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::triFace>&
|
||||
Foam::polygonTriangulate::triPoints() const
|
||||
{
|
||||
return triPoints_;
|
||||
}
|
||||
|
||||
|
||||
inline Foam::List<Foam::triFace> Foam::polygonTriangulate::triPoints
|
||||
(
|
||||
const UList<label>& polyPoints
|
||||
) const
|
||||
{
|
||||
return renumberPolyToGlobal(triPoints_, polyPoints);
|
||||
}
|
||||
|
||||
|
||||
inline Foam::triFace Foam::polygonTriangulate::triPoints
|
||||
(
|
||||
const label trii,
|
||||
const UList<label>& polyPoints
|
||||
) const
|
||||
{
|
||||
return renumberPolyToGlobal(triPoints_[trii], polyPoints);
|
||||
}
|
||||
|
||||
|
||||
inline const Foam::UList<Foam::FixedList<Foam::label, 3>>&
|
||||
Foam::polygonTriangulate::triEdges() const
|
||||
{
|
||||
return triEdges_;
|
||||
}
|
||||
|
||||
|
||||
inline Foam::List<Foam::FixedList<Foam::label, 3>>
|
||||
Foam::polygonTriangulate::triEdges
|
||||
(
|
||||
const UList<label>& polyEdges
|
||||
) const
|
||||
{
|
||||
return renumberPolyToGlobal(triEdges_, polyEdges);
|
||||
}
|
||||
|
||||
|
||||
inline Foam::FixedList<Foam::label, 3> Foam::polygonTriangulate::triEdges
|
||||
(
|
||||
const label trii,
|
||||
const UList<label>& polyEdges
|
||||
) const
|
||||
{
|
||||
return renumberPolyToGlobal(triEdges_[trii], polyEdges);
|
||||
}
|
||||
|
||||
|
||||
inline const Foam::UList<Foam::triFace>& Foam::polygonTriangulate::triangulate
|
||||
(
|
||||
const UIndirectList<point>& points,
|
||||
const vector& normal,
|
||||
const bool simple,
|
||||
const bool optimal
|
||||
)
|
||||
{
|
||||
points_ = points;
|
||||
|
||||
return triangulate(points_, normal, simple, optimal);
|
||||
}
|
||||
|
||||
|
||||
inline const Foam::UList<Foam::triFace>& Foam::polygonTriangulate::triangulate
|
||||
(
|
||||
const UIndirectList<point>& points,
|
||||
const bool simple,
|
||||
const bool optimal
|
||||
)
|
||||
{
|
||||
points_ = points;
|
||||
|
||||
return triangulate(points_, simple, optimal);
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -35,256 +35,6 @@ License
|
||||
const char* const Foam::face::typeName = "face";
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
|
||||
|
||||
Foam::tmp<Foam::vectorField>
|
||||
Foam::face::calcEdges(const pointField& points) const
|
||||
{
|
||||
tmp<vectorField> 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();
|
||||
|
||||
@ -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<vectorField> 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<unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
|
||||
label triangles
|
||||
(
|
||||
const pointField& points,
|
||||
DynamicList<face, SizeInc, SizeMult, SizeDiv>& 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
|
||||
|
||||
@ -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()
|
||||
|
||||
@ -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<unsigned SizeInc, unsigned SizeMult, unsigned SizeDiv>
|
||||
Foam::label Foam::face::triangles
|
||||
(
|
||||
const pointField& points,
|
||||
DynamicList<face, SizeInc, SizeMult, SizeDiv>& 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<class Type>
|
||||
Type Foam::face::average
|
||||
(
|
||||
|
||||
@ -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<CloudType>::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<CloudType>::initPolygons
|
||||
face f(identity(polyPoints.size()) + pointOffset);
|
||||
UIndirectList<point>(points_, f) = polyPoints;
|
||||
area_[facei] = f.mag(points_);
|
||||
|
||||
DynamicList<face> tris;
|
||||
f.triangles(points_, tris);
|
||||
faceTris_[facei].transfer(tris);
|
||||
|
||||
faces_[facei].transfer(f);
|
||||
|
||||
pointOffset += polyPoints.size();
|
||||
@ -527,7 +521,6 @@ Foam::ParticleCollector<CloudType>::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<CloudType>::ParticleCollector
|
||||
removeCollected_(pc.removeCollected_),
|
||||
points_(pc.points_),
|
||||
faces_(pc.faces_),
|
||||
faceTris_(pc.faceTris_),
|
||||
nSector_(pc.nSector_),
|
||||
radius_(pc.radius_),
|
||||
coordSys_(pc.coordSys_),
|
||||
|
||||
@ -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<face> faces_;
|
||||
|
||||
|
||||
// Polygon collector
|
||||
|
||||
//- Triangulation of faces
|
||||
List<List<face>> faceTris_;
|
||||
|
||||
// Concentric circles collector
|
||||
|
||||
//- Number of sectors per circle
|
||||
|
||||
@ -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<label> triToFace(2*patch.size());
|
||||
DynamicList<scalar> triMagSf(2*patch.size());
|
||||
DynamicList<face> triFace(2*patch.size());
|
||||
DynamicList<face> tris(5);
|
||||
DynamicList<triFace> 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<point>(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));
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -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_;
|
||||
|
||||
@ -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::scalarField> Foam::AMIInterpolation::patchMagSf
|
||||
|
||||
const pointField& patchPoints = patch.localPoints();
|
||||
|
||||
faceList patchFaceTris;
|
||||
triFaceList patchFaceTris;
|
||||
|
||||
forAll(result, patchFacei)
|
||||
{
|
||||
@ -162,12 +162,7 @@ Foam::tmp<Foam::scalarField> 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();
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -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<point, 3>
|
||||
tgtTri =
|
||||
@ -414,7 +414,7 @@ Foam::scalar Foam::sweptFaceAreaWeightAMI::interArea
|
||||
};
|
||||
|
||||
// Loop the source triangles
|
||||
forAllConstIter(faceList, srcFaceTris, srcIter)
|
||||
forAllConstIter(triFaceList, srcFaceTris, srcIter)
|
||||
{
|
||||
FixedList<point, 4>
|
||||
srcTri =
|
||||
|
||||
@ -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<point>(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);
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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
|
||||
{
|
||||
|
||||
@ -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<point>(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)
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -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)
|
||||
|
||||
@ -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<point>(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<point>(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++;
|
||||
}
|
||||
|
||||
@ -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<label> 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<point>(cutPoints, f)
|
||||
);
|
||||
forAll(triEngine.triPoints(), i)
|
||||
{
|
||||
dynCutFaces.append(triEngine.triPoints(i, f));
|
||||
dynCutCells.append(celli);
|
||||
}
|
||||
}
|
||||
|
||||
@ -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<point>(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<point>(origPoints, f)
|
||||
);
|
||||
forAll(triEngine.triPoints(), i)
|
||||
{
|
||||
surfFaces.append(triEngine.triPoints(i, f));
|
||||
surfCells.append(cellId);
|
||||
}
|
||||
}
|
||||
|
||||
@ -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<point>(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);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -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<Face>::triangulate
|
||||
else
|
||||
{
|
||||
// triangulate with points
|
||||
List<face> tmpTri(maxTri);
|
||||
polygonTriangulate triEngine;
|
||||
|
||||
label newFacei = 0;
|
||||
forAll(faceLst, facei)
|
||||
{
|
||||
// 'face' not '<Face>'
|
||||
const face& f = faceLst[facei];
|
||||
|
||||
label nTmp = 0;
|
||||
f.triangles(this->points(), nTmp, tmpTri);
|
||||
for (label triI = 0; triI < nTmp; triI++)
|
||||
triEngine.triangulate(UIndirectList<point>(this->points(), f));
|
||||
|
||||
forAll(triEngine.triPoints(), triI)
|
||||
{
|
||||
newFaces[newFacei] = Face
|
||||
(
|
||||
static_cast<labelUList&>(tmpTri[triI])
|
||||
);
|
||||
newFaces[newFacei] = triEngine.triPoints(triI, f);
|
||||
faceMap[newFacei] = facei;
|
||||
newFacei++;
|
||||
}
|
||||
|
||||
@ -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<Face>::read
|
||||
|
||||
readHeader(is, "PROSTAR_CELL");
|
||||
|
||||
// Create a triangulation engine
|
||||
polygonTriangulate triEngine;
|
||||
|
||||
DynamicList<Face> dynFaces;
|
||||
DynamicList<label> dynZones;
|
||||
DynamicList<word> dynNames;
|
||||
@ -201,24 +205,14 @@ bool Foam::fileFormats::STARCDsurfaceFormat<Face>::read
|
||||
SubList<label> vertices(vertexLabels, vertexLabels.size());
|
||||
if (mustTriangulate && nLabels > 3)
|
||||
{
|
||||
face f(vertices);
|
||||
triEngine.triangulate
|
||||
(
|
||||
UIndirectList<point>(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<labelUList&>(triFaces[facei])
|
||||
)
|
||||
);
|
||||
dynZones.append(zoneI);
|
||||
dynSizes[zoneI]++;
|
||||
dynFaces.append(triEngine.triPoints(trii, vertices));
|
||||
}
|
||||
}
|
||||
else
|
||||
|
||||
@ -1,7 +1,6 @@
|
||||
triSurfaceTools = triSurface/triSurfaceTools
|
||||
geometricSurfacePatch = triSurface/geometricSurfacePatch
|
||||
|
||||
faceTriangulation/faceTriangulation.C
|
||||
meshTriangulation/meshTriangulation.C
|
||||
|
||||
triSurface/triSurface.C
|
||||
|
||||
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#include "faceTriangulation.H"
|
||||
#include "plane.H"
|
||||
|
||||
|
||||
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
|
||||
|
||||
const Foam::scalar Foam::faceTriangulation::edgeRelTol = 1e-6;
|
||||
|
||||
|
||||
// Edge to the right of face vertex i
|
||||
Foam::label Foam::faceTriangulation::right(const label, label i)
|
||||
{
|
||||
return i;
|
||||
}
|
||||
|
||||
|
||||
// Edge to the left of face vertex i
|
||||
Foam::label Foam::faceTriangulation::left(const label size, label i)
|
||||
{
|
||||
return i ? i-1 : size-1;
|
||||
}
|
||||
|
||||
|
||||
// Calculate (normalised) edge vectors.
|
||||
// edges[i] gives edge between point i+1 and i.
|
||||
Foam::tmp<Foam::vectorField> Foam::faceTriangulation::calcEdges
|
||||
(
|
||||
const face& f,
|
||||
const pointField& points
|
||||
)
|
||||
{
|
||||
tmp<vectorField> 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<point>(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<vectorField> 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<point>(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<point>(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<point>(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)
|
||||
{}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -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 <http://www.gnu.org/licenses/>.
|
||||
|
||||
Class
|
||||
Foam::faceTriangulation
|
||||
|
||||
Description
|
||||
Triangulation of faces. Handles concave polygons as well
|
||||
(inefficiently)
|
||||
|
||||
Works by trying to subdivide the face at the vertex with 'flattest'
|
||||
internal angle (i.e. closest to 180 deg).
|
||||
|
||||
Based on routine 'Diagonal' in
|
||||
\verbatim
|
||||
"Efficient Triangulation of Simple Polygons"
|
||||
Godfried Toussaint, McGill University.
|
||||
\endverbatim
|
||||
|
||||
After construction is the list of triangles the face is decomposed into.
|
||||
(Or empty list if no valid triangulation could be found).
|
||||
|
||||
|
||||
SourceFiles
|
||||
faceTriangulation.C
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef faceTriangulation_H
|
||||
#define faceTriangulation_H
|
||||
|
||||
#include "triFaceList.H"
|
||||
#include "pointField.H"
|
||||
|
||||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||||
|
||||
namespace Foam
|
||||
{
|
||||
|
||||
// Forward declaration of classes
|
||||
|
||||
/*---------------------------------------------------------------------------*\
|
||||
Class faceTriangulation Declaration
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
class faceTriangulation
|
||||
:
|
||||
public triFaceList
|
||||
{
|
||||
// Static Data
|
||||
|
||||
//- Relative tolerance on edge.
|
||||
static const scalar edgeRelTol;
|
||||
|
||||
|
||||
// Static Member Functions
|
||||
|
||||
//- Edge to the right of face vertex i
|
||||
static label right(const label size, label i);
|
||||
|
||||
//- Edge to the left of face vertex i
|
||||
static label left(const label size, label i);
|
||||
|
||||
//- Calculate normalised edge vectors
|
||||
static tmp<vectorField> 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
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -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<point>(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<point>(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<point>(points, faces[facei])
|
||||
);
|
||||
|
||||
if (faceTris.empty())
|
||||
{
|
||||
WarningInFunction
|
||||
<< "Could not find triangulation for face " << facei
|
||||
<< " vertices " << faces[facei] << " coords "
|
||||
<< IndirectList<point>(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
|
||||
);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@ -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<labelledTri> tris(nElems);
|
||||
|
||||
@ -107,17 +111,11 @@ bool Foam::triSurface::readOFF(const fileName& OFFfileName)
|
||||
}
|
||||
else
|
||||
{
|
||||
faceList triFaces(f.nTriangles(points));
|
||||
triEngine.triangulate(UIndirectList<point>(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));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user