diff --git a/applications/test/vectorTools/Make/files b/applications/test/vectorTools/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..0b30b98f8f3b0c0c1048846c8e90a62733a12dda
--- /dev/null
+++ b/applications/test/vectorTools/Make/files
@@ -0,0 +1,3 @@
+Test-vectorTools.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-vectorTools
diff --git a/applications/test/vectorTools/Make/options b/applications/test/vectorTools/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..9e015e6078873ee32ae6d53ac1590010819dcb69
--- /dev/null
+++ b/applications/test/vectorTools/Make/options
@@ -0,0 +1 @@
+EXE_INC = -I$(FOAM_APP)/utilities/mesh/generation/cvMesh/vectorTools
diff --git a/applications/test/vectorTools/Test-vectorTools.C b/applications/test/vectorTools/Test-vectorTools.C
new file mode 100644
index 0000000000000000000000000000000000000000..0c85ac935bcfbe7aaa654be653ecb9635cdaa356
--- /dev/null
+++ b/applications/test/vectorTools/Test-vectorTools.C
@@ -0,0 +1,71 @@
+#include "vector.H"
+#include "IOstreams.H"
+#include "vectorTools.H"
+#include "unitConversion.H"
+
+using namespace Foam;
+
+
+void test(const vector& a, const vector& b, const scalar tolerance)
+{
+ Info<< "Vectors " << a << " and " << b
+ << " are (to tolerance of " << tolerance << "): ";
+
+ if (vectorTools::areParallel(a, b, tolerance))
+ Info<< " parallel ";
+
+ if (vectorTools::areOrthogonal(a, b, tolerance))
+ Info<< " orthogonal ";
+
+ if (vectorTools::areAcute(a, b))
+ Info<< " acute ";
+
+ if (vectorTools::areObtuse(a, b))
+ Info<< " obtuse ";
+
+ Info<< ", angle = " << vectorTools::degAngleBetween(a, b);
+
+ Info<< endl;
+}
+
+
+int main()
+{
+ vector a(1.0, 1.0, 1.0);
+ vector b(2.0, 2.0, 2.0);
+
+ test(a, b, 0.0);
+ test(a, b, VSMALL);
+ test(a, b, SMALL);
+ test(a, b, 1e-3);
+ test(a, b, 1e-1);
+
+ a = vector(1,0,0);
+ b = vector(0,2,0);
+
+ test(a, b, 0.0);
+ test(a, b, VSMALL);
+ test(a, b, SMALL);
+ test(a, b, 1e-3);
+ test(a, b, 1e-1);
+
+ a = vector(1,0,0);
+ b = vector(-1,0,0);
+
+ test(a, b, 0.0);
+ test(a, b, VSMALL);
+ test(a, b, SMALL);
+ test(a, b, 1e-3);
+ test(a, b, 1e-1);
+
+ a = vector(1,0,0);
+ b = vector(-1,2,0);
+
+ test(a, b, 0.0);
+ test(a, b, VSMALL);
+ test(a, b, SMALL);
+ test(a, b, 1e-3);
+ test(a, b, 1e-1);
+
+ return 0;
+}
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/files b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/files
index 29305c3320e3620372a63a9bea4c126aea8bed4b..217ec9d12e939f51c72fc6c3423b738f3d9c004a 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/files
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/files
@@ -4,6 +4,7 @@ conformalVoronoiMesh/conformalVoronoiMesh.C
conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
conformalVoronoiMesh/conformalVoronoiMeshIO.C
+conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
cvControls/cvControls.C
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options
index 711e4ae7738e8eacbef824ea378241d33bebcd97..1a1902b559ddda2a3ff6c37acca1316f89b6b7ea 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/Make/options
@@ -14,7 +14,8 @@ EXE_INC = \
-I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
-I$(LIB_SRC)/edgeMesh/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude \
- -I$(LIB_SRC)/triSurface/lnInclude
+ -I$(LIB_SRC)/triSurface/lnInclude \
+ -I../vectorTools
EXE_LIBS = \
-lmeshTools \
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3Ddefs.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3Ddefs.H
index 8efb9ec71126115c8f79fff2b2142775092201b7..3176e4dec42543fc4fdf9d49b9ef1178ec4e81a2 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3Ddefs.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/CGALTriangulation3Ddefs.H
@@ -72,6 +72,48 @@ typedef Delaunay::Cell_handle Cell_handle;
typedef Delaunay::Point Point;
+//- Spatial sort traits to use with a pair of point pointers and an integer.
+// Taken from a post on the CGAL lists: 2010-01/msg00004.html by
+// Sebastien Loriot (Geometry Factory).
+template<class Triangulation>
+struct Traits_for_spatial_sort
+:
+ public Triangulation::Geom_traits
+{
+ typedef typename Triangulation::Geom_traits Gt;
+
+ typedef std::pair<const typename Triangulation::Point*, int> Point_3;
+
+ struct Less_x_3
+ {
+ bool operator()(const Point_3& p, const Point_3& q) const
+ {
+ return typename Gt::Less_x_3()(*(p.first), *(q.first));
+ }
+ };
+
+ struct Less_y_3
+ {
+ bool operator()(const Point_3& p, const Point_3& q) const
+ {
+ return typename Gt::Less_y_3()(*(p.first), *(q.first));
+ }
+ };
+
+ struct Less_z_3
+ {
+ bool operator()(const Point_3& p, const Point_3& q) const
+ {
+ return typename Gt::Less_z_3()(*(p.first), *(q.first));
+ }
+ };
+
+ Less_x_3 less_x_3_object () const {return Less_x_3();}
+ Less_y_3 less_y_3_object () const {return Less_y_3();}
+ Less_z_3 less_z_3_object () const {return Less_z_3();}
+};
+
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index dd58fc5141707e9f78bad63d01d6bb5ecf5bd1b6..4e97af7a7a7583dea1535bae83092571ea12bbb0 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -39,6 +39,8 @@ defineTypeNameAndDebug(conformalVoronoiMesh, 0);
}
+const Foam::scalar Foam::conformalVoronoiMesh::tolParallel = 1e-3;
+
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
@@ -78,17 +80,17 @@ Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment
norm
);
- vector np = norm[0];
+ const vector np = norm[0];
// Generate equally spaced 'spokes' in a circle normal to the
// direction from the vertex to the closest point on the surface
// and look for a secondary intersection.
- vector d = surfHit.hitPoint() - pt;
+ const vector d = surfHit.hitPoint() - pt;
- tensor Rp = rotationTensor(vector(0,0,1), np);
+ const tensor Rp = rotationTensor(vector(0,0,1), np);
- label s = cvMeshControls().alignmentSearchSpokes();
+ const label s = cvMeshControls().alignmentSearchSpokes();
scalar closestSpokeHitDistance = GREAT;
@@ -96,7 +98,7 @@ Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment
label closestSpokeSurface = -1;
- scalar spanMag = geometryToConformTo_.globalBounds().mag();
+ const scalar spanMag = geometryToConformTo_.globalBounds().mag();
for (label i = 0; i < s; i++)
{
@@ -176,8 +178,7 @@ Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment
"("
"const Foam::point& pt"
") const"
- )
- << "No secondary surface hit found in spoke search "
+ ) << "No secondary surface hit found in spoke search "
<< "using " << s
<< " spokes, try increasing alignmentSearchSpokes."
<< endl;
@@ -198,7 +199,7 @@ Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment
// Secondary alignment
vector ns = np ^ na;
- if (mag(ns) < SMALL)
+ if (mag(ns) < SMALL)
{
FatalErrorIn("conformalVoronoiMesh::requiredAlignment")
<< "Parallel normals detected in spoke search." << nl
@@ -220,7 +221,7 @@ Foam::tensor Foam::conformalVoronoiMesh::requiredAlignment
void Foam::conformalVoronoiMesh::insertPoints
(
- std::list<Point>& points,
+ List<Point>& points,
bool distribute
)
{
@@ -241,11 +242,13 @@ void Foam::conformalVoronoiMesh::insertPoints
DynamicList<Foam::point> transferPoints;
+ List<Point> pointsOnProcessor;
+
for
(
- std::list<Point>::iterator pit=points.begin();
+ List<Point>::iterator pit = points.begin();
pit != points.end();
- // No action
+ ++pit
)
{
Foam::point p(topoint(*pit));
@@ -253,15 +256,20 @@ void Foam::conformalVoronoiMesh::insertPoints
if (!positionOnThisProc(p))
{
transferPoints.append(p);
-
- pit = points.erase(pit);
}
else
{
- ++pit;
+ pointsOnProcessor.append(*pit);
}
}
+ points.setSize(pointsOnProcessor.size());
+ forAll(pointsOnProcessor, pI)
+ {
+ points[pI] = pointsOnProcessor[pI];
+ }
+ pointsOnProcessor.clear();
+
// Send the points that are not on this processor to the appropriate
// place
Foam::autoPtr<Foam::mapDistribute> map
@@ -271,7 +279,7 @@ void Foam::conformalVoronoiMesh::insertPoints
forAll(transferPoints, tPI)
{
- points.push_back(toPoint(transferPoints[tPI]));
+ points.append(toPoint(transferPoints[tPI]));
}
label sizeChange = preDistributionSize - label(points.size());
@@ -307,16 +315,16 @@ void Foam::conformalVoronoiMesh::insertPoints
// using the range insert (faster than inserting points one by one)
insert(points.begin(), points.end());
- // Info<< "USING INDIVIDUAL INSERTION TO DETECT FAILURE" << endl;
- // for
- // (
- // std::list<Point>::iterator pit=points.begin();
- // pit != points.end();
- // ++pit
- // )
- // {
- // insertPoint(topoint(*pit), Vb::vtInternal);
- // }
+// Info<< "USING INDIVIDUAL INSERTION TO DETECT FAILURE" << endl;
+// for
+// (
+// List<Point>::iterator pit=points.begin();
+// pit != points.end();
+// ++pit
+// )
+// {
+// insertPoint(topoint(*pit), Vb::vtInternal);
+// }
label nInserted(number_of_vertices() - preInsertionSize);
@@ -386,31 +394,14 @@ void Foam::conformalVoronoiMesh::insertPoints
// << " points in total" << endl;
}
- // Using index is actually pointless, it is always zero. Keep for clarity
- // of code.
-
- forAll(pts, pI)
- {
- // creation of points and indices is done assuming that it will be
- // relative to the instantaneous number_of_vertices() at insertion.
-
- label type = types[pI];
-
- if (type > Vb::vtFar)
- {
- // This is a member of a point pair, don't use the type directly
- // (note that this routine never gets called for referredPoints
- // so type will never be -procI)
- type += number_of_vertices();
- }
-
- insertPoint
- (
- pts[pI],
- indices[pI] + number_of_vertices(),
- type
- );
- }
+ rangeInsertWithInfo
+ (
+ pts.begin(),
+ pts.end(),
+ *this,
+ indices,
+ types
+ );
}
@@ -507,734 +498,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>&
Foam::conformalVoronoiMesh::featurePointTree() const
{
@@ -1277,33 +540,17 @@ bool Foam::conformalVoronoiMesh::nearFeaturePt(const Foam::point& pt) const
}
-void Foam::conformalVoronoiMesh::reset()
+void Foam::conformalVoronoiMesh::reset(const bool distribute)
{
this->clear();
- insertBoundingPoints();
-}
-
-
-void Foam::conformalVoronoiMesh::insertBoundingPoints()
-{
- pointField farPts = geometryToConformTo_.globalBounds().points();
+ reinsertBoundingPoints();
- // Shift corners of bounds relative to origin
- farPts -= geometryToConformTo_.globalBounds().midpoint();
+ // Reinsert feature points, distributing them as necessary.
+ reinsertFeaturePoints(distribute);
+ //insertFeaturePoints();
- // 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);
- }
+ startOfInternalPoints_ = number_of_vertices();
}
@@ -1313,7 +560,7 @@ void Foam::conformalVoronoiMesh::insertInitialPoints()
timeCheck("Before initial points call");
- std::list<Point> initPts = initialPointsMethod_->initialPoints();
+ List<Point> initPts = initialPointsMethod_->initialPoints();
timeCheck("After initial points call");
@@ -1466,13 +713,8 @@ bool Foam::conformalVoronoiMesh::distributeBackground()
cellVertexTypes
);
- // Remove the entire tessellation
- reset();
-
- // Reinsert feature points, distributing them as necessary.
- reinsertFeaturePoints(true);
-
- startOfInternalPoints_ = number_of_vertices();
+ // Reset the entire tessellation
+ reset(true);
timeCheck("Distribution performed");
@@ -1533,7 +775,7 @@ bool Foam::conformalVoronoiMesh::distributeBackground()
void Foam::conformalVoronoiMesh::storeSizesAndAlignments()
{
- std::list<Point> storePts;
+ DynamicList<Point> storePts(number_of_vertices());
for
(
@@ -1544,17 +786,19 @@ void Foam::conformalVoronoiMesh::storeSizesAndAlignments()
{
if (vit->internalPoint())
{
- storePts.push_back(vit->point());
+ storePts.append(vit->point());
}
}
+ storePts.shrink();
+
storeSizesAndAlignments(storePts);
}
void Foam::conformalVoronoiMesh::storeSizesAndAlignments
(
- const std::list<Point>& storePts
+ const List<Point>& storePts
)
{
timeCheck("Start of storeSizesAndAlignments");
@@ -1571,7 +815,7 @@ void Foam::conformalVoronoiMesh::storeSizesAndAlignments
for
(
- std::list<Point>::const_iterator pit=storePts.begin();
+ List<Point>::const_iterator pit = storePts.begin();
pit != storePts.end();
++pit
)
@@ -1598,7 +842,7 @@ void Foam::conformalVoronoiMesh::storeSizesAndAlignments
void Foam::conformalVoronoiMesh::updateSizesAndAlignments
(
- const std::list<Point>& storePts
+ const List<Point>& storePts
)
{
// This function is only used in serial, the background redistribution
@@ -1909,7 +1153,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
Delaunay(),
runTime_(runTime),
rndGen_(64293*Pstream::myProcNo()),
- cvMeshControls_(*this, cvMeshDict),
+ cvMeshControls_(cvMeshDict),
allGeometry_
(
IOobject
@@ -1940,6 +1184,8 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
featureVertices_(),
featurePointLocations_(),
featurePointTreePtr_(),
+ edgeLocationTreePtr_(),
+ surfacePtLocationTreePtr_(),
sizeAndAlignmentLocations_(),
storedSizes_(),
storedAlignments_(),
@@ -1965,8 +1211,7 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
(
faceAreaWeightModel::New
(
- cvMeshDict.subDict("motionControl"),
- *this
+ cvMeshDict.subDict("motionControl")
)
),
decomposition_()
@@ -2100,7 +1345,7 @@ void Foam::conformalVoronoiMesh::move()
true
);
- std::list<Point> pointsToInsert;
+ DynamicList<Point> pointsToInsert(number_of_vertices());
for
(
@@ -2137,15 +1382,15 @@ void Foam::conformalVoronoiMesh::move()
forAll(alignmentDirsA, aA)
{
- const vector& a(alignmentDirsA[aA]);
+ const vector& a = alignmentDirsA[aA];
scalar maxDotProduct = 0.0;
forAll(alignmentDirsB, aB)
{
- const vector& b(alignmentDirsB[aB]);
+ const vector& b = alignmentDirsB[aB];
- scalar dotProduct = a & b;
+ const scalar dotProduct = a & b;
if (mag(dotProduct) > maxDotProduct)
{
@@ -2179,7 +1424,7 @@ void Foam::conformalVoronoiMesh::move()
&& pointToBeRetained[vB->index()] == true
)
{
- pointsToInsert.push_back
+ pointsToInsert.append
(
toPoint(0.5*(dVA + dVB))
);
@@ -2212,7 +1457,7 @@ void Foam::conformalVoronoiMesh::move()
alignmentDir *= -1;
}
- scalar alignmentDotProd = ((rAB/rABMag) & alignmentDir);
+ const scalar alignmentDotProd = ((rAB/rABMag) & alignmentDir);
if
(
@@ -2220,9 +1465,9 @@ void Foam::conformalVoronoiMesh::move()
> cvMeshControls().cosAlignmentAcceptanceAngle()
)
{
- scalar targetCellSize = averageCellSize(vA, vB);
+ const scalar targetCellSize = averageCellSize(vA, vB);
- scalar targetFaceArea = sqr(targetCellSize);
+ const scalar targetFaceArea = sqr(targetCellSize);
alignmentDir *= 0.5*targetCellSize;
@@ -2231,7 +1476,7 @@ void Foam::conformalVoronoiMesh::move()
// directions.
vector delta = alignmentDir - 0.5*rAB;
- scalar faceArea = dualFace.mag(dualVertices);
+ const scalar faceArea = dualFace.mag(dualVertices);
delta *= faceAreaWeightModel_->faceAreaWeight
(
@@ -2262,8 +1507,7 @@ void Foam::conformalVoronoiMesh::move()
)
{
// Prevent insertions spanning surfaces
-
- pointsToInsert.push_back
+ pointsToInsert.append
(
toPoint(0.5*(dVA + dVB))
);
@@ -2285,14 +1529,13 @@ void Foam::conformalVoronoiMesh::move()
// the short edge if neither attached
// point has already been identified to be
// removed.
-
if
(
pointToBeRetained[vA->index()] == true
&& pointToBeRetained[vB->index()] == true
)
{
- pointsToInsert.push_back
+ pointsToInsert.append
(
toPoint(0.5*(dVA + dVB))
);
@@ -2371,7 +1614,7 @@ void Foam::conformalVoronoiMesh::move()
// Only necessary if using an exact constructions kernel
// (extended precision)
- pointsToInsert.push_back
+ pointsToInsert.append
(
toPoint
(
@@ -2383,12 +1626,44 @@ void Foam::conformalVoronoiMesh::move()
}
}
+ pointsToInsert.shrink();
+
+ // Save displacements to file.
+ if (cvMeshControls().objOutput() && runTime_.outputTime())
+ {
+ Pout<< "Writing point displacement vectors to file." << endl;
+ OFstream str("displacements_" + runTime_.timeName() + ".obj");
+
+ for
+ (
+ Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+ vit != finite_vertices_end();
+ ++vit
+ )
+ {
+ if (vit->internalPoint())
+ {
+ if (pointToBeRetained[vit->index()] == true)
+ {
+ meshTools::writeOBJ(str, topoint(vit->point()));
+
+ str << "vn "
+ << displacementAccumulator[vit->index()][0] << " "
+ << displacementAccumulator[vit->index()][1] << " "
+ << displacementAccumulator[vit->index()][2] << " "
+ << endl;
+ }
+ }
+ }
+ }
+
// Remove the entire tessellation
reset();
- reinsertFeaturePoints();
-
- startOfInternalPoints_ = number_of_vertices();
+ if (cvMeshControls().objOutput() && runTime_.outputTime())
+ {
+ writePoints("featurePoints_" + runTime_.timeName() + ".obj", false);
+ }
timeCheck("Displacement calculated");
@@ -2396,13 +1671,21 @@ void Foam::conformalVoronoiMesh::move()
insertPoints(pointsToInsert);
+ if (cvMeshControls().objOutput() && runTime_.outputTime())
+ {
+ writePoints("points_" + runTime_.timeName() + ".obj", false);
+ }
+
timeCheck("Internal points inserted");
conformToSurface();
- timeCheck("After conformToSurface");
+ if (cvMeshControls().objOutput() && runTime_.outputTime())
+ {
+ writeBoundaryPoints("boundaryPoints_" + runTime_.timeName() + ".obj");
+ }
- updateSizesAndAlignments(pointsToInsert);
+ timeCheck("After conformToSurface");
// Write the intermediate mesh, do not filter the dual faces.
if (runTime_.outputTime())
@@ -2410,6 +1693,8 @@ void Foam::conformalVoronoiMesh::move()
writeMesh(runTime_.timeName(), false);
}
+ updateSizesAndAlignments(pointsToInsert);
+
Info<< nl
<< "Total displacement = " << totalDisp << nl
<< "Total distance = " << totalDist << nl
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index 53b56da7b20d734760541d8d9b892f0dc2fafdfa..4ee0b2a42d6ab720f3a3b81d46fecf4119cb6140 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -35,6 +35,10 @@ SourceFiles
conformalVoronoiMeshI.H
conformalVoronoiMesh.C
conformalVoronoiMeshIO.C
+ conformalVoronoiMeshConformToSurface.C
+ conformalVoronoiMeshFeaturePoints.C
+ conformalVoronoiMeshFeaturePointSpecialisations.C
+ conformalVoronoiMeshCalcDualMesh.C
\*---------------------------------------------------------------------------*/
@@ -44,6 +48,7 @@ SourceFiles
#define CGAL_INEXACT
#include "CGALTriangulation3Ddefs.H"
+#include <CGAL/Spatial_sort_traits_adapter_3.h>
#include "uint.H"
#include "ulong.H"
#include "searchableSurfaces.H"
@@ -57,6 +62,8 @@ SourceFiles
#include "plane.H"
#include "SortableList.H"
#include "meshTools.H"
+#include "dynamicIndexedOctree.H"
+#include "dynamicTreeDataPoint.H"
#include "indexedOctree.H"
#include "treeDataPoint.H"
#include "unitConversion.H"
@@ -73,6 +80,7 @@ SourceFiles
#include "processorPolyPatch.H"
#include "zeroGradientFvPatchFields.H"
#include "globalIndex.H"
+#include "pointFeatureEdgesTypes.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -116,6 +124,15 @@ public:
private:
+ // Static data
+
+ static const scalar tolParallel;
+
+ static const scalar searchConeAngle;
+
+ static const scalar searchAngleOppositeSurface;
+
+
// Private data
//- The time registry of the application
@@ -146,7 +163,7 @@ private:
//- Store the feature constraining points to be reinserted after a
// triangulation clear. Maintained with relative types and indices.
- std::list<Vb> featureVertices_;
+ List<Vb> featureVertices_;
//- Storing the locations of all of the features to be conformed to.
// Single pointField required by the featurePointTree.
@@ -155,6 +172,14 @@ private:
//- Search tree for feature point locations
mutable autoPtr<indexedOctree<treeDataPoint> > featurePointTreePtr_;
+ //- Search tree for edge point locations
+ mutable autoPtr<dynamicIndexedOctree<dynamicTreeDataPoint> >
+ edgeLocationTreePtr_;
+
+ //- Search tree for surface point locations
+ mutable autoPtr<dynamicIndexedOctree<dynamicTreeDataPoint> >
+ surfacePtLocationTreePtr_;
+
//- Store locations where the cell size and alignments will be
// pre-calculated and looked up
pointField sizeAndAlignmentLocations_;
@@ -170,7 +195,7 @@ private:
//- Store the surface and feature edge conformation locations to be
// reinserted
- std::list<Vb> surfaceConformationVertices_;
+ List<Vb> surfaceConformationVertices_;
//- Method for inserting initial points. Runtime selectable.
autoPtr<initialPointsMethod> initialPointsMethod_;
@@ -235,6 +260,12 @@ private:
const Foam::point& pt
) const;
+ //- Return the square of the local surface point exclusion distance
+ inline scalar surfacePtExclusionDistanceSqr
+ (
+ const Foam::point& pt
+ ) const;
+
//- Return the square of the local surface search distance
inline scalar surfaceSearchDistanceSqr(const Foam::point& pt) const;
@@ -279,7 +310,7 @@ private:
// processors
void insertPoints
(
- std::list<Point>& points,
+ List<Point>& points,
bool distribute = true
);
@@ -396,6 +427,8 @@ private:
//- Determine and insert point groups at the feature points
void insertFeaturePoints();
+ bool edgesShareNormal(const label e1, const label e2) const;
+
//- Create point groups at convex feature points
void createConvexFeaturePoints
(
@@ -420,12 +453,24 @@ private:
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
// structure
bool createSpecialisedFeaturePoint
(
const extendedFeatureEdgeMesh& feMesh,
- label ptI,
+ const labelList& pEds,
+ const pointFeatureEdgesTypes& pFEdgeTypes,
+ const List<extendedFeatureEdgeMesh::edgeStatus>& allEdStat,
+ const label ptI,
DynamicList<Foam::point>& pts,
DynamicList<label>& indices,
DynamicList<label>& types
@@ -434,6 +479,11 @@ private:
//- Store the locations of all of the features to be conformed to
void constructFeaturePointLocations();
+ List<pointIndexHit> findSurfacePtLocationsNearFeaturePoint
+ (
+ const Foam::point& featurePoint
+ ) const;
+
//- Reinsert stored feature point defining points
void reinsertFeaturePoints(bool distribute = false);
@@ -443,13 +493,18 @@ private:
//- Check if a location is in exclusion range around a feature point
bool nearFeaturePt(const Foam::point& pt) const;
- //- Clear the entire tesselation and insert bounding points
- void reset();
+ //- Clear the entire tesselation
+ // 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
// spanning huge distances
void insertBoundingPoints();
+ //- Reinsert the bounding points
+ void reinsertBoundingPoints();
+
//- Insert the initial points into the triangulation, based on the
// initialPointsMethod
void insertInitialPoints();
@@ -468,10 +523,10 @@ private:
//- Store data for sizeAndAlignmentLocations_, storedSizes_ and
// storedAlignments_ and initialise the sizeAndAlignmentTreePtr_
- void storeSizesAndAlignments(const std::list<Point>& storePts);
+ void storeSizesAndAlignments(const List<Point>& storePts);
//- Restore the sizes and alignments if required
- void updateSizesAndAlignments(const std::list<Point>& storePts);
+ void updateSizesAndAlignments(const List<Point>& storePts);
//- Demand driven construction of octree for and alignment points
const indexedOctree<treeDataPoint>& sizeAndAlignmentTree() const;
@@ -526,6 +581,14 @@ private:
const Delaunay::Finite_vertices_iterator& vit
) const;
+ //- Return all intersections
+ bool dualCellSurfaceAllIntersections
+ (
+ const Delaunay::Finite_vertices_iterator& vit,
+ DynamicList<pointIndexHit>& info,
+ DynamicList<label>& hitSurface
+ ) const;
+
//- Return false if the line is entirely outside the current processor
// domain, true is either point is inside, or the processor domain
// bounadry is intersected (i.e. the points are box outside but the
@@ -533,6 +596,7 @@ private:
// intersect.
bool clipLineToProc
(
+ const Foam::point& pt,
Foam::point& a,
Foam::point& b
) const;
@@ -624,21 +688,67 @@ private:
label callCount = 0
) const;
+ //- Find angle between the normals of two close surface points.
+ scalar angleBetweenSurfacePoints(Foam::point pA, Foam::point pB) const;
+
+ //- Check if a surface point is near another.
+ bool nearSurfacePoint
+ (
+ pointIndexHit& pHit,
+ label& surfaceHit,
+ DynamicList<Foam::point>& existingSurfacePtLocations
+ ) const;
+
+ //- Append a point to the surface point tree and the existing list
+ bool appendToSurfacePtTree
+ (
+ const Foam::point& pt,
+ DynamicList<Foam::point>& existingSurfacePtLocations
+ ) const;
+
+ //- Append a point to the edge location tree and the existing list
+ bool appendToEdgeLocationTree
+ (
+ const Foam::point& pt,
+ DynamicList<Foam::point>& existingEdgeLocations
+ ) const;
+
+ //- Check if a point is near any feature edge points.
+ bool pointIsNearFeatureEdgeLocation(const Foam::point& pt) const;
+
+ bool pointIsNearFeatureEdgeLocation
+ (
+ const Foam::point& pt,
+ pointIndexHit& info
+ ) const;
+
+ //- Check if a point is near any surface conformation points.
+ bool pointIsNearSurfaceLocation(const Foam::point& pt) const;
+
+ bool pointIsNearSurfaceLocation
+ (
+ const Foam::point& pt,
+ pointIndexHit& info
+ ) const;
+
//- Check if a location is in the exclusion range of an existing feature
//- edge conformation location
bool nearFeatureEdgeLocation
(
- const Foam::point& pt,
- DynamicList<Foam::point>& newEdgeLocations,
- pointField& existingEdgeLocations,
- autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree
+ pointIndexHit& pHit,
+ DynamicList<Foam::point>& existingEdgeLocations
) const;
- //- Build or rebuild the edgeLocationTree
+ //- Build or rebuild the edge location tree
void buildEdgeLocationTree
(
- autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree,
- const pointField& existingEdgeLocations
+ const DynamicList<Foam::point>& existingEdgeLocations
+ ) const;
+
+ //- Build or rebuild the surface point location tree
+ void buildSurfacePtLocationTree
+ (
+ const DynamicList<Foam::point>& existingSurfacePtLocations
) const;
//- Build or rebuild the sizeAndAlignmentTree
@@ -650,8 +760,8 @@ private:
(
const Delaunay::Finite_vertices_iterator& vit,
const Foam::point& vert,
- const pointIndexHit& surfHit,
- label hitSurface,
+ const DynamicList<pointIndexHit>& surfHit,
+ const DynamicList<label>& hitSurface,
scalar surfacePtReplaceDistCoeffSqr,
scalar edgeSearchDistCoeffSqr,
DynamicList<pointIndexHit>& surfaceHits,
@@ -659,8 +769,8 @@ private:
DynamicList<pointIndexHit>& featureEdgeHits,
DynamicList<label>& featureEdgeFeaturesHit,
DynamicList<Foam::point>& newEdgeLocations,
- pointField& existingEdgeLocations,
- autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree
+ DynamicList<Foam::point>& existingEdgeLocations,
+ DynamicList<Foam::point>& existingSurfacePtLocations
) const;
//- Store the surface conformation with the indices offset to be
@@ -996,6 +1106,9 @@ public:
//- Write Delaunay points to .obj file
void writePoints(const fileName& fName, bool internalOnly) const;
+ //- Write the boundary Delaunay points to .obj file
+ void writeBoundaryPoints(const fileName& fName) const;
+
//- Write list of points to file
void writePoints
(
@@ -1045,9 +1158,93 @@ public:
// actual cell size (cbrt(actual cell volume)).
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
// protrude out of the surface beyond a tolerance.
void findRemainingProtrusionSet(const fvMesh& mesh) const;
+
+
+ //- Function inserting points into a triangulation and setting the
+ // index and type data of the point in the correct order. This is
+ // faster than inserting points individually.
+ //
+ // Adapted from a post on the CGAL lists: 2010-01/msg00004.html by
+ // Sebastien Loriot (Geometry Factory).
+ //
+ // @todo problems putting declaration in the .C file. Function
+ // prototype is not recognised.
+ template<class Triangulation, class Point_iterator>
+ void rangeInsertWithInfo
+ (
+ Point_iterator begin,
+ Point_iterator end,
+ Triangulation& T,
+ DynamicList<label>& indices,
+ DynamicList<label>& types
+ )
+ {
+ typedef std::vector
+ <
+ std::pair<const typename Triangulation::Point*, label>
+ > vectorPairPointIndex;
+
+ vectorPairPointIndex points;
+ label index = 0;
+
+ for (Point_iterator it = begin; it != end; ++it)
+ {
+ points.push_back
+ (
+ std::make_pair(&(toPoint(*it)), index++)
+ );
+ }
+
+ std::random_shuffle(points.begin(), points.end());
+
+ spatial_sort
+ (
+ points.begin(),
+ points.end(),
+ Traits_for_spatial_sort<Triangulation>()
+ );
+
+ typename Triangulation::Cell_handle hint;
+
+ for
+ (
+ typename vectorPairPointIndex::const_iterator
+ p = points.begin();
+ p != points.end();
+ ++p
+ )
+ {
+ typename Triangulation::Locate_type lt;
+ typename Triangulation::Cell_handle c;
+ label li, lj;
+
+ c = T.locate(*(p->first), lt, li, lj, hint);
+
+ typename Triangulation::Vertex_handle v
+ = T.insert(*(p->first), lt, c, li, lj);
+
+ label oldIndex = p->second;
+
+ label type = types[oldIndex];
+
+ if (type > Vb::vtFar)
+ {
+ // This is a member of a point pair, don't use the type
+ // directly (note that this routine never gets called
+ // for referredPoints so type will never be -procI)
+ type += T.number_of_vertices() - 1;
+ }
+
+ v->index() = indices[oldIndex] + T.number_of_vertices() - 1;
+ v->type() = type;
+ }
+ }
};
@@ -1059,6 +1256,10 @@ public:
#include "conformalVoronoiMeshI.H"
+//#ifdef NoRepository
+//# include "conformalVoronoiMeshTemplates.C"
+//#endif
+
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index e764923699bf3ebfa50f76032db0893d6336128a..d2c23eef51c8c066db1cf1ac04ff4007489a1445 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -79,53 +79,59 @@ void Foam::conformalVoronoiMesh::calcDualMesh
}
}
- for
- (
- Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
- vit != finite_vertices_end();
- vit++
- )
- {
- std::list<Cell_handle> cells;
- incident_cells(vit, std::back_inserter(cells));
-
- bool hasProcPt = false;
-
- for
- (
- std::list<Cell_handle>::iterator cit=cells.begin();
- cit != cells.end();
- ++cit
- )
- {
- if
- (
- !(*cit)->vertex(0)->real()
- || !(*cit)->vertex(1)->real()
- || !(*cit)->vertex(2)->real()
- || !(*cit)->vertex(3)->real()
- )
- {
- hasProcPt = true;
-
- break;
- }
- }
-
- if (hasProcPt)
- {
- for
- (
- std::list<Cell_handle>::iterator cit=cells.begin();
- cit != cells.end();
- ++cit
- )
- {
- (*cit)->filterCount() =
- cvMeshControls().filterCountSkipThreshold() + 1;
- }
- }
- }
+// REMOVED BECAUSE THIS CODE STOPS ALL FACES NEAR ANY BOUNDARY (PROC OR REAL)
+// FROM BEING FILTERED.
+//
+// for
+// (
+// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+// vit != finite_vertices_end();
+// vit++
+// )
+// {
+// std::list<Cell_handle> cells;
+// incident_cells(vit, std::back_inserter(cells));
+//
+// bool hasProcPt = false;
+//
+// for
+// (
+// std::list<Cell_handle>::iterator cit = cells.begin();
+// cit != cells.end();
+// ++cit
+// )
+// {
+// // Allow filtering if any vertices are far points. Otherwise faces
+// // with boundary points attached to a cell with a far point will not
+// // be filtered.
+// if
+// (
+// ( !(*cit)->vertex(0)->real() && !(*cit)->vertex(0)->farPoint() )
+// || ( !(*cit)->vertex(1)->real() && !(*cit)->vertex(1)->farPoint() )
+// || ( !(*cit)->vertex(2)->real() && !(*cit)->vertex(2)->farPoint() )
+// || ( !(*cit)->vertex(3)->real() && !(*cit)->vertex(3)->farPoint() )
+// )
+// {
+// hasProcPt = true;
+//
+// break;
+// }
+// }
+//
+// if (hasProcPt)
+// {
+// for
+// (
+// std::list<Cell_handle>::iterator cit = cells.begin();
+// cit != cells.end();
+// ++cit
+// )
+// {
+// (*cit)->filterCount() =
+// cvMeshControls().filterCountSkipThreshold() + 1;
+// }
+// }
+// }
PackedBoolList boundaryPts(number_of_cells(), false);
@@ -250,7 +256,7 @@ void Foam::conformalVoronoiMesh::calcDualMesh
timeCheck("End of filtering iteration");
- } while (nBadQualityFaces > nInitialBadQualityFaces);
+ } while (nBadQualityFaces > 0); //nInitialBadQualityFaces);
}
}
@@ -987,7 +993,6 @@ Foam::label Foam::conformalVoronoiMesh::collapseFaces
if (dualFace.size() < 3)
{
- // This face has been collapsed already
continue;
}
@@ -995,8 +1000,6 @@ Foam::label Foam::conformalVoronoiMesh::collapseFaces
if (maxFC > cvMeshControls().filterCountSkipThreshold())
{
- // A vertex on this face has been limited too many
- // times, skip
continue;
}
@@ -1900,7 +1903,6 @@ void Foam::conformalVoronoiMesh::indexDualVertices
)
{
// This is a boundary dual vertex
-
boundaryPts[dualVertI] = true;
}
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index e65a5beefa2daa974addc6fea9a0c5faff76c9c6..925a19c099c057fa710b77bca13264838251740f 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -25,6 +25,15 @@ License
#include "conformalVoronoiMesh.H"
#include "backgroundMeshDecomposition.H"
+#include "vectorTools.H"
+
+using namespace Foam::vectorTools;
+
+const Foam::scalar Foam::conformalVoronoiMesh::searchConeAngle
+ = Foam::cos(degToRad(30));
+
+const Foam::scalar Foam::conformalVoronoiMesh::searchAngleOppositeSurface
+ = Foam::cos(degToRad(150));
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
@@ -125,12 +134,15 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
// Initialise containers to store the edge conformation locations
DynamicList<Foam::point> newEdgeLocations;
- pointField existingEdgeLocations(0);
+ DynamicList<Foam::point> existingEdgeLocations;
+ DynamicList<Foam::point> existingSurfacePtLocations;
+
+ buildEdgeLocationTree(existingEdgeLocations);
+ buildSurfacePtLocationTree(existingSurfacePtLocations);
- autoPtr<indexedOctree<treeDataPoint> > edgeLocationTree;
// Initialise the edgeLocationTree
- buildEdgeLocationTree(edgeLocationTree, existingEdgeLocations);
+ //buildEdgeLocationTree(edgeLocationTree, existingEdgeLocations);
label initialTotalHits = 0;
@@ -167,10 +179,10 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
// Initial surface protrusion conformation - nearest surface point
{
- scalar edgeSearchDistCoeffSqr =
+ const scalar edgeSearchDistCoeffSqr =
cvMeshControls().edgeSearchDistCoeffSqrInitial(reconfMode);
- scalar surfacePtReplaceDistCoeffSqr =
+ const scalar surfacePtReplaceDistCoeffSqr =
cvMeshControls().surfacePtReplaceDistCoeffSqrInitial(reconfMode);
DynamicList<pointIndexHit> surfaceHits;
@@ -188,46 +200,51 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
{
if (vit->internalPoint())
{
- Foam::point vert(topoint(vit->point()));
- scalar searchDistanceSqr = surfaceSearchDistanceSqr(vert);
- pointIndexHit surfHit;
- label hitSurface;
+ const Foam::point vert = topoint(vit->point());
- geometryToConformTo_.findSurfaceNearest
- (
- vert,
- searchDistanceSqr,
- surfHit,
- hitSurface
- );
+ if (!positionOnThisProc(vert))
+ {
+ continue;
+ }
- if (surfHit.hit())
+ DynamicList<pointIndexHit> surfHitList;
+ DynamicList<label> hitSurfaceList;
+
+ if
+ (
+ dualCellSurfaceAllIntersections
+ (
+ vit,
+ surfHitList,
+ hitSurfaceList
+ )
+ )
{
+ // 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();
- if (dualCellSurfaceAnyIntersection(vit))
- {
- // meshTools::writeOBJ(Pout, vert);
- // meshTools::writeOBJ(Pout, surfHit.hitPoint());
- // Pout<< "l cr0 cr1" << endl;
+ // meshTools::writeOBJ(Pout, vert);
+ // meshTools::writeOBJ(Pout, surfHit.hitPoint());
+ // Pout<< "l cr0 cr1" << endl;
- addSurfaceAndEdgeHits
- (
- vit,
- vert,
- surfHit,
- hitSurface,
- surfacePtReplaceDistCoeffSqr,
- edgeSearchDistCoeffSqr,
- surfaceHits,
- hitSurfaces,
- featureEdgeHits,
- featureEdgeFeaturesHit,
- newEdgeLocations,
- existingEdgeLocations,
- edgeLocationTree
- );
- }
+ addSurfaceAndEdgeHits
+ (
+ vit,
+ vert,
+ surfHitList,
+ hitSurfaceList,
+ surfacePtReplaceDistCoeffSqr,
+ edgeSearchDistCoeffSqr,
+ surfaceHits,
+ hitSurfaces,
+ featureEdgeHits,
+ featureEdgeFeaturesHit,
+ newEdgeLocations,
+ existingEdgeLocations,
+ existingSurfacePtLocations
+ );
}
}
}
@@ -339,7 +356,8 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
|| vit->referredInternalOrBoundaryPoint()
)
{
- Foam::point vert(topoint(vit->point()));
+ const Foam::point vert = topoint(vit->point());
+
pointIndexHit surfHit;
label hitSurface;
@@ -349,12 +367,18 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
if (surfHit.hit())
{
+ DynamicList<pointIndexHit> tmpPIH;
+ tmpPIH.append(surfHit);
+
+ DynamicList<label> tmpHS;
+ tmpHS.append(hitSurface);
+
addSurfaceAndEdgeHits
(
vit,
vert,
- surfHit,
- hitSurface,
+ tmpPIH,
+ tmpHS,
surfacePtReplaceDistCoeffSqr,
edgeSearchDistCoeffSqr,
surfaceHits,
@@ -363,13 +387,14 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
featureEdgeFeaturesHit,
newEdgeLocations,
existingEdgeLocations,
- edgeLocationTree
+ existingSurfacePtLocations
);
}
}
else if (vit->ppSlave() || vit->referredExternal())
{
- Foam::point vert(topoint(vit->point()));
+ const Foam::point vert = topoint(vit->point());
+
pointIndexHit surfHit;
label hitSurface;
@@ -379,12 +404,18 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
if (surfHit.hit())
{
+ DynamicList<pointIndexHit> tmpPIH;
+ tmpPIH.append(surfHit);
+
+ DynamicList<label> tmpHS;
+ tmpHS.append(hitSurface);
+
addSurfaceAndEdgeHits
(
vit,
vert,
- surfHit,
- hitSurface,
+ tmpPIH,
+ tmpHS,
surfacePtReplaceDistCoeffSqr,
edgeSearchDistCoeffSqr,
surfaceHits,
@@ -393,7 +424,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
featureEdgeFeaturesHit,
newEdgeLocations,
existingEdgeLocations,
- edgeLocationTree
+ existingSurfacePtLocations
);
}
}
@@ -497,10 +528,10 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
if (Pstream::parRun())
{
- Foam::point a = dE0;
- Foam::point b = dE1;
+ Foam::point& a = dE0;
+ Foam::point& b = dE1;
- bool inProc = clipLineToProc(a, b);
+ bool inProc = clipLineToProc(topoint(vit->point()), a, b);
// Check for the edge passing through a surface
if
@@ -530,74 +561,185 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
}
-bool Foam::conformalVoronoiMesh::clipLineToProc
+bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
(
- Foam::point& a,
- Foam::point& b
+ const Delaunay::Finite_vertices_iterator& vit,
+ DynamicList<pointIndexHit>& infoList,
+ DynamicList<label>& hitSurfaceList
) const
{
- bool intersects = false;
+ bool flagIntersection = false;
- if (decomposition_().positionOnThisProcessor(a))
- {
- // a is inside this processor
+ std::list<Facet> facets;
+ incident_facets(vit, std::back_inserter(facets));
- if (decomposition_().positionOnThisProcessor(b))
+ for
+ (
+ std::list<Facet>::iterator fit = facets.begin();
+ fit != facets.end();
+ ++fit
+ )
+ {
+ if
+ (
+ is_infinite(fit->first)
+ || is_infinite(fit->first->neighbor(fit->second))
+ )
{
- // both a and b inside, no clip required
- intersects = true;
+ continue;
}
- else
- {
- // b is outside, clip b to bounding box.
- pointIndexHit info = decomposition_().findLine(b, a);
+ // Construct the dual edge and search for intersections of the edge
+ // with the surface
+ Foam::point dE0 = topoint(dual(fit->first));
+ Foam::point dE1 = topoint(dual(fit->first->neighbor(fit->second)));
- if (info.hit())
- {
- intersects = true;
- b = info.hitPoint();
- }
- }
- }
- else
- {
- // a is outside this processor
+ pointIndexHit infoIntersection;
+ label hitSurfaceIntersection = -1;
- if (decomposition_().positionOnThisProcessor(b))
+ if (Pstream::parRun())
{
- // b is inside this processor, clip a to processor
+ bool inProc = clipLineToProc(topoint(vit->point()), dE0, dE1);
- pointIndexHit info = decomposition_().findLine(a, b);
-
- if (info.hit())
+ if (!inProc)
{
- intersects = true;
- a = info.hitPoint();
+ continue;
}
}
- else
+
+ geometryToConformTo_.findSurfaceNearestIntersection
+ (
+ dE0,
+ dE1,
+ infoIntersection,
+ hitSurfaceIntersection
+ );
+
+ if (infoIntersection.hit())
{
- // both a and b outside, but they can still intersect the processor
+ vectorField norm(1);
+
+ allGeometry_[hitSurfaceIntersection].getNormal
+ (
+ List<pointIndexHit>(1, infoIntersection),
+ norm
+ );
+
+ const vector& n = norm[0];
+
+ const Foam::point vertex = topoint(vit->point());
+
+ const plane p(infoIntersection.hitPoint(), n);
+
+ const plane::ray r(vertex, n);
+
+ const scalar d = p.normalIntersect(r);
- pointIndexHit info = decomposition_().findLine(a, b);
+ const Foam::point newPoint = vertex + d*n;
+
+ pointIndexHit info;
+ label hitSurface = -1;
+
+ geometryToConformTo_.findSurfaceNearest
+ (
+ newPoint,
+ surfaceSearchDistanceSqr(newPoint),
+ info,
+ hitSurface
+ );
+
+ bool rejectPoint = false;
if (info.hit())
{
- intersects = true;
- a = info.hitPoint();
+ if (!infoList.empty())
+ {
+ forAll(infoList, hitI)
+ {
+ // Reject point if the point is already added
+ if (infoList[hitI].index() == info.index())
+ {
+ rejectPoint = true;
+ break;
+ }
- info = decomposition_().findLine(b, a);
+ const Foam::point& p = infoList[hitI].hitPoint();
- if (info.hit())
- {
- b = info.hitPoint();
+ const scalar separationDistance
+ = mag(p - info.hitPoint());
+
+ const scalar minSepDist
+ = sqr(cvMeshControls().removalDistCoeff()
+ *targetCellSize(p));
+
+ // Reject the point if it is too close to another
+ // surface point.
+ // Could merge the points?
+ if (separationDistance < minSepDist)
+ {
+ rejectPoint = true;
+ break;
+ }
+ }
}
}
+
+ // The normal ray from the vertex will not always result in a hit
+ // because another surface may be in the way.
+ if (!rejectPoint && info.hit())
+ {
+ flagIntersection = true;
+ infoList.append(info);
+ hitSurfaceList.append(hitSurface);
+ }
}
}
- return intersects;
+ return flagIntersection;
+}
+
+
+bool Foam::conformalVoronoiMesh::clipLineToProc
+(
+ const Foam::point& pt,
+ Foam::point& a,
+ Foam::point& b
+) const
+{
+ bool inProc = false;
+
+ pointIndexHit findAnyIntersection = decomposition_().findLine(a, b);
+
+ if (!findAnyIntersection.hit())
+ {
+ pointIndexHit info = decomposition_().findLine(a, pt);
+
+ if (!info.hit())
+ {
+ inProc = true;
+ }
+ else
+ {
+ inProc = false;
+ }
+ }
+ else
+ {
+ pointIndexHit info = decomposition_().findLine(a, pt);
+
+ if (!info.hit())
+ {
+ inProc = true;
+ b = findAnyIntersection.hitPoint();
+ }
+ else
+ {
+ inProc = true;
+ a = findAnyIntersection.hitPoint();
+ }
+ }
+
+ return inProc;
}
@@ -1213,13 +1355,13 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
std::list<Facet> facets;
incident_facets(vit, std::back_inserter(facets));
- Foam::point vert(topoint(vit->point()));
+ const Foam::point vert = topoint(vit->point());
scalar maxProtrusionDistance = maxSurfaceProtrusion(vert);
for
(
- std::list<Facet>::iterator fit=facets.begin();
+ std::list<Facet>::iterator fit = facets.begin();
fit != facets.end();
++fit
)
@@ -1230,7 +1372,7 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
&& !is_infinite(fit->first->neighbor(fit->second))
)
{
- Foam::point edgeMid =
+ const Foam::point edgeMid =
0.5
*(
topoint(dual(fit->first))
@@ -1260,7 +1402,7 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
const vector& n = norm[0];
- scalar normalProtrusionDistance =
+ const scalar normalProtrusionDistance =
(edgeMid - surfHit.hitPoint()) & n;
if (normalProtrusionDistance > maxProtrusionDistance)
@@ -1311,7 +1453,7 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
for
(
- std::list<Facet>::iterator fit=facets.begin();
+ std::list<Facet>::iterator fit = facets.begin();
fit != facets.end();
++fit
)
@@ -1592,56 +1734,339 @@ void Foam::conformalVoronoiMesh::limitDisplacement
}
-bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
+Foam::scalar Foam::conformalVoronoiMesh::angleBetweenSurfacePoints
(
- const Foam::point& pt,
- DynamicList<Foam::point>& newEdgeLocations,
- pointField& existingEdgeLocations,
- autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree
+ Foam::point pA,
+ Foam::point pB
) const
{
- scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
+ pointIndexHit pAhit;
+ label pAsurfaceHit = -1;
- // 0.01 and 1000 determined from speed tests, varying the indexedOctree
- // rebuild frequency and balance of additions to queries.
+ const scalar searchDist = 5.0*targetCellSize(pA);
- if
+ geometryToConformTo_.findSurfaceNearest
(
- newEdgeLocations.size()
- >= max(0.01*existingEdgeLocations.size(), 1000)
- )
+ pA,
+ searchDist,
+ pAhit,
+ pAsurfaceHit
+ );
+
+ vectorField norm(1);
+
+ allGeometry_[pAsurfaceHit].getNormal
+ (
+ List<pointIndexHit>(1, pAhit),
+ norm
+ );
+
+ const vector nA = norm[0];
+
+ pointIndexHit pBhit;
+ label pBsurfaceHit = -1;
+
+ geometryToConformTo_.findSurfaceNearest
+ (
+ pB,
+ searchDist,
+ pBhit,
+ pBsurfaceHit
+ );
+
+ allGeometry_[pBsurfaceHit].getNormal
+ (
+ List<pointIndexHit>(1, pBhit),
+ norm
+ );
+
+ const vector nB = norm[0];
+
+ return vectorTools::cosPhi(nA, nB);
+}
+
+
+bool Foam::conformalVoronoiMesh::nearSurfacePoint
+(
+ pointIndexHit& pHit,
+ label& surfaceHit,
+ DynamicList<Foam::point>& existingSurfacePtLocations
+) const
+{
+ const Foam::point& pt = pHit.hitPoint();
+
+ pointIndexHit closePoint;
+
+ const bool closeToSurfacePt = pointIsNearSurfaceLocation(pt, closePoint);
+
+ if (closeToSurfacePt)
+ {
+ const scalar cosAngle
+ = angleBetweenSurfacePoints(pt, closePoint.hitPoint());
+
+ // @todo make this tolerance run-time selectable?
+ if (cosAngle < searchAngleOppositeSurface)
+ {
+ pointIndexHit pCloseHit;
+ label pCloseSurfaceHit = -1;
+
+ const scalar searchDist = targetCellSize(closePoint.hitPoint());
+
+ geometryToConformTo_.findSurfaceNearest
+ (
+ closePoint.hitPoint(),
+ searchDist,
+ pCloseHit,
+ pCloseSurfaceHit
+ );
+
+ vectorField norm(1);
+
+ allGeometry_[pCloseSurfaceHit].getNormal
+ (
+ List<pointIndexHit>(1, pCloseHit),
+ norm
+ );
+
+ const vector nA = norm[0];
+
+ pointIndexHit oppositeHit;
+ label oppositeSurfaceHit = -1;
+
+ geometryToConformTo_.findSurfaceNearestIntersection
+ (
+ closePoint.hitPoint() + SMALL*nA,
+ closePoint.hitPoint() + mag(pt - closePoint.hitPoint())*nA,
+ oppositeHit,
+ oppositeSurfaceHit
+ );
+
+ if (oppositeHit.hit())
+ {
+ // Replace point
+ pHit = oppositeHit;
+ surfaceHit = oppositeSurfaceHit;
+
+ Foam::point newPt = oppositeHit.hitPoint();
+
+ appendToSurfacePtTree(newPt, existingSurfacePtLocations);
+
+ return !closeToSurfacePt;
+
+ if (debug)
+ {
+ Info<< "Point " << pt
+ << " is close to " << closePoint.hitPoint()
+ << " so will be moved to " << oppositeHit.hitPoint()
+ << endl;
+ }
+ }
+ }
+ }
+ else
{
- existingEdgeLocations.append(newEdgeLocations);
+ appendToSurfacePtTree(pt, existingSurfacePtLocations);
+ }
+
+ return closeToSurfacePt;
+}
+
+
+bool Foam::conformalVoronoiMesh::appendToSurfacePtTree
+(
+ const Foam::point& pt,
+ DynamicList<Foam::point>& existingSurfacePtLocations
+) const
+{
+ label startIndex = existingSurfacePtLocations.size();
+
+ existingSurfacePtLocations.append(pt);
+
+ label endIndex = existingSurfacePtLocations.size();
+
+ return surfacePtLocationTreePtr_().insert(startIndex, endIndex);
+}
+
+
+bool Foam::conformalVoronoiMesh::appendToEdgeLocationTree
+(
+ const Foam::point& pt,
+ DynamicList<Foam::point>& existingEdgeLocations
+) const
+{
+ label startIndex = existingEdgeLocations.size();
+
+ existingEdgeLocations.append(pt);
+
+ label endIndex = existingEdgeLocations.size();
+
+ return edgeLocationTreePtr_().insert(startIndex, endIndex);
+}
+
+
+bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdgeLocation
+(
+ const Foam::point& pt
+) const
+{
+ const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
+
+ pointIndexHit info
+ = edgeLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
+
+ return info.hit();
+}
+
+
+bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdgeLocation
+(
+ const Foam::point& pt,
+ pointIndexHit& info
+) const
+{
+ const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
+
+ info = edgeLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
+
+ return info.hit();
+}
+
+
+bool Foam::conformalVoronoiMesh::pointIsNearSurfaceLocation
+(
+ const Foam::point& pt
+) const
+{
+ pointIndexHit info;
+
+ pointIsNearSurfaceLocation(pt, info);
+
+ return info.hit();
+}
+
+
+bool Foam::conformalVoronoiMesh::pointIsNearSurfaceLocation
+(
+ const Foam::point& pt,
+ pointIndexHit& info
+) const
+{
+ const scalar exclusionRangeSqr = surfacePtExclusionDistanceSqr(pt);
+
+ info = surfacePtLocationTreePtr_().findNearest(pt, exclusionRangeSqr);
+
+ return info.hit();
+}
+
+
+bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
+(
+ pointIndexHit& pHit,
+ DynamicList<Foam::point>& existingEdgeLocations
+) const
+{
+ const Foam::point pt = pHit.hitPoint();
+
+ const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
- buildEdgeLocationTree(edgeLocationTree, existingEdgeLocations);
+ pointIndexHit info;
- newEdgeLocations.clear();
+ bool closeToFeatureEdge = pointIsNearFeatureEdgeLocation(pt, info);
+
+ if (!closeToFeatureEdge)
+ {
+ appendToEdgeLocationTree(pt, existingEdgeLocations);
}
else
{
- // Search for the nearest point in newEdgeLocations.
- // Searching here first, because the intention is that the value of
- // newEdgeLocationsSizeLimit should make this faster by design.
+ // Check if the edge location that the new edge location is near to
+ // "might" be on a different edge. If so, add it anyway.
+ pointIndexHit edgeHit;
+ label featureHit = -1;
+
+ geometryToConformTo_.findEdgeNearest
+ (
+ pt,
+ exclusionRangeSqr,
+ edgeHit,
+ featureHit
+ );
+
+ const extendedFeatureEdgeMesh& eMesh
+ = geometryToConformTo_.features()[featureHit];
- if (min(magSqr(newEdgeLocations - pt)) <= exclusionRangeSqr)
+ const vector& edgeDir = eMesh.edgeDirections()[edgeHit.index()];
+
+ const vector lineBetweenPoints = pt - info.hitPoint();
+
+ const scalar cosAngle = vectorTools::cosPhi(edgeDir, lineBetweenPoints);
+
+ // Allow the point to be added if it is almost at right angles to the
+ // other point. Also check it is not the same point.
+// Info<< cosAngle<< " "
+// << radToDeg(acos(cosAngle)) << " "
+// << searchConeAngle << " "
+// << radToDeg(acos(searchConeAngle)) << endl;
+ if
+ (
+ mag(cosAngle) < searchConeAngle
+ && mag(lineBetweenPoints) > SMALL
+ )
{
- return true;
+ closeToFeatureEdge = false;
+ appendToEdgeLocationTree(pt, existingEdgeLocations);
}
}
+ return closeToFeatureEdge;
+
// Searching for the nearest point in existingEdgeLocations using the
// indexedOctree
- pointIndexHit info = edgeLocationTree().findNearest(pt, exclusionRangeSqr);
+ // Average the points...
+// if (info.hit())
+// {
+// Foam::point newPt = 0.5*(info.hitPoint() + pt);
+//
+// pHit.setPoint(newPt);
+//
+// //boolList toRemove(existingEdgeLocations.size(), false);
+//
+// forAll(existingEdgeLocations, pI)
+// {
+// if (pI == info.index())
+// {
+// //toRemove[pI] = true;
+// edgeLocationTree.remove(pI);
+// }
+// }
+////
+//// pointField newExistingEdgeLocations(existingEdgeLocations.size());
+////
+//// label count = 0;
+//// forAll(existingEdgeLocations, pI)
+//// {
+//// if (toRemove[pI] == false)
+//// {
+//// newExistingEdgeLocations[count++] = existingEdgeLocations[pI];
+//// }
+//// }
+////
+//// newExistingEdgeLocations.resize(count);
+////
+//// existingEdgeLocations = newExistingEdgeLocations;
+////
+//// existingEdgeLocations.append(newPt);
+//
+// return !info.hit();
+// }
- return info.hit();
}
void Foam::conformalVoronoiMesh::buildEdgeLocationTree
(
- autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree,
- const pointField& existingEdgeLocations
+ const DynamicList<Foam::point>& existingEdgeLocations
) const
{
treeBoundBox overallBb
@@ -1652,14 +2077,41 @@ void Foam::conformalVoronoiMesh::buildEdgeLocationTree
overallBb.min() -= Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
overallBb.max() += Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
- edgeLocationTree.reset
+ edgeLocationTreePtr_.reset
(
- new indexedOctree<treeDataPoint>
+ new dynamicIndexedOctree<dynamicTreeDataPoint>
(
- treeDataPoint(existingEdgeLocations),
+ dynamicTreeDataPoint(existingEdgeLocations),
overallBb, // overall search domain
- 10, // max levels
- 10.0, // maximum ratio of cubes v.s. cells
+ 10, // max levels, n/a
+ 20.0, // maximum ratio of cubes v.s. cells
+ 100.0 // max. duplicity; n/a since no bounding boxes.
+ )
+ );
+}
+
+
+void Foam::conformalVoronoiMesh::buildSurfacePtLocationTree
+(
+ const DynamicList<Foam::point>& existingSurfacePtLocations
+) const
+{
+ treeBoundBox overallBb
+ (
+ geometryToConformTo_.globalBounds().extend(rndGen_, 1e-4)
+ );
+
+ overallBb.min() -= Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+ overallBb.max() += Foam::point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+ surfacePtLocationTreePtr_.reset
+ (
+ new dynamicIndexedOctree<dynamicTreeDataPoint>
+ (
+ dynamicTreeDataPoint(existingSurfacePtLocations),
+ overallBb, // overall search domain
+ 10, // max levels, n/a
+ 20.0, // maximum ratio of cubes v.s. cells
100.0 // max. duplicity; n/a since no bounding boxes.
)
);
@@ -1702,8 +2154,8 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
(
const Delaunay::Finite_vertices_iterator& vit,
const Foam::point& vert,
- const pointIndexHit& surfHit,
- label hitSurface,
+ const DynamicList<pointIndexHit>& surfHit,
+ const DynamicList<label>& hitSurface,
scalar surfacePtReplaceDistCoeffSqr,
scalar edgeSearchDistCoeffSqr,
DynamicList<pointIndexHit>& surfaceHits,
@@ -1711,104 +2163,123 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
DynamicList<pointIndexHit>& featureEdgeHits,
DynamicList<label>& featureEdgeFeaturesHit,
DynamicList<Foam::point>& newEdgeLocations,
- pointField& existingEdgeLocations,
- autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree
+ DynamicList<Foam::point>& existingEdgeLocations,
+ DynamicList<Foam::point>& existingSurfacePtLocations
) const
{
bool keepSurfacePoint = true;
- if (nearFeaturePt(surfHit.hitPoint()))
+ const scalar cellSize = targetCellSize(vert);
+ const scalar cellSizeSqr = sqr(cellSize);
+
+ forAll(surfHit, sI)
{
- keepSurfacePoint = false;
+ pointIndexHit surfHitI = surfHit[sI];
+ label hitSurfaceI = hitSurface[sI];
- // If the triggering Vertex is part of a feature point, allow it to
- // conform to the surface
- if (vit->index() < startOfInternalPoints_)
+ if (!surfHitI.hit())
{
- surfaceHits.append(surfHit);
-
- hitSurfaces.append(hitSurface);
+ continue;
}
- }
- List<pointIndexHit> edHits;
+ bool isNearFeaturePt = nearFeaturePt(surfHitI.hitPoint());
- labelList featuresHit;
+ bool isNearSurfacePt = nearSurfacePoint
+ (
+ surfHitI,
+ hitSurfaceI,
+ existingSurfacePtLocations
+ );
- scalar targetCellSizeSqr = sqr(targetCellSize(vert));
+ if (isNearFeaturePt || isNearSurfacePt)
+ {
+ keepSurfacePoint = false;
- geometryToConformTo_.findEdgeNearestByType
- (
- surfHit.hitPoint(),
- edgeSearchDistCoeffSqr*targetCellSizeSqr,
- edHits,
- featuresHit
- );
+ // If the triggering Vertex is part of a feature point, allow it to
+ // conform to the surface.
+ // @todo Is this needed?! Shouldn't be any feature points here...
+ if (vit->index() < startOfInternalPoints_)
+ {
+ surfaceHits.append(surfHitI);
- // Gather edge locations but do not add them to newEdgeLocations inside the
- // loop as they will prevent nearby edge locations of different types being
- // conformed to.
+ hitSurfaces.append(hitSurfaceI);
+ }
+ }
- DynamicList<Foam::point> currentEdgeLocations;
+ List<List<pointIndexHit> > edHitsByFeature;
- forAll(edHits, i)
- {
- const pointIndexHit& edHit(edHits[i]);
+ labelList featuresHit;
- label featureHit = featuresHit[i];
+ const scalar searchRadiusSqr = edgeSearchDistCoeffSqr*cellSizeSqr;
- if (edHit.hit())
+ geometryToConformTo_.findAllNearestEdges
+ (
+ surfHitI.hitPoint(),
+ searchRadiusSqr,
+ edHitsByFeature,
+ featuresHit
+ );
+
+ // Gather edge locations but do not add them to newEdgeLocations inside
+ // the loop as they will prevent nearby edge locations of different
+ // types being conformed to.
+
+ DynamicList<Foam::point> currentEdgeLocations;
+
+ forAll(edHitsByFeature, i)
{
- if (!positionOnThisProc(edHit.hitPoint()))
- {
- continue;
- }
+ const label featureHit = featuresHit[i];
+
+ List<pointIndexHit>& edHits = edHitsByFeature[i];
- if (!nearFeaturePt(edHit.hitPoint()))
+ forAll(edHits, eHitI)
{
- if
- (
- magSqr(edHit.hitPoint() - surfHit.hitPoint())
- < surfacePtReplaceDistCoeffSqr*targetCellSizeSqr
- )
+ pointIndexHit& edHit = edHits[eHitI];
+
+ if (!nearFeaturePt(edHit.hitPoint()) && keepSurfacePoint)
{
- // If the point is within a given distance of a feature
- // edge, give control to edge control points instead, this
- // will prevent "pits" forming.
+ if
+ (
+ magSqr(edHit.hitPoint() - surfHitI.hitPoint())
+ < surfacePtReplaceDistCoeffSqr*cellSizeSqr
+ )
+ {
+ // If the point is within a given distance of a feature
+ // edge, give control to edge control points instead,
+ // this will prevent "pits" forming.
- keepSurfacePoint = false;
- }
+ keepSurfacePoint = false;
- if
- (
- !nearFeatureEdgeLocation
+ // NEED TO REMOVE FROM THE SURFACE TREE...
+ }
+
+ if
(
- edHit.hitPoint(),
- newEdgeLocations,
- existingEdgeLocations,
- edgeLocationTree
+ !nearFeatureEdgeLocation
+ (
+ edHit,
+ existingEdgeLocations
+ )
)
- )
- {
- // Do not place edge control points too close to a feature
- // point or existing edge control points
+ {
+ // Do not place edge control points too close to a
+ // feature point or existing edge control points
- featureEdgeHits.append(edHit);
- featureEdgeFeaturesHit.append(featureHit);
+ featureEdgeHits.append(edHit);
+ featureEdgeFeaturesHit.append(featureHit);
- currentEdgeLocations.append(edHit.hitPoint());
+ currentEdgeLocations.append(edHit.hitPoint());
+ }
}
}
}
- }
-
- newEdgeLocations.append(currentEdgeLocations);
- if (keepSurfacePoint)
- {
- surfaceHits.append(surfHit);
+ if (keepSurfacePoint)
+ {
+ surfaceHits.append(surfHitI);
- hitSurfaces.append(hitSurface);
+ hitSurfaces.append(hitSurfaceI);
+ }
}
}
@@ -1819,6 +2290,9 @@ void Foam::conformalVoronoiMesh::storeSurfaceConformation()
surfaceConformationVertices_.clear();
+ // Use a temporary dynamic list to speed up insertion.
+ DynamicList<Vb> tempSurfaceVertices(number_of_vertices()/10);
+
for
(
Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
@@ -1835,7 +2309,7 @@ void Foam::conformalVoronoiMesh::storeSurfaceConformation()
&& vit->index() >= startOfInternalPoints_
)
{
- surfaceConformationVertices_.push_back
+ tempSurfaceVertices.append
(
Vb
(
@@ -1847,6 +2321,10 @@ void Foam::conformalVoronoiMesh::storeSurfaceConformation()
}
}
+ tempSurfaceVertices.shrink();
+
+ surfaceConformationVertices_.transfer(tempSurfaceVertices);
+
Info<< " Stored "
<< returnReduce
(
@@ -1868,7 +2346,7 @@ void Foam::conformalVoronoiMesh::reinsertSurfaceConformation()
for
(
- std::list<Vb>::iterator vit=surfaceConformationVertices_.begin();
+ List<Vb>::iterator vit=surfaceConformationVertices_.begin();
vit != surfaceConformationVertices_.end();
++vit
)
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
index 95040bf865bd13ca748580a50147f96c4cb4c3a3..644d82c96130f880433b385389a4b3f3bfeba823 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C
@@ -24,33 +24,20 @@ License
\*---------------------------------------------------------------------------*/
#include "conformalVoronoiMesh.H"
+#include "vectorTools.H"
+
+using namespace Foam::vectorTools;
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
-bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
+Foam::List<Foam::extendedFeatureEdgeMesh::edgeStatus>
+Foam::conformalVoronoiMesh::calcPointFeatureEdgesTypes
(
const extendedFeatureEdgeMesh& feMesh,
- label ptI,
- DynamicList<Foam::point>& pts,
- DynamicList<label>& indices,
- DynamicList<label>& types
+ const labelList& pEds,
+ pointFeatureEdgesTypes& pFEdgeTypes
)
{
- 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());
forAll(pEds, i)
@@ -65,293 +52,388 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
{
case extendedFeatureEdgeMesh::EXTERNAL:
{
- nExternal++;
+ pFEdgeTypes.nExternal++;
break;
}
case extendedFeatureEdgeMesh::INTERNAL:
{
- nInternal++;
+ pFEdgeTypes.nInternal++;
break;
}
case extendedFeatureEdgeMesh::FLAT:
{
- nFlat++;
+ pFEdgeTypes.nFlat++;
break;
}
case extendedFeatureEdgeMesh::OPEN:
{
- nOpen++;
+ pFEdgeTypes.nOpen++;
break;
}
case extendedFeatureEdgeMesh::MULTIPLE:
{
- nMultiple++;
+ pFEdgeTypes.nMultiple++;
break;
}
case extendedFeatureEdgeMesh::NONE:
{
+ pFEdgeTypes.nNonFeature++;
break;
}
}
}
- if (nExternal == 2 && nInternal == 1)
+ if (debug)
{
- // if (!positionOnThisProc(pt))
- // {
- // continue;
- // }
+ Info<< pFEdgeTypes << endl;
+ }
- // Info<< "nExternal == 2 && nInternal == 1" << endl;
+ return allEdStat;
+}
- // const Foam::point& featPt = feMesh.points()[ptI];
- // scalar ppDist = pointPairDistance(featPt);
+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 vectorField& normals = feMesh.normals();
+ const Foam::point& featPt = feMesh.points()[ptI];
- // const labelListList& edgeNormals = feMesh.edgeNormals();
+ if (!positionOnThisProc(featPt))
+ {
+ return false;
+ }
- // label concaveEdgeI = pEds
- // [
- // findIndex(allEdStat, extendedFeatureEdgeMesh::INTERNAL)
- // ];
+ const label initialNumOfPoints = pts.size();
- // // // Find which planes are joined to the concave edge
+ const scalar ppDist = pointPairDistance(featPt);
- // // List<label> concaveEdgePlanes(2,label(-1));
+ const vectorField& normals = feMesh.normals();
- // // label concaveEdgeI = concaveEdges[0];
+ const labelListList& edgeNormals = feMesh.edgeNormals();
- // // // Pick up the two faces adjacent to the concave feature edge
- // // const labelList& eFaces = qSurf_.edgeFaces()[concaveEdgeI];
+ label concaveEdgeI = -1;
+ labelList convexEdgesI(2, -1);
+ label nConvex = 0;
- // // label faceA = eFaces[0];
+ forAll(pEds, i)
+ {
+ const extendedFeatureEdgeMesh::edgeStatus& eS = allEdStat[i];
- // // vector nA = qSurf_.faceNormals()[faceA];
+ 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);
+ }
+ }
- // // scalar maxNormalDotProduct = -SMALL;
+ const vector& concaveEdgePlaneANormal =
+ normals[edgeNormals[concaveEdgeI][0]];
- // // forAll(uniquePlaneNormals, uPN)
- // // {
- // // scalar normalDotProduct = nA & uniquePlaneNormals[uPN];
+ const vector& concaveEdgePlaneBNormal =
+ normals[edgeNormals[concaveEdgeI][1]];
- // // if (normalDotProduct > maxNormalDotProduct)
- // // {
- // // maxNormalDotProduct = normalDotProduct;
+ // 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
+ );
- // // concaveEdgePlanes[0] = uPN;
- // // }
- // // }
+ plane planeB
+ (
+ featPt + ppDist*concaveEdgePlaneBNormal,
+ concaveEdgePlaneBNormal
+ );
- // // label faceB = eFaces[1];
- // // vector nB = qSurf_.faceNormals()[faceB];
+ const vector& concaveEdgeDir = feMesh.edgeDirection
+ (
+ concaveEdgeI,
+ ptI
+ );
- // // maxNormalDotProduct = -SMALL;
+ // Todo,needed later but want to get rid of this.
+ const Foam::point concaveEdgeLocalFeatPt = featPt + ppDist*concaveEdgeDir;
- // // forAll(uniquePlaneNormals, uPN)
- // // {
- // // scalar normalDotProduct = nB & uniquePlaneNormals[uPN];
+ // Finding the nearest point on the intersecting line to the edge
+ // point. Floating point errors often occur using planePlaneIntersect
- // // if (normalDotProduct > maxNormalDotProduct)
- // // {
- // // maxNormalDotProduct = normalDotProduct;
+ plane planeF(concaveEdgeLocalFeatPt, concaveEdgeDir);
- // // concaveEdgePlanes[1] = uPN;
- // // }
- // // }
+ const Foam::point concaveEdgeExternalPt = planeF.planePlaneIntersect
+ (
+ planeA,
+ planeB
+ );
- // // const vector& concaveEdgePlaneANormal =
- // // uniquePlaneNormals[concaveEdgePlanes[0]];
+ // Redefine planes to be on the feature surfaces to project through
- // // const vector& concaveEdgePlaneBNormal =
- // // uniquePlaneNormals[concaveEdgePlanes[1]];
+ planeA = plane(featPt, concaveEdgePlaneANormal);
- // // label convexEdgesPlaneI;
+ planeB = plane(featPt, concaveEdgePlaneBNormal);
- // // if (findIndex(concaveEdgePlanes, 0) == -1)
- // // {
- // // convexEdgesPlaneI = 0;
- // // }
- // // else if (findIndex(concaveEdgePlanes, 1) == -1)
- // // {
- // // convexEdgesPlaneI = 1;
- // // }
- // // else
- // // {
- // // convexEdgesPlaneI = 2;
- // // }
+ const Foam::point internalPtA =
+ concaveEdgeExternalPt
+ - 2.0*planeA.distance(concaveEdgeExternalPt)
+ *concaveEdgePlaneANormal;
- // // const vector& convexEdgesPlaneNormal =
- // // uniquePlaneNormals[convexEdgesPlaneI];
+ pts.append(internalPtA);
+ indices.append(0);
+ types.append(4);
- // // const edge& concaveEdge = edges[concaveEdgeI];
+ const Foam::point internalPtB =
+ concaveEdgeExternalPt
+ - 2.0*planeB.distance(concaveEdgeExternalPt)
+ *concaveEdgePlaneBNormal;
- // // // Check direction of edge, if the feature point is at the end()
- // // // the reverse direction.
+ pts.append(internalPtB);
+ indices.append(0);
+ types.append(3);
- // // scalar edgeDirection = 1.0;
+ // Add the external points
- // // if (ptI == concaveEdge.end())
- // // {
- // // edgeDirection *= -1.0;
- // // }
+ Foam::point externalPtD;
+ Foam::point externalPtE;
+ vector convexEdgePlaneCNormal(0,0,0);
+ vector convexEdgePlaneDNormal(0,0,0);
- // const vector& concaveEdgePlaneANormal =
- // normals[edgeNormals[concaveEdgeI][0]];
+ const labelList& concaveEdgeNormals = edgeNormals[concaveEdgeI];
+ const labelList& convexEdgeANormals = edgeNormals[convexEdgesI[0]];
+ const labelList& convexEdgeBNormals = edgeNormals[convexEdgesI[1]];
- // const vector& concaveEdgePlaneBNormal =
- // normals[edgeNormals[concaveEdgeI][1]];
+ forAll(concaveEdgeNormals, edgeNormalI)
+ {
+ bool convexEdgeA = false;
+ bool convexEdgeB = false;
- // // 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
- // );
+ forAll(convexEdgeANormals, edgeAnormalI)
+ {
+ const vector& concaveNormal
+ = normals[concaveEdgeNormals[edgeNormalI]];
+ const vector& convexNormal
+ = normals[convexEdgeANormals[edgeAnormalI]];
+
+ Info<< "Angle between vectors = "
+ << degAngleBetween(concaveNormal, convexNormal) << endl;
+
+ // Need a looser tolerance, because sometimes adjacent triangles
+ // on the same surface will be slightly out of alignment.
+ if (areParallel(concaveNormal, convexNormal, tolParallel))
+ {
+ convexEdgeA = true;
+ }
+ }
- // plane planeB
- // (
- // featPt + ppDist*concaveEdgePlaneBNormal,
- // concaveEdgePlaneBNormal
- // );
+ forAll(convexEdgeBNormals, edgeBnormalI)
+ {
+ const vector& concaveNormal
+ = normals[concaveEdgeNormals[edgeNormalI]];
+ const vector& convexNormal
+ = normals[convexEdgeBNormals[edgeBnormalI]];
+
+ Info<< "Angle between vectors = "
+ << degAngleBetween(concaveNormal, convexNormal) << endl;
+
+ // Need a looser tolerance, because sometimes adjacent triangles
+ // on the same surface will be slightly out of alignment.
+ if (areParallel(concaveNormal, convexNormal, tolParallel))
+ {
+ convexEdgeB = true;
+ }
+ }
- // const vector& concaveEdgeDir = feMesh.edgeDirection
- // (
- // concaveEdgeI,
- // ptI
- // );
+ if ((convexEdgeA && convexEdgeB) || (!convexEdgeA && !convexEdgeB))
+ {
+ WarningIn
+ (
+ "Foam::conformalVoronoiMesh"
+ "::createSpecialisedFeaturePoint"
+ )
+ << "Both or neither of the convex edges share the concave "
+ << "edge's normal."
+ << " convexEdgeA = " << convexEdgeA
+ << " convexEdgeB = " << convexEdgeB
+ << endl;
+
+ // Remove points that have just been added before returning
+ for (label i = 0; i < 2; ++i)
+ {
+ pts.remove();
+ indices.remove();
+ types.remove();
+ }
+
+ return false;
+ }
- // Foam::point concaveEdgeLocalFeatPt = featPt + ppDist*concaveEdgeDir;
+ if (convexEdgeA)
+ {
+ forAll(convexEdgeANormals, edgeAnormalI)
+ {
+ const vector& concaveNormal
+ = normals[concaveEdgeNormals[edgeNormalI]];
+ const vector& convexNormal
+ = normals[convexEdgeANormals[edgeAnormalI]];
+
+ if
+ (
+ !areParallel(concaveNormal, convexNormal, tolParallel)
+ )
+ {
+ convexEdgePlaneCNormal = convexNormal;
+
+ plane planeC(featPt, convexEdgePlaneCNormal);
+
+ externalPtD =
+ internalPtA
+ + 2.0*planeC.distance(internalPtA)
+ *convexEdgePlaneCNormal;
+
+ pts.append(externalPtD);
+ indices.append(0);
+ types.append(-2);
+ }
+ }
+ }
- // // 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);
+ if (convexEdgeB)
+ {
+ forAll(convexEdgeBNormals, edgeBnormalI)
+ {
+ const vector& concaveNormal
+ = normals[concaveEdgeNormals[edgeNormalI]];
+ const vector& convexNormal
+ = normals[convexEdgeBNormals[edgeBnormalI]];
+
+ if
+ (
+ !areParallel(concaveNormal, convexNormal, tolParallel)
+ )
+ {
+ convexEdgePlaneDNormal = convexNormal;
+
+ plane planeD(featPt, convexEdgePlaneDNormal);
+
+ externalPtE =
+ internalPtB
+ + 2.0*planeD.distance(internalPtB)
+ *convexEdgePlaneDNormal;
+
+ pts.append(externalPtE);
+ indices.append(0);
+ types.append(-2);
+ }
+ }
+ }
+ }
- // Foam::point externalPtE =
- // internalPtB
- // + 2*planeC.distance(internalPtB)*convexEdgesPlaneNormal;
+ pts.append(concaveEdgeExternalPt);
+ indices.append(0);
+ types.append(-4);
- // 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);
+ const scalar totalAngle = radToDeg
+ (
+ constant::mathematical::pi
+ + radAngleBetween(concaveEdgePlaneANormal, concaveEdgePlaneBNormal)
+ );
- // 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;
+ if (totalAngle > cvMeshControls().maxQuadAngle())
+ {
+ // Add additional mitering points
+ //scalar angleSign = 1.0;
+
+ Info<<convexEdgePlaneCNormal << " " <<convexEdgePlaneDNormal << endl;
+
+ vector convexEdgesPlaneNormal = 0.5*(convexEdgePlaneCNormal + convexEdgePlaneDNormal);
+ 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;
+
+ pts.append(externalPtG);
+ indices.append(0);
+ types.append(-1);
+ }
- return false;
+ if (debug)
+ {
+ for (label ptI = initialNumOfPoints; ptI < pts.size(); ++ptI)
+ {
+ Info<< "Point " << ptI << " : ";
+ meshTools::writeOBJ(Info, pts[ptI]);
+ }
+ }
+
+ return true;
}
- else if (nExternal == 1 && nInternal == 2)
+ else if (pFEdgesTypes.nExternal == 1 && pFEdgesTypes.nInternal == 2)
{
// Info<< "nExternal == 1 && nInternal == 2" << endl;
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
new file mode 100644
index 0000000000000000000000000000000000000000..cb2a0c9340c3bd5c88c13c794c28f2c459400ff2
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C
@@ -0,0 +1,923 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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"
+#include "vectorTools.H"
+
+using namespace Foam::vectorTools;
+
+// * * * * * * * * * * * * 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 (areParallel(nA, nB))
+ {
+ // 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 (areParallel(nA, nB))
+ {
+ // 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 + radAngleBetween(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
+ const 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
+ (
+ 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();
+ featureVertices_.setSize(pointsToInsert.size());
+
+ forAll(pointsToInsert, pI)
+ {
+ featureVertices_[pI] =
+ Vb(toPoint(pointsToInsert[pI]), indices[pI], types[pI]);
+ }
+ }
+ else
+ {
+ for
+ (
+ 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;
+}
+
+
+bool Foam::conformalVoronoiMesh::edgesShareNormal
+(
+ const label e1,
+ const label e2
+) const
+{
+ const PtrList<extendedFeatureEdgeMesh>& feMeshes
+ (
+ geometryToConformTo_.features()
+ );
+
+ forAll(feMeshes, i)
+ {
+ const extendedFeatureEdgeMesh& feMesh(feMeshes[i]);
+
+ const vectorField& e1normals = feMesh.edgeNormals(e1);
+ const vectorField& e2normals = feMesh.edgeNormals(e2);
+
+ forAll(e1normals, nI1)
+ {
+ forAll(e2normals, nI2)
+ {
+ if (degAngleBetween(e1normals[nI1], e2normals[nI2]) < 1)
+ {
+ return true;
+ }
+ }
+ }
+ }
+
+ return false;
+}
+
+
+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;
+ }
+
+ const vectorField& featPtNormals = feMesh.featurePointNormals(ptI);
+
+ const scalar ppDist = - pointPairDistance(featPt);
+
+ vector cornerNormal = sum(featPtNormals);
+ cornerNormal /= mag(cornerNormal);
+
+ Foam::point internalPt = featPt + ppDist*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;
+ }
+
+ const vectorField& featPtNormals = feMesh.featurePointNormals(ptI);
+ const scalar ppDist = pointPairDistance(featPt);
+
+ vector cornerNormal = sum(featPtNormals);
+ cornerNormal /= mag(cornerNormal);
+
+ Foam::point externalPt = featPt + ppDist*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];
+
+ const plane planeN = plane(featPt, n);
+
+ const 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::createFeaturePoints
+//(
+// 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.mixedStart();
+// ++ptI
+// )
+// {
+// const Foam::point& featPt = feMesh.points()[ptI];
+//
+// if (!positionOnThisProc(featPt))
+// {
+// continue;
+// }
+//
+// const scalar searchRadiusSqr = 5.0*targetCellSize(featPt);
+//
+// labelList indices =
+// surfacePtLocationTreePtr_().findSphere(featPt, searchRadiusSqr);
+//
+// pointField nearestSurfacePoints(indices.size());
+//
+// forAll(indices, pI)
+// {
+// nearestSurfacePoints[pI] =
+// surfaceConformationVertices_[indices[pI]];
+// }
+//
+// forAll()
+//
+// // Now find the nearest points within the edge cones.
+//
+// // Calculate preliminary surface point locations
+//
+//
+// }
+// }
+//}
+
+
+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);
+
+// createFeaturePoints(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();
+ featureVertices_.setSize(pts.size());
+
+ forAll(pts, pI)
+ {
+ featureVertices_[pI] = 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);
+}
+
+
+Foam::List<Foam::pointIndexHit>
+Foam::conformalVoronoiMesh::findSurfacePtLocationsNearFeaturePoint
+(
+ const Foam::point& featurePoint
+) const
+{
+ DynamicList<pointIndexHit> dynPointList;
+
+ const scalar searchRadiusSqr = 3*targetCellSize(featurePoint);
+
+ labelList surfacePtList = surfacePtLocationTreePtr_().findSphere
+ (
+ featurePoint,
+ searchRadiusSqr
+ );
+
+ forAll(surfacePtList, elemI)
+ {
+ label index = surfacePtList[elemI];
+
+ const Foam::point& p
+ = surfacePtLocationTreePtr_().shapes().shapePoints()[index];
+
+ pointIndexHit nearHit(true, p, index);
+
+ dynPointList.append(nearHit);
+ }
+
+ return dynPointList.shrink();
+}
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
index 0cad8731f7f4c1480452ab34d3dbcf1c3e2bc4d0..a18b173ecd3c64890c1d9e24d65da8787ff0d297 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshI.H
@@ -186,6 +186,20 @@ inline Foam::scalar Foam::conformalVoronoiMesh::featureEdgeExclusionDistanceSqr
}
+inline Foam::scalar Foam::conformalVoronoiMesh::surfacePtExclusionDistanceSqr
+(
+ const Foam::point& pt
+) const
+{
+ return
+ sqr
+ (
+ targetCellSize(pt)
+ *cvMeshControls().surfacePtExclusionDistanceCoeff()
+ );
+}
+
+
inline Foam::scalar Foam::conformalVoronoiMesh::surfaceSearchDistanceSqr
(
const Foam::point& pt
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index 6430eab8ce21a2864899b96aba3c99b047cee3fe..bc9d04dfb1e79ca32216f53185c7443267f9996d 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -175,6 +175,30 @@ void Foam::conformalVoronoiMesh::writePoints
}
+void Foam::conformalVoronoiMesh::writeBoundaryPoints
+(
+ const fileName& fName
+) const
+{
+ OFstream str(runTime_.path()/fName);
+
+ Pout<< nl << "Writing boundary points to " << str.name() << endl;
+
+ for
+ (
+ Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+ vit != finite_vertices_end();
+ ++vit
+ )
+ {
+ if (!vit->internalPoint())
+ {
+ meshTools::writeOBJ(str, topoint(vit->point()));
+ }
+ }
+}
+
+
void Foam::conformalVoronoiMesh::writePoints
(
const fileName& fName,
@@ -597,29 +621,10 @@ void Foam::conformalVoronoiMesh::writeMesh
<< exit(FatalError);
}
- // pointIOField cellCs
- // (
- // IOobject
- // (
- // "cellCentres",
- // mesh.pointsInstance(),
- // polyMesh::meshSubDir,
- // mesh,
- // IOobject::NO_READ,
- // IOobject::AUTO_WRITE
- // ),
- // cellCentres
- // );
-
- // Info<< nl
- // << "Writing " << cellCs.name()
- // << " to " << cellCs.instance()
- // << endl;
-
- // cellCs.write();
-
writeCellSizes(mesh);
+ writeCellCentres(mesh);
+
findRemainingProtrusionSet(mesh);
}
@@ -805,6 +810,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
(
const fvMesh& mesh
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
index 7e5ea9112e416db3819440ad603034d86cbaf463..7ab4f10853427619285ce87430f1ccf14cd31e42 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
@@ -33,6 +33,10 @@ Description
An indexed form of CGAL::Triangulation_cell_base_3<K> used to keep
track of the Delaunay cells (tets) in the tessellation.
+SourceFiles
+ indexedCellI.H
+ indexedCell.C
+
\*---------------------------------------------------------------------------*/
#ifndef indexedCell_H
@@ -119,248 +123,72 @@ public:
};
- indexedCell()
- :
- Cb(),
- index_(ctFar),
- filterCount_(0)
- {}
-
-
- 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_;
- }
+ // Constructors
- inline bool farCell() const
- {
- return index_ == ctFar;
- }
-
- inline int& filterCount()
- {
- return filterCount_;
- }
+ inline indexedCell();
-
- inline int filterCount() const
- {
- return filterCount_;
- }
-
-
- //- Is the Delaunay cell real, i.e. any real vertex
- inline bool real() const
- {
- return
+ inline indexedCell
(
- this->vertex(0)->real()
- || this->vertex(1)->real()
- || this->vertex(2)->real()
- || this->vertex(3)->real()
+ Vertex_handle v0,
+ Vertex_handle v1,
+ Vertex_handle v2,
+ Vertex_handle v3
);
- }
-
- //- Does the Dual vertex form part of a processor patch
- inline bool parallelDualVertex() const
- {
- return
+ inline indexedCell
(
- this->vertex(0)->referred()
- || this->vertex(1)->referred()
- || this->vertex(2)->referred()
- || this->vertex(3)->referred()
+ Vertex_handle v0,
+ Vertex_handle v1,
+ Vertex_handle v2,
+ Vertex_handle v3,
+ Cell_handle n0,
+ Cell_handle n1,
+ Cell_handle n2,
+ Cell_handle n3
);
- }
- //- 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.
+ // Member Functions
- int masterProc = Foam::Pstream::nProcs() + 1;
+ inline int& cellIndex();
- 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());
- }
- }
+ inline int cellIndex() const;
- return masterProc;
- }
+ inline bool farCell() const;
+ inline int& filterCount();
- inline Foam::FixedList<Foam::label, 4> processorsAttached() const
- {
- if (!parallelDualVertex())
- {
- return Foam::FixedList<Foam::label, 4>(Foam::Pstream::myProcNo());
- }
+ inline int filterCount() const;
- 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;
- }
-
-
- // 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()
- );
- }
- }
+ //- Is the Delaunay cell real, i.e. any real vertex
+ inline bool real() const;
- // 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]);
- }
- }
- }
+ //- Does the Dual vertex form part of a processor patch
+ inline bool parallelDualVertex() const;
- return tVGI;
- }
+ //- Does the Dual vertex form part of a processor patch
+ inline Foam::label dualVertexMasterProc() const;
+ inline Foam::FixedList<Foam::label, 4> processorsAttached() const;
- // 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
+ //- 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)->internalOrBoundaryPoint()
- || this->vertex(1)->internalOrBoundaryPoint()
- || this->vertex(2)->internalOrBoundaryPoint()
- || this->vertex(3)->internalOrBoundaryPoint()
- );
- }
+ const Foam::globalIndex& globalDelaunayVertexIndices
+ ) const;
+ //- 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;
- // Is the Delaunay cell real or referred (or mixed), i.e. all vertices form
- // part of the real or referred internal or boundary definition
- inline bool anyInternalOrBoundaryDualVertex() const
- {
- return
- (
- this->vertex(0)->anyInternalOrBoundaryPoint()
- || this->vertex(1)->anyInternalOrBoundaryPoint()
- || this->vertex(2)->anyInternalOrBoundaryPoint()
- || this->vertex(3)->anyInternalOrBoundaryPoint()
- );
- }
-
+ //- Is the Delaunay cell real or referred (or mixed), i.e. all vertices
+ // form part of the real or referred internal or boundary definition
+ inline bool anyInternalOrBoundaryDualVertex() const;
- // 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
- {
- 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()
- )
- );
- }
+ //- 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;
// Info
@@ -387,6 +215,8 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#include "indexedCellI.H"
+
#ifdef NoRepository
# include "indexedCell.C"
#endif
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H
new file mode 100644
index 0000000000000000000000000000000000000000..5fbd378f9d9a220845245127070d9624f2ac9bba
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H
@@ -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 * * * * * * * * * * * * * //
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
index 0dae13edb3cad456c640f67406c0533eaed96741..246020fd172bbeb127089dc9bcc085c0aa7908e2 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
@@ -33,6 +33,10 @@ Description
An indexed form of CGAL::Triangulation_vertex_base_3<K> used to keep
track of the Delaunay vertices in the tessellation.
+SourceFiles
+ indexedVertexI.H
+ indexedVertex.C
+
\*---------------------------------------------------------------------------*/
#ifndef indexedVertex_H
@@ -88,12 +92,15 @@ class indexedVertex
// Lowest numbered is inside one (master).
int type_;
- // Required alignment of the dual cell of this vertex
+ //- Required alignment of the dual cell of this vertex
Foam::tensor alignment_;
- // Target size of the dual cell of this vertex
+ //- Target size of the dual cell of this vertex
Foam::scalar targetCellSize_;
+ //- Specify whether the vertex is fixed or movable.
+ bool vertexFixed_;
+
public:
@@ -104,6 +111,12 @@ public:
vtFar = INT_MIN + 2
};
+ enum vertexMotion
+ {
+ fixed = 0,
+ movable = 1
+ };
+
typedef typename Vb::Vertex_handle Vertex_handle;
typedef typename Vb::Cell_handle Cell_handle;
typedef typename Vb::Point Point;
@@ -116,268 +129,105 @@ public:
};
- indexedVertex()
- :
- Vb(),
- index_(vtInternal),
- type_(vtInternal),
- alignment_(),
- targetCellSize_(0.0)
- {}
-
-
- indexedVertex(const Point& p)
- :
- Vb(p),
- index_(vtInternal),
- type_(vtInternal),
- alignment_(),
- targetCellSize_(0.0)
- {}
-
-
- indexedVertex(const Point& p, int index, int type)
- :
- Vb(p),
- index_(index),
- type_(type),
- alignment_(),
- targetCellSize_(0.0)
- {}
-
-
- indexedVertex(const Point& p, Cell_handle f)
- :
- Vb(f, p),
- index_(vtInternal),
- type_(vtInternal),
- alignment_(),
- targetCellSize_(0.0)
- {}
-
-
- indexedVertex(Cell_handle f)
- :
- Vb(f),
- index_(vtInternal),
- type_(vtInternal),
- alignment_(),
- targetCellSize_(0.0)
- {}
-
-
- int& index()
- {
- return index_;
- }
+ // Constructors
+ inline indexedVertex();
- int index() const
- {
- return index_;
- }
+ inline indexedVertex(const Point& p);
+ inline indexedVertex(const Point& p, int index, int type);
- int& type()
- {
- return type_;
- }
+ inline indexedVertex(const Point& p, Cell_handle f);
+ inline indexedVertex(Cell_handle f);
- int type() const
- {
- return type_;
- }
+ // Member Functions
- inline Foam::tensor& alignment()
- {
- return alignment_;
- }
+ inline int& index();
+ inline int index() const;
- inline const Foam::tensor& alignment() const
- {
- return alignment_;
- }
+ inline int& type();
+ inline int type() const;
- inline Foam::scalar& targetCellSize()
- {
- return targetCellSize_;
- }
+ inline Foam::tensor& alignment();
+ inline const Foam::tensor& alignment() const;
- inline Foam::scalar targetCellSize() const
- {
- return targetCellSize_;
- }
+ inline Foam::scalar& targetCellSize();
+ inline Foam::scalar targetCellSize() const;
- inline bool uninitialised() const
- {
- return type_ == vtInternal && index_ == vtInternal;
- }
+ inline bool uninitialised() const;
+ //- Is point a far-point
+ inline bool farPoint() const;
- //- Is point a far-point
- inline bool farPoint() const
- {
- return type_ == vtFar;
- }
+ //- Is point internal, i.e. not on boundary
+ inline bool internalPoint() const;
+ //- Is point internal, i.e. not on boundary, external query.
+ inline static bool internalPoint(int type);
- //- Is point internal, i.e. not on boundary
- inline bool internalPoint() const
- {
- return internalPoint(type_);
- }
+ // is this a referred vertex
+ inline bool referred() const;
+ // is this a referred internal or boundary vertex
+ inline bool referredInternalOrBoundaryPoint() const;
- //- Is point internal, i.e. not on boundary, external query.
- inline static bool internalPoint(int type)
- {
- return type <= vtInternal;
- }
+ // is this a referred external (pair slave) vertex
+ inline bool referredExternal() const;
+ // 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;
- // is this a referred vertex
- inline bool referred() const
- {
- return (type_ < 0 && type_ > vtFar);
- }
-
-
- // is this a referred internal or boundary vertex
- inline bool referredInternalOrBoundaryPoint() const
- {
- return referred() && index_ >= 0;
- }
+ // For referred vertices, what is the original processor index
+ inline int procIndex() const;
+ inline static int encodeProcIndex(int procI);
- // is this a referred external (pair slave) vertex
- inline bool referredExternal() const
- {
- return referred() && index_ < 0;
- }
-
+ //- Set the point to be internal
+ inline void setInternal();
- // 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();
- }
+ //- Is point internal and near the boundary
+ inline bool nearBoundary() const;
+ //- Set the point to be near the boundary
+ inline void setNearBoundary();
- // For referred vertices, what is the original processor index
- inline int procIndex() const
- {
- if (referred())
- {
- return -(type_ + 1);
- }
- else
- {
- return -1;
- }
- }
+ //- Either master or slave of pointPair.
+ inline bool pairPoint() const;
+ //- Master of a pointPair is the lowest numbered one.
+ inline bool ppMaster() const;
- inline static int encodeProcIndex(int procI)
- {
- return -(procI + 1);
- }
+ //- Master of a pointPair is the lowest numbered one, external query.
+ inline static bool ppMaster(int index, int type);
+ //- Slave of a pointPair is the highest numbered one.
+ inline bool ppSlave() const;
- //- Set the point to be internal
- inline void setInternal()
- {
- type_ = vtInternal;
- }
+ //- Either original internal point or master of pointPair.
+ inline bool internalOrBoundaryPoint() const;
+ //- Either original internal point or master of pointPair.
+ // External query.
+ inline static bool internalOrBoundaryPoint(int index, int type);
- //- Is point internal and near the boundary
- inline bool nearBoundary() const
- {
- return type_ == vtNearBoundary;
- }
+ //- Is point near the boundary or part of the boundary definition
+ inline bool nearOrOnBoundary() const;
+ //- Either a real or referred internal or boundary point
+ inline bool anyInternalOrBoundaryPoint() const;
- //- 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();
- }
+ //- Is the vertex fixed or movable
+ inline bool isVertexFixed() const;
+ //- Fix the vertex so that it can't be moved
+ inline void setVertexFixed() const;
// inline void operator=(const Delaunay::Finite_vertices_iterator vit)
// {
@@ -414,6 +264,8 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#include "indexedVertexI.H"
+
#ifdef NoRepository
# include "indexedVertex.C"
#endif
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H
new file mode 100644
index 0000000000000000000000000000000000000000..c02ec4df2009489efa9fe91221a666b9f9ae0050
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H
@@ -0,0 +1,340 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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),
+ vertexFixed_(false)
+{}
+
+
+template<class Gt, class Vb>
+inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(const Point& p)
+:
+ Vb(p),
+ index_(vtInternal),
+ type_(vtInternal),
+ alignment_(),
+ targetCellSize_(0.0),
+ vertexFixed_(false)
+{}
+
+
+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),
+ vertexFixed_(false)
+{}
+
+
+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),
+ vertexFixed_(false)
+{}
+
+
+template<class Gt, class Vb>
+inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(Cell_handle f)
+:
+ Vb(f),
+ index_(vtInternal),
+ type_(vtInternal),
+ alignment_(),
+ targetCellSize_(0.0),
+ vertexFixed_(false)
+{}
+
+
+// * * * * * * * * * * * * * * * 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();
+}
+
+
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::isVertexFixed() const
+{
+ return vertexFixed_;
+}
+
+
+template<class Gt, class Vb>
+inline void CGAL::indexedVertex<Gt, Vb>::setVertexFixed() const
+{
+ vertexFixed_ = true;
+}
+
+
+// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/pointFeatureEdgesTypes.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/pointFeatureEdgesTypes.H
new file mode 100644
index 0000000000000000000000000000000000000000..3009f385b256556a202e55b03ff4cc286ba912a7
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/pointFeatureEdgesTypes.H
@@ -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
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
index 5e8a8d729746c8bf0c29a6203f56ca15a1e92521..b91d40121d48f23d8a1ff077c9ce679983ed6218 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.C
@@ -730,6 +730,49 @@ void Foam::conformationSurfaces::findEdgeNearestByType
}
+void Foam::conformationSurfaces::findAllNearestEdges
+(
+ const point& sample,
+ const scalar searchRadiusSqr,
+ List<List<pointIndexHit> >& edgeHitsByFeature,
+ List<label>& featuresHit
+) const
+{
+ // Initialise
+ //featuresHit.setSize(features_.size());
+ //featuresHit = -1;
+ //edgeHitsByFeature.setSize(features_.size());
+
+ // Work arrays
+ List<pointIndexHit> hitInfo(extendedFeatureEdgeMesh::nEdgeTypes);
+
+ forAll(features_, testI)
+ {
+ features_[testI].allNearestFeatureEdges
+ (
+ sample,
+ searchRadiusSqr,
+ hitInfo
+ );
+
+ bool anyHit = false;
+ forAll(hitInfo, hitI)
+ {
+ if (hitInfo[hitI].hit())
+ {
+ anyHit = true;
+ }
+ }
+
+ if (anyHit)
+ {
+ edgeHitsByFeature.append(hitInfo);
+ featuresHit.append(testI);
+ }
+ }
+}
+
+
void Foam::conformationSurfaces::writeFeatureObj(const fileName& prefix) const
{
OFstream ftStr(runTime_.time().path()/prefix + "_allFeatures.obj");
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
index b541bed0032027034eaf4a4d77a73df8af49fdf6..fe883fbaa5cae003313dc9a381ea16a5734a535b 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
@@ -166,7 +166,7 @@ public:
//- Check if point is closer to the surfaces to conform to than
// testDistSqr, in which case return false, otherwise assess in or
- // outside and erturn a result depending on the testForInside flag
+ // outside and return a result depending on the testForInside flag
Field<bool> wellInOutSide
(
const pointField& samplePts,
@@ -228,7 +228,7 @@ public:
label& hitSurface
) const;
- //- Find the nearest point to the sample and return it to the
+ //- Find the nearest point to the sample and return it to the
// pointIndexHit
void findSurfaceNearest
(
@@ -270,7 +270,17 @@ public:
const point& sample,
scalar nearestDistSqr,
List<pointIndexHit>& edgeHit,
- List<label>& featureHit
+ List<label>& featuresHit
+ ) const;
+
+ //- Find the nearest points on each feature edge that is within
+ // a given distance from the sample point
+ void findAllNearestEdges
+ (
+ const point& sample,
+ const scalar searchRadiusSqr,
+ List<List<pointIndexHit> >& edgeHitsByFeature,
+ List<label>& featuresHit
) const;
//- Find which patch is intersected by the line from one point to
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
index ce1994e7bd5014a50be06a2a96b3e71d91ed5577..9b59ce9b9b5d0b3a024b0e87445585e46c37ff9c 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.C
@@ -30,11 +30,9 @@ License
Foam::cvControls::cvControls
(
- const conformalVoronoiMesh& cvMesh,
const dictionary& cvMeshDict
)
:
- cvMesh_(cvMesh),
cvMeshDict_(cvMeshDict)
{
// Surface conformation controls
@@ -61,6 +59,11 @@ Foam::cvControls::cvControls
surfDict.lookup("featureEdgeExclusionDistanceCoeff")
);
+ specialiseFeaturePoints_ = Switch
+ (
+ surfDict.lookupOrDefault<Switch>("specialiseFeaturePoints", false)
+ );
+
surfaceSearchDistanceCoeff_ = readScalar
(
surfDict.lookup("surfaceSearchDistanceCoeff")
@@ -96,6 +99,11 @@ Foam::cvControls::cvControls
coarseDict.subDict("iteration")
);
+ surfacePtExclusionDistanceCoeff_ = readScalar
+ (
+ coarseInitialDict.lookup("surfacePtExclusionDistanceCoeff")
+ );
+
edgeSearchDistCoeffSqr_coarse_initial_ = sqr
(
readScalar
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
index 8b00b9a163297a944da282cf663e16d40f30e0f6..f6bdb755592a316533cef66e2bcc94710818a0a6 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControls.H
@@ -54,9 +54,6 @@ class cvControls
{
// Private data
- //- Reference to the conformalVoronoiMesh holding this cvControls object
- const conformalVoronoiMesh& cvMesh_;
-
//- Reference to the cvMeshDict
const dictionary& cvMeshDict_;
@@ -80,6 +77,9 @@ class cvControls
// fraction of the local target cell size
scalar featureEdgeExclusionDistanceCoeff_;
+ //- Switch for using specialised feature points
+ Switch specialiseFeaturePoints_;
+
//- Surface search distance coefficient - fraction of the local
// target cell size
scalar surfaceSearchDistanceCoeff_;
@@ -99,6 +99,11 @@ class cvControls
// Controls for coarse surface conformation
+ //- Distance to an existing surface conformation point location
+ // within which other surface point locations are excluded
+ // - fraction of the local target cell size
+ scalar surfacePtExclusionDistanceCoeff_;
+
//- Distance to search for feature edges near to
// surface protrusions - fraction of the local target
// cell size. Coarse conformation, initial protrusion
@@ -291,7 +296,6 @@ public:
//- Construct from references to conformalVoronoiMesh and dictionary
cvControls
(
- const conformalVoronoiMesh& cvMesh,
const dictionary& cvMeshDict
);
@@ -319,6 +323,12 @@ public:
//- Return the featureEdgeExclusionDistanceCoeff
inline scalar featureEdgeExclusionDistanceCoeff() const;
+ //- Return the surfacePtExclusionDistanceCoeff
+ inline scalar surfacePtExclusionDistanceCoeff() const;
+
+ //- Return whether to use specialised feature points
+ inline Switch specialiseFeaturePoints() const;
+
//- Return the surfaceSearchDistanceCoeff
inline scalar surfaceSearchDistanceCoeff() const;
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
index 5fd33640195e6e8cdccb0731a178387b280b88c3..dc3f391100bd0ee1a28a8007041e8a7812ffa4f5 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cvControls/cvControlsI.H
@@ -54,6 +54,15 @@ inline Foam::scalar Foam::cvControls::featureEdgeExclusionDistanceCoeff() const
return featureEdgeExclusionDistanceCoeff_;
}
+inline Foam::scalar Foam::cvControls::surfacePtExclusionDistanceCoeff() const
+{
+ return surfacePtExclusionDistanceCoeff_;
+}
+
+inline Foam::Switch Foam::cvControls::specialiseFeaturePoints() const
+{
+ return specialiseFeaturePoints_;
+}
inline Foam::scalar Foam::cvControls::surfaceSearchDistanceCoeff() const
{
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.C
index ec7a8ed29d8e0c7c30d8684d17dab152e7b5ca61..76d68aa073f753f848bb98628d47ca6cff78d8c7 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.C
@@ -42,12 +42,10 @@ defineRunTimeSelectionTable(faceAreaWeightModel, dictionary);
faceAreaWeightModel::faceAreaWeightModel
(
const word& type,
- const dictionary& relaxationDict,
- const conformalVoronoiMesh& cvMesh
+ const dictionary& relaxationDict
)
:
dictionary(relaxationDict),
- cvMesh_(cvMesh),
coeffDict_(subDict(type + "Coeffs"))
{}
@@ -56,8 +54,7 @@ faceAreaWeightModel::faceAreaWeightModel
autoPtr<faceAreaWeightModel> faceAreaWeightModel::New
(
- const dictionary& relaxationDict,
- const conformalVoronoiMesh& cvMesh
+ const dictionary& relaxationDict
)
{
word faceAreaWeightModelTypeName
@@ -85,7 +82,7 @@ autoPtr<faceAreaWeightModel> faceAreaWeightModel::New
<< exit(FatalError);
}
- return autoPtr<faceAreaWeightModel>(cstrIter()(relaxationDict, cvMesh));
+ return autoPtr<faceAreaWeightModel>(cstrIter()(relaxationDict));
}
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.H
index 0eaf7b14f6650ea66e05c57960d6c2068deac2fa..3e88f5bb74f64768956caf551750a2f0b78712f7 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/faceAreaWeightModel/faceAreaWeightModel.H
@@ -39,7 +39,6 @@ SourceFiles
#define faceAreaWeightModel_H
#include "point.H"
-#include "conformalVoronoiMesh.H"
#include "dictionary.H"
#include "autoPtr.H"
#include "runTimeSelectionTables.H"
@@ -62,9 +61,6 @@ protected:
// Protected data
- //- Reference to the conformalVoronoiMesh holding this cvControls object
- const conformalVoronoiMesh& cvMesh_;
-
//- Method coeffs dictionary
dictionary coeffDict_;
@@ -93,10 +89,9 @@ public:
faceAreaWeightModel,
dictionary,
(
- const dictionary& faceAreaWeightDict,
- const conformalVoronoiMesh& cvMesh
+ const dictionary& faceAreaWeightDict
),
- (faceAreaWeightDict, cvMesh)
+ (faceAreaWeightDict)
);
@@ -106,8 +101,7 @@ public:
faceAreaWeightModel
(
const word& type,
- const dictionary& faceAreaWeightDict,
- const conformalVoronoiMesh& cvMesh
+ const dictionary& faceAreaWeightDict
);
@@ -116,8 +110,7 @@ public:
//- Return a reference to the selected faceAreaWeightModel
static autoPtr<faceAreaWeightModel> New
(
- const dictionary& faceAreaWeightDict,
- const conformalVoronoiMesh& cvMesh
+ const dictionary& faceAreaWeightDict
);
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.C
index 2f7df3d3532da09d9bf65cf875a729214071620b..0fa6025c998b6e69a77b03bce8cbd75240a77afe 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.C
@@ -45,11 +45,10 @@ addToRunTimeSelectionTable
piecewiseLinearRamp::piecewiseLinearRamp
(
- const dictionary& faceAreaWeightDict,
- const conformalVoronoiMesh& cvMesh
+ const dictionary& faceAreaWeightDict
)
:
- faceAreaWeightModel(typeName, faceAreaWeightDict, cvMesh),
+ faceAreaWeightModel(typeName, faceAreaWeightDict),
lAF_(readScalar(coeffDict().lookup("lowerAreaFraction"))),
uAF_(readScalar(coeffDict().lookup("upperAreaFraction")))
{}
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.H
index 0c0ba587f32212b6dde84cb8c24a77d9681468e7..14f6a929bb9e9591830084f8c4e6d238f1e48868 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/faceAreaWeightModel/piecewiseLinearRamp/piecewiseLinearRamp.H
@@ -71,8 +71,7 @@ public:
//- Construct from components
piecewiseLinearRamp
(
- const dictionary& faceAreaWeightDict,
- const conformalVoronoiMesh& cvMesh
+ const dictionary& faceAreaWeightDict
);
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C
index bdc314ce64862453aaaf52fbaa1eda7f1dd2e344..84926910dfd50a9a09da5ebf8805711a3dea1749 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.C
@@ -177,7 +177,7 @@ bool Foam::autoDensity::combinedWellInside
void Foam::autoDensity::recurseAndFill
(
- std::list<Vb::Point>& initialPoints,
+ DynamicList<Vb::Point>& initialPoints,
const treeBoundBox& bb,
label levelLimit,
word recursionName
@@ -272,7 +272,7 @@ void Foam::autoDensity::recurseAndFill
bool Foam::autoDensity::fillBox
(
- std::list<Vb::Point>& initialPoints,
+ DynamicList<Vb::Point>& initialPoints,
const treeBoundBox& bb,
bool overlapping
) const
@@ -677,7 +677,7 @@ bool Foam::autoDensity::fillBox
// Pout<< "Final volume probability break accept"
// << endl;
- initialPoints.push_back
+ initialPoints.append
(
Vb::Point(p.x(), p.y(), p.z())
);
@@ -688,7 +688,7 @@ bool Foam::autoDensity::fillBox
break;
}
- initialPoints.push_back(Vb::Point(p.x(), p.y(), p.z()));
+ initialPoints.append(Vb::Point(p.x(), p.y(), p.z()));
volumeAdded += localVolume;
}
@@ -800,7 +800,7 @@ bool Foam::autoDensity::fillBox
// Pout<< "Final volume probability break accept"
// << endl;
- initialPoints.push_back
+ initialPoints.append
(
Vb::Point(p.x(), p.y(), p.z())
);
@@ -811,7 +811,7 @@ bool Foam::autoDensity::fillBox
break;
}
- initialPoints.push_back(Vb::Point(p.x(), p.y(), p.z()));
+ initialPoints.append(Vb::Point(p.x(), p.y(), p.z()));
volumeAdded += localVolume;
}
@@ -849,10 +849,7 @@ autoDensity::autoDensity
:
initialPointsMethod(typeName, initialPointsDict, cvMesh),
globalTrialPoints_(0),
- minCellSizeLimit_
- (
- detailsDict().lookupOrDefault<scalar>("minCellSizeLimit", 0.0)
- ),
+ minCellSizeLimit_(readScalar(detailsDict().lookup("minCellSizeLimit"))),
minLevels_(readLabel(detailsDict().lookup("minLevels"))),
maxSizeRatio_(readScalar(detailsDict().lookup("maxSizeRatio"))),
volRes_(readLabel(detailsDict().lookup("sampleResolution"))),
@@ -882,7 +879,7 @@ autoDensity::autoDensity
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-std::list<Vb::Point> autoDensity::initialPoints() const
+List<Vb::Point> autoDensity::initialPoints() const
{
treeBoundBox hierBB;
@@ -902,7 +899,18 @@ std::list<Vb::Point> autoDensity::initialPoints() const
);
}
- std::list<Vb::Point> initialPoints;
+ // Initialise size of points list.
+ const scalar volumeBoundBox = Foam::pow3(hierBB.typDim());
+ const scalar volumeSmallestCell = Foam::pow3(minCellSizeLimit_);
+
+ const int initialPointEstimate
+ = min
+ (
+ static_cast<int>(volumeBoundBox/(volumeSmallestCell + SMALL)/10),
+ 1e6
+ );
+
+ DynamicList<Vb::Point> initialPoints(initialPointEstimate);
Info<< nl << " " << typeName << endl;
@@ -919,6 +927,8 @@ std::list<Vb::Point> autoDensity::initialPoints() const
"recursionBox"
);
+ initialPoints.shrink();
+
label nInitialPoints = initialPoints.size();
if (Pstream::parRun())
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.H
index 905a3acdf6e87d1e7b25827479c2c6b0cf018c77..4b9063c5cc518aabe479cddc4a46fbc1e3d9d2d6 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/autoDensity/autoDensity.H
@@ -116,7 +116,7 @@ private:
//- Descend into octants of the supplied bound box
void recurseAndFill
(
- std::list<Vb::Point>& initialPoints,
+ DynamicList<Vb::Point>& initialPoints,
const treeBoundBox& bb,
label levelLimit,
word recursionName
@@ -127,7 +127,7 @@ private:
// in favour of further recursion.
bool fillBox
(
- std::list<Vb::Point>& initialPoints,
+ DynamicList<Vb::Point>& initialPoints,
const treeBoundBox& bb,
bool overlapping
) const;
@@ -156,7 +156,7 @@ public:
// Member Functions
//- Return the initial points for the conformalVoronoiMesh
- virtual std::list<Vb::Point> initialPoints() const;
+ virtual List<Vb::Point> initialPoints() const;
};
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/bodyCentredCubic/bodyCentredCubic.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/bodyCentredCubic/bodyCentredCubic.C
index 0376f1389ac6f8e56b3f6b52f5269423159691a5..f41760311f35c36fad1390f785afab990052522d 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/bodyCentredCubic/bodyCentredCubic.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/bodyCentredCubic/bodyCentredCubic.C
@@ -56,7 +56,7 @@ bodyCentredCubic::bodyCentredCubic
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-std::list<Vb::Point> bodyCentredCubic::initialPoints() const
+List<Vb::Point> bodyCentredCubic::initialPoints() const
{
boundBox bb;
@@ -91,7 +91,7 @@ std::list<Vb::Point> bodyCentredCubic::initialPoints() const
scalar pert = randomPerturbationCoeff_*cmptMin(delta);
- std::list<Vb::Point> initialPoints;
+ DynamicList<Vb::Point> initialPoints(ni*nj*nk/10);
for (label i = 0; i < ni; i++)
{
@@ -177,13 +177,13 @@ std::list<Vb::Point> bodyCentredCubic::initialPoints() const
{
const point& p(points[i]);
- initialPoints.push_back(Vb::Point(p.x(), p.y(), p.z()));
+ initialPoints.append(Vb::Point(p.x(), p.y(), p.z()));
}
}
}
}
- return initialPoints;
+ return initialPoints.shrink();
}
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/bodyCentredCubic/bodyCentredCubic.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/bodyCentredCubic/bodyCentredCubic.H
index 434c3280c3c2d5f2d29e1b25455dbdf953f27650..c3d1217d6bd58d942263bfd9cd3baaa4067c89dd 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/bodyCentredCubic/bodyCentredCubic.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/bodyCentredCubic/bodyCentredCubic.H
@@ -89,7 +89,7 @@ public:
// Member Functions
//- Return the initial points for the conformalVoronoiMesh
- virtual std::list<Vb::Point> initialPoints() const;
+ virtual List<Vb::Point> initialPoints() const;
};
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/faceCentredCubic/faceCentredCubic.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/faceCentredCubic/faceCentredCubic.C
index c9ade503cb13832b3f8dcd49975b8eb4b7e551c9..cafccb55747b6dd52c77a6b4d9c9c9a292e0e9d0 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/faceCentredCubic/faceCentredCubic.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/faceCentredCubic/faceCentredCubic.C
@@ -56,7 +56,7 @@ faceCentredCubic::faceCentredCubic
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-std::list<Vb::Point> faceCentredCubic::initialPoints() const
+List<Vb::Point> faceCentredCubic::initialPoints() const
{
boundBox bb;
@@ -91,7 +91,7 @@ std::list<Vb::Point> faceCentredCubic::initialPoints() const
scalar pert = randomPerturbationCoeff_*cmptMin(delta);
- std::list<Vb::Point> initialPoints;
+ DynamicList<Vb::Point> initialPoints(ni*nj*nk/10);
for (label i = 0; i < ni; i++)
{
@@ -238,13 +238,13 @@ std::list<Vb::Point> faceCentredCubic::initialPoints() const
{
const point& p(points[i]);
- initialPoints.push_back(Vb::Point(p.x(), p.y(), p.z()));
+ initialPoints.append(Vb::Point(p.x(), p.y(), p.z()));
}
}
}
}
- return initialPoints;
+ return initialPoints.shrink();
}
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/faceCentredCubic/faceCentredCubic.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/faceCentredCubic/faceCentredCubic.H
index c9194b99b65bb1398e0ce68160d94e7e60e02f67..b7660d023bb403b91fd3eb642dbe3d2011a262b6 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/faceCentredCubic/faceCentredCubic.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/faceCentredCubic/faceCentredCubic.H
@@ -89,7 +89,7 @@ public:
// Member Functions
//- Return the initial points for the conformalVoronoiMesh
- virtual std::list<Vb::Point> initialPoints() const;
+ virtual List<Vb::Point> initialPoints() const;
};
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/initialPointsMethod/initialPointsMethod.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/initialPointsMethod/initialPointsMethod.H
index 6feb9ac79d93ec2a03fa23674c31313727528dff..d4334ddf98a5a1c227cbe39b37175b5bba45e6c8 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/initialPointsMethod/initialPointsMethod.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/initialPointsMethod/initialPointsMethod.H
@@ -139,7 +139,7 @@ public:
}
//- Return the initial points for the conformalVoronoiMesh
- virtual std::list<Vb::Point> initialPoints() const = 0;
+ virtual List<Vb::Point> initialPoints() const = 0;
};
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
index 0c9792710952ede38abb3a1f6757a9cab41abb35..f037435af1d1641e0a8f0335958cfb8e4399decd 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.C
@@ -51,7 +51,7 @@ pointFile::pointFile
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-std::list<Vb::Point> pointFile::initialPoints() const
+List<Vb::Point> pointFile::initialPoints() const
{
pointIOField points
(
@@ -69,7 +69,7 @@ std::list<Vb::Point> pointFile::initialPoints() const
if (points.empty())
{
- FatalErrorIn("std::list<Vb::Point> pointFile::initialPoints() const")
+ FatalErrorIn("List<Vb::Point> pointFile::initialPoints() const")
<< "Point file contain no points"
<< exit(FatalError) << endl;
}
@@ -126,8 +126,6 @@ std::list<Vb::Point> pointFile::initialPoints() const
}
}
- std::list<Vb::Point> initialPoints;
-
Field<bool> insidePoints = cvMesh_.geometryToConformTo().wellInside
(
points,
@@ -138,16 +136,20 @@ std::list<Vb::Point> pointFile::initialPoints() const
)
);
+ DynamicList<Vb::Point> initialPoints(insidePoints.size()/10);
+
forAll(insidePoints, i)
{
if (insidePoints[i])
{
const point& p(points[i]);
- initialPoints.push_back(Vb::Point(p.x(), p.y(), p.z()));
+ initialPoints.append(Vb::Point(p.x(), p.y(), p.z()));
}
}
+ initialPoints.shrink();
+
label nPointsRejected = points.size() - initialPoints.size();
if (Pstream::parRun())
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.H
index 74c0775d405bd58f13929327b84821be665f1c40..a60e23e67a2ede9db8cd0bd730fdb0e11b6104ff 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/pointFile/pointFile.H
@@ -85,7 +85,7 @@ public:
// Member Functions
//- Return the initial points for the conformalVoronoiMesh
- virtual std::list<Vb::Point> initialPoints() const;
+ virtual List<Vb::Point> initialPoints() const;
};
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C
index 9413b38501f9b0fc08f7ca7cb30aaca482a65320..c440f4cdc95a252ba38942954344b816c3b94374 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.C
@@ -56,7 +56,7 @@ uniformGrid::uniformGrid
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
-std::list<Vb::Point> uniformGrid::initialPoints() const
+List<Vb::Point> uniformGrid::initialPoints() const
{
boundBox bb;
@@ -91,7 +91,8 @@ std::list<Vb::Point> uniformGrid::initialPoints() const
scalar pert = randomPerturbationCoeff_*cmptMin(delta);
- std::list<Vb::Point> initialPoints;
+ // Initialise points list
+ DynamicList<Vb::Point> initialPoints(ni*nj*nk/10);
for (label i = 0; i < ni; i++)
{
@@ -153,13 +154,13 @@ std::list<Vb::Point> uniformGrid::initialPoints() const
{
const point& p(points[i]);
- initialPoints.push_back(Vb::Point(p.x(), p.y(), p.z()));
+ initialPoints.append(Vb::Point(p.x(), p.y(), p.z()));
}
}
}
}
- return initialPoints;
+ return initialPoints.shrink();
}
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.H
index 655efae4702f422525b402501bb2114d01c729a4..158a49b59e9372be1b8e22d30925a4b0b76205ba 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/initialPointsMethod/uniformGrid/uniformGrid.H
@@ -89,7 +89,7 @@ public:
// Member Functions
//- Return the initial points for the conformalVoronoiMesh
- virtual std::list<Vb::Point> initialPoints() const;
+ virtual List<Vb::Point> initialPoints() const;
};
diff --git a/applications/utilities/mesh/generation/cvMesh/cvMesh.C b/applications/utilities/mesh/generation/cvMesh/cvMesh.C
index fa95c852c11c5ccaa307ab80e7f63c4e2ff0bb19..daf03aa3766ab1bb75e0df22ea441f84579a3dc2 100644
--- a/applications/utilities/mesh/generation/cvMesh/cvMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/cvMesh.C
@@ -59,6 +59,8 @@ int main(int argc, char *argv[])
)
);
+ conformalVoronoiMesh::debug = true;
+
conformalVoronoiMesh mesh(runTime, cvMeshDict);
while (runTime.loop())
@@ -72,7 +74,7 @@ int main(int argc, char *argv[])
<< nl << endl;
}
- mesh.writeMesh(runTime.constant());
+ mesh.writeMesh(runTime.constant(), true);
Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
diff --git a/applications/utilities/mesh/generation/cvMesh/vectorTools/vectorTools.H b/applications/utilities/mesh/generation/cvMesh/vectorTools/vectorTools.H
new file mode 100644
index 0000000000000000000000000000000000000000..9fb39a159bafbfb8de882cb4b469fbf6e7d14af8
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/vectorTools/vectorTools.H
@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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/>.
+
+Class
+ Foam::vectorTools
+
+Description
+ Functions for analysing the relationships between vectors
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef vectorTools_H
+#define vectorTools_H
+
+#include "vector.H"
+#include "unitConversion.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Namespace vectorTools Declaration
+\*---------------------------------------------------------------------------*/
+
+//- Collection of functions for testing relationships between two vectors.
+namespace vectorTools
+{
+ //- Test if a and b are parallel: a.b = 1
+ // Uses the cross product, so the tolerance is proportional to
+ // the sine of the angle between a and b in radians
+ template <typename T>
+ bool areParallel
+ (
+ const Vector<T>& a,
+ const Vector<T>& b,
+ const T& tolerance = SMALL
+ )
+ {
+ return (mag(a ^ b) < tolerance) ? true : false;
+// return ( mag( mag(a & b)/(mag(a)*mag(b)) - 1.0 ) < tolerance )
+// ? true
+// : false;
+ }
+
+ //- Test if a and b are orthogonal: a.b = 0
+ // Uses the dot product, so the tolerance is proportional to
+ // the cosine of the angle between a and b in radians
+ template <typename T>
+ bool areOrthogonal
+ (
+ const Vector<T>& a,
+ const Vector<T>& b,
+ const T& tolerance = SMALL
+ )
+ {
+ return (mag(a & b) < tolerance) ? true : false;
+ }
+
+ //- Test if angle between a and b is acute: a.b > 0
+ template <typename T>
+ bool areAcute
+ (
+ const Vector<T>& a,
+ const Vector<T>& b
+ )
+ {
+ return ((a & b) > 0) ? true : false;
+ }
+
+ //- Test if angle between a and b is obtuse: a.b < 0
+ template <typename T>
+ bool areObtuse
+ (
+ const Vector<T>& a,
+ const Vector<T>& b
+ )
+ {
+ return ((a & b) < 0) ? true : false;
+ }
+
+ //- Calculate angle between a and b in radians
+ template <typename T>
+ T cosPhi
+ (
+ const Vector<T>& a,
+ const Vector<T>& b,
+ const T& tolerance = SMALL
+ )
+ {
+ scalar cosPhi = (a & b)/(mag(a)*mag(b) + tolerance);
+
+ // Enforce bounding between -1 and 1
+ return min(max(cosPhi, -1), 1);
+ }
+
+ //- Calculate angle between a and b in radians
+ template <typename T>
+ T radAngleBetween
+ (
+ const Vector<T>& a,
+ const Vector<T>& b,
+ const T& tolerance = SMALL
+ )
+ {
+ scalar cosPhi = (a & b)/(mag(a)*mag(b) + tolerance);
+
+ // Enforce bounding between -1 and 1
+ return acos( min(max(cosPhi, -1), 1) );
+ }
+
+ //- Calculate angle between a and b in degrees
+ template <typename T>
+ T degAngleBetween
+ (
+ const Vector<T>& a,
+ const Vector<T>& b,
+ const T& tolerance = SMALL
+ )
+ {
+ return radToDeg(radAngleBetween(a, b, tolerance));
+ }
+
+} // End namespace vectorTools
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C b/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C
index 2179e67b0c0a3085b405abd2e795c193dcd98b25..2951cab590a4c70903d73f1a3a5d719774751dcb 100644
--- a/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C
+++ b/applications/utilities/mesh/manipulation/checkMesh/printMeshStats.C
@@ -8,6 +8,7 @@
#include "pyrMatcher.H"
#include "tetWedgeMatcher.H"
#include "tetMatcher.H"
+#include "IOmanip.H"
void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
@@ -62,22 +63,19 @@ void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
<< nInternalEdges-nInternal1Edges << nl;
}
- Info<< " faces: "
- << returnReduce(mesh.faces().size(), sumOp<label>()) << nl
- << " internal faces: "
- << returnReduce(mesh.faceNeighbour().size(), sumOp<label>()) << nl
- << " cells: "
- << returnReduce(mesh.cells().size(), sumOp<label>()) << nl
- << " boundary patches: "
- << mesh.boundaryMesh().size() << nl
- << " point zones: "
- << mesh.pointZones().size() << nl
- << " face zones: "
- << mesh.faceZones().size() << nl
- << " cell zones: "
- << mesh.cellZones().size() << nl
- << endl;
+ label nFaces = returnReduce(mesh.faces().size(), sumOp<label>());
+ label nIntFaces = returnReduce(mesh.faceNeighbour().size(), sumOp<label>());
+ label nCells = returnReduce(mesh.cells().size(), sumOp<label>());
+ Info<< " faces: " << nFaces << nl
+ << " internal faces: " << nIntFaces << nl
+ << " cells: " << nCells << nl
+ << " faces per cell: " << scalar(nFaces + nIntFaces)/nCells << nl
+ << " boundary patches: " << mesh.boundaryMesh().size() << nl
+ << " point zones: " << mesh.pointZones().size() << nl
+ << " face zones: " << mesh.faceZones().size() << nl
+ << " cell zones: " << mesh.cellZones().size() << nl
+ << endl;
// Construct shape recognizers
hexMatcher hex;
@@ -96,6 +94,8 @@ void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
label nTetWedge = 0;
label nUnknown = 0;
+ Map<label> polyhedralFaces;
+
for (label cellI = 0; cellI < mesh.nCells(); cellI++)
{
if (hex.isA(mesh, cellI))
@@ -125,6 +125,7 @@ void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
else
{
nUnknown++;
+ polyhedralFaces(mesh.cells()[cellI].size())++;
}
}
@@ -144,5 +145,21 @@ void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
<< " tet wedges: " << nTetWedge << nl
<< " tetrahedra: " << nTet << nl
<< " polyhedra: " << nUnknown
- << nl << endl;
+ << endl;
+
+ if (nUnknown > 0)
+ {
+ Pstream::mapCombineGather(polyhedralFaces, plusEqOp<label>());
+
+ Info<< " Breakdown of polyhedra by number of faces:" << endl;
+ Info<< " faces" << " number of cells" << endl;
+
+ forAllConstIter(Map<label>, polyhedralFaces, iter)
+ {
+ Info<< setf(std::ios::right) << setw(13)
+ << iter.key() << " " << iter() << nl;
+ }
+ }
+
+ Info<< endl;
}