BUG: cvMesh: added tolerance for deciding whether vectors are parallel.

There was a bug in feature point handling. normals that were effectively
parallel were not being picked up as being parallel, so a tolerance
has been added as a static const.

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 \

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  * * * * * * * * * * * * //
 

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMesh.H

@@ -121,6 +121,8 @@ public:
 
 private:
 
+        static const scalar tolParallel;
+
     // Private data
 
         //- The time registry of the application
@@ -658,7 +660,7 @@ private:
         //- edge conformation location
         bool nearFeatureEdgeLocation
         (
-            const Foam::point& pt,
+            pointIndexHit& pHit,
             DynamicList<Foam::point>& newEdgeLocations,
             pointField& existingEdgeLocations,
             autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C

@@ -1821,40 +1821,45 @@ void Foam::conformalVoronoiMesh::limitDisplacement
 
 bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
 (
-    const Foam::point& pt,
+    pointIndexHit& pHit,
     DynamicList<Foam::point>& newEdgeLocations,
     pointField& existingEdgeLocations,
     autoPtr<indexedOctree<treeDataPoint> >& edgeLocationTree
 ) const
 {
+    const Foam::point pt = pHit.hitPoint();
+
     scalar exclusionRangeSqr = featureEdgeExclusionDistanceSqr(pt);
 
     // 0.01 and 1000 determined from speed tests, varying the indexedOctree
     // rebuild frequency and balance of additions to queries.
 
-    if
+//    if
-    (
+//    (
-        newEdgeLocations.size()
+//        newEdgeLocations.size()
-     >= max(0.01*existingEdgeLocations.size(), 1000)
+//     >= max(0.01*existingEdgeLocations.size(), 1000)
-    )
+//    )
-    {
+//    {
+        const label originalSize = existingEdgeLocations.size();
+
         existingEdgeLocations.append(newEdgeLocations);
 
         buildEdgeLocationTree(edgeLocationTree, existingEdgeLocations);
 
         newEdgeLocations.clear();
-    }
+//    }
-    else
+//    else
-    {
+//    {
-        // Search for the nearest point in newEdgeLocations.
+//        // Search for the nearest point in newEdgeLocations.
-        // Searching here first, because the intention is that the value of
+//        // Searching here first, because the intention is that the value of
-        // newEdgeLocationsSizeLimit should make this faster by design.
+//        // newEdgeLocationsSizeLimit should make this faster by design.
-
+//
-        if (min(magSqr(newEdgeLocations - pt)) <= exclusionRangeSqr)
+//        if (min(magSqr(newEdgeLocations - pt)) <= exclusionRangeSqr)
-        {
+//        {
-            return true;
+//            // Average the points...
-        }
+//            return true;
-    }
+//        }
+//    }
 
     // Searching for the nearest point in existingEdgeLocations using the
     // indexedOctree
@@ -1862,6 +1867,44 @@ bool Foam::conformalVoronoiMesh::nearFeatureEdgeLocation
     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)
+//        {
+//            const Foam::point& existingPoint = existingEdgeLocations[pI];
+//
+//            if (existingPoint == info.hitPoint())
+//            {
+//                toRemove[pI] = true;
+//                //Info<< "Replace " << pt << " and " << info.hitPoint() << " with " << newPt << endl;
+//            }
+//        }
+//
+//        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();
 }
@@ -2048,7 +2091,7 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
                     (
                         !nearFeatureEdgeLocation
                         (
-                            edHit.hitPoint(),
+                            edHit,
                             newEdgeLocations,
                             existingEdgeLocations,
                             edgeLocationTree
@@ -2114,9 +2157,9 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
                 continue;
             }
 
-            const pointIndexHit& edHit = edHits[i];
+            pointIndexHit& edHit = edHits[i];
 
-            const label featureHit = featuresHit[i];
+            label featureHit = featuresHit[i];
 
             if (edHit.hit())
             {
@@ -2144,7 +2187,7 @@ void Foam::conformalVoronoiMesh::addSurfaceAndEdgeHits
                     (
                         !nearFeatureEdgeLocation
                         (
-                            edHit.hitPoint(),
+                            edHit,
                             newEdgeLocations,
                             existingEdgeLocations,
                             edgeLocationTree

applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePointSpecialisations.C

@@ -248,7 +248,7 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
                 const vector& convexNormal
                     = normals[convexEdgeANormals[edgeAnormalI]];
 
-                if (areParallel(concaveNormal, convexNormal))
+                if (areParallel(concaveNormal, convexNormal, tolParallel))
                 {
                     convexEdgeA = true;
                 }
@@ -263,7 +263,7 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
 
                 // Need a looser tolerance, because sometimes adjacent triangles
                 // on the same surface will be slightly out of alignment.
-                if (areParallel(concaveNormal, convexNormal))
+                if (areParallel(concaveNormal, convexNormal, tolParallel))
                 {
                     convexEdgeB = true;
                 }
@@ -291,8 +291,7 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
 
                     if
                     (
-                        areObtuse(concaveNormal, convexNormal)
+                        !areParallel(concaveNormal, convexNormal, tolParallel)
-                     || areOrthogonal(concaveNormal, convexNormal)
                     )
                     {
                         convexEdgePlaneCNormal = convexNormal;
@@ -322,8 +321,7 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
 
                     if
                     (
-                        areObtuse(concaveNormal, convexNormal)
+                        !areParallel(concaveNormal, convexNormal, tolParallel)
-                     || areOrthogonal(concaveNormal, convexNormal)
                     )
                     {
                         convexEdgePlaneDNormal = convexNormal;
@@ -400,6 +398,7 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
             pts.append(externalPtG);
             indices.append(0);
             types.append(-1);
+        }
 
         if (debug)
         {
@@ -409,7 +408,6 @@ bool Foam::conformalVoronoiMesh::createSpecialisedFeaturePoint
                 meshTools::writeOBJ(Info, pts[ptI]);
             }
         }
-        }
 
         return true;
     }

applications/utilities/mesh/generation/cvMesh/vectorTools/vectorTools.H

@@ -21,11 +21,6 @@ License
     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
     Foam::vectorTools
 
@@ -54,8 +49,9 @@ namespace Foam
 //- 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
     (
@@ -71,6 +67,8 @@ namespace vectorTools
     }
 
     //- 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
     (
@@ -127,7 +125,6 @@ namespace vectorTools
     {
         return Foam::radToDeg(radAngleBetween(a, b, tolerance));
     }
-
 } // End namespace vectorTools
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //