zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: cvMesh: 3D CV Mesh Generator

Browse Source 
- Remove cvMesh from constructor arguments for faceAreaWeightModel and cvControls
    - Add I.H files for indexedCell and indexedVertex classes
    - Separate feature point functions; put into conformalVoronoiMeshFeaturePoints.C
    - Add feature point specialisations for 2 external and 1 internal edge
    - Add a struct for feature point types
    - Add a writeCellCentres function
This commit is contained in:
laurence
2011-12-13 17:50:22 +00:00
parent f29a53ecbd
commit f07935e6c6
20 changed files with 2028 additions and 1534 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/files
View File

@ -4,6 +4,7 @@ conformalVoronoiMesh/conformalVoronoiMesh.C
conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
conformalVoronoiMesh/conformalVoronoiMeshIO.C conformalVoronoiMesh/conformalVoronoiMeshIO.C
conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
cvControls/cvControls.C cvControls/cvControls.C

784
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
View File

@ -507,734 +507,6 @@ void Foam::conformalVoronoiMesh::insertEdgePointGroups
} }
void Foam::conformalVoronoiMesh::createEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
label edgeI = edHit.index();
extendedFeatureEdgeMesh::edgeStatus edStatus = feMesh.getEdgeStatus(edgeI);
switch (edStatus)
{
case extendedFeatureEdgeMesh::EXTERNAL:
{
createExternalEdgePointGroup(feMesh, edHit, pts, indices, types);
break;
}
case extendedFeatureEdgeMesh::INTERNAL:
{
createInternalEdgePointGroup(feMesh, edHit, pts, indices, types);
break;
}
case extendedFeatureEdgeMesh::FLAT:
{
createFlatEdgePointGroup(feMesh, edHit, pts, indices, types);
break;
}
case extendedFeatureEdgeMesh::OPEN:
{
createOpenEdgePointGroup(feMesh, edHit, pts, indices, types);
break;
}
case extendedFeatureEdgeMesh::MULTIPLE:
{
createMultipleEdgePointGroup(feMesh, edHit, pts, indices, types);
break;
}
case extendedFeatureEdgeMesh::NONE:
{
break;
}
}
}
void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
const Foam::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]];
if (mag(nA & nB) > 1 - SMALL)
{
// The normals are nearly parallel, so this is too sharp a feature to
// conform to.
return;
}
// Normalised distance of reference point from edge point
vector refVec((nA + nB)/(1 + (nA & nB)));
if (magSqr(refVec) > sqr(5.0))
{
// Limit the size of the conformation
ppDist *= 5.0/mag(refVec);
// Pout<< nl << "createExternalEdgePointGroup limit "
// << "edgePt " << edgePt << nl
// << "refVec " << refVec << nl
// << "mag(refVec) " << mag(refVec) << nl
// << "ppDist " << ppDist << nl
// << "nA " << nA << nl
// << "nB " << nB << nl
// << "(nA & nB) " << (nA & nB) << nl
// << endl;
}
// Convex. So refPt will be inside domain and hence a master point
Foam::point refPt = edgePt - ppDist*refVec;
// Result when the points are eventually inserted.
// Add number_of_vertices() at insertion of first vertex to all numbers:
// pt index type
// refPt 0 1
// reflectedA 1 0
// reflectedB 2 0
// Insert the master point pairing the the first slave
pts.append(refPt);
indices.append(0);
types.append(1);
// Insert the slave points by reflecting refPt in both faces.
// with each slave refering to the master
Foam::point reflectedA = refPt + 2*ppDist*nA;
pts.append(reflectedA);
indices.append(0);
types.append(-1);
Foam::point reflectedB = refPt + 2*ppDist*nB;
pts.append(reflectedB);
indices.append(0);
types.append(-2);
}
void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
const Foam::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]];
if (mag(nA & nB) > 1 - SMALL)
{
// The normals are nearly parallel, so this is too sharp a feature to
// conform to.
return;
}
// Normalised distance of reference point from edge point
vector refVec((nA + nB)/(1 + (nA & nB)));
if (magSqr(refVec) > sqr(5.0))
{
// Limit the size of the conformation
ppDist *= 5.0/mag(refVec);
// Pout<< nl << "createInternalEdgePointGroup limit "
// << "edgePt " << edgePt << nl
// << "refVec " << refVec << nl
// << "mag(refVec) " << mag(refVec) << nl
// << "ppDist " << ppDist << nl
// << "nA " << nA << nl
// << "nB " << nB << nl
// << "(nA & nB) " << (nA & nB) << nl
// << endl;
}
// Concave. master and reflected points inside the domain.
Foam::point refPt = edgePt - ppDist*refVec;
// Generate reflected master to be outside.
Foam::point reflMasterPt = refPt + 2*(edgePt - refPt);
// Reflect reflMasterPt in both faces.
Foam::point reflectedA = reflMasterPt - 2*ppDist*nA;
Foam::point reflectedB = reflMasterPt - 2*ppDist*nB;
scalar totalAngle =
radToDeg(constant::mathematical::pi + acos(nA & nB));
// 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 = 2 + nAddPoints;
// Add number_of_vertices() at insertion of first vertex to all numbers:
// Result for nAddPoints 1 when the points are eventually inserted
// pt index type
// reflectedA 0 2
// reflectedB 1 2
// reflMasterPt 2 0
// Result for nAddPoints 1 when the points are eventually inserted
// pt index type
// reflectedA 0 3
// reflectedB 1 3
// refPt 2 3
// reflMasterPt 3 0
// Result for nAddPoints 2 when the points are eventually inserted
// pt index type
// reflectedA 0 4
// reflectedB 1 4
// reflectedAa 2 4
// reflectedBb 3 4
// reflMasterPt 4 0
// Master A is inside.
pts.append(reflectedA);
indices.append(0);
types.append(reflectedMaster--);
// Master B is inside.
pts.append(reflectedB);
indices.append(0);
types.append(reflectedMaster--);
if (nAddPoints == 1)
{
// One additinal point is the reflection of the slave point,
// i.e. the original reference point
pts.append(refPt);
indices.append(0);
types.append(reflectedMaster--);
}
else if (nAddPoints == 2)
{
Foam::point reflectedAa = refPt + ppDist*nB;
pts.append(reflectedAa);
indices.append(0);
types.append(reflectedMaster--);
Foam::point reflectedBb = refPt + ppDist*nA;
pts.append(reflectedBb);
indices.append(0);
types.append(reflectedMaster--);
}
// Slave is outside.
pts.append(reflMasterPt);
indices.append(0);
types.append(-(nAddPoints + 2));
}
void Foam::conformalVoronoiMesh::createFlatEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
const Foam::point& edgePt = edHit.hitPoint();
scalar ppDist = pointPairDistance(edgePt);
const vectorField& feNormals = feMesh.normals();
const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
// As this is a flat edge, there are two normals by definition
const vector& nA = feNormals[edNormalIs[0]];
const vector& nB = feNormals[edNormalIs[1]];
// Average normal to remove any bias to one side, although as this
// is a flat edge, the normals should be essentially the same
vector n = 0.5*(nA + nB);
// Direction along the surface to the control point, sense of edge
// direction not important, as +s and -s can be used because this
// is a flat edge
vector s = 2.0*ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
createPointPair(ppDist, edgePt + s, n, pts, indices, types);
createPointPair(ppDist, edgePt - s, n, pts, indices, types);
}
void Foam::conformalVoronoiMesh::createOpenEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
Info<< "NOT INSERTING OPEN EDGE POINT GROUP, NOT IMPLEMENTED" << endl;
}
void Foam::conformalVoronoiMesh::createMultipleEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
Info<< "NOT INSERTING MULTIPLE EDGE POINT GROUP, NOT IMPLEMENTED"
<< endl;
}
void Foam::conformalVoronoiMesh::insertFeaturePoints()
{
Info<< nl << "Conforming to feature points" << endl;
DynamicList<Foam::point> pts;
DynamicList<label> indices;
DynamicList<label> types;
label preFeaturePointSize = number_of_vertices();
createConvexFeaturePoints(pts, indices, types);
createConcaveFeaturePoints(pts, indices, types);
createMixedFeaturePoints(pts, indices, types);
// Insert the created points, distributing to the appropriate processor
insertPoints(pts, indices, types, true);
if (cvMeshControls().objOutput())
{
writePoints("featureVertices.obj", pts);
}
label nFeatureVertices = number_of_vertices() - preFeaturePointSize;
if (Pstream::parRun())
{
reduce(nFeatureVertices, sumOp<label>());
}
if (nFeatureVertices > 0)
{
Info<< " Inserted " << nFeatureVertices
<< " feature vertices" << endl;
}
featureVertices_.clear();
forAll(pts, pI)
{
featureVertices_.push_back
(
Vb(toPoint(pts[pI]), indices[pI], types[pI])
);
}
constructFeaturePointLocations();
}
void Foam::conformalVoronoiMesh::createConvexFeaturePoints
(
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
const PtrList<extendedFeatureEdgeMesh>& feMeshes
(
geometryToConformTo_.features()
);
forAll(feMeshes, i)
{
const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
for
(
label ptI = feMesh.convexStart();
ptI < feMesh.concaveStart();
ptI++
)
{
const Foam::point& featPt = feMesh.points()[ptI];
if (!positionOnThisProc(featPt))
{
continue;
}
vectorField featPtNormals = feMesh.featurePointNormals(ptI);
vector cornerNormal = sum(featPtNormals);
cornerNormal /= mag(cornerNormal);
Foam::point internalPt =
featPt - pointPairDistance(featPt)*cornerNormal;
// Result when the points are eventually inserted (example n = 4)
// Add number_of_vertices() at insertion of first vertex to all
// numbers:
// pt index type
// internalPt 0 1
// externalPt0 1 0
// externalPt1 2 0
// externalPt2 3 0
// externalPt3 4 0
pts.append(internalPt);
indices.append(0);
types.append(1);
label internalPtIndex = -1;
forAll(featPtNormals, nI)
{
const vector& n = featPtNormals[nI];
plane planeN = plane(featPt, n);
Foam::point externalPt =
internalPt + 2.0 * planeN.distance(internalPt) * n;
pts.append(externalPt);
indices.append(0);
types.append(internalPtIndex--);
}
}
}
}
void Foam::conformalVoronoiMesh::createConcaveFeaturePoints
(
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
const PtrList<extendedFeatureEdgeMesh>& feMeshes
(
geometryToConformTo_.features()
);
forAll(feMeshes, i)
{
const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
for
(
label ptI = feMesh.concaveStart();
ptI < feMesh.mixedStart();
ptI++
)
{
const Foam::point& featPt = feMesh.points()[ptI];
if (!positionOnThisProc(featPt))
{
continue;
}
vectorField featPtNormals = feMesh.featurePointNormals(ptI);
vector cornerNormal = sum(featPtNormals);
cornerNormal /= mag(cornerNormal);
Foam::point externalPt =
featPt + pointPairDistance(featPt)*cornerNormal;
label externalPtIndex = featPtNormals.size();
// Result when the points are eventually inserted (example n = 5)
// Add number_of_vertices() at insertion of first vertex to all
// numbers:
// pt index type
// internalPt0 0 5
// internalPt1 1 5
// internalPt2 2 5
// internalPt3 3 5
// internalPt4 4 5
// externalPt 5 4
forAll(featPtNormals, nI)
{
const vector& n = featPtNormals[nI];
plane planeN = plane(featPt, n);
Foam::point internalPt =
externalPt - 2.0 * planeN.distance(externalPt) * n;
pts.append(internalPt);
indices.append(0);
types.append(externalPtIndex--);
}
pts.append(externalPt);
indices.append(0);
types.append(-1);
}
}
}
void Foam::conformalVoronoiMesh::createMixedFeaturePoints
(
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
if (cvMeshControls().mixedFeaturePointPPDistanceCoeff() < 0)
{
Info<<nl << "Skipping specialised handling for mixed feature points"
<< endl;
return;
}
const PtrList<extendedFeatureEdgeMesh>& feMeshes
(
geometryToConformTo_.features()
);
forAll(feMeshes, i)
{
const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
for
(
label ptI = feMesh.mixedStart();
ptI < feMesh.nonFeatureStart();
ptI++
)
{
if
(
!createSpecialisedFeaturePoint(feMesh, ptI, pts, indices, types)
)
{
// Specialisations available for some mixed feature points. For
// non-specialised feature points, inserting mixed internal and
// external edge groups at feature point.
const labelList& pEds(feMesh.pointEdges()[ptI]);
// Skipping unsupported mixed feature point types
bool skipEdge = false;
forAll(pEds, e)
{
label edgeI = pEds[e];
extendedFeatureEdgeMesh::edgeStatus edStatus =
feMesh.getEdgeStatus(edgeI);
if
(
edStatus == extendedFeatureEdgeMesh::OPEN
|| edStatus == extendedFeatureEdgeMesh::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 Foam::point& pt(feMesh.points()[ptI]);
if (!positionOnThisProc(pt))
{
continue;
}
scalar edgeGroupDistance = mixedFeaturePointDistance(pt);
forAll(pEds, e)
{
label edgeI = pEds[e];
Foam::point edgePt =
pt + edgeGroupDistance*feMesh.edgeDirection(edgeI, ptI);
pointIndexHit edgeHit(true, edgePt, edgeI);
createEdgePointGroup(feMesh, edgeHit, pts, indices, types);
}
}
}
}
}
void Foam::conformalVoronoiMesh::constructFeaturePointLocations()
{
DynamicList<Foam::point> ftPtLocs;
const PtrList<extendedFeatureEdgeMesh>& feMeshes
(
geometryToConformTo_.features()
);
forAll(feMeshes, i)
{
const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
if (cvMeshControls().mixedFeaturePointPPDistanceCoeff() < 0)
{
// Ignoring mixed feature points
for
(
label ptI = feMesh.convexStart();
ptI < feMesh.mixedStart();
ptI++
)
{
ftPtLocs.append(feMesh.points()[ptI]);
}
}
else
{
for
(
label ptI = feMesh.convexStart();
ptI < feMesh.nonFeatureStart();
ptI++
)
{
ftPtLocs.append(feMesh.points()[ptI]);
}
}
}
featurePointLocations_.transfer(ftPtLocs);
}
void Foam::conformalVoronoiMesh::reinsertFeaturePoints(bool distribute)
{
Info<< nl << "Reinserting stored feature points" << endl;
label preReinsertionSize(number_of_vertices());
if (distribute)
{
DynamicList<Foam::point> pointsToInsert;
DynamicList<label> indices;
DynamicList<label> types;
for
(
std::list<Vb>::iterator vit=featureVertices_.begin();
vit != featureVertices_.end();
++vit
)
{
pointsToInsert.append(topoint(vit->point()));
indices.append(vit->index());
types.append(vit->type());
}
// Insert distributed points
insertPoints(pointsToInsert, indices, types, true);
// Save points in new distribution
featureVertices_.clear();
forAll(pointsToInsert, pI)
{
featureVertices_.push_back
(
Vb(toPoint(pointsToInsert[pI]), indices[pI], types[pI])
);
}
}
else
{
for
(
std::list<Vb>::iterator vit=featureVertices_.begin();
vit != featureVertices_.end();
++vit
)
{
// Assuming that all of the reinsertions are pair points, and that
// the index and type are relative, i.e. index 0 and type relative
// to it.
insertPoint
(
vit->point(),
vit->index() + number_of_vertices(),
vit->type() + number_of_vertices()
);
}
}
Info<< " Reinserted "
<< returnReduce
(
label(number_of_vertices()) - preReinsertionSize,
sumOp<label>()
)
<< " vertices" << endl;
}
const Foam::indexedOctree<Foam::treeDataPoint>& const Foam::indexedOctree<Foam::treeDataPoint>&
Foam::conformalVoronoiMesh::featurePointTree() const Foam::conformalVoronoiMesh::featurePointTree() const
{ {
@ -1277,33 +549,17 @@ bool Foam::conformalVoronoiMesh::nearFeaturePt(const Foam::point& pt) const
} }
void Foam::conformalVoronoiMesh::reset() void Foam::conformalVoronoiMesh::reset(const bool distribute)
{ {
this->clear(); this->clear();
insertBoundingPoints(); reinsertBoundingPoints();
}
// Reinsert feature points, distributing them as necessary.
reinsertFeaturePoints(distribute);
//insertFeaturePoints();
void Foam::conformalVoronoiMesh::insertBoundingPoints() startOfInternalPoints_ = number_of_vertices();
{
pointField farPts = geometryToConformTo_.globalBounds().points();
// Shift corners of bounds relative to origin
farPts -= geometryToConformTo_.globalBounds().midpoint();
// Scale the box up
farPts *= 10.0;
// Shift corners of bounds back to be relative to midpoint
farPts += geometryToConformTo_.globalBounds().midpoint();
limitBounds_ = treeBoundBox(farPts);
forAll(farPts, fPI)
{
insertPoint(farPts[fPI], Vb::vtFar);
}
} }
@ -1466,13 +722,8 @@ bool Foam::conformalVoronoiMesh::distributeBackground()
cellVertexTypes cellVertexTypes
); );
// Remove the entire tessellation // Reset the entire tessellation
reset(); reset(true);
// Reinsert feature points, distributing them as necessary.
reinsertFeaturePoints(true);
startOfInternalPoints_ = number_of_vertices();
timeCheck("Distribution performed"); timeCheck("Distribution performed");
@ -1909,7 +1160,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
Delaunay(), Delaunay(),
runTime_(runTime), runTime_(runTime),
rndGen_(64293*Pstream::myProcNo()), rndGen_(64293*Pstream::myProcNo()),
cvMeshControls_(*this, cvMeshDict), cvMeshControls_(cvMeshDict),
allGeometry_ allGeometry_
( (
IOobject IOobject
@ -1965,8 +1216,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
( (
faceAreaWeightModel::New faceAreaWeightModel::New
( (
cvMeshDict.subDict("motionControl"), cvMeshDict.subDict("motionControl")
*this
) )
), ),
decomposition_() decomposition_()
@ -2003,6 +1253,8 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
buildSurfaceConformation(rmCoarse); buildSurfaceConformation(rmCoarse);
is_valid(true);
// The introduction of the surface conformation may have distorted the // The introduction of the surface conformation may have distorted the
// balance of vertices, distribute if necessary. // balance of vertices, distribute if necessary.
if (distributeBackground()) if (distributeBackground())
@ -2386,9 +1638,10 @@ void Foam::conformalVoronoiMesh::move()
// Remove the entire tessellation // Remove the entire tessellation
reset(); reset();
reinsertFeaturePoints(); if (cvMeshControls().objOutput() && runTime_.outputTime())
{
startOfInternalPoints_ = number_of_vertices(); writePoints("featurePoints_" + runTime_.timeName() + ".obj", false);
}
timeCheck("Displacement calculated"); timeCheck("Displacement calculated");
@ -2400,6 +1653,11 @@ void Foam::conformalVoronoiMesh::move()
conformToSurface(); conformToSurface();
if (cvMeshControls().objOutput() && runTime_.outputTime())
{
writePoints("points_" + runTime_.timeName() + ".obj", false);
}
timeCheck("After conformToSurface"); timeCheck("After conformToSurface");
updateSizesAndAlignments(pointsToInsert); updateSizesAndAlignments(pointsToInsert);

31
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
View File

@ -35,6 +35,10 @@ SourceFiles
conformalVoronoiMeshI.H conformalVoronoiMeshI.H
conformalVoronoiMesh.C conformalVoronoiMesh.C
conformalVoronoiMeshIO.C conformalVoronoiMeshIO.C
conformalVoronoiMeshConformToSurface.C
conformalVoronoiMeshFeaturePoints.C
conformalVoronoiMeshFeaturePointSpecialisations.C
conformalVoronoiMeshCalcDualMesh.C
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
@ -73,6 +77,7 @@ SourceFiles
#include "processorPolyPatch.H" #include "processorPolyPatch.H"
#include "zeroGradientFvPatchFields.H" #include "zeroGradientFvPatchFields.H"
#include "globalIndex.H" #include "globalIndex.H"
#include "pointFeatureEdgesTypes.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -420,12 +425,24 @@ private:
DynamicList<label>& types DynamicList<label>& types
); );
//- Fill the pointFeatureEdgesType struct with the types of feature
// edges that are attached to the point.
List<extendedFeatureEdgeMesh::edgeStatus> calcPointFeatureEdgesTypes
(
const extendedFeatureEdgeMesh& feMesh,
const labelList& pEds,
pointFeatureEdgesTypes& pFEdgeTypes
);
//- Create feature point groups if a specialisation exists for the //- Create feature point groups if a specialisation exists for the
// structure // structure
bool createSpecialisedFeaturePoint bool createSpecialisedFeaturePoint
( (
const extendedFeatureEdgeMesh& feMesh, const extendedFeatureEdgeMesh& feMesh,
label ptI, const labelList& pEds,
const pointFeatureEdgesTypes& pFEdgeTypes,
const List<extendedFeatureEdgeMesh::edgeStatus>& allEdStat,
const label ptI,
DynamicList<Foam::point>& pts, DynamicList<Foam::point>& pts,
DynamicList<label>& indices, DynamicList<label>& indices,
DynamicList<label>& types DynamicList<label>& types
@ -443,13 +460,18 @@ private:
//- Check if a location is in exclusion range around a feature point //- Check if a location is in exclusion range around a feature point
bool nearFeaturePt(const Foam::point& pt) const; bool nearFeaturePt(const Foam::point& pt) const;
//- Clear the entire tesselation and insert bounding points //- Clear the entire tesselation
void reset(); // Reinsert bounding points, feature points and recalculate
// startOfInternalPoints_
void reset(const bool distribute = false);
//- Insert far points in a large bounding box to avoid dual edges //- Insert far points in a large bounding box to avoid dual edges
// spanning huge distances // spanning huge distances
void insertBoundingPoints(); void insertBoundingPoints();
//- Reinsert the bounding points
void reinsertBoundingPoints();
//- Insert the initial points into the triangulation, based on the //- Insert the initial points into the triangulation, based on the
// initialPointsMethod // initialPointsMethod
void insertInitialPoints(); void insertInitialPoints();
@ -1045,6 +1067,9 @@ public:
// actual cell size (cbrt(actual cell volume)). // actual cell size (cbrt(actual cell volume)).
void writeCellSizes(const fvMesh& mesh) const; void writeCellSizes(const fvMesh& mesh) const;
//- Calculate and write the cell centres.
void writeCellCentres(const fvMesh& mesh) const;
//- Find the cellSet of the boundary cells which have points that //- Find the cellSet of the boundary cells which have points that
// protrude out of the surface beyond a tolerance. // protrude out of the surface beyond a tolerance.
void findRemainingProtrusionSet(const fvMesh& mesh) const; void findRemainingProtrusionSet(const fvMesh& mesh) const;

11
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
View File

@ -203,10 +203,13 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
if (surfHit.hit()) if (surfHit.hit())
{ {
vit->setNearBoundary();
if (dualCellSurfaceAnyIntersection(vit)) if (dualCellSurfaceAnyIntersection(vit))
{ {
// This used to be just before this if statement.
// Moved because a point is only near the boundary if
// the dual cell intersects the surface.
vit->setNearBoundary();
// meshTools::writeOBJ(Pout, vert); // meshTools::writeOBJ(Pout, vert);
// meshTools::writeOBJ(Pout, surfHit.hitPoint()); // meshTools::writeOBJ(Pout, surfHit.hitPoint());
// Pout<< "l cr0 cr1" << endl; // Pout<< "l cr0 cr1" << endl;
@ -497,8 +500,8 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
if (Pstream::parRun()) if (Pstream::parRun())
{ {
Foam::point a = dE0; Foam::point& a = dE0;
Foam::point b = dE1; Foam::point& b = dE1;
bool inProc = clipLineToProc(a, b); bool inProc = clipLineToProc(a, b);

573
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
View File

@ -27,30 +27,14 @@ License
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint Foam::List<Foam::extendedFeatureEdgeMesh::edgeStatus>
Foam::conformalVoronoiMesh::calcPointFeatureEdgesTypes
( (
const extendedFeatureEdgeMesh& feMesh, const extendedFeatureEdgeMesh& feMesh,
label ptI, const labelList& pEds,
DynamicList<Foam::point>& pts, pointFeatureEdgesTypes& pFEdgeTypes
DynamicList<label>& indices,
DynamicList<label>& types
) )
{ {
const labelList& pEds(feMesh.pointEdges()[ptI]);
if (pEds.size() != 3)
{
// Only three edge specialisations available
return false;
}
label nExternal = 0;
label nInternal = 0;
label nFlat = 0;
label nOpen = 0;
label nMultiple = 0;
List<extendedFeatureEdgeMesh::edgeStatus> allEdStat(pEds.size()); List<extendedFeatureEdgeMesh::edgeStatus> allEdStat(pEds.size());
forAll(pEds, i) forAll(pEds, i)
@ -65,293 +49,322 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
{ {
case extendedFeatureEdgeMesh::EXTERNAL: case extendedFeatureEdgeMesh::EXTERNAL:
{ {
nExternal++; pFEdgeTypes.nExternal++;
break; break;
} }
case extendedFeatureEdgeMesh::INTERNAL: case extendedFeatureEdgeMesh::INTERNAL:
{ {
nInternal++; pFEdgeTypes.nInternal++;
break; break;
} }
case extendedFeatureEdgeMesh::FLAT: case extendedFeatureEdgeMesh::FLAT:
{ {
nFlat++; pFEdgeTypes.nFlat++;
break; break;
} }
case extendedFeatureEdgeMesh::OPEN: case extendedFeatureEdgeMesh::OPEN:
{ {
nOpen++; pFEdgeTypes.nOpen++;
break; break;
} }
case extendedFeatureEdgeMesh::MULTIPLE: case extendedFeatureEdgeMesh::MULTIPLE:
{ {
nMultiple++; pFEdgeTypes.nMultiple++;
break; break;
} }
case extendedFeatureEdgeMesh::NONE: case extendedFeatureEdgeMesh::NONE:
{ {
pFEdgeTypes.nNonFeature++;
break; break;
} }
} }
} }
if (nExternal == 2 && nInternal == 1) if (debug)
{ {
// if (!positionOnThisProc(pt)) Info<< pFEdgeTypes << endl;
// {
// continue;
// }
// Info<< "nExternal == 2 && nInternal == 1" << endl;
// const Foam::point& featPt = feMesh.points()[ptI];
// scalar ppDist = pointPairDistance(featPt);
// const vectorField& normals = feMesh.normals();
// const labelListList& edgeNormals = feMesh.edgeNormals();
// label concaveEdgeI = pEds
// [
// findIndex(allEdStat, extendedFeatureEdgeMesh::INTERNAL)
// ];
// // // Find which planes are joined to the concave edge
// // List<label> concaveEdgePlanes(2,label(-1));
// // label concaveEdgeI = concaveEdges[0];
// // // Pick up the two faces adjacent to the concave feature edge
// // const labelList& eFaces = qSurf_.edgeFaces()[concaveEdgeI];
// // label faceA = eFaces[0];
// // vector nA = qSurf_.faceNormals()[faceA];
// // scalar maxNormalDotProduct = -SMALL;
// // forAll(uniquePlaneNormals, uPN)
// // {
// // scalar normalDotProduct = nA & uniquePlaneNormals[uPN];
// // if (normalDotProduct > maxNormalDotProduct)
// // {
// // maxNormalDotProduct = normalDotProduct;
// // concaveEdgePlanes[0] = uPN;
// // }
// // }
// // label faceB = eFaces[1];
// // vector nB = qSurf_.faceNormals()[faceB];
// // maxNormalDotProduct = -SMALL;
// // forAll(uniquePlaneNormals, uPN)
// // {
// // scalar normalDotProduct = nB & uniquePlaneNormals[uPN];
// // if (normalDotProduct > maxNormalDotProduct)
// // {
// // maxNormalDotProduct = normalDotProduct;
// // concaveEdgePlanes[1] = uPN;
// // }
// // }
// // const vector& concaveEdgePlaneANormal =
// // uniquePlaneNormals[concaveEdgePlanes[0]];
// // const vector& concaveEdgePlaneBNormal =
// // uniquePlaneNormals[concaveEdgePlanes[1]];
// // label convexEdgesPlaneI;
// // if (findIndex(concaveEdgePlanes, 0) == -1)
// // {
// // convexEdgesPlaneI = 0;
// // }
// // else if (findIndex(concaveEdgePlanes, 1) == -1)
// // {
// // convexEdgesPlaneI = 1;
// // }
// // else
// // {
// // convexEdgesPlaneI = 2;
// // }
// // const vector& convexEdgesPlaneNormal =
// // uniquePlaneNormals[convexEdgesPlaneI];
// // const edge& concaveEdge = edges[concaveEdgeI];
// // // Check direction of edge, if the feature point is at the end()
// // // the reverse direction.
// // scalar edgeDirection = 1.0;
// // if (ptI == concaveEdge.end())
// // {
// // edgeDirection *= -1.0;
// // }
// const vector& concaveEdgePlaneANormal =
// normals[edgeNormals[concaveEdgeI][0]];
// const vector& concaveEdgePlaneBNormal =
// normals[edgeNormals[concaveEdgeI][1]];
// // Intersect planes parallel to the concave edge planes offset
// // by ppDist and the plane defined by featPt and the edge
// // vector.
// plane planeA
// (
// featPt + ppDist*concaveEdgePlaneANormal,
// concaveEdgePlaneANormal
// );
// plane planeB
// (
// featPt + ppDist*concaveEdgePlaneBNormal,
// concaveEdgePlaneBNormal
// );
// const vector& concaveEdgeDir = feMesh.edgeDirection
// (
// concaveEdgeI,
// ptI
// );
// Foam::point concaveEdgeLocalFeatPt = featPt + ppDist*concaveEdgeDir;
// // Finding the nearest point on the intersecting line to the
// // edge point. Floating point errors often encountered using
// // planePlaneIntersect
// plane planeF(concaveEdgeLocalFeatPt, concaveEdgeDir);
// Foam::point concaveEdgeExternalPt = planeF.planePlaneIntersect
// (
// planeA,
// planeB
// );
// label concaveEdgeExternalPtI = number_of_vertices() + 4;
// // Redefine planes to be on the feature surfaces to project
// // through
// planeA = plane(featPt, concaveEdgePlaneANormal);
// planeB = plane(featPt, concaveEdgePlaneBNormal);
// Foam::point internalPtA =
// concaveEdgeExternalPt
// - 2*planeA.distance(concaveEdgeExternalPt)
// *concaveEdgePlaneANormal;
// label internalPtAI = insertPoint
// (
// internalPtA,
// concaveEdgeExternalPtI
// );
// Foam::point internalPtB =
// concaveEdgeExternalPt
// - 2*planeB.distance(concaveEdgeExternalPt)
// *concaveEdgePlaneBNormal;
// label internalPtBI = insertPoint
// (
// internalPtB,
// concaveEdgeExternalPtI
// );
// // TEMPORARY VARIABLE TO TEST
// vector convexEdgesPlaneNormal = -concaveEdgeDir;
// plane planeC(featPt, convexEdgesPlaneNormal);
// Foam::point externalPtD =
// internalPtA
// + 2*planeC.distance(internalPtA)*convexEdgesPlaneNormal;
// insertPoint(externalPtD, internalPtAI);
// Foam::point externalPtE =
// internalPtB
// + 2*planeC.distance(internalPtB)*convexEdgesPlaneNormal;
// insertPoint(externalPtE, internalPtBI);
// insertPoint(concaveEdgeExternalPt, internalPtAI);
// Info<< nl << "# featPt " << endl;
// meshTools::writeOBJ(Info, featPt);
// Info<< "# internalPtA" << endl;
// meshTools::writeOBJ(Info, internalPtA);
// Info<< "# internalPtB" << endl;
// meshTools::writeOBJ(Info, internalPtB);
// Info<< "# externalPtD" << endl;
// meshTools::writeOBJ(Info, externalPtD);
// Info<< "# externalPtE" << endl;
// meshTools::writeOBJ(Info, externalPtE);
// scalar totalAngle = radToDeg
// (
// constant::mathematical::pi +
// acos(mag(concaveEdgePlaneANormal & concaveEdgePlaneBNormal))
// );
// if (totalAngle > cvMeshControls().maxQuadAngle())
// {
// // Add additional mitering points
// scalar angleSign = 1.0;
// if
// (
// geometryToConformTo_.outside
// (
// featPt - convexEdgesPlaneNormal*ppDist
// )
// )
// {
// angleSign = -1.0;
// }
// scalar phi =
// angleSign*acos(concaveEdgeDir & -convexEdgesPlaneNormal);
// scalar guard =
// (
// 1 + sin(phi)*ppDist/mag
// (
// concaveEdgeLocalFeatPt - concaveEdgeExternalPt
// )
// )/cos(phi) - 1;
// Foam::point internalPtF =
// concaveEdgeExternalPt
// + (2 + guard)*(concaveEdgeLocalFeatPt - concaveEdgeExternalPt);
// label internalPtFI =
// insertPoint(internalPtF, number_of_vertices() + 1);
// Foam::point externalPtG =
// internalPtF
// + 2*planeC.distance(internalPtF) * convexEdgesPlaneNormal;
// insertPoint(externalPtG, internalPtFI);
// }
// return true;
return false;
} }
else if (nExternal == 1 && nInternal == 2)
return allEdStat;
}
bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
(
const extendedFeatureEdgeMesh& feMesh,
const labelList& pEds,
const pointFeatureEdgesTypes& pFEdgesTypes,
const List<extendedFeatureEdgeMesh::edgeStatus>& allEdStat,
const label ptI,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
if
(
pFEdgesTypes.nExternal == 2
&& pFEdgesTypes.nInternal == 1
&& pEds.size() == 3
)
{
Info<< "nExternal == 2 && nInternal == 1" << endl;
const Foam::point& featPt = feMesh.points()[ptI];
if (!positionOnThisProc(featPt))
{
return false;
}
const scalar ppDist = pointPairDistance(featPt);
const vectorField& normals = feMesh.normals();
const labelListList& edgeNormals = feMesh.edgeNormals();
label concaveEdgeI = -1;
labelList convexEdgesI(2, -1);
label nConvex = 0;
forAll(pEds, i)
{
const extendedFeatureEdgeMesh::edgeStatus& eS = allEdStat[i];
if (eS == extendedFeatureEdgeMesh::INTERNAL)
{
concaveEdgeI = pEds[i];
}
else if (eS == extendedFeatureEdgeMesh::EXTERNAL)
{
convexEdgesI[nConvex++] = pEds[i];
}
else if (eS == extendedFeatureEdgeMesh::FLAT)
{
WarningIn("Foam::conformalVoronoiMesh::"
"createSpecialisedFeaturePoint")
<< "Edge " << eS << " is flat"
<< endl;
}
else
{
FatalErrorIn("Foam::conformalVoronoiMesh::"
"createSpecialisedFeaturePoint")
<< "Edge " << eS << " not concave/convex"
<< exit(FatalError);
}
}
const vector& concaveEdgePlaneANormal =
normals[edgeNormals[concaveEdgeI][0]];
const vector& concaveEdgePlaneBNormal =
normals[edgeNormals[concaveEdgeI][1]];
// Intersect planes parallel to the concave edge planes offset
// by ppDist and the plane defined by featPt and the edge vector.
plane planeA
(
featPt + ppDist*concaveEdgePlaneANormal,
concaveEdgePlaneANormal
);
plane planeB
(
featPt + ppDist*concaveEdgePlaneBNormal,
concaveEdgePlaneBNormal
);
const vector& concaveEdgeDir = feMesh.edgeDirection
(
concaveEdgeI,
ptI
);
// Todo,needed later but want to get rid of this.
const Foam::point concaveEdgeLocalFeatPt = featPt + ppDist*concaveEdgeDir;
// Finding the nearest point on the intersecting line to the edge
// point. Floating point errors often occur using planePlaneIntersect
plane planeF(concaveEdgeLocalFeatPt, concaveEdgeDir);
const Foam::point concaveEdgeExternalPt = planeF.planePlaneIntersect
(
planeA,
planeB
);
// Redefine planes to be on the feature surfaces to project through
planeA = plane(featPt, concaveEdgePlaneANormal);
planeB = plane(featPt, concaveEdgePlaneBNormal);
const Foam::point internalPtA =
concaveEdgeExternalPt
- 2*planeA.distance(concaveEdgeExternalPt)
*concaveEdgePlaneANormal;
pts.append(internalPtA);
indices.append(0);
types.append(4);
const Foam::point internalPtB =
concaveEdgeExternalPt
- 2*planeB.distance(concaveEdgeExternalPt)
*concaveEdgePlaneBNormal;
pts.append(internalPtB);
indices.append(0);
types.append(3);
// Add the external points
Foam::point externalPtD;
Foam::point externalPtE;
vector convexEdgePlaneCNormal;
vector convexEdgePlaneDNormal;
const labelList& concaveEdgeNormals = edgeNormals[concaveEdgeI];
const labelList& convexEdgeANormals = edgeNormals[convexEdgesI[0]];
const labelList& convexEdgeBNormals = edgeNormals[convexEdgesI[1]];
forAll(concaveEdgeNormals, edgeNormalI)
{
forAll(convexEdgeANormals, edgeAnormalI)
{
const vector& concaveNormal
= normals[concaveEdgeNormals[edgeNormalI]];
const vector& convexNormal
= normals[convexEdgeANormals[edgeAnormalI]];
const scalar orientation = concaveNormal & convexNormal;
if (orientation <= 0.0)
{
convexEdgePlaneCNormal = convexNormal;
plane planeC(featPt, convexEdgePlaneCNormal);
externalPtD =
internalPtA
+ 2.0*planeC.distance(internalPtA)*convexEdgePlaneCNormal;
pts.append(externalPtD);
indices.append(0);
types.append(-2);
}
}
forAll(convexEdgeBNormals, edgeBnormalI)
{
const vector& concaveNormal
= normals[concaveEdgeNormals[edgeNormalI]];
const vector& convexNormal
= normals[convexEdgeBNormals[edgeBnormalI]];
const scalar orientation = concaveNormal & convexNormal;
if (orientation <= 0.0)
{
convexEdgePlaneDNormal = convexNormal;
plane planeD(featPt, convexEdgePlaneDNormal);
externalPtE =
internalPtB
+ 2.0*planeD.distance(internalPtB)*convexEdgePlaneDNormal;
pts.append(externalPtE);
indices.append(0);
types.append(-2);
}
}
}
pts.append(concaveEdgeExternalPt);
indices.append(0);
types.append(-4);
const scalar totalAngle = radToDeg
(
constant::mathematical::pi +
acos(mag(concaveEdgePlaneANormal & concaveEdgePlaneBNormal))
);
if (totalAngle > cvMeshControls().maxQuadAngle())
{
// Add additional mitering points
//scalar angleSign = 1.0;
vector convexEdgesPlaneNormal = convexEdgePlaneCNormal;
plane planeM(featPt, convexEdgesPlaneNormal);
// if
// (
// geometryToConformTo_.outside
// (
// featPt - convexEdgesPlaneNormal*ppDist
// )
// )
// {
// angleSign = -1.0;
// }
// scalar phi =
// angleSign*acos(concaveEdgeDir & -convexEdgesPlaneNormal);
//
// scalar guard =
// (
// 1.0 + sin(phi)*ppDist/mag
// (
// concaveEdgeLocalFeatPt - concaveEdgeExternalPt
// )
// )/cos(phi) - 1.0;
const Foam::point internalPtF =
concaveEdgeExternalPt
//+ (2.0 + guard)*(concaveEdgeLocalFeatPt - concaveEdgeExternalPt);
+ 2.0*(concaveEdgeLocalFeatPt - concaveEdgeExternalPt);
pts.append(internalPtF);
indices.append(0);
types.append(1);
const Foam::point externalPtG =
internalPtF
//+ 2.0*planeM.distance(internalPtF) * convexEdgesPlaneNormal;
+ 2.0*planeM.distance(internalPtF)*convexEdgesPlaneNormal;
pts.append(externalPtG);
indices.append(0);
types.append(-1);
if (debug)
{
Info<< nl << "# featPt " << endl;
meshTools::writeOBJ(Info, featPt);
Info<< "# internalPtA" << endl;
meshTools::writeOBJ(Info, internalPtA);
Info<< "# internalPtB" << endl;
meshTools::writeOBJ(Info, internalPtB);
Info<< "# externalPtD" << endl;
meshTools::writeOBJ(Info, externalPtD);
Info<< "# externalPtE" << endl;
meshTools::writeOBJ(Info, externalPtE);
Info<< "# internalPtF" << endl;
meshTools::writeOBJ(Info, internalPtF);
Info<< "# externalPtG" << endl;
meshTools::writeOBJ(Info, externalPtG);
}
}
return true;
}
else if (pFEdgesTypes.nExternal == 1 && pFEdgesTypes.nInternal == 2)
{ {
// Info<< "nExternal == 1 && nInternal == 2" << endl; // Info<< "nExternal == 1 && nInternal == 2" << endl;

801
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C Normal file
View File

@ -0,0 +1,801 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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 "conformalVoronoiMesh.H"
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
void Foam::conformalVoronoiMesh::insertBoundingPoints()
{
pointField farPts = geometryToConformTo_.globalBounds().points();
// Shift corners of bounds relative to origin
farPts -= geometryToConformTo_.globalBounds().midpoint();
// Scale the box up
farPts *= 10.0;
// Shift corners of bounds back to be relative to midpoint
farPts += geometryToConformTo_.globalBounds().midpoint();
limitBounds_ = treeBoundBox(farPts);
forAll(farPts, fPI)
{
insertPoint(farPts[fPI], Vb::vtFar);
}
}
void Foam::conformalVoronoiMesh::createEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
label edgeI = edHit.index();
extendedFeatureEdgeMesh::edgeStatus edStatus = feMesh.getEdgeStatus(edgeI);
switch (edStatus)
{
case extendedFeatureEdgeMesh::EXTERNAL:
{
createExternalEdgePointGroup(feMesh, edHit, pts, indices, types);
break;
}
case extendedFeatureEdgeMesh::INTERNAL:
{
createInternalEdgePointGroup(feMesh, edHit, pts, indices, types);
break;
}
case extendedFeatureEdgeMesh::FLAT:
{
createFlatEdgePointGroup(feMesh, edHit, pts, indices, types);
break;
}
case extendedFeatureEdgeMesh::OPEN:
{
createOpenEdgePointGroup(feMesh, edHit, pts, indices, types);
break;
}
case extendedFeatureEdgeMesh::MULTIPLE:
{
createMultipleEdgePointGroup(feMesh, edHit, pts, indices, types);
break;
}
case extendedFeatureEdgeMesh::NONE:
{
break;
}
}
}
void Foam::conformalVoronoiMesh::createExternalEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
const Foam::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]];
if (mag(nA & nB) > 1 - SMALL)
{
// The normals are nearly parallel, so this is too sharp a feature to
// conform to.
return;
}
// Normalised distance of reference point from edge point
vector refVec((nA + nB)/(1 + (nA & nB)));
if (magSqr(refVec) > sqr(5.0))
{
// Limit the size of the conformation
ppDist *= 5.0/mag(refVec);
// Pout<< nl << "createExternalEdgePointGroup limit "
// << "edgePt " << edgePt << nl
// << "refVec " << refVec << nl
// << "mag(refVec) " << mag(refVec) << nl
// << "ppDist " << ppDist << nl
// << "nA " << nA << nl
// << "nB " << nB << nl
// << "(nA & nB) " << (nA & nB) << nl
// << endl;
}
// Convex. So refPt will be inside domain and hence a master point
Foam::point refPt = edgePt - ppDist*refVec;
// Result when the points are eventually inserted.
// Add number_of_vertices() at insertion of first vertex to all numbers:
// pt index type
// refPt 0 1
// reflectedA 1 0
// reflectedB 2 0
// Insert the master point pairing the the first slave
pts.append(refPt);
indices.append(0);
types.append(1);
// Insert the slave points by reflecting refPt in both faces.
// with each slave refering to the master
Foam::point reflectedA = refPt + 2*ppDist*nA;
pts.append(reflectedA);
indices.append(0);
types.append(-1);
Foam::point reflectedB = refPt + 2*ppDist*nB;
pts.append(reflectedB);
indices.append(0);
types.append(-2);
}
void Foam::conformalVoronoiMesh::createInternalEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
const Foam::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]];
if (mag(nA & nB) > 1 - SMALL)
{
// The normals are nearly parallel, so this is too sharp a feature to
// conform to.
return;
}
// Normalised distance of reference point from edge point
vector refVec((nA + nB)/(1 + (nA & nB)));
if (magSqr(refVec) > sqr(5.0))
{
// Limit the size of the conformation
ppDist *= 5.0/mag(refVec);
// Pout<< nl << "createInternalEdgePointGroup limit "
// << "edgePt " << edgePt << nl
// << "refVec " << refVec << nl
// << "mag(refVec) " << mag(refVec) << nl
// << "ppDist " << ppDist << nl
// << "nA " << nA << nl
// << "nB " << nB << nl
// << "(nA & nB) " << (nA & nB) << nl
// << endl;
}
// Concave. master and reflected points inside the domain.
Foam::point refPt = edgePt - ppDist*refVec;
// Generate reflected master to be outside.
Foam::point reflMasterPt = refPt + 2*(edgePt - refPt);
// Reflect reflMasterPt in both faces.
Foam::point reflectedA = reflMasterPt - 2*ppDist*nA;
Foam::point reflectedB = reflMasterPt - 2*ppDist*nB;
scalar totalAngle =
radToDeg(constant::mathematical::pi + acos(nA & nB));
// 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 = 2 + nAddPoints;
// Add number_of_vertices() at insertion of first vertex to all numbers:
// Result for nAddPoints 1 when the points are eventually inserted
// pt index type
// reflectedA 0 2
// reflectedB 1 2
// reflMasterPt 2 0
// Result for nAddPoints 1 when the points are eventually inserted
// pt index type
// reflectedA 0 3
// reflectedB 1 3
// refPt 2 3
// reflMasterPt 3 0
// Result for nAddPoints 2 when the points are eventually inserted
// pt index type
// reflectedA 0 4
// reflectedB 1 4
// reflectedAa 2 4
// reflectedBb 3 4
// reflMasterPt 4 0
// Master A is inside.
pts.append(reflectedA);
indices.append(0);
types.append(reflectedMaster--);
// Master B is inside.
pts.append(reflectedB);
indices.append(0);
types.append(reflectedMaster--);
if (nAddPoints == 1)
{
// One additinal point is the reflection of the slave point,
// i.e. the original reference point
pts.append(refPt);
indices.append(0);
types.append(reflectedMaster--);
}
else if (nAddPoints == 2)
{
Foam::point reflectedAa = refPt + ppDist*nB;
pts.append(reflectedAa);
indices.append(0);
types.append(reflectedMaster--);
Foam::point reflectedBb = refPt + ppDist*nA;
pts.append(reflectedBb);
indices.append(0);
types.append(reflectedMaster--);
}
// Slave is outside.
pts.append(reflMasterPt);
indices.append(0);
types.append(-(nAddPoints + 2));
}
void Foam::conformalVoronoiMesh::createFlatEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
const Foam::point& edgePt = edHit.hitPoint();
const scalar ppDist = pointPairDistance(edgePt);
const vectorField& feNormals = feMesh.normals();
const labelList& edNormalIs = feMesh.edgeNormals()[edHit.index()];
// As this is a flat edge, there are two normals by definition
const vector& nA = feNormals[edNormalIs[0]];
const vector& nB = feNormals[edNormalIs[1]];
// Average normal to remove any bias to one side, although as this
// is a flat edge, the normals should be essentially the same
vector n = 0.5*(nA + nB);
// Direction along the surface to the control point, sense of edge
// direction not important, as +s and -s can be used because this
// is a flat edge
vector s = ppDist*(feMesh.edgeDirections()[edHit.index()] ^ n);
createPointPair(ppDist, edgePt + s, n, pts, indices, types);
createPointPair(ppDist, edgePt - s, n, pts, indices, types);
}
void Foam::conformalVoronoiMesh::createOpenEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
Info<< "NOT INSERTING OPEN EDGE POINT GROUP, NOT IMPLEMENTED" << endl;
}
void Foam::conformalVoronoiMesh::createMultipleEdgePointGroup
(
const extendedFeatureEdgeMesh& feMesh,
const pointIndexHit& edHit,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
Info<< "NOT INSERTING MULTIPLE EDGE POINT GROUP, NOT IMPLEMENTED" << endl;
}
void Foam::conformalVoronoiMesh::reinsertBoundingPoints()
{
tmp<pointField> farPts = limitBounds_.points();
forAll(farPts(), fPI)
{
insertPoint(farPts()[fPI], Vb::vtFar);
}
}
void Foam::conformalVoronoiMesh::reinsertFeaturePoints(bool distribute)
{
Info<< nl << "Reinserting stored feature points" << endl;
label preReinsertionSize(number_of_vertices());
if (distribute)
{
DynamicList<Foam::point> pointsToInsert;
DynamicList<label> indices;
DynamicList<label> types;
for
(
std::list<Vb>::iterator vit=featureVertices_.begin();
vit != featureVertices_.end();
++vit
)
{
pointsToInsert.append(topoint(vit->point()));
indices.append(vit->index());
types.append(vit->type());
}
// Insert distributed points
insertPoints(pointsToInsert, indices, types, true);
// Save points in new distribution
featureVertices_.clear();
forAll(pointsToInsert, pI)
{
featureVertices_.push_back
(
Vb(toPoint(pointsToInsert[pI]), indices[pI], types[pI])
);
}
}
else
{
for
(
std::list<Vb>::iterator vit=featureVertices_.begin();
vit != featureVertices_.end();
++vit
)
{
// Assuming that all of the reinsertions are pair points, and that
// the index and type are relative, i.e. index 0 and type relative
// to it.
insertPoint
(
vit->point(),
vit->index() + number_of_vertices(),
vit->type() + number_of_vertices()
);
}
}
Info<< " Reinserted "
<< returnReduce
(
label(number_of_vertices()) - preReinsertionSize,
sumOp<label>()
)
<< " vertices" << endl;
}
void Foam::conformalVoronoiMesh::createConvexFeaturePoints
(
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
const PtrList<extendedFeatureEdgeMesh>& feMeshes
(
geometryToConformTo_.features()
);
forAll(feMeshes, i)
{
const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
for
(
label ptI = feMesh.convexStart();
ptI < feMesh.concaveStart();
ptI++
)
{
const Foam::point& featPt = feMesh.points()[ptI];
if (!positionOnThisProc(featPt))
{
continue;
}
vectorField featPtNormals = feMesh.featurePointNormals(ptI);
vector cornerNormal = sum(featPtNormals);
cornerNormal /= mag(cornerNormal);
Foam::point internalPt =
featPt - pointPairDistance(featPt)*cornerNormal;
// Result when the points are eventually inserted (example n = 4)
// Add number_of_vertices() at insertion of first vertex to all
// numbers:
// pt index type
// internalPt 0 1
// externalPt0 1 0
// externalPt1 2 0
// externalPt2 3 0
// externalPt3 4 0
pts.append(internalPt);
indices.append(0);
types.append(1);
label internalPtIndex = -1;
forAll(featPtNormals, nI)
{
const vector& n = featPtNormals[nI];
plane planeN = plane(featPt, n);
Foam::point externalPt =
internalPt + 2.0 * planeN.distance(internalPt) * n;
pts.append(externalPt);
indices.append(0);
types.append(internalPtIndex--);
}
}
}
}
void Foam::conformalVoronoiMesh::createConcaveFeaturePoints
(
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
const PtrList<extendedFeatureEdgeMesh>& feMeshes
(
geometryToConformTo_.features()
);
forAll(feMeshes, i)
{
const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
for
(
label ptI = feMesh.concaveStart();
ptI < feMesh.mixedStart();
ptI++
)
{
const Foam::point& featPt = feMesh.points()[ptI];
if (!positionOnThisProc(featPt))
{
continue;
}
vectorField featPtNormals = feMesh.featurePointNormals(ptI);
vector cornerNormal = sum(featPtNormals);
cornerNormal /= mag(cornerNormal);
Foam::point externalPt =
featPt + pointPairDistance(featPt)*cornerNormal;
label externalPtIndex = featPtNormals.size();
// Result when the points are eventually inserted (example n = 5)
// Add number_of_vertices() at insertion of first vertex to all
// numbers:
// pt index type
// internalPt0 0 5
// internalPt1 1 5
// internalPt2 2 5
// internalPt3 3 5
// internalPt4 4 5
// externalPt 5 4
forAll(featPtNormals, nI)
{
const vector& n = featPtNormals[nI];
plane planeN = plane(featPt, n);
Foam::point internalPt =
externalPt - 2.0 * planeN.distance(externalPt) * n;
pts.append(internalPt);
indices.append(0);
types.append(externalPtIndex--);
}
pts.append(externalPt);
indices.append(0);
types.append(-1);
}
}
}
void Foam::conformalVoronoiMesh::createMixedFeaturePoints
(
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
)
{
if (cvMeshControls().mixedFeaturePointPPDistanceCoeff() < 0)
{
Info<<nl << "Skipping specialised handling for mixed feature points"
<< endl;
return;
}
const PtrList<extendedFeatureEdgeMesh>& feMeshes
(
geometryToConformTo_.features()
);
forAll(feMeshes, i)
{
const extendedFeatureEdgeMesh& feMesh = feMeshes[i];
const labelListList& pointsEdges = feMesh.pointEdges();
const pointField& points = feMesh.points();
for
(
label ptI = feMesh.mixedStart();
ptI < feMesh.nonFeatureStart();
ptI++
)
{
const labelList& pEds = pointsEdges[ptI];
pointFeatureEdgesTypes pFEdgeTypes(ptI);
const List<extendedFeatureEdgeMesh::edgeStatus> allEdStat
= calcPointFeatureEdgesTypes(feMesh, pEds, pFEdgeTypes);
bool specialisedSuccess = false;
if (cvMeshControls().specialiseFeaturePoints())
{
specialisedSuccess = createSpecialisedFeaturePoint
(
feMesh, pEds, pFEdgeTypes, allEdStat, ptI,
pts, indices, types
);
}
if (!specialisedSuccess)
{
// Specialisations available for some mixed feature points. For
// non-specialised feature points, inserting mixed internal and
// external edge groups at feature point.
// Skipping unsupported mixed feature point types
bool skipEdge = false;
forAll(pEds, e)
{
const label edgeI = pEds[e];
const extendedFeatureEdgeMesh::edgeStatus edStatus =
feMesh.getEdgeStatus(edgeI);
if
(
edStatus == extendedFeatureEdgeMesh::OPEN
|| edStatus == extendedFeatureEdgeMesh::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 Foam::point& pt = points[ptI];
if (!positionOnThisProc(pt))
{
continue;
}
const scalar edgeGroupDistance = mixedFeaturePointDistance(pt);
forAll(pEds, e)
{
const label edgeI = pEds[e];
const Foam::point edgePt =
pt + edgeGroupDistance*feMesh.edgeDirection(edgeI, ptI);
const pointIndexHit edgeHit(true, edgePt, edgeI);
createEdgePointGroup(feMesh, edgeHit, pts, indices, types);
}
}
}
}
}
void Foam::conformalVoronoiMesh::insertFeaturePoints()
{
Info<< nl << "Conforming to feature points" << endl;
DynamicList<Foam::point> pts;
DynamicList<label> indices;
DynamicList<label> types;
const label preFeaturePointSize = number_of_vertices();
createConvexFeaturePoints(pts, indices, types);
createConcaveFeaturePoints(pts, indices, types);
createMixedFeaturePoints(pts, indices, types);
// Insert the created points, distributing to the appropriate processor
insertPoints(pts, indices, types, true);
if (cvMeshControls().objOutput())
{
writePoints("featureVertices.obj", pts);
}
label nFeatureVertices = number_of_vertices() - preFeaturePointSize;
if (Pstream::parRun())
{
reduce(nFeatureVertices, sumOp<label>());
}
if (nFeatureVertices > 0)
{
Info<< " Inserted " << nFeatureVertices
<< " feature vertices" << endl;
}
featureVertices_.clear();
forAll(pts, pI)
{
featureVertices_.push_back
(
Vb(toPoint(pts[pI]), indices[pI], types[pI])
);
}
constructFeaturePointLocations();
}
void Foam::conformalVoronoiMesh::constructFeaturePointLocations()
{
DynamicList<Foam::point> ftPtLocs;
const PtrList<extendedFeatureEdgeMesh>& feMeshes
(
geometryToConformTo_.features()
);
forAll(feMeshes, i)
{
const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
if (cvMeshControls().mixedFeaturePointPPDistanceCoeff() < 0)
{
// Ignoring mixed feature points
for
(
label ptI = feMesh.convexStart();
ptI < feMesh.mixedStart();
ptI++
)
{
ftPtLocs.append(feMesh.points()[ptI]);
}
}
else
{
for
(
label ptI = feMesh.convexStart();
ptI < feMesh.nonFeatureStart();
ptI++
)
{
ftPtLocs.append(feMesh.points()[ptI]);
}
}
}
featurePointLocations_.transfer(ftPtLocs);
}

30
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
View File

@ -620,6 +620,8 @@ void Foam::conformalVoronoiMesh::writeMesh
writeCellSizes(mesh); writeCellSizes(mesh);
writeCellCentres(mesh);
findRemainingProtrusionSet(mesh); findRemainingProtrusionSet(mesh);
} }
@ -805,6 +807,34 @@ void Foam::conformalVoronoiMesh::writeCellSizes
} }
void Foam::conformalVoronoiMesh::writeCellCentres
(
const fvMesh& mesh
) const
{
Info<< "Writing components of cellCentre positions to volScalarFields"
<< " ccx, ccy, ccz in " << runTime_.timeName() << endl;
for (direction i=0; i<vector::nComponents; i++)
{
volScalarField cci
(
IOobject
(
"cc" + word(vector::componentNames[i]),
runTime_.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh.C().component(i)
);
cci.write();
}
}
void Foam::conformalVoronoiMesh::findRemainingProtrusionSet void Foam::conformalVoronoiMesh::findRemainingProtrusionSet
( (
const fvMesh& mesh const fvMesh& mesh

290
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
View File

@ -33,6 +33,10 @@ Description
An indexed form of CGAL::Triangulation_cell_base_3<K> used to keep An indexed form of CGAL::Triangulation_cell_base_3<K> used to keep
track of the Delaunay cells (tets) in the tessellation. track of the Delaunay cells (tets) in the tessellation.
SourceFiles
indexedCellI.H
indexedCell.C
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#ifndef indexedCell_H #ifndef indexedCell_H
@ -119,248 +123,72 @@ public:
}; };
indexedCell() // Constructors
:
Cb(),
index_(ctFar),
filterCount_(0)
{}
inline indexedCell();
indexedCell inline indexedCell
(
Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3
)
:
Cb(v0, v1, v2, v3),
index_(ctFar),
filterCount_(0)
{}
indexedCell
(
Vertex_handle v0,
Vertex_handle v1,
Vertex_handle v2,
Vertex_handle v3,
Cell_handle n0,
Cell_handle n1,
Cell_handle n2,
Cell_handle n3
)
:
Cb(v0, v1, v2, v3, n0, n1, n2, n3),
index_(ctFar),
filterCount_(0)
{}
int& cellIndex()
{
return index_;
}
int cellIndex() const
{
return index_;
}
inline bool farCell() const
{
return index_ == ctFar;
}
inline int& filterCount()
{
return filterCount_;
}
inline int filterCount() const
{
return filterCount_;
}
//- Is the Delaunay cell real, i.e. any real vertex
inline bool real() const
{
return
( (
this->vertex(0)->real() Vertex_handle v0,
|| this->vertex(1)->real() Vertex_handle v1,
|| this->vertex(2)->real() Vertex_handle v2,
|| this->vertex(3)->real() Vertex_handle v3
);
}
//- Does the Dual vertex form part of a processor patch
inline bool parallelDualVertex() const
{
return
(
this->vertex(0)->referred()
|| this->vertex(1)->referred()
|| this->vertex(2)->referred()
|| this->vertex(3)->referred()
);
}
//- Does the Dual vertex form part of a processor patch
inline Foam::label dualVertexMasterProc() const
{
if (!parallelDualVertex())
{
return -1;
}
// The master processor is the lowest numbered of the four on this tet.
int masterProc = Foam::Pstream::nProcs() + 1;
for (int i = 0; i < 4; i++)
{
if (this->vertex(i)->referred())
{
masterProc = min(masterProc, this->vertex(i)->procIndex());
}
else
{
masterProc = min(masterProc, Foam::Pstream::myProcNo());
}
}
return masterProc;
}
inline Foam::FixedList<Foam::label, 4> processorsAttached() const
{
if (!parallelDualVertex())
{
return Foam::FixedList<Foam::label, 4>(Foam::Pstream::myProcNo());
}
Foam::FixedList<Foam::label, 4> procsAttached
(
Foam::Pstream::myProcNo()
); );
for (int i = 0; i < 4; i++) inline indexedCell
{
if (this->vertex(i)->referred())
{
procsAttached[i] = this->vertex(i)->procIndex();
}
}
return procsAttached;
}
// Using the globalIndex object, return a list of four (sorted) global
// Delaunay vertex indices that uniquely identify this tet in parallel
inline Foam::tetCell vertexGlobalIndices
(
const Foam::globalIndex& globalDelaunayVertexIndices
) const
{
// tetVertexGlobalIndices
Foam::tetCell tVGI;
for (int i = 0; i < 4; i++)
{
Vertex_handle v = this->vertex(i);
// Finding the global index of each Delaunay vertex
if (v->referred())
{
// Referred vertices may have negative indices
tVGI[i] = globalDelaunayVertexIndices.toGlobal
(
v->procIndex(),
Foam::mag(v->index())
);
}
else
{
tVGI[i] = globalDelaunayVertexIndices.toGlobal
(
Foam::Pstream::myProcNo(),
v->index()
);
}
}
// bubble sort
for (int i = 0; i < tVGI.size(); i++)
{
for (int j = tVGI.size() - 1 ; j > i; j--)
{
if (tVGI[j - 1] > tVGI[j])
{
Foam::Swap(tVGI[j - 1], tVGI[j]);
}
}
}
return tVGI;
}
// Is the Delaunay cell part of the final dual mesh, i.e. any vertex form
// part of the internal or boundary definition
inline bool internalOrBoundaryDualVertex() const
{
return
( (
this->vertex(0)->internalOrBoundaryPoint() Vertex_handle v0,
|| this->vertex(1)->internalOrBoundaryPoint() Vertex_handle v1,
|| this->vertex(2)->internalOrBoundaryPoint() Vertex_handle v2,
|| this->vertex(3)->internalOrBoundaryPoint() Vertex_handle v3,
Cell_handle n0,
Cell_handle n1,
Cell_handle n2,
Cell_handle n3
); );
}
// Is the Delaunay cell real or referred (or mixed), i.e. all vertices form // Member Functions
// part of the real or referred internal or boundary definition
inline bool anyInternalOrBoundaryDualVertex() const inline int& cellIndex();
{
return inline int cellIndex() const;
inline bool farCell() const;
inline int& filterCount();
inline int filterCount() const;
//- Is the Delaunay cell real, i.e. any real vertex
inline bool real() const;
//- Does the Dual vertex form part of a processor patch
inline bool parallelDualVertex() const;
//- Does the Dual vertex form part of a processor patch
inline Foam::label dualVertexMasterProc() const;
inline Foam::FixedList<Foam::label, 4> processorsAttached() const;
//- Using the globalIndex object, return a list of four (sorted) global
// Delaunay vertex indices that uniquely identify this tet in parallel
inline Foam::tetCell vertexGlobalIndices
( (
this->vertex(0)->anyInternalOrBoundaryPoint() const Foam::globalIndex& globalDelaunayVertexIndices
|| this->vertex(1)->anyInternalOrBoundaryPoint() ) const;
|| this->vertex(2)->anyInternalOrBoundaryPoint()
|| this->vertex(3)->anyInternalOrBoundaryPoint()
);
}
//- Is the Delaunay cell part of the final dual mesh, i.e. any vertex
// form part of the internal or boundary definition
inline bool internalOrBoundaryDualVertex() const;
// A dual vertex on the boundary will result from a Delaunay cell with //- Is the Delaunay cell real or referred (or mixed), i.e. all vertices
// least one Delaunay vertex outside and at least one inside // form part of the real or referred internal or boundary definition
inline bool boundaryDualVertex() const inline bool anyInternalOrBoundaryDualVertex() const;
{
return //- A dual vertex on the boundary will result from a Delaunay cell with
( // least one Delaunay vertex outside and at least one inside
( inline bool boundaryDualVertex() const;
this->vertex(0)->internalOrBoundaryPoint()
|| this->vertex(1)->internalOrBoundaryPoint()
|| this->vertex(2)->internalOrBoundaryPoint()
|| this->vertex(3)->internalOrBoundaryPoint()
)
&& (
!this->vertex(0)->internalOrBoundaryPoint()
|| !this->vertex(1)->internalOrBoundaryPoint()
|| !this->vertex(2)->internalOrBoundaryPoint()
|| !this->vertex(3)->internalOrBoundaryPoint()
)
);
}
// Info // Info
@ -387,6 +215,8 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "indexedCellI.H"
#ifdef NoRepository #ifdef NoRepository
# include "indexedCell.C" # include "indexedCell.C"
#endif #endif

288
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H Normal file
View File

@ -0,0 +1,288 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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/>.
As a special exception, you have permission to link this program with the
CGAL library and distribute executables, as long as you follow the
requirements of the GNU GPL in regard to all of the software in the
executable aside from CGAL.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Gt, class Cb>
CGAL::indexedCell<Gt, Cb>::indexedCell()
:
Cb(),
index_(ctFar),
filterCount_(0)
{}
template<class Gt, class Cb>
CGAL::indexedCell<Gt, Cb>::indexedCell
(
Vertex_handle v0, Vertex_handle v1, Vertex_handle v2, Vertex_handle v3
)
:
Cb(v0, v1, v2, v3),
index_(ctFar),
filterCount_(0)
{}
template<class Gt, class Cb>
CGAL::indexedCell<Gt, Cb>::indexedCell
(
Vertex_handle v0,
Vertex_handle v1,
Vertex_handle v2,
Vertex_handle v3,
Cell_handle n0,
Cell_handle n1,
Cell_handle n2,
Cell_handle n3
)
:
Cb(v0, v1, v2, v3, n0, n1, n2, n3),
index_(ctFar),
filterCount_(0)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Gt, class Cb>
int& CGAL::indexedCell<Gt, Cb>::cellIndex()
{
return index_;
}
template<class Gt, class Cb>
int CGAL::indexedCell<Gt, Cb>::cellIndex() const
{
return index_;
}
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::farCell() const
{
return index_ == ctFar;
}
template<class Gt, class Cb>
inline int& CGAL::indexedCell<Gt, Cb>::filterCount()
{
return filterCount_;
}
template<class Gt, class Cb>
inline int CGAL::indexedCell<Gt, Cb>::filterCount() const
{
return filterCount_;
}
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::real() const
{
return
(
this->vertex(0)->real()
|| this->vertex(1)->real()
|| this->vertex(2)->real()
|| this->vertex(3)->real()
);
}
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::parallelDualVertex() const
{
return
(
this->vertex(0)->referred()
|| this->vertex(1)->referred()
|| this->vertex(2)->referred()
|| this->vertex(3)->referred()
);
}
template<class Gt, class Cb>
inline Foam::label CGAL::indexedCell<Gt, Cb>::dualVertexMasterProc() const
{
if (!parallelDualVertex())
{
return -1;
}
// The master processor is the lowest numbered of the four on this tet.
int masterProc = Foam::Pstream::nProcs() + 1;
for (int i = 0; i < 4; i++)
{
if (this->vertex(i)->referred())
{
masterProc = min(masterProc, this->vertex(i)->procIndex());
}
else
{
masterProc = min(masterProc, Foam::Pstream::myProcNo());
}
}
return masterProc;
}
template<class Gt, class Cb>
inline Foam::FixedList<Foam::label, 4>
CGAL::indexedCell<Gt, Cb>::processorsAttached() const
{
if (!parallelDualVertex())
{
return Foam::FixedList<Foam::label, 4>(Foam::Pstream::myProcNo());
}
Foam::FixedList<Foam::label, 4> procsAttached
(
Foam::Pstream::myProcNo()
);
for (int i = 0; i < 4; i++)
{
if (this->vertex(i)->referred())
{
procsAttached[i] = this->vertex(i)->procIndex();
}
}
return procsAttached;
}
template<class Gt, class Cb>
inline Foam::tetCell CGAL::indexedCell<Gt, Cb>::vertexGlobalIndices
(
const Foam::globalIndex& globalDelaunayVertexIndices
) const
{
// tetVertexGlobalIndices
Foam::tetCell tVGI;
for (int i = 0; i < 4; i++)
{
Vertex_handle v = this->vertex(i);
// Finding the global index of each Delaunay vertex
if (v->referred())
{
// Referred vertices may have negative indices
tVGI[i] = globalDelaunayVertexIndices.toGlobal
(
v->procIndex(),
Foam::mag(v->index())
);
}
else
{
tVGI[i] = globalDelaunayVertexIndices.toGlobal
(
Foam::Pstream::myProcNo(),
v->index()
);
}
}
// bubble sort
for (int i = 0; i < tVGI.size(); i++)
{
for (int j = tVGI.size() - 1 ; j > i; j--)
{
if (tVGI[j - 1] > tVGI[j])
{
Foam::Swap(tVGI[j - 1], tVGI[j]);
}
}
}
return tVGI;
}
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::internalOrBoundaryDualVertex() const
{
return
(
this->vertex(0)->internalOrBoundaryPoint()
|| this->vertex(1)->internalOrBoundaryPoint()
|| this->vertex(2)->internalOrBoundaryPoint()
|| this->vertex(3)->internalOrBoundaryPoint()
);
}
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::anyInternalOrBoundaryDualVertex() const
{
return
(
this->vertex(0)->anyInternalOrBoundaryPoint()
|| this->vertex(1)->anyInternalOrBoundaryPoint()
|| this->vertex(2)->anyInternalOrBoundaryPoint()
|| this->vertex(3)->anyInternalOrBoundaryPoint()
);
}
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::boundaryDualVertex() const
{
return
(
(
this->vertex(0)->internalOrBoundaryPoint()
|| this->vertex(1)->internalOrBoundaryPoint()
|| this->vertex(2)->internalOrBoundaryPoint()
|| this->vertex(3)->internalOrBoundaryPoint()
)
&& (
!this->vertex(0)->internalOrBoundaryPoint()
|| !this->vertex(1)->internalOrBoundaryPoint()
|| !this->vertex(2)->internalOrBoundaryPoint()
|| !this->vertex(3)->internalOrBoundaryPoint()
)
);
}
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //

288
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
View File

@ -33,6 +33,10 @@ Description
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep
track of the Delaunay vertices in the tessellation. track of the Delaunay vertices in the tessellation.
SourceFiles
indexedVertexI.H
indexedVertex.C
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#ifndef indexedVertex_H #ifndef indexedVertex_H
@ -116,267 +120,99 @@ public:
}; };
indexedVertex() // Constructors
:
Vb(),
index_(vtInternal),
type_(vtInternal),
alignment_(),
targetCellSize_(0.0)
{}
inline indexedVertex();
indexedVertex(const Point& p) inline indexedVertex(const Point& p);
:
Vb(p),
index_(vtInternal),
type_(vtInternal),
alignment_(),
targetCellSize_(0.0)
{}
inline indexedVertex(const Point& p, int index, int type);
indexedVertex(const Point& p, int index, int type) inline indexedVertex(const Point& p, Cell_handle f);
:
Vb(p),
index_(index),
type_(type),
alignment_(),
targetCellSize_(0.0)
{}
inline indexedVertex(Cell_handle f);
indexedVertex(const Point& p, Cell_handle f)
:
Vb(f, p),
index_(vtInternal),
type_(vtInternal),
alignment_(),
targetCellSize_(0.0)
{}
// Member Functions
indexedVertex(Cell_handle f) inline int& index();
:
Vb(f),
index_(vtInternal),
type_(vtInternal),
alignment_(),
targetCellSize_(0.0)
{}
inline int index() const;
int& index() inline int& type();
{
return index_;
}
inline int type() const;
int index() const inline Foam::tensor& alignment();
{
return index_;
}
inline const Foam::tensor& alignment() const;
int& type() inline Foam::scalar& targetCellSize();
{
return type_;
}
inline Foam::scalar targetCellSize() const;
int type() const inline bool uninitialised() const;
{
return type_;
}
//- Is point a far-point
inline bool farPoint() const;
inline Foam::tensor& alignment() //- Is point internal, i.e. not on boundary
{ inline bool internalPoint() const;
return alignment_;
}
//- Is point internal, i.e. not on boundary, external query.
inline static bool internalPoint(int type);
inline const Foam::tensor& alignment() const // is this a referred vertex
{ inline bool referred() const;
return alignment_;
}
// is this a referred internal or boundary vertex
inline bool referredInternalOrBoundaryPoint() const;
inline Foam::scalar& targetCellSize() // is this a referred external (pair slave) vertex
{ inline bool referredExternal() const;
return targetCellSize_;
}
// Is this a "real" point on this processor, i.e. is it internal or part
// of the boundary description, and not a "far" or "referred" point
inline bool real() const;
inline Foam::scalar targetCellSize() const // For referred vertices, what is the original processor index
{ inline int procIndex() const;
return targetCellSize_;
}
inline static int encodeProcIndex(int procI);
inline bool uninitialised() const //- Set the point to be internal
{ inline void setInternal();
return type_ == vtInternal && index_ == vtInternal;
}
//- Is point internal and near the boundary
inline bool nearBoundary() const;
//- Is point a far-point //- Set the point to be near the boundary
inline bool farPoint() const inline void setNearBoundary();
{
return type_ == vtFar;
}
//- Either master or slave of pointPair.
inline bool pairPoint() const;
//- Is point internal, i.e. not on boundary //- Master of a pointPair is the lowest numbered one.
inline bool internalPoint() const inline bool ppMaster() const;
{
return internalPoint(type_);
}
//- Master of a pointPair is the lowest numbered one, external query.
inline static bool ppMaster(int index, int type);
//- Is point internal, i.e. not on boundary, external query. //- Slave of a pointPair is the highest numbered one.
inline static bool internalPoint(int type) inline bool ppSlave() const;
{
return type <= vtInternal;
}
//- Either original internal point or master of pointPair.
inline bool internalOrBoundaryPoint() const;
// is this a referred vertex //- Either original internal point or master of pointPair.
inline bool referred() const // External query.
{ inline static bool internalOrBoundaryPoint(int index, int type);
return (type_ < 0 && type_ > vtFar);
}
//- Is point near the boundary or part of the boundary definition
inline bool nearOrOnBoundary() const;
// is this a referred internal or boundary vertex //- Either a real or referred internal or boundary point
inline bool referredInternalOrBoundaryPoint() const inline bool anyInternalOrBoundaryPoint() const;
{
return referred() && index_ >= 0;
}
// is this a referred external (pair slave) vertex
inline bool referredExternal() const
{
return referred() && index_ < 0;
}
// is this a "real" point on this processor, i.e. is it internal or part of
// the boundary description, and not a "far" or "referred" point
inline bool real() const
{
return internalPoint() || pairPoint();
}
// For referred vertices, what is the original processor index
inline int procIndex() const
{
if (referred())
{
return -(type_ + 1);
}
else
{
return -1;
}
}
inline static int encodeProcIndex(int procI)
{
return -(procI + 1);
}
//- Set the point to be internal
inline void setInternal()
{
type_ = vtInternal;
}
//- Is point internal and near the boundary
inline bool nearBoundary() const
{
return type_ == vtNearBoundary;
}
//- Set the point to be near the boundary
inline void setNearBoundary()
{
type_ = vtNearBoundary;
}
//- Either master or slave of pointPair.
inline bool pairPoint() const
{
return type_ >= 0;
}
//- Master of a pointPair is the lowest numbered one.
inline bool ppMaster() const
{
return ppMaster(index_, type_);
}
//- Master of a pointPair is the lowest numbered one, external query.
inline static bool ppMaster(int index, int type)
{
if (index >= 0 && type > index)
{
return true;
}
return false;
}
//- Slave of a pointPair is the highest numbered one.
inline bool ppSlave() const
{
if (type_ >= 0 && type_ < index_)
{
return true;
}
else
{
return false;
}
}
//- Either original internal point or master of pointPair.
inline bool internalOrBoundaryPoint() const
{
return internalOrBoundaryPoint(index_, type_);
}
//- Either original internal point or master of pointPair, external query.
inline static bool internalOrBoundaryPoint(int index, int type)
{
return internalPoint(type) || ppMaster(index, type);
}
//- Is point near the boundary or part of the boundary definition
inline bool nearOrOnBoundary() const
{
return pairPoint() || nearBoundary();
}
//- Either a real or referred internal or boundary point
inline bool anyInternalOrBoundaryPoint() const
{
return internalOrBoundaryPoint() || referredInternalOrBoundaryPoint();
}
// inline void operator=(const Delaunay::Finite_vertices_iterator vit) // inline void operator=(const Delaunay::Finite_vertices_iterator vit)
@ -414,6 +250,8 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "indexedVertexI.H"
#ifdef NoRepository #ifdef NoRepository
# include "indexedVertex.C" # include "indexedVertex.C"
#endif #endif

321
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H Normal file
View File

@ -0,0 +1,321 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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/>.
As a special exception, you have permission to link this program with the
CGAL library and distribute executables, as long as you follow the
requirements of the GNU GPL in regard to all of the software in the
executable aside from CGAL.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Gt, class Vb>
inline CGAL::indexedVertex<Gt, Vb>::indexedVertex()
:
Vb(),
index_(vtInternal),
type_(vtInternal),
alignment_(),
targetCellSize_(0.0)
{}
template<class Gt, class Vb>
inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(const Point& p)
:
Vb(p),
index_(vtInternal),
type_(vtInternal),
alignment_(),
targetCellSize_(0.0)
{}
template<class Gt, class Vb>
inline CGAL::indexedVertex<Gt, Vb>::indexedVertex
(
const Point& p,
int index,
int type
)
:
Vb(p),
index_(index),
type_(type),
alignment_(),
targetCellSize_(0.0)
{}
template<class Gt, class Vb>
inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(const Point& p, Cell_handle f)
:
Vb(f, p),
index_(vtInternal),
type_(vtInternal),
alignment_(),
targetCellSize_(0.0)
{}
template<class Gt, class Vb>
inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(Cell_handle f)
:
Vb(f),
index_(vtInternal),
type_(vtInternal),
alignment_(),
targetCellSize_(0.0)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Gt, class Vb>
inline int& CGAL::indexedVertex<Gt, Vb>::index()
{
return index_;
}
template<class Gt, class Vb>
inline int CGAL::indexedVertex<Gt, Vb>::index() const
{
return index_;
}
template<class Gt, class Vb>
inline int& CGAL::indexedVertex<Gt, Vb>::type()
{
return type_;
}
template<class Gt, class Vb>
inline int CGAL::indexedVertex<Gt, Vb>::type() const
{
return type_;
}
template<class Gt, class Vb>
inline Foam::tensor& CGAL::indexedVertex<Gt, Vb>::alignment()
{
return alignment_;
}
template<class Gt, class Vb>
inline const Foam::tensor& CGAL::indexedVertex<Gt, Vb>::alignment() const
{
return alignment_;
}
template<class Gt, class Vb>
inline Foam::scalar& CGAL::indexedVertex<Gt, Vb>::targetCellSize()
{
return targetCellSize_;
}
template<class Gt, class Vb>
inline Foam::scalar CGAL::indexedVertex<Gt, Vb>::targetCellSize() const
{
return targetCellSize_;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::uninitialised() const
{
return type_ == vtInternal && index_ == vtInternal;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::farPoint() const
{
return type_ == vtFar;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::internalPoint() const
{
return internalPoint(type_);
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::internalPoint(int type)
{
return type <= vtInternal;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::referred() const
{
return (type_ < 0 && type_ > vtFar);
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::referredInternalOrBoundaryPoint() const
{
return referred() && index_ >= 0;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::referredExternal() const
{
return referred() && index_ < 0;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::real() const
{
return internalPoint() || pairPoint();
}
template<class Gt, class Vb>
inline int CGAL::indexedVertex<Gt, Vb>::procIndex() const
{
if (referred())
{
return -(type_ + 1);
}
else
{
return -1;
}
}
template<class Gt, class Vb>
inline int CGAL::indexedVertex<Gt, Vb>::encodeProcIndex(int procI)
{
return -(procI + 1);
}
template<class Gt, class Vb>
inline void CGAL::indexedVertex<Gt, Vb>::setInternal()
{
type_ = vtInternal;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::nearBoundary() const
{
return type_ == vtNearBoundary;
}
template<class Gt, class Vb>
inline void CGAL::indexedVertex<Gt, Vb>::setNearBoundary()
{
type_ = vtNearBoundary;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::pairPoint() const
{
return type_ >= 0;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::ppMaster() const
{
return ppMaster(index_, type_);
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::ppMaster(int index, int type)
{
if (index >= 0 && type > index)
{
return true;
}
return false;
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::ppSlave() const
{
if (type_ >= 0 && type_ < index_)
{
return true;
}
else
{
return false;
}
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::internalOrBoundaryPoint() const
{
return internalOrBoundaryPoint(index_, type_);
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::internalOrBoundaryPoint
(
int index,
int type
)
{
return internalPoint(type) || ppMaster(index, type);
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::nearOrOnBoundary() const
{
return pairPoint() || nearBoundary();
}
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::anyInternalOrBoundaryPoint() const
{
return internalOrBoundaryPoint() || referredInternalOrBoundaryPoint();
}
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //

87
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/pointFeatureEdgesTypes.H Normal file
View File

@ -0,0 +1,87 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 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/>.
As a special exception, you have permission to link this program with the
CGAL library and distribute executables, as long as you follow the
requirements of the GNU GPL in regard to all of the software in the
executable aside from CGAL.
Class
pointFeatureEdgesTypes
Description
struct for holding information on the types of feature edges attached to
feature points
\*---------------------------------------------------------------------------*/
#ifndef pointFeatureEdgesTypes_H
#define pointFeatureEdgesTypes_H
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class pointFeatureEdgesTypes Declaration
\*---------------------------------------------------------------------------*/
//- Hold the types of feature edges attached to the point.
struct pointFeatureEdgesTypes
{
label ptI;
label nExternal, nInternal;
label nFlat, nOpen, nMultiple, nNonFeature;
pointFeatureEdgesTypes(const label pointLabel)
{
ptI = pointLabel;
nExternal = nInternal = 0;
nFlat = nOpen = nMultiple = nNonFeature = 0;
}
friend Ostream& operator<<(Ostream& os, const pointFeatureEdgesTypes& p)
{
os << "E " << "I " << "F " << "O " << "M " << "N "
<< "Pt" << nl
<< p.nExternal << " " << p.nInternal << " "
<< p.nFlat << " " << p.nOpen << " "
<< p.nMultiple << " " << p.nNonFeature << " "
<< p.ptI
<< endl;
return os;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

7
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
View File

@ -30,11 +30,9 @@ License
Foam::cvControls::cvControls Foam::cvControls::cvControls
( (
const conformalVoronoiMesh& cvMesh,
const dictionary& cvMeshDict const dictionary& cvMeshDict
) )
: :
cvMesh_(cvMesh),
cvMeshDict_(cvMeshDict) cvMeshDict_(cvMeshDict)
{ {
// Surface conformation controls // Surface conformation controls
@ -61,6 +59,11 @@ Foam::cvControls::cvControls
surfDict.lookup("featureEdgeExclusionDistanceCoeff") surfDict.lookup("featureEdgeExclusionDistanceCoeff")
); );
specialiseFeaturePoints_ = Switch
(
surfDict.lookupOrDefault<Switch>("specialiseFeaturePoints", false)
);
surfaceSearchDistanceCoeff_ = readScalar surfaceSearchDistanceCoeff_ = readScalar
( (
surfDict.lookup("surfaceSearchDistanceCoeff") surfDict.lookup("surfaceSearchDistanceCoeff")

10
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
View File

@ -54,9 +54,6 @@ class cvControls
{ {
// Private data // Private data
//- Reference to the conformalVoronoiMesh holding this cvControls object
const conformalVoronoiMesh& cvMesh_;
//- Reference to the cvMeshDict //- Reference to the cvMeshDict
const dictionary& cvMeshDict_; const dictionary& cvMeshDict_;
@ -80,6 +77,9 @@ class cvControls
// fraction of the local target cell size // fraction of the local target cell size
scalar featureEdgeExclusionDistanceCoeff_; scalar featureEdgeExclusionDistanceCoeff_;
//- Switch for using specialised feature points
Switch specialiseFeaturePoints_;
//- Surface search distance coefficient - fraction of the local //- Surface search distance coefficient - fraction of the local
// target cell size // target cell size
scalar surfaceSearchDistanceCoeff_; scalar surfaceSearchDistanceCoeff_;
@ -291,7 +291,6 @@ public:
//- Construct from references to conformalVoronoiMesh and dictionary //- Construct from references to conformalVoronoiMesh and dictionary
cvControls cvControls
( (
const conformalVoronoiMesh& cvMesh,
const dictionary& cvMeshDict const dictionary& cvMeshDict
); );
@ -319,6 +318,9 @@ public:
//- Return the featureEdgeExclusionDistanceCoeff //- Return the featureEdgeExclusionDistanceCoeff
inline scalar featureEdgeExclusionDistanceCoeff() const; inline scalar featureEdgeExclusionDistanceCoeff() const;
//- Return whether to use specialised feature points
inline Switch specialiseFeaturePoints() const;
//- Return the surfaceSearchDistanceCoeff //- Return the surfaceSearchDistanceCoeff
inline scalar surfaceSearchDistanceCoeff() const; inline scalar surfaceSearchDistanceCoeff() const;

4
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
View File

@ -54,6 +54,10 @@ inline Foam::scalar Foam::cvControls::featureEdgeExclusionDistanceCoeff() const
return featureEdgeExclusionDistanceCoeff_; return featureEdgeExclusionDistanceCoeff_;
} }
inline Foam::Switch Foam::cvControls::specialiseFeaturePoints() const
{
return specialiseFeaturePoints_;
}
inline Foam::scalar Foam::cvControls::surfaceSearchDistanceCoeff() const inline Foam::scalar Foam::cvControls::surfaceSearchDistanceCoeff() const
{ {

9
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.C
View File

@ -42,12 +42,10 @@ defineRunTimeSelectionTable(faceAreaWeightModel, dictionary);
faceAreaWeightModel::faceAreaWeightModel faceAreaWeightModel::faceAreaWeightModel
( (
const word& type, const word& type,
const dictionary& relaxationDict, const dictionary& relaxationDict
const conformalVoronoiMesh& cvMesh
) )
: :
dictionary(relaxationDict), dictionary(relaxationDict),
cvMesh_(cvMesh),
coeffDict_(subDict(type + "Coeffs")) coeffDict_(subDict(type + "Coeffs"))
{} {}
@ -56,8 +54,7 @@ faceAreaWeightModel::faceAreaWeightModel
autoPtr<faceAreaWeightModel> faceAreaWeightModel::New autoPtr<faceAreaWeightModel> faceAreaWeightModel::New
( (
const dictionary& relaxationDict, const dictionary& relaxationDict
const conformalVoronoiMesh& cvMesh
) )
{ {
word faceAreaWeightModelTypeName word faceAreaWeightModelTypeName
@ -85,7 +82,7 @@ autoPtr<faceAreaWeightModel> faceAreaWeightModel::New
<< exit(FatalError); << exit(FatalError);
} }
return autoPtr<faceAreaWeightModel>(cstrIter()(relaxationDict, cvMesh)); return autoPtr<faceAreaWeightModel>(cstrIter()(relaxationDict));
} }

15
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.H
View File

@ -39,7 +39,6 @@ SourceFiles
#define faceAreaWeightModel_H #define faceAreaWeightModel_H
#include "point.H" #include "point.H"
#include "conformalVoronoiMesh.H"
#include "dictionary.H" #include "dictionary.H"
#include "autoPtr.H" #include "autoPtr.H"
#include "runTimeSelectionTables.H" #include "runTimeSelectionTables.H"
@ -62,9 +61,6 @@ protected:
// Protected data // Protected data
//- Reference to the conformalVoronoiMesh holding this cvControls object
const conformalVoronoiMesh& cvMesh_;
//- Method coeffs dictionary //- Method coeffs dictionary
dictionary coeffDict_; dictionary coeffDict_;
@ -93,10 +89,9 @@ public:
faceAreaWeightModel, faceAreaWeightModel,
dictionary, dictionary,
( (
const dictionary& faceAreaWeightDict, const dictionary& faceAreaWeightDict
const conformalVoronoiMesh& cvMesh
), ),
(faceAreaWeightDict, cvMesh) (faceAreaWeightDict)
); );
@ -106,8 +101,7 @@ public:
faceAreaWeightModel faceAreaWeightModel
( (
const word& type, const word& type,
const dictionary& faceAreaWeightDict, const dictionary& faceAreaWeightDict
const conformalVoronoiMesh& cvMesh
); );
@ -116,8 +110,7 @@ public:
//- Return a reference to the selected faceAreaWeightModel //- Return a reference to the selected faceAreaWeightModel
static autoPtr<faceAreaWeightModel> New static autoPtr<faceAreaWeightModel> New
( (
const dictionary& faceAreaWeightDict, const dictionary& faceAreaWeightDict
const conformalVoronoiMesh& cvMesh
); );

5
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.C
View File

@ -45,11 +45,10 @@ addToRunTimeSelectionTable
piecewiseLinearRamp::piecewiseLinearRamp piecewiseLinearRamp::piecewiseLinearRamp
( (
const dictionary& faceAreaWeightDict, const dictionary& faceAreaWeightDict
const conformalVoronoiMesh& cvMesh
) )
: :
faceAreaWeightModel(typeName, faceAreaWeightDict, cvMesh), faceAreaWeightModel(typeName, faceAreaWeightDict),
lAF_(readScalar(coeffDict().lookup("lowerAreaFraction"))), lAF_(readScalar(coeffDict().lookup("lowerAreaFraction"))),
uAF_(readScalar(coeffDict().lookup("upperAreaFraction"))) uAF_(readScalar(coeffDict().lookup("upperAreaFraction")))
{} {}

3
applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.H
View File

@ -71,8 +71,7 @@ public:
//- Construct from components //- Construct from components
piecewiseLinearRamp piecewiseLinearRamp
( (
const dictionary& faceAreaWeightDict, const dictionary& faceAreaWeightDict
const conformalVoronoiMesh& cvMesh
); );

4
applications/utilities/mesh/generation/cvMesh/cvMesh.C
View File

@ -59,6 +59,8 @@ int main(int argc, char *argv[])
) )
); );
conformalVoronoiMesh::debug = true;
conformalVoronoiMesh mesh(runTime, cvMeshDict); conformalVoronoiMesh mesh(runTime, cvMeshDict);
while (runTime.loop()) while (runTime.loop())
@ -72,7 +74,7 @@ int main(int argc, char *argv[])
<< nl << endl; << nl << endl;
} }
mesh.writeMesh(runTime.constant()); mesh.writeMesh(runTime.constant(), true);
Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s"