Commit 8eef91c5 authored by Mark Olesen's avatar Mark Olesen
Browse files

ENH: improve edge access for face/triFace

- additional rcEdge(), rcEdges() methods for reverse order walk

- accept generic edge() method as alternative to faceEdge() for
  single edge retrieval.

- edge() method with points -> returns the vector

- reduce the number of operations in edgeDirection methods

DEFEATURE: remove longestEdge global function

- deprecated and replaced by face::longestEdge() method (2017-04)
parent 5eb48c44
......@@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016 OpenCFD Ltd.
Copyright (C) 2016-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -33,6 +33,8 @@ Description
#include "argList.H"
#include "labelledTri.H"
#include "faceList.H"
#include "triFaceList.H"
#include "pointList.H"
#include "ListOps.H"
......@@ -63,6 +65,35 @@ void testSign
}
template<class Face>
void testEdges(const Face& f)
{
const label nEdges = f.nEdges();
Info<< "face: " << f << nl
<< "flip: " << f.reverseFace() << nl
<< " fc edges:" << flatOutput(f.edges()) << nl
<< " rc edges:" << flatOutput(f.rcEdges()) << nl;
Info<< " forward edges" << nl;
for (label edgei = 0; edgei < nEdges; ++edgei)
{
Info<< " " << edgei << " : " << f.edge(edgei) << nl;
}
Info<< " reverse edges" << nl;
for (label edgei = 0; edgei < nEdges; ++edgei)
{
Info<< " " << edgei << " : " << f.rcEdge(edgei) << nl;
}
}
void testCompare(const triFace& a, const triFace& b)
{
Info<< "compare: " << a << " with " << b
<< " == " << triFace::compare(a, b) << nl;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
......@@ -92,7 +123,7 @@ int main(int argc, char *argv[])
{ 2, 2, 2}
});
face f1{ 1, 2, 3, 4 };
face f1({1, 2, 3, 4});
Info<< "face:"; faceInfo(f1, points1); Info << nl;
testSign(f1, points1, testPoints);
......@@ -100,7 +131,7 @@ int main(int argc, char *argv[])
testSign(f1, points2, testPoints);
Info<< nl;
triFace t1{ 1, 2, 3 };
triFace t1({1, 2, 3});
Info<< "triFace:"; faceInfo(t1, points1); Info << nl;
testSign(t1, points1, testPoints);
......@@ -115,11 +146,12 @@ int main(int argc, char *argv[])
f1 = t1.triFaceFace();
Info<< "face:" << f1 << nl;
// expect these to fail
#if 0
// Expect failure, but triggers abort which cannot be caught
const bool throwingError = FatalError.throwExceptions();
try
{
labelledTri l1{ 1, 2, 3, 10, 24 };
labelledTri l1({1, 2, 3, 10, 24});
Info<< "labelled:" << l1 << nl;
}
catch (const Foam::error& err)
......@@ -128,11 +160,12 @@ int main(int argc, char *argv[])
<< "Caught FatalError " << err << nl << endl;
}
FatalError.throwExceptions(throwingError);
#endif
labelledTri l2{ 1, 2, 3 };
labelledTri l2({1, 2, 3});
Info<< "labelled:" << l2 << nl;
labelledTri l3{ 1, 2, 3, 10 };
labelledTri l3({1, 2, 3, 10});
Info<< "labelled:" << l3 << nl;
t1.flip();
......@@ -141,6 +174,62 @@ int main(int argc, char *argv[])
Info<< "flip:" << t1 << nl;
Info<< "flip:" << l3 << nl;
{
triFaceList faceList1
({
triFace{1, 2, 3},
triFace{4, 2, 100},
triFace{1, 3, 2},
});
Info<< nl << "Test edges" << nl;
for (const auto& f : faceList1)
{
testEdges(f);
Info<< nl;
}
}
{
faceList faceList1
({
face{1, 2, 3, 4},
face{1, 4, 3, 2},
face{4, 2, 100, 8, 35},
});
Info<< nl << "Test edges" << nl;
for (const auto& f : faceList1)
{
testEdges(f);
Info<< nl;
}
}
{
triFaceList faceList1
({
triFace{1, 2, 3},
triFace{1, 3, 2},
triFace{3, 1, 2},
triFace{4, 5, 1},
});
Info<< nl << "Test triFace compare" << nl;
for (const triFace& a : faceList1)
{
for (const triFace& b : faceList1)
{
testCompare(a, b);
}
}
Info<< nl;
}
Info<< "\nEnd\n" << endl;
return 0;
}
......
......@@ -31,6 +31,7 @@ License
#include "triPointRef.H"
#include "mathematicalConstants.H"
#include "ConstCirculator.H"
#include <algorithm>
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
......@@ -40,25 +41,18 @@ const char* const Foam::face::typeName = "face";
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::tmp<Foam::vectorField>
Foam::face::calcEdges(const UList<point>& points) const
Foam::face::calcEdgeVectors(const UList<point>& points) const
{
tmp<vectorField> tedges(new vectorField(size()));
vectorField& edges = tedges.ref();
auto tedgeVecs = tmp<vectorField>::New(size());
auto& edgeVecs = tedgeVecs.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;
edgeVecs[i] = vector(points[nextLabel(i)] - points[thisLabel(i)]);
edgeVecs[i].normalise();
}
return tedges;
return tedgeVecs;
}
......@@ -68,11 +62,11 @@ Foam::scalar Foam::face::edgeCos
const label index
) const
{
label leftEdgeI = left(index);
label rightEdgeI = right(index);
const vector& leftEdge = edges[rcIndex(index)];
const vector& rightEdge = edges[index];
// Note negate on left edge to get correct left-pointing edge.
return -(edges[leftEdgeI] & edges[rightEdgeI]);
return -(leftEdge & rightEdge);
}
......@@ -90,12 +84,12 @@ Foam::label Foam::face::mostConcaveAngle
forAll(edges, i)
{
label leftEdgeI = left(i);
label rightEdgeI = right(i);
const vector& leftEdge = edges[rcIndex(i)];
const vector& rightEdge = edges[i];
vector edgeNormal = edges[rightEdgeI] ^ edges[leftEdgeI];
vector edgeNormal = (rightEdge ^ leftEdge);
scalar edgeCos = edges[leftEdgeI] & edges[rightEdgeI];
scalar edgeCos = (leftEdge & rightEdge);
scalar edgeAngle = acos(max(-1.0, min(1.0, edgeCos)));
scalar angle;
......@@ -133,15 +127,15 @@ Foam::label Foam::face::split
faceList& quadFaces
) const
{
label oldIndices = (triI + quadI);
const label oldIndices = (triI + quadI);
if (size() <= 2)
if (size() < 3)
{
FatalErrorInFunction
<< "Serious problem: asked to split a face with < 3 vertices"
<< abort(FatalError);
}
if (size() == 3)
else if (size() == 3)
{
// Triangle. Just copy.
if (mode == COUNTTRIANGLE || mode == COUNTQUAD)
......@@ -166,7 +160,7 @@ Foam::label Foam::face::split
else if (mode == SPLITTRIANGLE)
{
// Start at point with largest internal angle.
const vectorField edges(calcEdges(points));
const vectorField edges(calcEdgeVectors(points));
scalar minAngle;
label startIndex = mostConcaveAngle(points, edges, minAngle);
......@@ -197,13 +191,13 @@ Foam::label Foam::face::split
{
// General case. Like quad: search for largest internal angle.
const vectorField edges(calcEdges(points));
const vectorField edges(calcEdgeVectors(points));
scalar minAngle = 1;
label startIndex = mostConcaveAngle(points, edges, minAngle);
scalar bisectAngle = minAngle/2;
vector rightEdge = edges[right(startIndex)];
const vector& rightEdge = edges[startIndex];
//
// Look for opposite point which as close as possible bisects angle
......@@ -222,7 +216,7 @@ Foam::label Foam::face::split
points[operator[](index)]
- points[operator[](startIndex)]
);
splitEdge /= Foam::mag(splitEdge) + VSMALL;
splitEdge.normalise();
const scalar splitCos = splitEdge & rightEdge;
const scalar splitAngle = acos(max(-1.0, min(1.0, splitCos)));
......@@ -786,62 +780,57 @@ Foam::tensor Foam::face::inertia
Foam::edgeList Foam::face::edges() const
{
const labelList& points = *this;
const labelList& verts = *this;
const label nVerts = verts.size();
edgeList theEdges(nVerts);
edgeList e(points.size());
// Last edge closes the polygon
theEdges.last().first() = verts.last();
theEdges.last().second() = verts[0];
for (label pointi = 0; pointi < points.size() - 1; ++pointi)
for (label verti = 0; verti < nVerts - 1; ++verti)
{
e[pointi] = edge(points[pointi], points[pointi + 1]);
theEdges[verti].first() = verts[verti];
theEdges[verti].second() = verts[verti + 1];
}
// Add last edge
e.last() = edge(points.last(), points[0]);
return e;
return theEdges;
}
int Foam::face::edgeDirection(const edge& e) const
Foam::edgeList Foam::face::rcEdges() const
{
forAll(*this, i)
const labelList& verts = *this;
const label nVerts = verts.size();
edgeList theEdges(nVerts);
// First edge closes the polygon
theEdges.first().first() = verts[0];
theEdges.first().second() = verts.last();
for (label verti = 1; verti < nVerts; ++verti)
{
if (operator[](i) == e.first())
{
if (operator[](rcIndex(i)) == e.second())
{
// Reverse direction
return -1;
}
else if (operator[](fcIndex(i)) == e.second())
{
// Forward direction
return 1;
}
theEdges[verti].first() = verts[nVerts - verti];
theEdges[verti].second() = verts[nVerts - verti - 1];
}
// No match
return 0;
}
else if (operator[](i) == e.second())
{
if (operator[](rcIndex(i)) == e.first())
{
// Forward direction
return 1;
}
else if (operator[](fcIndex(i)) == e.first())
{
// Reverse direction
return -1;
}
return theEdges;
}
// No match
return 0;
}
int Foam::face::edgeDirection(const Foam::edge& e) const
{
const label idx = find(e.first());
if (idx != -1)
{
if (e.second() == nextLabel(idx)) return 1; // Forward
if (e.second() == prevLabel(idx)) return -1; // Reverse
}
// Not found
return 0;
return 0; // Not found
}
......@@ -894,29 +883,26 @@ Foam::label Foam::face::trianglesQuads
Foam::label Foam::face::longestEdge(const UList<point>& pts) const
{
const edgeList& eds = this->edges();
const labelList& verts = *this;
const label nVerts = verts.size();
label longestEdgeI = -1;
scalar longestEdgeLength = -SMALL;
// Last edge closes the polygon. Use it to initialize loop
label longest = nVerts - 1;
scalar longestLen = Foam::edge(verts.first(), verts.last()).mag(pts);
forAll(eds, edI)
// Examine other edges
for (label edgei = 0; edgei < nVerts - 1; ++edgei)
{
scalar edgeLength = eds[edI].mag(pts);
scalar edgeLen = Foam::edge(verts[edgei], verts[edgei + 1]).mag(pts);
if (edgeLength > longestEdgeLength)
if (longestLen < edgeLen)
{
longestEdgeI = edI;
longestEdgeLength = edgeLength;
longest = edgei;
longestLen = edgeLen;
}
}
return longestEdgeI;
}
Foam::label Foam::longestEdge(const face& f, const UList<point>& pts)
{
return f.longestEdge(pts);
return longest;
}
......
......@@ -76,14 +76,8 @@ class face
{
// 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
tmp<vectorField> calcEdgeVectors
(
const UList<point>& points
) const;
......@@ -190,7 +184,7 @@ public:
void flip();
//- Return the points corresponding to this face
inline pointField points(const UList<point>& points) const;
inline pointField points(const UList<point>& pts) const;
//- Centre point of face
point centre(const UList<point>& points) const;
......@@ -227,13 +221,14 @@ public:
// Navigation through face vertices
//- Return true if the point label is found in face.
inline bool found(const label pointLabel) const;
//- Find local index on face for the point label,
//- Find local index on face for the point label, same as find()
// \return position in face (0,1,2,...) or -1 if not found.
inline label which(const label pointLabel) const;
//- The vertex on face - identical to operator[], but with naming
//- similar to nextLabel(), prevLabel()
inline label thisLabel(const label i) const;
//- Next vertex on face
inline label nextLabel(const label i) const;
......@@ -342,25 +337,45 @@ public:
) const;
//- Return number of edges
inline label nEdges() const;
inline label nEdges() const noexcept;
//- Return i-th face edge in forward walk order.
//- The faceEdge(0) is the edge between [0] and [1]
inline Foam::edge faceEdge(const label edgei) const;
//- Return i-th face edge in forward walk order.
//- Identical to faceEdge() but with generic name
inline Foam::edge edge(const label edgei) const;
//- Return edges in face point ordering,
//- Return vector of i-th face edge in forward walk order.
inline vector edge(const label edgei, const UList<point>& pts) const;
//- Return i-th face edge in reverse walk order.
//- The rcEdge(0) is the edge between [0] and [n-1]
inline Foam::edge rcEdge(const label edgei) const;
//- Return vector of i-th face edge in reverse walk order.
inline vector rcEdge(const label edgei, const UList<point>& pts) const;
//- Return list of edges in forward walk order.
// i.e. edges()[0] is edge between [0] and [1]
edgeList edges() const;
//- Return n-th face edge
inline edge faceEdge(const label n) const;
//- Return list of edges in reverse walk order.
// i.e. rcEdges()[0] is edge between [0] and [n-1]
edgeList rcEdges() const;
//- The edge direction on the face
//- Test the edge direction on the face
// \return
// - 0: edge not found on the face
// - +1: forward (counter-clockwise) on the face
// - -1: reverse (clockwise) on the face
int edgeDirection(const edge& e) const;
int edgeDirection(const Foam::edge& e) const;
//- Find the longest edge on a face.
label longestEdge(const UList<point>& pts) const;
// Face splitting utilities
//- Number of triangles after splitting
......@@ -493,13 +508,6 @@ struct offsetOp<face>
};
//- Deprecated(2017-04) find the longest edge on a face.
//- Face point labels index into pts.
// \deprecated(2017-04) use class method instead
FOAM_DEPRECATED_FOR(2017-04, "use face::longestEdge() method")
label longestEdge(const face& f, const UList<point>& pts);
// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * //
inline bool operator==(const face& a, const face& b);
......
......@@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011 OpenFOAM Foundation
Copyright (C) 2017-2020 OpenCFD Ltd.
Copyright (C) 2017-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
......@@ -26,20 +26,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(const label sz)
......@@ -98,23 +84,19 @@ inline Foam::face::face(Istream& is)
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline Foam::pointField Foam::face::points
(
const UList<point>& meshPoints
) const
inline Foam::pointField Foam::face::points(const UList<point>& pts) const
{
// There are as many points as there are labels for them
pointField p(size());
// For each point in list, set it to the point in 'pnts' addressed
// by 'labs'
label i = 0;
auto iter = p.begin();
for (const label pointi : *this)
{
p[i++] = meshPoints[pointi];
*iter = pts[pointi];
++iter;
}
// Return list
return p;
}
......@@ -133,22 +115,54 @@ inline Foam::scalar Foam::face::mag(const UList<point>& p) const
}
inline Foam::label Foam::face::nEdges() const
inline Foam::label Foam::face::nEdges() const noexcept
{
// for a closed polygon a number of edges is the same as number of points
return size();
}
inline Foam::edge Foam::face::faceEdge(const label n) const
inline Foam::edge Foam::face::faceEdge(const label edgei) const
{
return Foam::edge(thisLabel(edgei), nextLabel(edgei));
}
inline Foam::edge Foam::face::edge(const label edgei) const