zhyang-dev/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

fixed checkOrientation method - triSurface (and PrimitivePatchExtra)

Browse Source 
- the old code used the edge information and examined the next face edges
    to find the orientation. This fails since the direction of the edge
    itself is missing.

  - simpler: find the edge start on both faces, check the next face point.
    If they are the same, the edge goes in the same direction on both faces
    and thus the orientation is incorrect.
This commit is contained in:
Mark Olesen
2008-08-10 19:28:30 +02:00
parent 120e32a58d
commit b87f997fdf
4 changed files with 93 additions and 157 deletions
Download Patch File Download Diff File Expand all files Collapse all files

77
src/OpenFOAM/meshes/primitiveMesh/PrimitivePatchExtra/PrimitivePatchExtra.H
View File

@ -79,8 +79,8 @@ private:
//- Edge-face addressing (sorted)
mutable labelListList* sortedEdgeFacesPtr_;
//- Label of face that 'owns' edge (i.e. e.vec() is righthanded walk
// along face)
//- Label of face that 'owns' edge
// i.e. e.vec() is righthanded walk along face
mutable labelList* edgeOwnerPtr_;
@ -92,42 +92,6 @@ private:
//- Calculate owner
void calcEdgeOwner() const;
protected:
// Protected Member Functions
// Edit
//- Fill faceZone with currentZone for every face reachable
// from faceI without crossing edge marked in borderEdge.
// Note: faceZone has to be sized nFaces before calling this fun.
void markZone
(
const boolList& borderEdge,
const label faceI,
const label currentZone,
labelList& faceZone
) const;
//- (size and) fills faceZone with zone of face.
// Zone is area reachable by edge crossing without crossing borderEdge
// (bool for every edge in surface). Returns number of zones.
label markZones
(
const boolList& borderEdge,
labelList& faceZone
) const;
//- Determine the mapping for a sub mesh.
// Only include faces for which boolList entry is true
// Sets: pointMap: from new to old localPoints
// faceMap: new to old faces
void subsetMap
(
const boolList& include,
labelList& pointMap,
labelList& faceMap
) const;
public:
// Constructors
@ -183,14 +147,47 @@ public:
// Check
//- Check triply (or more) connected edges.
// Return list of faces sharing these edges.
//- Check multiply-connected edges.
void checkEdges(const bool verbose) const;
//- Check orientation (normals) and normals of neighbouring faces
boolList checkOrientation(const bool verbose) const;
// Edit
//- Fill faceZone with currentZone for every face reachable
// from faceI without crossing edge marked in borderEdge.
// Note: faceZone has to be sized nFaces before calling.
void markZone
(
const boolList& borderEdge,
const label faceI,
const label currentZone,
labelList& faceZone
) const;
//- (size and) fills faceZone with zone of face.
// Zone is area reachable by edge crossing without crossing borderEdge
// (bool for every edge in surface). Returns number of zones.
label markZones
(
const boolList& borderEdge,
labelList& faceZone
) const;
//- Determine the mapping for a sub-mesh.
// Only include faces for which boolList entry is true
// @param[out] pointMap mapping new to old localPoints
// @param[out] faceMap mapping new to old faces
void subsetMap
(
const boolList& include,
labelList& pointMap,
labelList& faceMap
) const;
};

84
src/OpenFOAM/meshes/primitiveMesh/PrimitivePatchExtra/PrimitivePatchExtraCleanup.C
View File

@ -28,14 +28,6 @@ License
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Check/fix edges with more than two faces
@ -52,16 +44,15 @@ checkEdges
const bool verbose
) const
{
const labelListList& eFaces =
PrimitivePatch<Face, ListType, PointField>::edgeFaces();
const labelListList& eFaces = TemplateType::edgeFaces();
const edgeList& edgeLst = TemplateType::edges();
const edgeList& edgeLst =
PrimitivePatch<Face, ListType, PointField>::edges();
forAll (eFaces, edgeI)
forAll(eFaces, edgeI)
{
const labelList& myFaces = eFaces[edgeI];
// boundary edges have one face
// interior edges have two faces
if (myFaces.size() == 0)
{
FatalErrorIn("PrimitivePatchExtra::checkEdges(bool verbose)")
@ -97,15 +88,15 @@ checkOrientation
const bool verbose
) const
{
const ListType<FaceType>& faceLst = *this;
const edgeList& edgeLst = TemplateType::edges();
const labelListList& faceEs = TemplateType::faceEdges();
const List<FaceType>& faceLst = TemplateType::faces();
const label numEdges = TemplateType::nEdges();
const pointField& pointLst = TemplateType::points();
const vectorField& normLst = TemplateType::faceNormals();
// Check edge normals, face normals, point normals.
forAll (faceEs, faceI)
forAll(faceEs, faceI)
{
const labelList& edgeLabels = faceEs[faceI];
@ -118,7 +109,7 @@ checkOrientation
}
bool valid = true;
forAll (edgeLabels, i)
forAll(edgeLabels, i)
{
if (edgeLabels[i] < 0 || edgeLabels[i] >= numEdges)
{
@ -126,8 +117,8 @@ checkOrientation
(
"PrimitivePatchExtra::checkOrientation(bool)"
) << "edge number " << edgeLabels[i] << " on face " << faceI
<< " out of range"
<< "\nThis usually means that the input surface has "
<< " out-of-range\n"
<< "This usually means that the input surface has "
<< "edges with more than 2 faces connected.\n"
<< endl;
valid = false;
@ -141,7 +132,7 @@ checkOrientation
//
//- Compute normal from 3 points, use the first as the origin
//
// minor warpage should not be a problem
const FaceType& f = faceLst[faceI];
const point p0(pointLst[f[0]]);
const point p1(pointLst[f[1]]);
@ -169,58 +160,34 @@ checkOrientation
// edge
boolList borderEdge(numEdges, false);
forAll (edgeLst, edgeI)
forAll(edgeLst, edgeI)
{
const edge& e = edgeLst[edgeI];
const labelList& neighbours = eFaces[edgeI];
if (neighbours.size() == 2)
{
// Two faces, A and B. Check if edge orientation is
// anticlockwise on both.
const labelList& fEdgesA = faceEs[neighbours[0]];
const labelList& fEdgesB = faceEs[neighbours[1]];
const FaceType& faceA = faceLst[neighbours[0]];
const FaceType& faceB = faceLst[neighbours[1]];
// Get next edge after edgeI
label nextEdgeA = fEdgesA.fcIndex(findIndex(fEdgesA, edgeI));
label nextEdgeB = fEdgesB.fcIndex(findIndex(fEdgesB, edgeI));
// Now check if nextEdgeA and nextEdgeB have any common points
// The edge cannot be going in the same direction if both faces
// are oriented counterclockwise.
// Thus the next face point *must* different between the faces.
if
(
edgeLst[nextEdgeA].start() == edgeLst[nextEdgeB].start()
|| edgeLst[nextEdgeA].start() == edgeLst[nextEdgeB].end()
|| edgeLst[nextEdgeA].end() == edgeLst[nextEdgeB].start()
|| edgeLst[nextEdgeA].end() == edgeLst[nextEdgeB].end()
faceA[faceA.fcIndex(findIndex(faceA, e.start()))]
== faceB[faceB.fcIndex(findIndex(faceB, e.start()))]
)
{
borderEdge[edgeI] = true;
if (verbose)
{
// just list first three points
// to simplify generating the message
WarningIn("PrimitivePatchExtra::checkOrientation(bool)")
<< "face orientation incorrect." << nl
<< "edge neighbours:" << neighbours << nl
<< "face " << neighbours[0] << " has edges "
<< fEdgesA << nl
<< " with points " << nl
<< " " << edgeLst[fEdgesA[0]].start() << ' '
<< edgeLst[fEdgesA[0]].end() << nl
<< " " << edgeLst[fEdgesA[1]].start() << ' '
<< edgeLst[fEdgesA[1]].end() << nl
<< " " << edgeLst[fEdgesA[2]].start() << ' '
<< edgeLst[fEdgesA[2]].end()
<< endl
<< "face " << neighbours[1] << " has edges "
<< fEdgesB << nl
<< " with points " << nl
<< " " << edgeLst[fEdgesB[0]].start() << ' '
<< edgeLst[fEdgesB[0]].end() << nl
<< " " << edgeLst[fEdgesB[1]].start() << ' '
<< edgeLst[fEdgesB[1]].end() << nl
<< " " << edgeLst[fEdgesB[2]].start() << ' '
<< edgeLst[fEdgesB[2]].end() << nl
<< endl;
<< "edge[" << edgeI << "] " << e
<< " between faces " << neighbours << ":" << nl
<< "face[" << neighbours[0] << "] " << faceA << nl
<< "face[" << neighbours[1] << "] " << faceB << endl;
}
}
}
@ -228,10 +195,9 @@ checkOrientation
{
if (verbose)
{
const edge& e = edgeLst[edgeI];
WarningIn("PrimitivePatchExtra::checkOrientation(bool)")
<< "Wrong number of edge neighbours." << endl
<< "Edge:" << e
<< "edge[" << edgeI << "] " << e
<< "with points:" << locPointsLst[e.start()]
<< ' ' << locPointsLst[e.end()]
<< " has neighbours:" << neighbours << endl;

10
src/OpenFOAM/meshes/primitiveMesh/PrimitivePatchExtra/PrimitivePatchExtraSearch.C
View File

@ -26,8 +26,6 @@ License
#include "PrimitivePatchExtra.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Finds area, starting at faceI, delimited by borderEdge. Marks all visited
@ -58,13 +56,13 @@ void Foam::PrimitivePatchExtra<Face, ListType, PointField, PointType>::markZone
// Pick up neighbours of changedFaces
DynamicList<label> newChangedFaces(2*changedFaces.size());
forAll (changedFaces, i)
forAll(changedFaces, i)
{
label faceI = changedFaces[i];
const labelList& fEdges = faceEs[faceI];
forAll(fEdges, i)
forAllfEdges, i)
{
label edgeI = fEdges[i];
@ -202,7 +200,7 @@ subsetMap
boolList pointHad(numPoints, false);
forAll (include, oldFaceI)
forAll(include, oldFaceI)
{
if (include[oldFaceI])
{
@ -212,7 +210,7 @@ subsetMap
// Renumber labels for face
const FaceType& f = locFaces[oldFaceI];
forAll (f, fp)
forAll(f, fp)
{
const label ptLabel = f[fp];
if (!pointHad[ptLabel])

79
src/triSurface/triSurface/triSurface.C
View File

@ -161,7 +161,7 @@ void triSurface::printTriangle
const pointField& points
)
{
os
os
<< pre.c_str() << "vertex numbers:"
<< f[0] << ' ' << f[1] << ' ' << f[2] << endl
<< pre.c_str() << "vertex coords :"
@ -274,13 +274,13 @@ void triSurface::checkTriangles(const bool verbose)
"triSurface::checkTriangles(bool verbose)"
) << "triangles share the same vertices:\n"
<< " face 1 :" << faceI << endl;
printTriangle(Warning, " ", f, points());
printTriangle(Warning, " ", f, points());
Warning
<< endl
<< " face 2 :"
<< neighbours[neighbourI] << endl;
printTriangle(Warning, " ", n, points());
printTriangle(Warning, " ", n, points());
}
break;
@ -375,8 +375,8 @@ boolList triSurface::checkOrientation(const bool verbose)
(
"triSurface::checkOrientation(bool)"
) << "edge number " << edgeLabels[i] << " on face " << facei
<< " out of range"
<< "\nThis usually means that the input surface has "
<< " out-of-range\n"
<< "This usually means that the input surface has "
<< "edges with more than 2 triangles connected.\n"
<< endl;
valid = false;
@ -413,63 +413,39 @@ boolList triSurface::checkOrientation(const bool verbose)
const labelListList& eFaces = edgeFaces();
// Storage for holding status of edge. True if normal flips across this
// edge
boolList borderEdge(nEdges(), false);
forAll(es, edgei)
forAll(es, edgeI)
{
const labelList& neighbours = eFaces[edgei];
const edge& e = es[edgeI];
const labelList& neighbours = eFaces[edgeI];
if (neighbours.size() == 2)
{
// Two triangles, A and B. Check if edge orientation is
// anticlockwise on both.
const labelList& fEdgesA = faceEdges()[neighbours[0]];
const labelledTri& faceA = (*this)[neighbours[0]];
const labelledTri& faceB = (*this)[neighbours[1]];
// Get next edge after edgei
label nextEdgeA = fEdgesA.fcIndex(findIndex(fEdgesA, edgei));
const labelList& fEdgesB = faceEdges()[neighbours[1]];
label nextEdgeB = fEdgesB.fcIndex(findIndex(fEdgesB, edgei));
// Now check if nextEdgeA and nextEdgeB have any common points
// The edge cannot be going in the same direction if both faces
// are oriented counterclockwise.
// Thus the next face point *must* different between the faces.
if
(
(es[nextEdgeA].start() == es[nextEdgeB].start())
|| (es[nextEdgeA].start() == es[nextEdgeB].end())
|| (es[nextEdgeA].end() == es[nextEdgeB].start())
|| (es[nextEdgeA].end() == es[nextEdgeB].end())
faceA[faceA.fcIndex(findIndex(faceA, e.start()))]
== faceB[faceB.fcIndex(findIndex(faceB, e.start()))]
)
{
borderEdge[edgei] = true;
borderEdge[edgeI] = true;
if (verbose)
{
WarningIn("triSurface::checkOrientation(bool)")
<< "Triangle orientation incorrect." << endl
<< "edge neighbours:" << neighbours << endl
<< "triangle " << neighbours[0] << " has edges "
<< fEdgesA << endl
<< " with points " << endl
<< " " << es[fEdgesA[0]].start() << ' '
<< es[fEdgesA[0]].end() << endl
<< " " << es[fEdgesA[1]].start() << ' '
<< es[fEdgesA[1]].end() << endl
<< " " << es[fEdgesA[2]].start() << ' '
<< es[fEdgesA[2]].end() << endl
<< "triangle " << neighbours[1] << " has edges "
<< fEdgesB << endl
<< " with points " << endl
<< " " << es[fEdgesB[0]].start() << ' '
<< es[fEdgesB[0]].end() << endl
<< " " << es[fEdgesB[1]].start() << ' '
<< es[fEdgesB[1]].end() << endl
<< " " << es[fEdgesB[2]].start() << ' '
<< es[fEdgesB[2]].end() << endl
<< endl;
WarningIn("PrimitivePatchExtra::checkOrientation(bool)")
<< "face orientation incorrect." << nl
<< "edge[" << edgeI << "] " << e
<< " between faces " << neighbours << ":" << nl
<< "face[" << neighbours[0] << "] " << faceA << nl
<< "face[" << neighbours[1] << "] " << faceB << endl;
}
}
}
@ -477,15 +453,14 @@ boolList triSurface::checkOrientation(const bool verbose)
{
if (verbose)
{
const edge& e = es[edgei];
WarningIn("triSurface::checkOrientation(bool)")
<< "Wrong number of edge neighbours." << endl
<< "Edge:" << e
<< "edge[" << edgeI << "] " << e
<< "with points:" << localPoints()[e.start()]
<< ' ' << localPoints()[e.end()]
<< " has neighbours:" << neighbours << endl;
}
borderEdge[edgei] = true;
borderEdge[edgeI] = true;
}
}
@ -1109,7 +1084,7 @@ void triSurface::subsetMeshMap
if (!pointHad[a])
{
pointHad[a] = true;
pointMap[pointI++] = a;
pointMap[pointI++] = a;
}
label b = tri[1];
@ -1147,7 +1122,7 @@ triSurface triSurface::subsetMesh
// Fill pointMap, faceMap
subsetMeshMap(include, pointMap, faceMap);
// Create compact coordinate list and forward mapping array
pointField newPoints(pointMap.size());