diff --git a/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/files b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/files
new file mode 100644
index 0000000000000000000000000000000000000000..73871a7532e1c273d0030d8c8bc4180e1ecef2a7
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/files
@@ -0,0 +1,3 @@
+checkCvMesh.C
+
+EXE = $(FOAM_APPBIN)/checkCvMesh
diff --git a/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/options b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/options
new file mode 100644
index 0000000000000000000000000000000000000000..ba68fd3819b53037f9711dfbb9691846a215f4d3
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/Make/options
@@ -0,0 +1,14 @@
+EXE_INC = \
+ -I$(LIB_SRC)/finiteVolume/lnInclude \
+ -I$(LIB_SRC)/triSurface/lnInclude \
+ -I$(LIB_SRC)/mesh/autoMesh/lnInclude \
+ -I$(LIB_SRC)/dynamicMesh/lnInclude \
+ -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+ -lfiniteVolume \
+ -ldynamicMesh \
+ -ltriSurface \
+ -lautoMesh \
+ -lmeshTools
+
diff --git a/applications/utilities/mesh/generation/cvMesh/checkCvMesh/checkCvMesh.C b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/checkCvMesh.C
new file mode 100644
index 0000000000000000000000000000000000000000..df2a3b90f81834625c22ab8fc3b490c4c0b73d99
--- /dev/null
+++ b/applications/utilities/mesh/generation/cvMesh/checkCvMesh/checkCvMesh.C
@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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/>.
+
+Application
+ checkCvMesh
+
+Description
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "fvMesh.H"
+#include "autoSnapDriver.H"
+#include "faceSet.H"
+#include "motionSmoother.H"
+#include "timeSelector.H"
+
+
+using namespace Foam;
+
+
+int main(int argc, char *argv[])
+{
+# include "addOverwriteOption.H"
+
+# include "setRootCase.H"
+# include "createTime.H"
+
+ instantList timeDirs = timeSelector::select0(runTime, args);
+
+# include "createMesh.H"
+
+ runTime.functionObjects().off();
+
+ forAll(timeDirs, timeI)
+ {
+ runTime.setTime(timeDirs[timeI], timeI);
+
+ Info<< "Create mesh for time = " << runTime.timeName()
+ << nl << endl;
+
+ mesh.readUpdate();
+
+ Info<< "Read mesh in = "
+ << runTime.cpuTimeIncrement() << " s" << endl;
+
+ // Check patches and faceZones are synchronised
+ mesh.boundaryMesh().checkParallelSync(true);
+ meshRefinement::checkCoupledFaceZones(mesh);
+
+ // Read meshing dictionary
+ IOdictionary cvMeshDict
+ (
+ IOobject
+ (
+ "cvMeshDict",
+ runTime.system(),
+ mesh,
+ IOobject::MUST_READ_IF_MODIFIED,
+ IOobject::NO_WRITE
+ )
+ );
+
+ // mesh motion and mesh quality parameters
+ const dictionary& meshQualityDict
+ = cvMeshDict.subDict("meshQualityControls");
+
+
+ Info<< "Checking initial mesh ..." << endl;
+ faceSet wrongFaces(mesh, "wrongFaces", mesh.nFaces()/100);
+ motionSmoother::checkMesh(false, mesh, meshQualityDict, wrongFaces);
+
+ const label nInitErrors = returnReduce
+ (
+ wrongFaces.size(),
+ sumOp<label>()
+ );
+
+ Info<< "Detected " << nInitErrors << " illegal faces"
+ << " (concave, zero area or negative cell pyramid volume)"
+ << endl;
+
+ if (nInitErrors > 0)
+ {
+ Info<< "Writing " << nInitErrors
+ << " faces in error to set "
+ << wrongFaces.name() << endl;
+
+ wrongFaces.instance() = mesh.pointsInstance();
+ wrongFaces.write();
+ }
+
+ Info<< nl << "End of time " << runTime.timeName() << nl << endl;
+ }
+
+ Info<< "End\n" << endl;
+
+ return 0;
+
+}
+
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
index edafdb50f912e02a89acc6133fc0d7cf7491cd5c..a2d44f0ee1ab506045eabdeac7ccb7300ad3d275 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.C
@@ -724,6 +724,8 @@ void Foam::backgroundMeshDecomposition::buildPatchAndTree()
globalBackgroundBounds_ = treeBoundBox(bbMin, bbMax);
+ octreeNearestDistances_ = bFTreePtr_().calcNearestDistance();
+
if (cvMesh_.cvMeshControls().objOutput())
{
OFstream fStr
@@ -795,6 +797,7 @@ Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
),
boundaryFacesPtr_(),
bFTreePtr_(),
+ octreeNearestDistances_(),
allBackgroundMeshBounds_(Pstream::nProcs()),
globalBackgroundBounds_(),
decomposeDict_
@@ -1150,6 +1153,8 @@ bool Foam::backgroundMeshDecomposition::positionOnThisProcessor
const point& pt
) const
{
+// return bFTreePtr_().findAnyOverlap(pt, 0.0);
+
return
bFTreePtr_().getVolumeType(pt)
== indexedOctree<treeDataBPatch>::INSIDE;
@@ -1176,6 +1181,7 @@ bool Foam::backgroundMeshDecomposition::overlapsThisProcessor
const treeBoundBox& box
) const
{
+// return !procBounds().contains(box);
return !bFTreePtr_().findBox(box).empty();
}
@@ -1183,9 +1189,11 @@ bool Foam::backgroundMeshDecomposition::overlapsThisProcessor
bool Foam::backgroundMeshDecomposition::overlapsThisProcessor
(
const point& centre,
- scalar radiusSqr
+ const scalar radiusSqr
) const
{
+ //return bFTreePtr_().findAnyOverlap(centre, radiusSqr);
+
return bFTreePtr_().findNearest(centre, radiusSqr).hit();
}
@@ -1645,6 +1653,7 @@ Foam::labelListList Foam::backgroundMeshDecomposition::overlapsProcessors
// If the sphere finds a nearest element of the patch, then it overlaps
sphereOverlapsCandidate[sI] = bFTreePtr_().findNearest(c, rSqr).hit();
+ //sphereOverlapsCandidate[sI] = bFTreePtr_().findAnyOverlap(c, rSqr);
}
map().reverseDistribute
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
index 87ab88af356873a7fc4f4512645ac46873806110..d8e6388cfb00e28d73c38928d9be482179689697 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecomposition.H
@@ -111,6 +111,8 @@ class backgroundMeshDecomposition
//- Search tree for the boundaryFaces_ patch
autoPtr<indexedOctree<treeDataBPatch> > bFTreePtr_;
+ List<scalar> octreeNearestDistances_;
+
//- The bounds of all background meshes on all processors
treeBoundBoxList allBackgroundMeshBounds_;
@@ -225,16 +227,16 @@ public:
//- Are the given positions inside the domain of this decomposition
boolList positionOnThisProcessor(const List<point>& pts) const;
- //- Does the given box overlap the faces of the bounday of this
+ //- Does the given box overlap the faces of the boundary of this
// processor
bool overlapsThisProcessor(const treeBoundBox& box) const;
- //- Does the given sphere overlap the faces of the bounday of this
+ //- Does the given sphere overlap the faces of the boundary of this
// processor
bool overlapsThisProcessor
(
const point& centre,
- scalar radiusSqr
+ const scalar radiusSqr
) const;
//- Find nearest intersection of line between start and end, (exposing
@@ -289,6 +291,12 @@ public:
//- Return access to the underlying mesh
inline const fvMesh& mesh() const;
+ //- Return access to the underlying tree
+ inline const indexedOctree<treeDataBPatch>& tree() const;
+
+ //- Return access to the nearest distance of the octree nodes
+ inline const List<scalar>& octreeNearestDistances() const;
+
//- Return the boundBox of this processor
inline const treeBoundBox& procBounds() const;
};
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecompositionI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecompositionI.H
index df659cfd9cc0fed0f996b2bb97ff0b68622a1dcf..2897d3fc5830e133ba305b23436df6f829ccff99 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecompositionI.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/backgroundMeshDecomposition/backgroundMeshDecompositionI.H
@@ -30,6 +30,18 @@ const Foam::fvMesh& Foam::backgroundMeshDecomposition::mesh() const
return mesh_;
}
+const Foam::indexedOctree<Foam::treeDataBPatch>&
+Foam::backgroundMeshDecomposition::tree() const
+{
+ return bFTreePtr_();
+}
+
+const Foam::List<Foam::scalar>&
+Foam::backgroundMeshDecomposition::octreeNearestDistances() const
+{
+ return octreeNearestDistances_;
+}
+
const Foam::treeBoundBox& Foam::backgroundMeshDecomposition::procBounds() const
{
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
index e8cedc49c2cbd5e05256b513c4be1250f8140e4f..958fcf6763ec79332bf98faee94882df8eb0c640 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.C
@@ -240,9 +240,9 @@ void Foam::conformalVoronoiMesh::insertPoints
{
label preDistributionSize(points.size());
- DynamicList<Foam::point> transferPoints(points.size()/2);
+ DynamicList<Foam::point> transferPoints;
- DynamicList<Point> pointsOnProcessor(points.size()/2);
+ DynamicList<Point> pointsOnProcessor;
for
(
@@ -277,12 +277,16 @@ void Foam::conformalVoronoiMesh::insertPoints
decomposition_().distributePoints(transferPoints)
);
+ const label oldSize = points.size();
+
+ points.setSize(oldSize + transferPoints.size());
+
forAll(transferPoints, tPI)
{
- points.append(toPoint(transferPoints[tPI]));
+ points[tPI + oldSize] = toPoint(transferPoints[tPI]);
}
- label sizeChange = preDistributionSize - label(points.size());
+ label sizeChange = preDistributionSize - points.size();
// if (mag(sizeChange) > 0)
// {
@@ -417,7 +421,6 @@ void Foam::conformalVoronoiMesh::insertPoints
// );
// }
-
rangeInsertWithInfo
(
pts.begin(),
@@ -1246,6 +1249,8 @@ Foam::conformalVoronoiMesh::conformalVoronoiMesh
// better balance the surface conformation load.
distributeBackground();
+// conformToSurface();
+
buildSurfaceConformation(rmCoarse);
// The introduction of the surface conformation may have distorted the
@@ -1313,7 +1318,7 @@ void Foam::conformalVoronoiMesh::move()
{
cit->cellIndex() = dualVertI;
- dualVertices[dualVertI] = topoint(dual(cit));
+ dualVertices[dualVertI] = cit->dual();
dualVertI++;
}
@@ -1511,7 +1516,6 @@ void Foam::conformalVoronoiMesh::move()
(
toPoint(0.5*(dVA + dVB))
);
-
}
}
else if
@@ -1671,9 +1675,57 @@ void Foam::conformalVoronoiMesh::move()
insertPoints(pointsToInsert);
+ // Remove internal points that have been inserted outside the surface.
+ label internalPtIsOutside = 0;
+
+ for
+ (
+ Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+ vit != finite_vertices_end();
+ ++vit
+ )
+ {
+ if (vit->internalPoint())
+ {
+ bool inside
+ = geometryToConformTo_.inside(topoint(vit->point()));
+
+ if (!inside)
+ {
+ remove(vit);
+ internalPtIsOutside++;
+ }
+ }
+ }
+
+ Info<< " " << internalPtIsOutside
+ << " internal points were inserted outside the domain. "
+ << "They have been removed." << endl;
+
+ // Fix points that have not been significantly displaced
+// for
+// (
+// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+// vit != finite_vertices_end();
+// ++vit
+// )
+// {
+// if (vit->internalPoint())
+// {
+// if
+// (
+// mag(displacementAccumulator[vit->index()])
+// < 0.1*targetCellSize(topoint(vit->point()))
+// )
+// {
+// vit->setVertexFixed();
+// }
+// }
+// }
+
if (cvMeshControls().objOutput() && runTime_.outputTime())
{
- writePoints("points_" + runTime_.timeName() + ".obj", false);
+ writePoints("points_" + runTime_.timeName() + ".obj", true);
}
timeCheck("Internal points inserted");
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
index b590eb31cd67383146216cfbdc4f211a38397137..a514c660a1d96350897478ca826dba7efe5af1d8 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H
@@ -48,7 +48,6 @@ SourceFiles
#define CGAL_INEXACT
#include "CGALTriangulation3Ddefs.H"
-#include <CGAL/Spatial_sort_traits_adapter_3.h>
#include "uint.H"
#include "ulong.H"
#include "searchableSurfaces.H"
@@ -595,6 +594,12 @@ private:
Foam::point& b
) const;
+ label removeProcessorBoundarySeeds(bool reinsertBoundPts);
+
+ void seedProcessorBoundarySurfaces(bool seedProcessors);
+
+ label numberOfUnusedReferredPoints() const;
+
//- Build the parallelInterfaces of the mesh
void buildParallelInterface
(
@@ -1240,7 +1245,7 @@ public:
Traits_for_spatial_sort<Triangulation>()
);
- typename Triangulation::Cell_handle hint;
+ typename Triangulation::Vertex_handle hint;
for
(
@@ -1250,16 +1255,9 @@ public:
++p
)
{
- typename Triangulation::Locate_type lt;
- typename Triangulation::Cell_handle c;
- label li, lj;
-
- c = T.locate(*(p->first), lt, li, lj, hint);
-
const size_t checkInsertion = T.number_of_vertices();
- typename Triangulation::Vertex_handle v
- = T.insert(*(p->first), lt, c, li, lj);
+ hint = T.insert(*(p->first), hint);
if (checkInsertion != T.number_of_vertices() - 1)
{
@@ -1278,12 +1276,11 @@ public:
// type directly (note that this routine never gets
// called for referredPoints so type will never be
// -procI
- type += T.number_of_vertices() - 1;
+ type += checkInsertion;
}
- v->index() = indices[oldIndex]
- + T.number_of_vertices() - 1;
- v->type() = type;
+ hint->index() = indices[oldIndex] + checkInsertion;
+ hint->type() = type;
}
}
}
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
index 04ef2af045d22226200bce8d1943ec4ed9438696..d103d3c0fcb1ff968452105319a3c304e5a65578 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C
@@ -1890,7 +1890,7 @@ void Foam::conformalVoronoiMesh::indexDualVertices
{
cit->cellIndex() = dualVertI;
- pts[dualVertI] = topoint(dual(cit));
+ pts[dualVertI] = cit->dual();
if
(
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
index 847974a2f4ed3c01bdc652fe52baafafcb8bc96c..7f4c7abb6ef2cd34a4513fda2e94a41cf5af0d9b 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C
@@ -41,6 +41,11 @@ void Foam::conformalVoronoiMesh::conformToSurface()
{
reconformationMode reconfMode = reconformationControl();
+ if (Pstream::parRun())
+ {
+ seedProcessorBoundarySurfaces(true);
+ }
+
if (reconfMode == rmNone)
{
// Reinsert stored surface conformation
@@ -69,6 +74,16 @@ void Foam::conformalVoronoiMesh::conformToSurface()
storeSurfaceConformation();
}
+ if (Pstream::parRun())
+ {
+ label nFarPoints = removeProcessorBoundarySeeds(true);
+
+ reduce(nFarPoints, sumOp<label>());
+
+ Info<< " Removed " << nFarPoints
+ << " far points from the mesh." << endl;
+ }
+
// reportSurfaceConformationQuality();
}
@@ -140,13 +155,8 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
buildEdgeLocationTree(existingEdgeLocations);
buildSurfacePtLocationTree(existingSurfacePtLocations);
-
- // Initialise the edgeLocationTree
- //buildEdgeLocationTree(edgeLocationTree, existingEdgeLocations);
-
label initialTotalHits = 0;
-
// Surface protrusion conformation is done in two steps.
// 1. the dual edges (of all internal vertices) can stretch to
// 'infinity' so any intersection would be badly behaved. So
@@ -382,7 +392,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
(
Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
- vit++
+ ++vit
)
{
// The initial surface conformation has already identified the
@@ -511,14 +521,18 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation
timeCheck("Conformation iteration " + name(iterationNo));
- // Update the parallel interface
- buildParallelInterface
- (
- referralVertices,
- receivedVertices,
- false,
- name(iterationNo)
- );
+ // Only need to update the interface if there are surface/edge hits
+ if (totalHits > 0)
+ {
+ // Update the parallel interface
+ buildParallelInterface
+ (
+ referralVertices,
+ receivedVertices,
+ false,
+ name(iterationNo)
+ );
+ }
iterationNo++;
@@ -563,8 +577,8 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAnyIntersection
continue;
}
- Foam::point dE0 = topoint(dual(fit->first));
- Foam::point dE1 = topoint(dual(fit->first->neighbor(fit->second)));
+ Foam::point dE0 = fit->first->dual();
+ Foam::point dE1 = fit->first->neighbor(fit->second)->dual();
if (Pstream::parRun())
{
@@ -631,8 +645,8 @@ bool Foam::conformalVoronoiMesh::dualCellSurfaceAllIntersections
// 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)));
+ Foam::point dE0 = fit->first->dual();
+ Foam::point dE1 = fit->first->neighbor(fit->second)->dual();
pointIndexHit infoIntersection;
label hitSurfaceIntersection = -1;
@@ -801,6 +815,144 @@ void Foam::conformalVoronoiMesh::buildParallelInterface
}
+Foam::label Foam::conformalVoronoiMesh::removeProcessorBoundarySeeds
+(
+ bool reinsertBoundPts
+)
+{
+ label nFarPoints = 0;
+
+
+
+ std::list<Vertex_handle> toRemove;
+
+ for
+ (
+ Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+ vit != finite_vertices_end();
+ ++vit
+ )
+ {
+ if (vit->farPoint())
+ {
+ //remove(vit);
+ toRemove.push_back(vit);
+ nFarPoints++;
+ }
+ }
+
+ // This function removes the points in the iterator range and then
+ // retriangulates.
+ timeCheck("Start Removing Seeded Points " + name(toRemove.size()));
+
+ remove_cluster(toRemove.begin(), toRemove.end());
+
+
+ // Need to do this to make sure the triangulation is well-behaved
+ if (reinsertBoundPts)
+ {
+ reinsertBoundingPoints();
+ }
+
+
+
+// DynamicList<Foam::point> toAdd;
+// DynamicList<label> indices;
+// DynamicList<label> types;
+//
+// for
+// (
+// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+// vit != finite_vertices_end();
+// ++vit
+// )
+// {
+// if (!vit->farPoint())
+// {
+// toAdd.append(topoint(vit->point()));
+// indices.append(vit->index());
+// types.append(vit->type());
+//
+// nFarPoints++;
+// }
+// }
+//
+// this->clear();
+//
+// // Need to do this to make sure the triangulation is well-behaved
+// if (reinsertBoundPts)
+// {
+// reinsertBoundingPoints();
+// }
+//
+// forAll(toAdd, pI)
+// {
+// insertPoint(toAdd[pI], indices[pI], types[pI]);
+// }
+
+
+ timeCheck("End Removing Seeded Points");
+
+ return nFarPoints;
+}
+
+
+void Foam::conformalVoronoiMesh::seedProcessorBoundarySurfaces
+(
+ bool seedProcessors
+)
+{
+ //removeProcessorBoundarySeeds(false);
+
+ // Loop over processor patch faces and insert a point in the centre of the
+ // face
+ const fvMesh& mesh = decomposition_().mesh();
+ const polyBoundaryMesh& bMesh = mesh.boundaryMesh();
+
+ DynamicList<Foam::point> pts;
+ DynamicList<label> indices;
+ DynamicList<label> types;
+
+ label nFarPoints = 0;
+
+ const scalar normalDistance = 5.0;
+ const scalar pert = 0.1*(rndGen_.scalar01() - 0.5);
+
+ forAll(bMesh, patchI)
+ {
+ const polyPatch& patch = bMesh[patchI];
+
+ if (!seedProcessors && isA<processorPolyPatch>(patch))
+ {
+ continue;
+ }
+
+ forAll(patch, faceI)
+ {
+ if (faceI % 1 == 0)
+ {
+ const face& f = patch[faceI];
+
+ pts.append
+ (
+ f.centre(mesh.points())
+ + pert*normalDistance*f.normal(mesh.points())
+ );
+ indices.append(nFarPoints++);
+ types.append(Vb::vtFar);
+ }
+ }
+ }
+
+ insertPoints(pts, indices, types, false);
+
+ reduce(nFarPoints, sumOp<label>());
+
+ Info<< " Inserted " << nFarPoints
+ << " far points into the mesh." << endl;
+}
+
+
void Foam::conformalVoronoiMesh::buildParallelInterface
(
List<labelHashSet>& referralVertices,
@@ -835,22 +987,88 @@ void Foam::conformalVoronoiMesh::buildParallelInterface
timeCheck("After buildParallelInterfaceAll");
}
- if (initialEdgeReferral)
- {
- // Used as an initial pass to localise the vertex referring - find
- // vertices whose dual edges pierce nearby processor volumes and refer
- // them to establish a sensible boundary interface region before
- // running a circumsphere assessment.
- buildParallelInterfaceIntersection
- (
- referralVertices,
- receivedVertices,
- outputName
- );
+ // Reject points that are not near the boundary from the subsequent
+ // searches
- timeCheck("After buildParallelInterfaceIntersection");
- }
+// label nearProcCount = 0;
+// label notNearProcCount = 0;
+// boundBox quickRejectionBox(decomposition_().procBounds());
+//
+//// Pout<< "Processor boundBox: " << quickRejectionBox << endl;
+//
+// quickRejectionBox.inflate(-0.1);
+//
+// OFstream str("rejectionBox_" + name(Pstream::myProcNo()) + ".obj");
+//
+// meshTools::writeOBJ
+// (
+// str,
+// boundBox::faces(),
+// quickRejectionBox.points()
+// );
+//
+// Pout<< "Inflated rejection boundBox: " << quickRejectionBox << endl;
+
+// for
+// (
+// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+// vit != finite_vertices_end();
+// vit++
+// )
+// {
+// if (vit->real() && !vit->nearProcBoundary())
+// {
+// const Foam::point& pt = topoint(vit->point());
+// const scalar range = targetCellSize(pt);
+//
+//// if
+//// (
+//// decomposition_().overlapsThisProcessor
+//// (
+//// pt,
+//// range
+//// )
+//// )
+// if (!quickRejectionBox.contains(pt))
+// {
+// vit->setNearProcBoundary();
+// nearProcCount++;
+// }
+// else
+// {
+// //vit->setNearProcBoundary();
+// notNearProcCount++;
+// }
+// }
+// }
+//
+// reduce(nearProcCount, sumOp<label>());
+// reduce(notNearProcCount, sumOp<label>());
+//
+// timeCheck
+// (
+// "End of potential intersection search. "
+// + name(nearProcCount) + " are near a processor boundary. "
+// + name(notNearProcCount) + " are not."
+// );
+
+// if (initialEdgeReferral)
+// {
+// // Used as an initial pass to localise the vertex referring - find
+// // vertices whose dual edges pierce nearby processor volumes and refer
+// // them to establish a sensible boundary interface region before
+// // running a circumsphere assessment.
+//
+// buildParallelInterfaceIntersection
+// (
+// referralVertices,
+// receivedVertices,
+// outputName
+// );
+//
+// timeCheck("After buildParallelInterfaceIntersection");
+// }
buildParallelInterfaceInfluence
(
@@ -860,6 +1078,56 @@ void Foam::conformalVoronoiMesh::buildParallelInterface
);
timeCheck("After buildParallelInterface");
+
+ // Check all referred vertices are actually used on the processor.
+ label nUnusedReferred = numberOfUnusedReferredPoints();
+
+ reduce(nUnusedReferred, sumOp<label>());
+
+ Info<< " Number of referred points that are not used : "
+ << nUnusedReferred << " (approximate)" << endl;
+}
+
+
+Foam::label Foam::conformalVoronoiMesh::numberOfUnusedReferredPoints() const
+{
+ label nUnusedPoints = 0;
+
+ for
+ (
+ Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+ vit != finite_vertices_end();
+ ++vit
+ )
+ {
+ if (vit->referred())
+ {
+ std::list<Vertex_handle> adjVertices;
+ finite_adjacent_vertices(vit, std::back_inserter(adjVertices));
+
+ bool isUsed = false;
+
+ for
+ (
+ std::list<Vertex_handle>::iterator adjVit = adjVertices.begin();
+ adjVit != adjVertices.end();
+ ++adjVit
+ )
+ {
+ if ((*adjVit)->real())
+ {
+ isUsed = true;
+ }
+ }
+
+ if (!isUsed)
+ {
+ nUnusedPoints++;
+ }
+ }
+ }
+
+ return nUnusedPoints;
}
@@ -880,7 +1148,7 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceAll
(
Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
vit != finite_vertices_end();
- vit++
+ ++vit
)
{
if (!vit->farPoint())
@@ -966,7 +1234,10 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceIntersection
// If either Delaunay cell at the end of the Dual edge is infinite,
// skip.
- if (!is_infinite(c1) && !is_infinite(c2))
+ if
+ (
+ !is_infinite(c1) && !is_infinite(c2)
+ )
{
// The Delaunauy cells at either end of the dual edge need to be
// real, i.e. all vertices form part of the internal or boundary
@@ -977,12 +1248,14 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceIntersection
&& c2->internalOrBoundaryDualVertex()
)
{
- Foam::point a = topoint(dual(c1));
- Foam::point b = topoint(dual(c2));
+ const Foam::point& a = c1->dual();
+ const Foam::point& b = c2->dual();
// Only if the dual edge cuts the boundary of this processor is
// it going to be counted.
- if (decomposition_().findLineAny(a, b).hit())
+ pointIndexHit info = decomposition_().findLineAny(a, b);
+
+ if (info.hit())
{
dE0.append(a);
dE1.append(b);
@@ -995,11 +1268,22 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceIntersection
fIOuter++;
}
+ reduce(fIOuter, sumOp<label>());
+
+ timeCheck
+ (
+ "End of actual intersection search over "
+ + name(fIOuter)
+ + " faces."
+ );
+
// Preform intersections in both directions, as there is no sense
// associated with the Dual edge
List<List<pointIndexHit> > intersectionForward(intersectsProc(dE0, dE1));
List<List<pointIndexHit> > intersectionReverse(intersectsProc(dE1, dE0));
+ timeCheck("End of find processor intersection");
+
// Reset counter
fIOuter = 0;
@@ -1175,13 +1459,58 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
DynamicList<Foam::point> circumcentre;
DynamicList<scalar> circumradiusSqr;
- PackedBoolList testCellInfluence(number_of_cells(), false);
// Index outer (all) Delaunauy cells for whether they are potential
// overlaps, index (inner) the list of tests an results.
label cIInner = 0;
label cIOuter = 0;
+
+ label cellIndexCount = 0;
+ for
+ (
+ Delaunay::Finite_cells_iterator cit = finite_cells_begin();
+ cit != finite_cells_end();
+ ++cit
+ )
+ {
+ cit->cellIndex() = cellIndexCount++;
+ }
+
+ timeCheck("End of cell Indexing");
+
+ labelList testCellInfluence(number_of_cells(), 0);
+
+ label nQuickRejections = 0;
+
+ for
+ (
+ Delaunay::Finite_cells_iterator cit = finite_cells_begin();
+ cit != finite_cells_end();
+ ++cit
+ )
+ {
+ const Foam::point& cc = cit->dual();
+
+ const scalar crSqr = magSqr(cc - topoint(cit->vertex(0)->point()));
+
+ if
+ (
+ decomposition_().tree().quickCircumsphereRejection
+ (
+ cc,
+ crSqr,
+ decomposition_().octreeNearestDistances()
+ )
+ )
+ {
+ nQuickRejections++;
+ testCellInfluence[cit->cellIndex()] = -1;
+ }
+ }
+
+ timeCheck("End of octreeNearestDistances calculation");
+
for
(
Delaunay::Finite_cells_iterator cit = finite_cells_begin();
@@ -1196,14 +1525,15 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
// The Delaunay cells to assess have to be real, i.e. all vertices form
// part of the internal or any part of the boundary definition
- if (cit->real())
+ if
+ (
+ (testCellInfluence[cit->cellIndex()] == 0)
+ && (cit->real() || cit->hasFarPoint())
+ )
{
- Foam::point cc(topoint(dual(cit)));
+ const Foam::point& cc = cit->dual();
- scalar crSqr
- (
- magSqr(cc - topoint(cit->vertex(0)->point()))
- );
+ const scalar crSqr = magSqr(cc - topoint(cit->vertex(0)->point()));
// Only if the circumsphere overlaps the boundary of this processor
// is there a chance of it overlapping others
@@ -1212,7 +1542,7 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
circumcentre.append(cc);
circumradiusSqr.append(crSqr);
- testCellInfluence[cIOuter] = true;
+ testCellInfluence[cit->cellIndex()] = 1;
}
}
@@ -1221,6 +1551,9 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
timeCheck("End of testing cell influence");
+ Pout<< "Number of quick rejections = " << nQuickRejections << endl;
+ Pout<< "Number of influences = " << circumcentre.size() << endl;
+
// Increasing the circumspheres to increase the overlaps and compensate for
// floating point errors missing some referrals
labelListList circumsphereOverlaps
@@ -1242,7 +1575,7 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
)
{
// Pre-tested circumsphere potential influence
- if (testCellInfluence[cIOuter])
+ if (testCellInfluence[cit->cellIndex()] == 1)
{
const labelList& citOverlaps = circumsphereOverlaps[cIInner];
@@ -1289,6 +1622,70 @@ void Foam::conformalVoronoiMesh::buildParallelInterfaceInfluence
cIOuter++;
}
+
+// label nFarPoints = removeProcessorBoundarySeeds(true);
+//
+// reduce(nFarPoints, sumOp<label>());
+//
+// Info<< " Removed " << nFarPoints
+// << " far points from the mesh." << endl;
+
+
+// seedProcessorBoundarySurfaces(false);
+
+// cIInner = 0;
+// cIOuter = 0;
+
+
+ // Relying on the order of iteration of cells being the same as before
+// for
+// (
+// Delaunay::Finite_cells_iterator cit = finite_cells_begin();
+// cit != finite_cells_end();
+// ++cit
+// )
+// {
+// // Pre-tested circumsphere potential influence
+// if (testCellInfluence[cIOuter])
+// {
+// const labelList& citOverlaps = circumsphereOverlaps[cIInner];
+//
+// forAll(citOverlaps, cOI)
+// {
+// label procI = citOverlaps[cOI];
+//
+// recursiveCircumsphereSearch
+// (
+// cit,
+// procI,
+// referralVertices,
+// checkedCells,
+// parallelInfluencePoints,
+// parallelInfluenceIndices,
+// targetProcessor
+// );
+// }
+//
+// cIInner++;
+// }
+//
+// cIOuter++;
+// }
+
+// for
+// (
+// Delaunay::Finite_vertices_iterator vit = finite_vertices_begin();
+// vit != finite_vertices_end();
+// ++vit
+// )
+// {
+// if (vit->referred())
+// {
+// //Pout << "REMOVE: " << topoint(vit->point()) << endl;
+// remove(vit);
+// }
+// }
+
referVertices
(
targetProcessor,
@@ -1344,6 +1741,8 @@ void Foam::conformalVoronoiMesh::referVertices
);
}
+ timeCheck("Start of referVertices " + stageName + " insertion.");
+
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
const labelList& constructMap = pointMap.constructMap()[procI];
@@ -1397,7 +1796,7 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
hitSurfaceLargest = -1;
std::list<Facet> facets;
- incident_facets(vit, std::back_inserter(facets));
+ finite_incident_facets(vit, std::back_inserter(facets));
const Foam::point vert = topoint(vit->point());
@@ -1410,52 +1809,45 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceProtrusion
++fit
)
{
- if
+ const Foam::point edgeMid =
+ 0.5
+ *(
+ fit->first->dual()
+ + fit->first->neighbor(fit->second)->dual()
+ );
+
+ pointIndexHit surfHit;
+ label hitSurface;
+
+ geometryToConformTo_.findSurfaceAnyIntersection
(
- !is_infinite(fit->first)
- && !is_infinite(fit->first->neighbor(fit->second))
- )
- {
- const Foam::point edgeMid =
- 0.5
- *(
- topoint(dual(fit->first))
- + topoint(dual(fit->first->neighbor(fit->second)))
- );
+ vert,
+ edgeMid,
+ surfHit,
+ hitSurface
+ );
- pointIndexHit surfHit;
- label hitSurface;
+ if (surfHit.hit())
+ {
+ vectorField norm(1);
- geometryToConformTo_.findSurfaceAnyIntersection
+ allGeometry_[hitSurface].getNormal
(
- vert,
- edgeMid,
- surfHit,
- hitSurface
+ List<pointIndexHit>(1, surfHit),
+ norm
);
- if (surfHit.hit())
- {
- vectorField norm(1);
-
- allGeometry_[hitSurface].getNormal
- (
- List<pointIndexHit>(1, surfHit),
- norm
- );
-
- const vector& n = norm[0];
+ const vector& n = norm[0];
- const scalar normalProtrusionDistance =
- (edgeMid - surfHit.hitPoint()) & n;
+ const scalar normalProtrusionDistance =
+ (edgeMid - surfHit.hitPoint()) & n;
- if (normalProtrusionDistance > maxProtrusionDistance)
- {
- surfHitLargest = surfHit;
- hitSurfaceLargest = hitSurface;
+ if (normalProtrusionDistance > maxProtrusionDistance)
+ {
+ surfHitLargest = surfHit;
+ hitSurfaceLargest = hitSurface;
- maxProtrusionDistance = normalProtrusionDistance;
- }
+ maxProtrusionDistance = normalProtrusionDistance;
}
}
}
@@ -1489,9 +1881,9 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
hitSurfaceLargest = -1;
std::list<Facet> facets;
- incident_facets(vit, std::back_inserter(facets));
+ finite_incident_facets(vit, std::back_inserter(facets));
- Foam::point vert(topoint(vit->point()));
+ const Foam::point vert = topoint(vit->point());
scalar minIncursionDistance = -maxSurfaceProtrusion(vert);
@@ -1502,57 +1894,50 @@ void Foam::conformalVoronoiMesh::dualCellLargestSurfaceIncursion
++fit
)
{
- if
+ const Foam::point edgeMid =
+ 0.5
+ *(
+ fit->first->dual()
+ + fit->first->neighbor(fit->second)->dual()
+ );
+
+ pointIndexHit surfHit;
+ label hitSurface;
+
+ geometryToConformTo_.findSurfaceAnyIntersection
(
- !is_infinite(fit->first)
- && !is_infinite(fit->first->neighbor(fit->second))
- )
- {
- Foam::point edgeMid =
- 0.5
- *(
- topoint(dual(fit->first))
- + topoint(dual(fit->first->neighbor(fit->second)))
- );
+ vert,
+ edgeMid,
+ surfHit,
+ hitSurface
+ );
- pointIndexHit surfHit;
- label hitSurface;
+ if (surfHit.hit())
+ {
+ vectorField norm(1);
- geometryToConformTo_.findSurfaceAnyIntersection
+ allGeometry_[hitSurface].getNormal
(
- vert,
- edgeMid,
- surfHit,
- hitSurface
+ List<pointIndexHit>(1, surfHit),
+ norm
);
- if (surfHit.hit())
- {
- vectorField norm(1);
-
- allGeometry_[hitSurface].getNormal
- (
- List<pointIndexHit>(1, surfHit),
- norm
- );
-
- const vector& n = norm[0];
+ const vector& n = norm[0];
- scalar normalIncursionDistance =
- (edgeMid - surfHit.hitPoint()) & n;
+ scalar normalIncursionDistance =
+ (edgeMid - surfHit.hitPoint()) & n;
- if (normalIncursionDistance < minIncursionDistance)
- {
- surfHitLargest = surfHit;
- hitSurfaceLargest = hitSurface;
+ if (normalIncursionDistance < minIncursionDistance)
+ {
+ surfHitLargest = surfHit;
+ hitSurfaceLargest = hitSurface;
- minIncursionDistance = normalIncursionDistance;
+ minIncursionDistance = normalIncursionDistance;
- // Info<< nl << "# Incursion: " << endl;
- // meshTools::writeOBJ(Info, vert);
- // meshTools::writeOBJ(Info, edgeMid);
- // Info<< "l Na Nb" << endl;
- }
+ // Info<< nl << "# Incursion: " << endl;
+ // meshTools::writeOBJ(Info, vert);
+ // meshTools::writeOBJ(Info, edgeMid);
+ // Info<< "l Na Nb" << endl;
}
}
}
@@ -1948,6 +2333,35 @@ bool Foam::conformalVoronoiMesh::appendToEdgeLocationTree
}
+Foam::List<Foam::pointIndexHit>
+Foam::conformalVoronoiMesh::nearestFeatureEdgeLocations
+(
+ const Foam::point& pt
+) const
+{
+ const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
+
+ labelList elems
+ = edgeLocationTreePtr_().findSphere(pt, exclusionRangeSqr);
+
+ DynamicList<pointIndexHit> dynPointHit;
+
+ forAll(elems, elemI)
+ {
+ label index = elems[elemI];
+
+ const Foam::point& pointI
+ = edgeLocationTreePtr_().shapes().shapePoints()[index];
+
+ pointIndexHit nearHit(true, pointI, index);
+
+ dynPointHit.append(nearHit);
+ }
+
+ return dynPointHit;
+}
+
+
bool Foam::conformalVoronoiMesh::pointIsNearFeatureEdgeLocation
(
const Foam::point& pt
@@ -2009,103 +2423,77 @@ bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
DynamicList<Foam::point>& existingEdgeLocations
) const
{
- const Foam::point pt = pHit.hitPoint();
+ Foam::point pt = pHit.hitPoint();
const scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
- pointIndexHit info;
-
- bool closeToFeatureEdge = pointIsNearFeatureEdgeLocation(pt, info);
+ bool closeToFeatureEdge = pointIsNearFeatureEdgeLocation(pt);
- if (!closeToFeatureEdge)
- {
- appendToEdgeLocationTree(pt, existingEdgeLocations);
- }
- else
+ if (closeToFeatureEdge)
{
- // 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;
+ List<pointIndexHit> nearHits = nearestFeatureEdgeLocations(pt);
- geometryToConformTo_.findEdgeNearest
- (
- pt,
- exclusionRangeSqr,
- edgeHit,
- featureHit
- );
+ forAll(nearHits, elemI)
+ {
+ pointIndexHit& info = nearHits[elemI];
- const extendedFeatureEdgeMesh& eMesh
- = geometryToConformTo_.features()[featureHit];
+ // 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;
- const vector& edgeDir = eMesh.edgeDirections()[edgeHit.index()];
+ geometryToConformTo_.findEdgeNearest
+ (
+ pt,
+ exclusionRangeSqr,
+ edgeHit,
+ featureHit
+ );
- const vector lineBetweenPoints = pt - info.hitPoint();
+ const extendedFeatureEdgeMesh& eMesh
+ = geometryToConformTo_.features()[featureHit];
- const scalar cosAngle = vectorTools::cosPhi(edgeDir, lineBetweenPoints);
+ const vector& edgeDir = eMesh.edgeDirections()[edgeHit.index()];
- // 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
- )
- {
- closeToFeatureEdge = false;
- appendToEdgeLocationTree(pt, existingEdgeLocations);
- }
- }
+ const vector lineBetweenPoints = pt - info.hitPoint();
- return closeToFeatureEdge;
+ const scalar cosAngle
+ = vectorTools::cosPhi(edgeDir, lineBetweenPoints);
- // Searching for the nearest point in existingEdgeLocations using the
- // indexedOctree
+ // 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;
- // 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();
-// }
+ if
+ (
+ mag(cosAngle) < searchConeAngle
+ && (
+ mag(lineBetweenPoints)
+ > cvMeshControls().pointPairDistanceCoeff()*targetCellSize(pt)
+ )
+ )
+ {
+ pt = edgeHit.hitPoint();
+ pHit.setPoint(pt);
+ closeToFeatureEdge = false;
+ }
+ else
+ {
+ closeToFeatureEdge = true;
+ break;
+ }
+ }
+ }
+
+ if (!closeToFeatureEdge)
+ {
+ appendToEdgeLocationTree(pt, existingEdgeLocations);
+ }
+ return closeToFeatureEdge;
}
@@ -2265,12 +2653,6 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
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)
{
const label featureHit = featuresHit[i];
@@ -2281,39 +2663,43 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
{
pointIndexHit& edHit = edHits[eHitI];
- if (!nearFeaturePt(edHit.hitPoint()) && keepSurfacePoint)
+ if (edHit.hit())
{
- if
- (
- magSqr(edHit.hitPoint() - surfHitI.hitPoint())
- < surfacePtReplaceDistCoeffSqr*cellSizeSqr
- )
+ if (!nearFeaturePt(edHit.hitPoint()))
{
- // 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;
- // NEED TO REMOVE FROM THE SURFACE TREE...
- }
+ // NEED TO REMOVE FROM THE SURFACE TREE...
+ surfacePtLocationTreePtr_().remove
+ (
+ existingSurfacePtLocations.size()
+ );
+ }
- if
- (
- !nearFeatureEdgeLocation
+ if
(
- edHit,
- existingEdgeLocations
+ !nearFeatureEdgeLocation
+ (
+ edHit,
+ existingEdgeLocations
+ )
)
- )
- {
- // Do not place edge control points too close to a
- // feature point or existing edge control points
-
- featureEdgeHits.append(edHit);
- featureEdgeFeaturesHit.append(featureHit);
-
- currentEdgeLocations.append(edHit.hitPoint());
+ {
+ // Do not place edge control points too close to a
+ // feature point or existing edge control points
+ featureEdgeHits.append(edHit);
+ featureEdgeFeaturesHit.append(featureHit);
+ }
}
}
}
@@ -2322,7 +2708,6 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
if (keepSurfacePoint)
{
surfaceHits.append(surfHitI);
-
hitSurfaces.append(hitSurfaceI);
}
}
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
index ab7ce33f8563760f144c6aff33f6893a1a7e0078..89afb87134a98ce7b7c9c5ec37afe4764ed4f4b2 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshIO.C
@@ -125,7 +125,7 @@ void Foam::conformalVoronoiMesh::drawDelaunayCell
os << "# cell index: " << label(c->cellIndex()) << endl;
os << "# circumradius "
- << mag(topoint(dual(c)) - topoint(c->vertex(0)->point()))
+ << mag(c->dual() - topoint(c->vertex(0)->point()))
<< endl;
for (int i = 0; i < 4; i++)
@@ -144,7 +144,7 @@ void Foam::conformalVoronoiMesh::drawDelaunayCell
os << "# cicumcentre " << endl;
- meshTools::writeOBJ(os, topoint(dual(c)));
+ meshTools::writeOBJ(os, c->dual());
os << "l " << 1 + offset << " " << 5 + offset << endl;
}
@@ -167,7 +167,7 @@ void Foam::conformalVoronoiMesh::writePoints
++vit
)
{
- if (!internalOnly || vit->internalOrBoundaryPoint())
+ if (!internalOnly || vit->internalPoint())
{
meshTools::writeOBJ(str, topoint(vit->point()));
}
@@ -241,7 +241,7 @@ void Foam::conformalVoronoiMesh::writeProcessorInterface
{
if (!cit->farCell())
{
- points[cit->cellIndex()] = topoint(dual(cit));
+ points[cit->cellIndex()] = cit->dual();
}
}
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
index 7ab4f10853427619285ce87430f1ccf14cd31e42..5fc7e970b1131c78607d13f71c91ff3b90aeb276 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell.H
@@ -51,6 +51,7 @@ SourceFiles
#include "Swap.H"
#include "InfoProxy.H"
#include "tetCell.H"
+#include "typeInfo.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -154,6 +155,8 @@ public:
inline int cellIndex() const;
+ inline const Foam::point& dual();
+
inline bool farCell() const;
inline int& filterCount();
@@ -163,6 +166,11 @@ public:
//- Is the Delaunay cell real, i.e. any real vertex
inline bool real() const;
+ //- Does the Delaunay cell have a far point
+ inline bool hasFarPoint() const;
+
+ inline bool hasInternalPoint() const;
+
//- Does the Dual vertex form part of a processor patch
inline bool parallelDualVertex() const;
@@ -190,6 +198,8 @@ public:
// least one Delaunay vertex outside and at least one inside
inline bool boundaryDualVertex() const;
+ inline bool nearProcBoundary() const;
+
// Info
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H
index 5fbd378f9d9a220845245127070d9624f2ac9bba..6450b01908a44149dbde6d8abc6d8e927f6083b1 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCellI.H
@@ -86,6 +86,24 @@ int CGAL::indexedCell<Gt, Cb>::cellIndex() const
}
+template<class Gt, class Cb>
+const Foam::point& CGAL::indexedCell<Gt, Cb>::dual()
+{
+#ifdef CGAL_INEXACT
+ return reinterpret_cast<const Foam::point&>(this->circumcenter());
+#else
+ const typename Gt::Point_3& P = this->circumcenter();
+
+ return
+ (
+ CGAL::to_double(P.x()),
+ CGAL::to_double(P.y()),
+ CGAL::to_double(P.z())
+ );
+#endif
+}
+
+
template<class Gt, class Cb>
inline bool CGAL::indexedCell<Gt, Cb>::farCell() const
{
@@ -112,10 +130,45 @@ 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()
+ (
+ this->vertex(0)->real()
+ || this->vertex(1)->real()
+ || this->vertex(2)->real()
+ || this->vertex(3)->real()
+ )
+ &&
+ !(
+ this->vertex(0)->farPoint()
+ || this->vertex(1)->farPoint()
+ || this->vertex(2)->farPoint()
+ || this->vertex(3)->farPoint()
+ )
+ );
+}
+
+
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::hasFarPoint() const
+{
+ return
+ (
+ this->vertex(0)->farPoint()
+ || this->vertex(1)->farPoint()
+ || this->vertex(2)->farPoint()
+ || this->vertex(3)->farPoint()
+ );
+}
+
+
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::hasInternalPoint() const
+{
+ return
+ (
+ this->vertex(0)->internalPoint()
+ || this->vertex(1)->internalPoint()
+ || this->vertex(2)->internalPoint()
+ || this->vertex(3)->internalPoint()
);
}
@@ -285,4 +338,17 @@ inline bool CGAL::indexedCell<Gt, Cb>::boundaryDualVertex() const
}
+template<class Gt, class Cb>
+inline bool CGAL::indexedCell<Gt, Cb>::nearProcBoundary() const
+{
+ return
+ (
+ this->vertex(0)->nearProcBoundary()
+ || this->vertex(1)->nearProcBoundary()
+ || this->vertex(2)->nearProcBoundary()
+ || this->vertex(3)->nearProcBoundary()
+ );
+}
+
+
// * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * //
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.C
index 9ad21f92feade712950e9aed5d5cb6d1a21a3eaf..2173e7f5fa7b6be6a0469b004a46e9148f4195f7 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.C
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.C
@@ -29,7 +29,6 @@ License
\*---------------------------------------------------------------------------*/
#include "indexedVertex.H"
-//#include "conformalVoronoiMesh.H"
#include "point.H"
// * * * * * * * * * * * * * * * * IOStream operators * * * * * * * * * * * //
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
index 246020fd172bbeb127089dc9bcc085c0aa7908e2..525c2156e82b60682ca93eb20d44fc4d34d051d7 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertex.H
@@ -45,6 +45,7 @@ SourceFiles
#include <CGAL/Triangulation_3.h>
#include "tensor.H"
#include "InfoProxy.H"
+#include "point.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -101,6 +102,8 @@ class indexedVertex
//- Specify whether the vertex is fixed or movable.
bool vertexFixed_;
+ bool nearProcBoundary_;
+
public:
@@ -198,6 +201,12 @@ public:
//- Set the point to be near the boundary
inline void setNearBoundary();
+ //- Is point internal and near a proc boundary
+ inline bool nearProcBoundary() const;
+
+ //- Set the point to be near a proc boundary
+ inline void setNearProcBoundary();
+
//- Either master or slave of pointPair.
inline bool pairPoint() const;
@@ -227,7 +236,7 @@ public:
inline bool isVertexFixed() const;
//- Fix the vertex so that it can't be moved
- inline void setVertexFixed() const;
+ inline void setVertexFixed();
// inline void operator=(const Delaunay::Finite_vertices_iterator vit)
// {
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H
index c02ec4df2009489efa9fe91221a666b9f9ae0050..1e225eaa2ed02a2674b57a713fbc142f05d451b4 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedVertexI.H
@@ -38,7 +38,8 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex()
type_(vtInternal),
alignment_(),
targetCellSize_(0.0),
- vertexFixed_(false)
+ vertexFixed_(false),
+ nearProcBoundary_(false)
{}
@@ -50,7 +51,8 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(const Point& p)
type_(vtInternal),
alignment_(),
targetCellSize_(0.0),
- vertexFixed_(false)
+ vertexFixed_(false),
+ nearProcBoundary_(false)
{}
@@ -67,7 +69,8 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex
type_(type),
alignment_(),
targetCellSize_(0.0),
- vertexFixed_(false)
+ vertexFixed_(false),
+ nearProcBoundary_(false)
{}
@@ -79,7 +82,8 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(const Point& p, Cell_handle f)
type_(vtInternal),
alignment_(),
targetCellSize_(0.0),
- vertexFixed_(false)
+ vertexFixed_(false),
+ nearProcBoundary_(false)
{}
@@ -91,7 +95,8 @@ inline CGAL::indexedVertex<Gt, Vb>::indexedVertex(Cell_handle f)
type_(vtInternal),
alignment_(),
targetCellSize_(0.0),
- vertexFixed_(false)
+ vertexFixed_(false),
+ nearProcBoundary_(false)
{}
@@ -251,6 +256,20 @@ inline void CGAL::indexedVertex<Gt, Vb>::setNearBoundary()
}
+template<class Gt, class Vb>
+inline bool CGAL::indexedVertex<Gt, Vb>::nearProcBoundary() const
+{
+ return nearProcBoundary_;
+}
+
+
+template<class Gt, class Vb>
+inline void CGAL::indexedVertex<Gt, Vb>::setNearProcBoundary()
+{
+ nearProcBoundary_ = true;
+}
+
+
template<class Gt, class Vb>
inline bool CGAL::indexedVertex<Gt, Vb>::pairPoint() const
{
@@ -331,7 +350,7 @@ inline bool CGAL::indexedVertex<Gt, Vb>::isVertexFixed() const
template<class Gt, class Vb>
-inline void CGAL::indexedVertex<Gt, Vb>::setVertexFixed() const
+inline void CGAL::indexedVertex<Gt, Vb>::setVertexFixed()
{
vertexFixed_ = true;
}
diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
index d8e5fd016616633358be21ac667efdac3eab22bf..ce56a977aa8e676c6f30f655a2730a90f2e4c3d8 100644
--- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
+++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformationSurfaces/conformationSurfaces.H
@@ -285,7 +285,9 @@ public:
) const;
//- Find the nearest points on each feature edge that is within
- // a given distance from the sample point
+ // a given distance from the sample point. Will need to check for
+ // a hit or a miss because near edges may not have a nearest point
+ // on them which is perpendicular to the sample point.
void findAllNearestEdges
(
const point& sample,
diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
index 17ef34fb6f710730e879bcdc4a2005d7b320bc07..024f7d27ecf4203c1684e3295bc7e20ad52692fc 100644
--- a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
+++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C
@@ -379,6 +379,124 @@ Foam::label Foam::indexedOctree<Type>::compactContents
}
+template <class Type>
+bool Foam::indexedOctree<Type>::quickCircumsphereRejection
+(
+ const label nodeI,
+ const point& cc,
+ const scalar crSqr,
+ const List<scalar>& nearestDistances
+) const
+{
+ const node& nod = nodes_[nodeI];
+
+ volumeType nodeType = volumeType(nodeTypes_.get(nodeI<<3));
+
+ //scalar boxDist = nearestDistances[nodeI] + 0.5*nod.bb_.mag();
+ scalar boxDist = crSqr + magSqr(cc - nod.bb_.midpoint());
+
+ if
+ (
+ nodeType == INSIDE
+ //&& (crSqr < sqr(boxDist))
+ && (boxDist < sqr(nearestDistances[nodeI]))
+ )
+ {
+ return true;
+ }
+ else
+ {
+ direction octant = nod.bb_.subOctant(cc);
+
+ labelBits index = nod.subNodes_[octant];
+
+ if (isNode(index))
+ {
+ return quickCircumsphereRejection
+ (
+ getNode(index),
+ cc,
+ crSqr,
+ nearestDistances
+ );
+ }
+ else
+ {
+ return false;
+ }
+ }
+}
+
+
+template <class Type>
+bool Foam::indexedOctree<Type>::quickCircumsphereRejection
+(
+ const point& cc,
+ const scalar crSqr,
+ const List<scalar>& nearestDistances
+) const
+{
+ if (nodes_.size())
+ {
+ return quickCircumsphereRejection
+ (
+ 0,
+ cc,
+ crSqr,
+ nearestDistances
+ );
+ }
+
+ return false;
+}
+
+
+template <class Type>
+Foam::scalar
+Foam::indexedOctree<Type>::calcNearestDistance
+(
+ const label nodeI
+) const
+{
+ const node& nod = nodes_[nodeI];
+
+ const point& nodeCentre = nod.bb_.midpoint();
+
+ scalar nearestDistance = 0.0;
+
+ pointIndexHit pHit = findNearest(nodeCentre, sqr(GREAT));
+
+ if (pHit.hit())
+ {
+ nearestDistance = mag(pHit.hitPoint() - nodeCentre);
+ }
+ else
+ {
+ WarningIn("Foam::indexedOctree<Type>::calcNearestDistance(const label)")
+ << "Cannot calculate distance of nearest point on surface from "
+ << "the midpoint of the octree node. Returning distance of zero."
+ << endl;
+ }
+
+ return nearestDistance;
+}
+
+
+template <class Type>
+Foam::List<Foam::scalar>
+Foam::indexedOctree<Type>::calcNearestDistance() const
+{
+ List<scalar> nearestDistances(nodes_.size());
+
+ forAll(nearestDistances, nodeI)
+ {
+ nearestDistances[nodeI] = calcNearestDistance(nodeI);
+ }
+
+ return nearestDistances;
+}
+
+
// Pre-calculates wherever possible the volume status per node/subnode.
// Recurses to determine status of lowest level boxes. Level above is
// combination of octants below.
@@ -540,6 +658,67 @@ Foam::indexedOctree<Type>::getSide
// ~~~~~~~~~~~~~~
//
+
+//template <class Type>
+//bool Foam::indexedOctree<Type>::findAnyOverlap
+//(
+// const label nodeI,
+// const point& sample,
+// const scalar nearestDistSqr
+//) const
+//{
+// const node& nod = nodes_[nodeI];
+//
+// // Determine order to walk through octants
+// FixedList<direction, 8> octantOrder;
+// nod.bb_.searchOrder(sample, octantOrder);
+//
+// // Go into all suboctants (one containing sample first) and update
+// // nearest.
+// for (direction i = 0; i < 8; i++)
+// {
+// direction octant = octantOrder[i];
+//
+// labelBits index = nod.subNodes_[octant];
+//
+// if (isNode(index))
+// {
+// label subNodeI = getNode(index);
+//
+// const treeBoundBox& subBb = nodes_[subNodeI].bb_;
+//
+// if (overlaps(subBb.min(), subBb.max(), nearestDistSqr, sample))
+// {
+// return findAnyOverlap
+// (
+// subNodeI,
+// sample,
+// nearestDistSqr
+// );
+// }
+// }
+// else if (isContent(index))
+// {
+// if
+// (
+// overlaps
+// (
+// nod.bb_,
+// octant,
+// nearestDistSqr,
+// sample
+// )
+// )
+// {
+// return true;
+// }
+// }
+// }
+//
+// return false;
+//}
+
+
// Find nearest point starting from nodeI
template <class Type>
void Foam::indexedOctree<Type>::findNearest
@@ -1614,7 +1793,6 @@ void Foam::indexedOctree<Type>::traverseNode
}
}
-
const node& nod = nodes_[nodeI];
labelBits index = nod.subNodes_[octant];
@@ -1781,6 +1959,11 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
label i = 0;
for (; i < 100000; i++)
{
+// if (isLineInsideOrOutside(nodeI, treeStart, treeEnd))
+// {
+// return hitInfo;
+// }
+
// Ray-trace to end of current node. Updates point (either on triangle
// in case of hit or on node bounding box in case of miss)
@@ -1935,6 +2118,38 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
}
+//template <class Type>
+//bool Foam::indexedOctree<Type>::isLineInsideOrOutside
+//(
+// const label nodeI,
+// const point& start,
+// const point& end
+//) const
+//{
+// const node& nod = nodes_[nodeI];
+//
+// direction startOctant = nod.bb_.subOctant(start);
+// direction endOctant = nod.bb_.subOctant(end);
+//
+// if (startOctant == endOctant)
+// {
+// volumeType startOctantType
+// = volumeType(nodeTypes_.get((nodeI<<3) + startOctant));
+//
+// if
+// (
+// startOctantType == INSIDE || startOctantType == OUTSIDE
+// )
+// {
+// //Info<< nodeI << " | " << start << " " << end << endl;
+// return true;
+// }
+// }
+//
+// return false;
+//}
+
+
// Find first intersection
template <class Type>
Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
@@ -2559,6 +2774,27 @@ Foam::scalar& Foam::indexedOctree<Type>::perturbTol()
}
+//template <class Type>
+//bool Foam::indexedOctree<Type>::findAnyOverlap
+//(
+// const point& sample,
+// const scalar startDistSqr
+//) const
+//{
+// if (nodes_.size())
+// {
+// return findAnyOverlap
+// (
+// 0,
+// sample,
+// startDistSqr
+// );
+// }
+//
+// return false;
+//}
+
+
template <class Type>
Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
(
diff --git a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H
index aa493f34e92c40b7adb599d4a15dac0d0d55af1c..40ee7a8a32cf54d0403a53fd477c4046cecf5cbc 100644
--- a/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H
+++ b/src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H
@@ -201,6 +201,16 @@ private:
label& compactI
);
+ scalar calcNearestDistance(const label nodeI) const;
+
+ bool quickCircumsphereRejection
+ (
+ const label nodeI,
+ const point& cc,
+ const scalar crSqr,
+ const List<scalar>& nearestDistances
+ ) const;
+
//- Determine inside/outside per node (mixed if cannot be
// determined). Only valid for closed shapes.
volumeType calcVolumeType(const label nodeI) const;
@@ -320,6 +330,13 @@ private:
const bool verbose = false
) const;
+// bool isLineInsideOrOutside
+// (
+// const label nodeI,
+// const point& start,
+// const point& end
+// ) const;
+
//- Find any or nearest intersection of line between start and end.
pointIndexHit findLine
(
@@ -532,6 +549,19 @@ public:
const scalar nearestDistSqr
) const;
+// bool findAnyOverlap
+// (
+// const point& sample,
+// const scalar nearestDistSqr
+// ) const;
+//
+// bool findAnyOverlap
+// (
+// const label nodeI,
+// const point& sample,
+// const scalar nearestDistSqr
+// ) const;
+
//- Low level: calculate nearest starting from subnode.
void findNearest
(
@@ -628,6 +658,16 @@ public:
CompareOp& cop
) const;
+ //- Return a list containing the nearest distance of nodes to any
+ // shapes
+ List<scalar> calcNearestDistance() const;
+
+ bool quickCircumsphereRejection
+ (
+ const point& cc,
+ const scalar crSqr,
+ const List<scalar>& nearestDistances
+ ) const;
// Write
diff --git a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
index 0ebb7b58ab6d0e54c7bc91092e02a51df9fadaa8..81f853aa417f1b3f2f0b264437bcc49ee47695d5 100644
--- a/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
+++ b/src/edgeMesh/extendedFeatureEdgeMesh/extendedFeatureEdgeMesh.C
@@ -546,6 +546,49 @@ Foam::extendedFeatureEdgeMesh::extendedFeatureEdgeMesh
}
featurePointNormals_ = featurePointNormals;
+
+
+ // Create featurePointEdges_
+ // For each feature point, stores a list of edges which are arranged in
+ // order according to their connectivity
+
+// featurePointEdges_.setSize(nonFeatureStart_);
+//
+// Info<< sFeatEds.size() << " " << surf.pointEdges().size() << " "
+// << edges().size() << endl;
+//
+// forAll(sFeatEds, eI)
+// {
+// Info<< "Edge " << eI << " " << sFeatEds[eI] << " "
+// << edges()[eI] << endl;
+// }
+//
+// const edgeList& edges = eds;
+//
+// for (label i = 0; i < nonFeatureStart_; i++)
+// {
+// const labelList& ptEds = surf.pointEdges()[i];
+//
+// DynamicList<label> tmpFtEdges;
+//
+// forAll(ptEds, eI)
+// {
+// const label edgeI = ptEds[eI];
+// const edge& e = sFeatEds[edgeI];
+//
+// Info<< "Edges: " << edgeI << " " << e << endl;
+//
+// forAll(sFeatEds, fEdgeI)
+// {
+// if (edges[fEdgeI] == e)
+// {
+// tmpFtEdges.append(fEdgeI);
+// }
+// }
+// }
+//
+// featurePointEdges_[i] = tmpFtEdges;
+// }
}
@@ -821,7 +864,7 @@ void Foam::extendedFeatureEdgeMesh::allNearestFeatureEdges
if (!hitPoint.hit())
{
- nearHit = pointIndexHit(true, hitPoint.missPoint(), hitIndex);
+ nearHit = pointIndexHit(false, hitPoint.missPoint(), hitIndex);
}
else
{