Commit 24b10b2e authored by graham's avatar graham
Browse files

WIP commit. Modification to triSurfaceTools for triangle edge/face finding

tolerance.  Adding edgeDirection function to featureEdge to automatically
reverse direction if necessary.  Feature point handling experiments.
parent f00b1562
......@@ -91,7 +91,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
{
timeCheck();
conformToFeaturePoints();
createFeaturePoints();
timeCheck();
insertInitialPoints();
......@@ -294,7 +294,7 @@ void Foam::conformalVoronoiMesh::insertInternalEdgePointGroup
const vector& nB = feNormals[edNormalIs[1]];
// Concave. master and reflected points inside the domain.
point refPt = edgePt - ppDist*(nA + nB)/(1 + (nA & nB) + VSMALL );
point refPt = edgePt - ppDist*(nA + nB)/(1 + (nA & nB) + VSMALL);
// Generate reflected master to be outside.
point reflMasterPt = refPt + 2*(edgePt - refPt);
......@@ -375,8 +375,30 @@ void Foam::conformalVoronoiMesh::insertMultipleEdgePointGroup
}
void Foam::conformalVoronoiMesh::conformToFeaturePoints()
void Foam::conformalVoronoiMesh::createFeaturePoints()
{
const boundBox& bb = geometryToConformTo().bounds();
scalar span
(
max(mag(bb.max().x()), mag(bb.min().x()))
+ max(mag(bb.max().y()), mag(bb.min().y()))
+ max(mag(bb.max().z()), mag(bb.min().z()))
);
Info<< nl << "Creating bounding points" << endl;
scalar bigSpan = 10*span;
insertPoint(point(-bigSpan, -bigSpan, -bigSpan), Vb::FAR_POINT);
insertPoint(point(-bigSpan, -bigSpan, bigSpan), Vb::FAR_POINT);
insertPoint(point(-bigSpan, bigSpan, -bigSpan), Vb::FAR_POINT);
insertPoint(point(-bigSpan, bigSpan, bigSpan), Vb::FAR_POINT);
insertPoint(point( bigSpan, -bigSpan, -bigSpan), Vb::FAR_POINT);
insertPoint(point( bigSpan, -bigSpan, bigSpan), Vb::FAR_POINT);
insertPoint(point( bigSpan, bigSpan, -bigSpan), Vb::FAR_POINT);
insertPoint(point( bigSpan, bigSpan , bigSpan), Vb::FAR_POINT);
Info<< nl << "Conforming to feature points" << endl;
insertConvexFeaturesPoints();
......@@ -421,14 +443,14 @@ void Foam::conformalVoronoiMesh::insertConvexFeaturesPoints()
for
(
label pI = feMesh.convexStart();
pI < feMesh.concaveStart();
pI++
label ptI = feMesh.convexStart();
ptI < feMesh.concaveStart();
ptI++
)
{
vectorField featPtNormals = feMesh.featurePointNormals(pI);
vectorField featPtNormals = feMesh.featurePointNormals(ptI);
const point& featPt = feMesh.points()[pI];
const point& featPt = feMesh.points()[ptI];
vector cornerNormal = sum(featPtNormals);
cornerNormal /= mag(cornerNormal);
......@@ -464,14 +486,14 @@ void Foam::conformalVoronoiMesh::insertConcaveFeaturePoints()
for
(
label pI = feMesh.concaveStart();
pI < feMesh.mixedStart();
pI++
label ptI = feMesh.concaveStart();
ptI < feMesh.mixedStart();
ptI++
)
{
vectorField featPtNormals = feMesh.featurePointNormals(pI);
vectorField featPtNormals = feMesh.featurePointNormals(ptI);
const point& featPt = feMesh.points()[pI];
const point& featPt = feMesh.points()[ptI];
vector cornerNormal = sum(featPtNormals);
cornerNormal /= mag(cornerNormal);
......@@ -494,7 +516,7 @@ void Foam::conformalVoronoiMesh::insertConcaveFeaturePoints()
internalPtIndex = insertPoint(internalPt, externalPtIndex);
}
insertPoint(externalPt,internalPtIndex);
insertPoint(externalPt, internalPtIndex);
}
}
}
......@@ -502,24 +524,167 @@ void Foam::conformalVoronoiMesh::insertConcaveFeaturePoints()
void Foam::conformalVoronoiMesh::insertMixedFeaturePoints()
{
Info<< " Mixed feature points not implemented." << endl;
const PtrList<featureEdgeMesh>& feMeshes(geometryToConformTo_.features());
// const PtrList<featureEdgeMesh>& feMeshes(geometryToConformTo_.features());
forAll(feMeshes, i)
{
const featureEdgeMesh& feMesh(feMeshes[i]);
// forAll(feMeshes, i)
// {
// const featureEdgeMesh& feMesh(feMeshes[i]);
for
(
label ptI = feMesh.mixedStart();
ptI < feMesh.nonFeatureStart();
ptI++
)
{
labelList pEds(feMesh.pointEdges()[ptI]);
// for
// (
// label pI = feMesh.mixedStart();
// pI < feMesh.nonFeatureStart();
// pI++
// )
// {
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// }
// }
// Hard coding mixed type decision making to be an integer sum, as
// only 3 edge types allowed:
// edgeTypeSum == 1 is 2 EXTERNAL, 1 INTERNAL
// edgeTypeSum == 2 is 1 EXTERNAL, 2 INTERNAL
label edgeTypeSum = 0;
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Skipping unsupported mixed feature point types
bool skipEdge = false;
// if (pEds.size() > 3)
// {
// Info<< " Mixed feature points of 3 edges only supported."
// << " Point " << ptI << " has " << pEds.size()
// << " edges." << endl;
// skipEdge = true;
// }
forAll(pEds, e)
{
label edgeI = pEds[e];
featureEdgeMesh::edgeStatus edStatus =
feMesh.getEdgeStatus(edgeI);
edgeTypeSum += edStatus;
if
(
edStatus == featureEdgeMesh::FLAT
|| edStatus == featureEdgeMesh::OPEN
|| edStatus == featureEdgeMesh::MULTIPLE
)
{
Info<< " Edge type " << edStatus
<< " found for mixed feature point " << ptI
<< ". Not supported."
<< endl;
skipEdge = true;
}
}
if(skipEdge)
{
Info<< " Skipping point " << ptI << nl << endl;
continue;
}
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
const point& pt(feMesh.points()[ptI]);
// if (edgeTypeSum == 1 || edgeTypeSum == 2)
// {
// Info<< " 2 EXTERNAL, 1 INTERNAL mixed feature point" << nl
// << " or 1 EXTERNAL, 2 INTERNAL mixed feature point" << endl;
scalar ppDist = pointPairDistance(pt);
forAll(pEds, e)
{
label edgeI = pEds[e];
const vectorField& feNormals = feMesh.normals();
const labelList& edNormalIs =
feMesh.edgeNormals()[edgeI];
// As this is an external or internal edge, there are two
// normals by definition
const vector& nA = feNormals[edNormalIs[0]];
const vector& nB = feNormals[edNormalIs[1]];
point edgePt =
pt + 8.0*ppDist*feMesh.edgeDirection(edgeI, ptI);
// Concave. master and reflected points inside the domain.
point refPt =
edgePt - ppDist*(nA + nB)/(1 + (nA & nB) + VSMALL);
featureEdgeMesh::edgeStatus edStatus =
feMesh.getEdgeStatus(edgeI);
if (edStatus == featureEdgeMesh::INTERNAL)
{
// Generate reflected master to be outside.
point reflMasterPt = refPt + 2*(edgePt - refPt);
// Reflect reflMasterPt in both faces.
point reflectedA = reflMasterPt - 2*ppDist*nA;
point reflectedB = reflMasterPt - 2*ppDist*nB;
// index of reflMaster
label reflectedMaster = number_of_vertices() + 2;
// Master A is inside.
label reflectedAI =
insertPoint(reflectedA, reflectedMaster);
// Master B is inside.
insertPoint(reflectedB, reflectedMaster);
// Slave is outside.
insertPoint(reflMasterPt, reflectedAI);
}
else if (edStatus == featureEdgeMesh::EXTERNAL)
{
// Insert the master point referring the the first slave
label masterPtIndex =
insertPoint(refPt, number_of_vertices() + 1);
// Insert the slave points by reflecting refPt in both
// faces. with each slave refering to the master
point reflectedA = refPt + 2*ppDist*nA;
insertPoint(reflectedA, masterPtIndex);
point reflectedB = refPt + 2*ppDist*nB;
insertPoint(reflectedB, masterPtIndex);
}
}
// }
// else if (edgeTypeSum == 2)
// {
// Info<< "1 EXTERNAL, 2 INTERNAL mixed feature point" << endl;
// }
// else
// {
// FatalErrorIn
// (
// "Foam::conformalVoronoiMesh::insertMixedFeaturePoints"
// )
// << "Invalid edgeTypeSum " << edgeTypeSum
// << " for point " << ptI
// << exit(FatalError);
// }
}
}
}
......
......@@ -214,8 +214,8 @@ class conformalVoronoiMesh
const pointIndexHit& edHit
);
//- Insert point groups at the feature points.
void conformToFeaturePoints();
//- Determine and insert point groups at the feature points.
void createFeaturePoints();
//- Insert point groups at convex feature points
void insertConvexFeaturesPoints();
......
......@@ -308,6 +308,9 @@ public:
//- Return the edgeDirection vectors
inline const vectorField& edgeDirections() const;
//- Return the direction of edgeI, pointing away from ptI
inline vector edgeDirection(label edgeI, label ptI) const;
//- Return the indices of the normals that are adjacent to the
// feature edges
inline const labelListList& edgeNormals() const;
......
......@@ -97,6 +97,35 @@ inline const Foam::vectorField& Foam::featureEdgeMesh::edgeDirections() const
}
inline Foam::vector Foam::featureEdgeMesh::edgeDirection
(
label edgeI,
label ptI
) const
{
const edge& e = edges()[edgeI];
if (ptI == e.start())
{
return edgeDirections()[edgeI];
}
else if (ptI == e.end())
{
return -edgeDirections()[edgeI];
}
else
{
FatalErrorIn("Foam::featureEdgeMesh::edgedirection")
<< "Requested ptI " << ptI << " is not a point on the requested "
<< "edgeI " << edgeI << ". edgeI start and end: "
<< e.start() << " " << e.end()
<< exit(FatalError);
return vector::zero;
}
}
inline const Foam::labelListList& Foam::featureEdgeMesh::edgeNormals() const
{
return edgeNormals_;
......
......@@ -158,7 +158,7 @@ void Foam::triSurfaceTools::greenRefine
}
// Refine all triangles marked for refinement.
// Refine all triangles marked for refinement.
Foam::triSurface Foam::triSurfaceTools::doRefine
(
const triSurface& surf,
......@@ -358,7 +358,7 @@ void Foam::triSurfaceTools::protectNeighbours
forAll(myFaces, myFaceI)
{
label faceI = myFaces[myFaceI];
if ((faceStatus[faceI] == ANYEDGE) || (faceStatus[faceI] >= 0))
{
faceStatus[faceI] = NOEDGE;
......@@ -392,10 +392,10 @@ Foam::labelHashSet Foam::triSurfaceTools::getCollapsedFaces
{
facesToBeCollapsed.insert(myFaces[myFaceI]);
}
// From faces using v1 check if they share an edge with faces
// using v2.
// - share edge: are part of 'splay' tree and will collapse if edge
// - share edge: are part of 'splay' tree and will collapse if edge
// collapses
const labelList& v1Faces = surf.pointFaces()[v1];
......@@ -795,7 +795,7 @@ bool Foam::triSurfaceTools::collapseCreatesFold
//
// // Now neighbours contains first layer of triangles outside of
// // collapseFaces
// // There should be
// // There should be
// // -two if edgeI is a boundary edge
// // since the outside 'edge' of collapseFaces should
// // form a triangle and the face connected to edgeI is not inserted.
......@@ -816,7 +816,7 @@ bool Foam::triSurfaceTools::collapseCreatesFold
// for (label j = i+1; j < neighbourList.size(); i++)
// {
// const labelList& faceJEdges = surf.faceEdges()[neighbourList[j]];
//
//
// // Check if faceI and faceJ share an edge
// forAll(faceIEdges, fI)
// {
......@@ -1194,7 +1194,7 @@ Foam::surfaceLocation Foam::triSurfaceTools::visitFaces
nearest.setHit();
nearest.triangle() = triI;
break;
}
}
else
{
// Which edge is cut.
......@@ -1913,7 +1913,7 @@ Foam::triSurface Foam::triSurfaceTools::greenRefine
// Storage for new faces
DynamicList<labelledTri> newFaces(0);
pointField newPoints(surf.localPoints());
newPoints.setSize(surf.nPoints() + surf.nEdges());
label newPointI = surf.nPoints();
......@@ -2225,6 +2225,24 @@ Foam::triSurfaceTools::sideType Foam::triSurfaceTools::surfaceSide
// Nearest to face interior. Use faceNormal to determine side
scalar c = (sample - nearestPoint) & surf.faceNormals()[nearestFaceI];
c /= (mag(sample - nearestPoint)*mag(surf.faceNormals()[nearestFaceI]));
if (mag(c) < 0.9)
{
FatalErrorIn("triSurfaceTools::surfaceSide")
<< "nearestPoint identified as being on triangle face "
<< "but vector from nearestPoint to sample is not "
<< "perpendicular to the normal." << nl
<< "sample: " << sample << nl
<< "nearestPoint: " << nearestPoint << nl
<< "normalised dot product: " << c << nl
<< "triangle vertices: " << nl
<< " " << points[f[0]] << nl
<< " " << points[f[1]] << nl
<< " " << points[f[2]] << nl
<< abort(FatalError);
}
if (c > 0)
{
return OUTSIDE;
......@@ -2635,7 +2653,7 @@ Foam::triSurface Foam::triSurfaceTools::delaunay2D(const List<vector2D>& pts)
// // ~~~~~~~~~~~~~~~
//
// label vertI = 0;
// for
// for
// (
// Vertex_iterator it = T.vertices_begin();
// it != T.vertices_end();
......@@ -2849,7 +2867,7 @@ Foam::surfaceLocation Foam::triSurfaceTools::classify
// Nearest point could be on point or edge. Retest.
label index, elemType;
//bool inside =
//bool inside =
triPointRef(s[triI].tri(s.points())).classify
(
trianglePoint,
......@@ -2865,7 +2883,7 @@ Foam::surfaceLocation Foam::triSurfaceTools::classify
nearest.setHit();
nearest.setIndex(triI);
nearest.elementType() = triPointRef::NONE;
}
}
else if (elemType == triPointRef::EDGE)
{
nearest.setMiss();
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment