Conformation to internal and external feature edges.

This commit is contained in:
graham
2009-04-30 20:21:50 +01:00
parent b487a13281
commit c5ee731267
7 changed files with 357 additions and 85 deletions

View File

@ -66,12 +66,12 @@ int main(int argc, char *argv[])
runTime++;
Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
<< endl;
}
Info << nl << "End\n" << endl;
Info<< nl << "End\n" << endl;
return 0;
}

View File

@ -134,7 +134,7 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
{
FatalErrorIn("Foam::conformalVoronoiMesh::insertPointPairs")
<< "surfacePpDist, surfacePoints and surfaceNormals are not "
<< "the same size. Sizes"
<< "the same size. Sizes "
<< surfacePpDist.size() << ' '
<< surfacePoints.size() << ' '
<< surfaceNormals.size()
@ -158,6 +158,189 @@ void Foam::conformalVoronoiMesh::insertSurfacePointPairs
}
void Foam::conformalVoronoiMesh::insertEdgePointGroups
(
const List<pointIndexHit>& edgeHits,
const labelList& featuresHit
)
{
if (edgeHits.size() != featuresHit.size())
{
FatalErrorIn("Foam::conformalVoronoiMesh::insertEdgePointGroups")
<< "edgeHits and featuresHit are not the same size. Sizes "
<< edgeHits.size() << ' '
<< featuresHit.size()
<< exit(FatalError);
}
forAll(edgeHits, i)
{
const featureEdgeMesh& feMesh
(
geometryToConformTo_.features()[featuresHit[i]]
);
label edgeI = edgeHits[i].index();
featureEdgeMesh::edgeStatus edStatus = feMesh.getEdgeStatus(edgeI);
if (edStatus == featureEdgeMesh::EXTERNAL)
{
insertExternalEdgePointGroup(feMesh, edgeHits[i]);
}
else if (edStatus == featureEdgeMesh::INTERNAL)
{
insertInternalEdgePointGroup(feMesh, edgeHits[i]);
}
else if (edStatus == featureEdgeMesh::FLAT)
{
insertFlatEdgePointGroup(feMesh, edgeHits[i]);
}
else if (edStatus == featureEdgeMesh::OPEN)
{
insertOpenEdgePointGroup(feMesh, edgeHits[i]);
}
else if (edStatus == featureEdgeMesh::MULTIPLE)
{
insertMultipleEdgePointGroup(feMesh, edgeHits[i]);
}
}
}
void Foam::conformalVoronoiMesh::insertExternalEdgePointGroup
(
const featureEdgeMesh& feMesh,
const pointIndexHit& edHit
)
{
const point& edgePt = edHit.hitPoint();
scalar ppDist = pointPairDistance(edgePt);
const vectorField& feNormals = feMesh.normals();
const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
// As this is an external edge, there are two normals by definition
const vector& nA = feNormals[edNormalIs[0]];
const vector& nB = feNormals[edNormalIs[1]];
// Convex. So refPt will be inside domain and hence a master point
point refPt = edgePt - ppDist*(nA + nB)/(1 + (nA & nB));
// 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);
}
void Foam::conformalVoronoiMesh::insertInternalEdgePointGroup
(
const featureEdgeMesh& feMesh,
const pointIndexHit& edHit
)
{
const point& edgePt = edHit.hitPoint();
scalar ppDist = pointPairDistance(edgePt);
const vectorField& feNormals = feMesh.normals();
const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
// As this is an external edge, there are two normals by definition
const vector& nA = feNormals[edNormalIs[0]];
const vector& nB = feNormals[edNormalIs[1]];
// Concave. master and reflected points inside the domain.
point refPt = edgePt - ppDist*(nA + nB)/(1 + (nA & nB));
// Generate reflected master to be outside.
point reflMasterPt = edgePt + 2*(edgePt - refPt);
// Reflect reflMasterPt in both faces.
point reflectedA = reflMasterPt - 2*ppDist*nA;
point reflectedB = reflMasterPt - 2*ppDist*nB;
scalar totalAngle =
180*(mathematicalConstant::pi + acos(mag(nA & nB)))
/mathematicalConstant::pi;
// Number of quadrants the angle should be split into
int nQuads = int(totalAngle/cvMeshControls_.maxQuadAngle()) + 1;
// The number of additional master points needed to obtain the
// required number of quadrants.
int nAddPoints = min(max(nQuads - 2, 0), 2);
// index of reflMaster
label reflectedMaster = number_of_vertices() + 2 + nAddPoints;
// Master A is inside.
label reflectedAI = insertPoint(reflectedA, reflectedMaster);
// Master B is inside.
insertPoint(reflectedB, reflectedMaster);
if (nAddPoints == 1)
{
// One additinal point is the reflection of the slave point,
// i.e. the original reference point
insertPoint(refPt, reflectedMaster);
}
else if (nAddPoints == 2)
{
point reflectedAa = refPt + ppDist*nB;
insertPoint(reflectedAa, reflectedMaster);
point reflectedBb = refPt + ppDist*nA;
insertPoint(reflectedBb, reflectedMaster);
}
// Slave is outside.
insertPoint(reflMasterPt, reflectedAI);
}
void Foam::conformalVoronoiMesh::insertFlatEdgePointGroup
(
const featureEdgeMesh& feMesh,
const pointIndexHit& edHit
)
{
Info<< " NOT INSERTING FLAT EDGE POINT GROUP, NOT IMPLEMENTED" << endl;
}
void Foam::conformalVoronoiMesh::insertOpenEdgePointGroup
(
const featureEdgeMesh& feMesh,
const pointIndexHit& edHit
)
{
Info<< " NOT INSERTING OPEN EDGE POINT GROUP, NOT IMPLEMENTED" << endl;
}
void Foam::conformalVoronoiMesh::insertMultipleEdgePointGroup
(
const featureEdgeMesh& feMesh,
const pointIndexHit& edHit
)
{
Info<< " NOT INSERTING MULTIPLE EDGE POINT GROUP, NOT IMPLEMENTED"
<< endl;
}
void Foam::conformalVoronoiMesh::conformToFeaturePoints()
{
Info<< nl << "Conforming to feature points" << endl;
@ -1020,9 +1203,13 @@ void Foam::conformalVoronoiMesh::conformToSurface()
startOfSurfacePointPairs_ = number_of_vertices();
Info<< " EDGE DISTANCE COEFFS HARD-CODED." << endl;
scalar edgeSearchDistCoeffSqr = sqr(1.25);
scalar surfacePtReplaceDistCoeffSqr = sqr(0.3);
// Surface protrusion conformation
label nIterations = 1;
label nIterations = 2;
for(label iterationNo = 0; iterationNo < nIterations; iterationNo++)
{
@ -1031,6 +1218,10 @@ void Foam::conformalVoronoiMesh::conformToSurface()
DynamicList<point> surfacePoints;
DynamicList<vector> surfaceNormals;
DynamicList<pointIndexHit> featureEdgeHits;
DynamicList<label> featureEdgeFeaturesHit;
DynamicList<point> tmpFeatureEdgePoints;
for
(
Triangulation::Finite_vertices_iterator vit =
@ -1043,18 +1234,18 @@ void Foam::conformalVoronoiMesh::conformToSurface()
{
point vert(topoint(vit->point()));
scalar searchDistanceSqr = surfaceSearchDistanceSqr(vert);
pointIndexHit pHit;
pointIndexHit surfHit;
vector normal;
geometryToConformTo_.findSurfaceNearestAndNormal
(
vert,
searchDistanceSqr,
pHit,
surfHit,
normal
);
if (pHit.hit())
if (surfHit.hit())
{
vit->setNearBoundary();
@ -1064,10 +1255,49 @@ void Foam::conformalVoronoiMesh::conformToSurface()
// edge, shift it to being an edge control point
// instead, this will prevent "pits" forming.
surfacePpDist.append(pointPairDistance(vert));
surfaceSearchDistSqr.append(searchDistanceSqr);
surfacePoints.append(pHit.hitPoint());
surfaceNormals.append(normal);
pointIndexHit edHit;
label featureHit;
scalar targetCellSizeSqr = sqr(targetCellSize(vert));
geometryToConformTo_.findEdgeNearest
(
surfHit.hitPoint(),
edgeSearchDistCoeffSqr*targetCellSizeSqr,
edHit,
featureHit
);
bool keepSurfacePoint = true;
if (edHit.hit())
{
// Note that edge type classification can be used
// here if necessary, i.e to avoid over-populating
// internal edges
featureEdgeHits.append(edHit);
featureEdgeFeaturesHit.append(featureHit);
tmpFeatureEdgePoints.append(edHit.hitPoint());
if
(
magSqr(edHit.hitPoint() - surfHit.hitPoint())
< surfacePtReplaceDistCoeffSqr*targetCellSizeSqr
)
{
keepSurfacePoint = false;
}
}
if (keepSurfacePoint)
{
surfacePpDist.append(pointPairDistance(vert));
surfaceSearchDistSqr.append(searchDistanceSqr);
surfacePoints.append(surfHit.hitPoint());
surfaceNormals.append(normal);
}
}
}
}
@ -1075,7 +1305,9 @@ void Foam::conformalVoronoiMesh::conformToSurface()
Info<< nl <<" iterationNo " << iterationNo << nl
<< " number_of_vertices " << number_of_vertices() << nl
<< " surfacePoints.size() " << surfacePoints.size() << endl;
<< " surfacePoints.size() " << surfacePoints.size() << nl
<< " featureEdgeHits.size() " << featureEdgeHits.size()
<< endl;
insertSurfacePointPairs
(
@ -1088,46 +1320,17 @@ void Foam::conformalVoronoiMesh::conformToSurface()
)
);
// Feature edge conformation.
geometryToConformTo_.findEdgeNearest
insertEdgePointGroups
(
pointField(surfacePoints),
scalarField(surfaceSearchDistSqr)
featureEdgeHits,
featureEdgeFeaturesHit
);
// After the surface conformation points are added, any points that are
// still protruding the surface may be protruding from edges, so
// identify the points and test if they are close to a feature edge
DynamicList<point> edgeGenerationPoints;
for
(
Triangulation::Finite_vertices_iterator vit =
finite_vertices_begin();
vit != finite_vertices_end();
vit++
)
{
if (vit->nearBoundary())
{
if (dualCellSurfaceIntersection(vit))
{
// Test to see if near to an edge, conform to the nearest
// point on that edge if so.
edgeGenerationPoints.append(topoint(vit->point()));
}
}
}
writePoints
(
"edgeGenerationLocations_" + name(iterationNo) + ".obj",
edgeGenerationPoints
tmpFeatureEdgePoints
);
}
}
@ -1137,6 +1340,8 @@ void Foam::conformalVoronoiMesh::move()
scalar relaxation = relaxationModel_->relaxation();
Info<< nl << " Relaxation = " << relaxation << endl;
timeCheck();
}

View File

@ -160,6 +160,50 @@ class conformalVoronoiMesh
const fileName fName = fileName::null
);
//- Insert groups of points to conform to an edge given a list of
// pointIndexHits specifying the location and edge index of the point
// to be conformed to on the corresponding entry in featureHit
void insertEdgePointGroups
(
const List<pointIndexHit>& edgeHits,
const labelList& featuresHit
);
//- Insert points to conform to an external edge
void insertExternalEdgePointGroup
(
const featureEdgeMesh& feMesh,
const pointIndexHit& edHit
);
//- Insert points to conform to an internal edge
void insertInternalEdgePointGroup
(
const featureEdgeMesh& feMesh,
const pointIndexHit& edHit
);
//- Insert points to conform to a flat edge
void insertFlatEdgePointGroup
(
const featureEdgeMesh& feMesh,
const pointIndexHit& edHit
);
//- Insert points to conform to an open edge
void insertOpenEdgePointGroup
(
const featureEdgeMesh& feMesh,
const pointIndexHit& edHit
);
//- Insert points to conform to multiply connected edge
void insertMultipleEdgePointGroup
(
const featureEdgeMesh& feMesh,
const pointIndexHit& edHit
);
//- Insert point groups at the feature points.
void conformToFeaturePoints();

View File

@ -370,7 +370,7 @@ void Foam::conformationSurfaces::findSurfaceNearestAndNormal
void Foam::conformationSurfaces::findSurfaceNearestAndNormal
(
const pointField& samples,
scalarField nearestDistSqr,
const scalarField& nearestDistSqr,
List<pointIndexHit>& hitInfo,
vectorField& normals
) const
@ -416,21 +416,46 @@ void Foam::conformationSurfaces::findSurfaceNearestAndNormal
void Foam::conformationSurfaces::findEdgeNearest
(
const pointField& samples,
const scalarField& nearestDistSqr
const point& sample,
scalar nearestDistSqr,
pointIndexHit& edgeHit,
label& featureHit
) const
{
pointField samples(1, sample);
scalarField nearestDistsSqr(1, nearestDistSqr);
labelList nearestSurfaces;
List<pointIndexHit> nearestInfo;
List<pointIndexHit> edgeHits;
labelList featuresHit;
findEdgeNearest
(
samples,
nearestDistsSqr,
edgeHits,
featuresHit
);
edgeHit = edgeHits[0];
featureHit = featuresHit[0];
}
void Foam::conformationSurfaces::findEdgeNearest
(
const pointField& samples,
const scalarField& nearestDistsSqr,
List<pointIndexHit>& edgeHits,
labelList& featuresHit
) const
{
// Initialise
nearestSurfaces.setSize(samples.size());
nearestSurfaces = -1;
nearestInfo.setSize(samples.size());
featuresHit.setSize(samples.size());
featuresHit = -1;
edgeHits.setSize(samples.size());
// Work arrays
scalarField minDistSqr(nearestDistSqr);
scalarField minDistSqr(nearestDistsSqr);
List<pointIndexHit> hitInfo(samples.size());
forAll(features_, testI)
@ -452,34 +477,11 @@ void Foam::conformationSurfaces::findEdgeNearest
hitInfo[pointI].hitPoint()
- samples[pointI]
);
nearestInfo[pointI] = hitInfo[pointI];
nearestSurfaces[pointI] = testI;
edgeHits[pointI] = hitInfo[pointI];
featuresHit[pointI] = testI;
}
}
}
OFstream edgeNearestStr("testFindEdgeNearest.obj");
Pout<< "Writing edge nearest points to " << edgeNearestStr.name() << endl;
forAll(nearestInfo, i)
{
if (nearestInfo[i].hit())
{
meshTools::writeOBJ(edgeNearestStr, nearestInfo[i].hitPoint());
}
}
// labelList& edgeLabel;
// pointField& edgePoint;
// List<vectorField>& adjacentNormals;
// edgeLabel.setSize(samples.size());
// edgePoint.setSize(samples.size());
// adjacentNormals.setSize(samples.size());
// edgeLabel[i] = pHit.index();
// edgePoint[i] = pHit.hitPoint();
// adjacentNormals[i] = edgeNormals(edgeLabel[i]);
}

View File

@ -166,17 +166,28 @@ public:
void findSurfaceNearestAndNormal
(
const pointField& sample,
scalarField nearestDistSqr,
List<pointIndexHit>& pHit,
vectorField& normal
const pointField& samples,
const scalarField& nearestDistSqr,
List<pointIndexHit>& hitInfo,
vectorField& normals
) const;
//- Find the nearest point on any feature edge
void findEdgeNearest
(
const point& sample,
scalar nearestDistSqr,
pointIndexHit& edgeHit,
label& featureHit
) const;
void findEdgeNearest
(
const pointField& samples,
const scalarField& nearestDistSqr
const scalarField& nearestDistsSqr,
List<pointIndexHit>& edgeHits,
labelList& featuresHit
) const;
// Write

View File

@ -130,6 +130,9 @@ public:
//- Return the surfaceSearchDistanceCoeff
inline scalar surfaceSearchDistanceCoeff() const;
//- Return the maxQuadAngle
inline scalar maxQuadAngle() const;
//- Return the minimumEdgeLengthCoeff
inline scalar minimumEdgeLengthCoeff() const;

View File

@ -51,6 +51,13 @@ inline Foam::scalar Foam::cvControls::surfaceSearchDistanceCoeff() const
return surfaceSearchDistanceCoeff_;
}
inline Foam::scalar Foam::cvControls::maxQuadAngle() const
{
return maxQuadAngle_;
}
inline Foam::scalar Foam::cvControls::minimumEdgeLengthCoeff() const
{
return minimumEdgeLengthCoeff_;