An intermediate version of edge capturing + a class for removing user-

git-svn-id: https://pl5.projectlocker.com/igui/meshGeneration/svn@22 fdcce57e-7e00-11e2-b579-49867b4cea03

meshLibrary/Make/files

@@ -10,7 +10,7 @@ meshSurfaceCheckEdgeTypes = utilities/surfaceTools/meshSurfaceCheckEdgeTypes
 meshSurfaceCutter  =  utilities/surfaceTools/meshSurfaceCutter
 meshSurfaceMapper  =  utilities/surfaceTools/meshSurfaceMapper
 edgeExtraction = utilities/surfaceTools/edgeExtraction
-edgeExtractionFunctions = $(edgeExtraction)/edgeExtractionFunctions
+edgeExtractor = $(edgeExtraction)/edgeExtractor
 meshSurfaceEdgeExtractor = utilities/surfaceTools/meshSurfaceEdgeExtractor
 meshSurfaceEdgeExtractorNonTopo = utilities/surfaceTools/meshSurfaceEdgeExtractorNonTopo
 meshSurfaceEdgeExtractorFUN = utilities/surfaceTools/meshSurfaceEdgeExtractorFUN
@@ -79,6 +79,7 @@ triSurfaceDetectPlanarRegions = utilities/triSurfaceTools/triSurfaceDetectPlanar
 triSurfaceDetectFeatureEdges = utilities/triSurfaceTools/triSurfaceDetectFeatureEdges
 triSurfaceClassifyEdges = utilities/triSurfaceTools/triSurfaceClassifyEdges
 triSurfacePatchManipulator = utilities/triSurfaceTools/triSurfacePatchManipulator
+triSurfaceRemoveFacets = utilities/triSurfaceTools/triSurfaceRemoveFacets
 
 polyMeshGen = utilities/meshes/polyMeshGen
 writePatch  = utilities/meshes/polyMeshGen/writePatch
@@ -290,7 +291,7 @@ $(meshSurfaceMapper)/meshSurfaceMapperMapVertices.C
 $(meshSurfaceMapper)/meshSurfaceMapperCornersAndEdges.C
 $(meshSurfaceMapper)/meshSurfaceMapperPremapVertices.C
 
-$(edgeExtractionFunctions)/edgeExtractionFunctions.C
+$(edgeExtractor)/edgeExtractor.C
 
 $(meshSurfaceEdgeExtractor)/meshSurfaceEdgeExtractor.C
 $(meshSurfaceEdgeExtractor)/meshSurfaceEdgeCreateEdgeVertices.C
@@ -425,6 +426,9 @@ $(triSurfaceClassifyEdges)/triSurfaceClassifyEdgesFunctions.C
 $(triSurfacePatchManipulator)/triSurfacePatchManipulator.C
 $(triSurfacePatchManipulator)/triSurfacePatchManipulatorFunctions.C
 
+$(triSurfaceRemoveFacets)/triSurfaceRemoveFacets.C
+$(triSurfaceRemoveFacets)/triSurfaceRemoveFacetsFunctions.C
+
 $(cartesianMeshExtractor)/cartesianMeshExtractor.C
 $(cartesianMeshExtractor)/cartesianMeshExtractorPointsAndAddressing.C
 $(cartesianMeshExtractor)/cartesianMeshExtractorPolyMesh.C

meshLibrary/cartesianMesh/cartesianMeshGenerator/cartesianMeshGenerator.C

@@ -351,6 +351,8 @@ cartesianMeshGenerator::cartesianMeshGenerator(const Time& time)
 
     surfacePtr_ = new triSurf(db_.path()/surfaceFile);
 
+    surfacePtr_->writeSurface("geom.fms");
+
     if( surfacePtr_->featureEdges().size() != 0 )
     {
         //- create surface patches based on the feature edges

meshLibrary/utilities/helperClasses/parallelHelpers/labelledScalar/labelledScalar.H

@@ -48,18 +48,18 @@ namespace Foam
 /*---------------------------------------------------------------------------*\
                            Class labelledScalar Declaration
 \*---------------------------------------------------------------------------*/
-    
+
 class labelledScalar
 {
     // Private data
         //- label
         label sLabel_;
-    
+
         //- value
         scalar value_;
-    
+
     public:
-        
+
     // Constructors
         //- Null construct
         labelledScalar()
@@ -67,52 +67,84 @@ class labelledScalar
             sLabel_(-1),
             value_(0.0)
         {}
-        
+
         //- Construct from label and value
         labelledScalar(const label sl, const scalar s)
         :
             sLabel_(sl),
             value_(s)
         {}
-            
+
     // Destructor
         ~labelledScalar()
         {}
-            
+
     // Member functions
         //- return scalar label
         inline label scalarLabel() const
         {
             return sLabel_;
         }
-        
+
         //- return the value
         inline const scalar& value() const
         {
             return value_;
         }
-        
+
     // Member operators
-        
+
         inline void operator=(const labelledScalar& ls)
         {
             sLabel_ = ls.sLabel_;
             value_ = ls.value_;
         }
-        
+
         inline bool operator==(const labelledScalar& ls) const
         {
             if( sLabel_ == ls.sLabel_ )
                 return true;
-            
+
             return false;
         }
-        
+
+        inline bool operator<(const labelledScalar& ls) const
+        {
+            if( value_ < ls.value_ )
+                return true;
+
+            return false;
+        }
+
+        inline bool operator<=(const labelledScalar& ls) const
+        {
+            if( value_ <= ls.value_ )
+                return true;
+
+            return false;
+        }
+
+        inline bool operator>(const labelledScalar& ls) const
+        {
+            if( value_ > ls.value_ )
+                return true;
+
+            return false;
+        }
+
+        inline bool operator>=(const labelledScalar& ls) const
+        {
+            if( value_ >= ls.value_ )
+                return true;
+
+            return false;
+        }
+
         inline bool operator!=(const labelledScalar& ls) const
         {
             return !this->operator==(ls);
         }
-        
+
     // Friend operators
         friend Ostream& operator<<(Ostream& os, const labelledScalar& ls)
         {
@@ -125,21 +157,21 @@ class labelledScalar
 
             return os;
         }
-        
+
         friend Istream& operator>>(Istream& is, labelledScalar& ls)
         {
             // Read beginning of labelledScalar
             is.readBegin("labelledScalar");
-        
+
             is >> ls.sLabel_;
             is >> ls.value_;
-            
+
             // Read end of labelledScalar
             is.readEnd("labelledScalar");
-        
+
             // Check state of Istream
             is.check("operator>>(Istream&, labelledScalar");
-            
+
             return is;
         }
 };

meshLibrary/utilities/octrees/meshOctree/meshOctree.H

@@ -171,7 +171,7 @@ public:
             const point& p
         ) const;
 
-        //- find nearest surface edges vertex to the given vertex
+        //- find nearest feature-edges vertex to a given vertex
         bool findNearestEdgePoint
         (
             const point& p,
@@ -188,6 +188,25 @@ public:
             scalar& distSq
         ) const;
 
+        //- find nearest corner point
+        bool findNearestCorner
+        (
+            point& nearest,
+            scalar& distSq,
+            const point& p,
+            const DynList<label>& patches
+        ) const;
+
+        //- find the nearest vertex to the given patches
+        bool findNearestVertexToPatches
+        (
+            point& nearest,
+            scalar& distSq,
+            const point& p,
+            const DynList<label>& patches,
+            const scalar tol = 1e-4
+        ) const;
+
         //- return a reference to the surface
         inline const triSurf& surface() const;
 

meshLibrary/utilities/octrees/meshOctree/meshOctreeFindNearestSurfacePoint.C

@@ -30,6 +30,7 @@ Description
 #include "triSurf.H"
 #include "demandDrivenData.H"
 #include "helperFunctions.H"
+#include "HashSet.H"
 
 // #define DEBUGSearch
 
@@ -352,6 +353,178 @@ bool meshOctree::findNearestVertexToTheEdge
     return found;
 }
 
+bool meshOctree::findNearestCorner
+(
+    point& nearest,
+    scalar& distSq,
+    const point& p,
+    const DynList<label>& patches
+) const
+{
+
+
+    const label cLabel = findLeafContainingVertex(p);
+    vector sizeVec;
+    if( cLabel < 0 )
+    {
+        sizeVec.x() = sizeVec.y() = sizeVec.z() = searchRange_;
+    }
+    else
+    {
+        const scalar s = 1.5 * leaves_[cLabel]->size(rootBox_);
+        sizeVec.x() = sizeVec.y() = sizeVec.z() = s;
+    }
+
+    //- find nearest surface vertex to the point p
+    bool found;
+    label iterationI(0);
+    DynList<const meshOctreeCube*, 256> neighbours;
+
+    const pointField& points = surface_.points();
+    const VRWGraph& pEdges = surface_.pointEdges();
+    const VRWGraph& eFacets = surface_.edgeFacets();
+
+    do
+    {
+        found = false;
+        boundBox bb(p - sizeVec, p + sizeVec);
+        distSq = Foam::sqr(sizeVec.x());
+
+        neighbours.clear();
+        findLeavesContainedInBox(bb, neighbours);
+        labelHashSet checkedPoint;
+
+        //- find nearest projection
+        forAll(neighbours, neiI)
+        {
+            if( !neighbours[neiI]->hasContainedElements() )
+                continue;
+
+            const VRWGraph& ct =
+                neighbours[neiI]->slotPtr()->containedTriangles_;
+            const constRow el = ct[neighbours[neiI]->containedElements()];
+            forAll(el, tI)
+            {
+                const labelledTri& tri = surface_[el[tI]];
+
+                forAll(tri, pI)
+                {
+                    const label spI = tri[pI];
+
+                    if( checkedPoint.found(spI) )
+                        continue;
+
+                    checkedPoint.insert(spI);
+
+                    DynList<label> nodePatches;
+                    label nEdges(0);
+
+                    forAllRow(pEdges, spI, i)
+                    {
+                        const label eI = pEdges(spI, i);
+
+                        if( pEdges.sizeOfRow(eI) < 2 )
+                            break;
+
+                        if(
+                            surface_[eFacets(eI, 0)].region() !=
+                            surface_[eFacets(eI, 1)].region()
+                        )
+                        {
+                            //- found an edge attached to this vertex
+                            ++nEdges;
+                            nodePatches.appendIfNotIn
+                            (
+                                surface_[eFacets(eI, 0)].region()
+                            );
+                            nodePatches.appendIfNotIn
+                            (
+                                surface_[eFacets(eI, 1)].region()
+                            );
+                        }
+                    }
+
+                    if( nEdges > 2 )
+                    {
+                        //- check if all required patches
+                        //- are present at this corner
+                        nEdges = 0;
+                        forAll(patches, i)
+                        {
+                            if( nodePatches.contains(patches[i]) )
+                                ++nEdges;
+                        }
+
+                        if( nEdges >= patches.size() )
+                        {
+                            //- all patches are present, check the distance
+                            const scalar dSq = magSqr(points[spI] - p);
+
+                            if( dSq < distSq )
+                            {
+                                distSq = dSq;
+                                found = true;
+                                nearest = points[spI];
+                            }
+                        }
+                    }
+                }
+            }
+        }
+
+        if( !found )
+            sizeVec *= 2.0;
+
+    } while( !found && (iterationI++ < 5) );
+
+    return found;
+}
+
+bool meshOctree::findNearestVertexToPatches
+(
+    point& nearest,
+    scalar& distSq,
+    const point& p,
+    const DynList<label>& patches,
+    const scalar tol
+) const
+{
+    if( patches.size() == 0 )
+        return false;
+
+    point mapPointApprox(p);
+    scalar distSqApprox;
+    label iter(0);
+    while( iter++ < 20 )
+    {
+        point newP(vector::zero);
+        forAll(patches, patchI)
+        {
+            point np;
+            this->findNearestSurfacePointInRegion
+            (
+                np,
+                distSqApprox,
+                patches[patchI],
+                mapPointApprox
+            );
+
+            newP += np;
+        }
+
+        newP /= patches.size();
+        if( magSqr(newP - mapPointApprox) < tol * searchRange_ )
+            break;
+
+        mapPointApprox = newP;
+    }
+
+    distSq = magSqr(mapPointApprox - p);
+    nearest = mapPointApprox;
+
+    return true;
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

meshLibrary/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizer.H

@@ -64,13 +64,13 @@ class meshSurfaceOptimizer
     // Private data
         //- reference to the mesh surface
         meshSurfaceEngine& surfaceEngine_;
-    
+
         //- reference to the octree
         const meshOctree& meshOctree_;
-    
+
         //- type of surface vertex
         List<direction> vertexType_;
-    
+
         //- form trihedra
         mutable LongList<triFace>* trianglesPtr_;
         mutable VRWGraph* pointTrianglesPtr_;
@@ -78,15 +78,15 @@ class meshSurfaceOptimizer
     // Private member functions
         //- classify surface vertices as PARTITION, EDGE, CORNER
         void classifySurfaceVertices();
-    
+
         //- calculate surface triangulation
         void calculateTrianglesAndAddressing() const;
         inline const LongList<triFace>& triangles() const;
         inline const VRWGraph& pointTriangles() const;
-        
+
         //- find vertices which are inverted
         label findInvertedVertices(boolList& smoothVertex) const;
-    
+
         //- transform into a 2D space in plane
         inline bool transformIntoPlane
         (
@@ -97,7 +97,7 @@ class meshSurfaceOptimizer
             DynList<point>& pts,
             DynList<triFace>& trias
         ) const;
-        
+
         //- transform into a 2D space in plane for parallel boundaries
         inline bool transformIntoPlaneParallel
         (
@@ -109,7 +109,7 @@ class meshSurfaceOptimizer
             DynList<point>& pts,
             DynList<triFace>& trias
         ) const;
-        
+
         //- new position of a node after laplacian smoothing
         //- the position is the average of neighbouring vertex positions
         inline point newPositionLaplacian
@@ -117,7 +117,7 @@ class meshSurfaceOptimizer
             const label bpI,
             const bool transformIntoPlane = true
         ) const;
-        
+
         //- new position of a node after laplacian smoothing
         //- the position is the average of neighbouring face centres
         inline point newPositionLaplacianFC
@@ -125,7 +125,7 @@ class meshSurfaceOptimizer
             const label bpI,
             const bool transformIntoPlane = true
         ) const;
-        
+
         //- new position of a node after laplacian smoothing
         //- the position is the weighted average of neighbouring face centres
         //- the weights are the magnitudes of neighbouring face area vectors
@@ -134,7 +134,7 @@ class meshSurfaceOptimizer
             const label bpI,
             const bool transformIntoPlane = true
         ) const;
-        
+
         //- new position of a node after using surfaceOptimizer
         inline point newPositionSurfaceOptimizer
         (
@@ -142,7 +142,7 @@ class meshSurfaceOptimizer
             const scalar tol = 0.001,
             const bool transformIntoPlane = true
         ) const;
-        
+
         //- new position of edge vertex
         //- the position is the average of neighbouring edge vertices
         inline point newEdgePositionLaplacian
@@ -150,7 +150,7 @@ class meshSurfaceOptimizer
             const label bpI,
             const bool transformIntoPlane = true
         ) const;
-    
+
         //- smooth the node using the laplacian smoother
         //- new position is the average of the neighbouring vertices
         void nodeDisplacementLaplacian
@@ -158,7 +158,7 @@ class meshSurfaceOptimizer
             const label bpI,
             const bool transformIntoPlane = true
         ) const;
-        
+
         //- smooth the node using the laplacian smoother
         //- new position is the average of the centres of faces attached
         //- to the vertex
@@ -167,26 +167,26 @@ class meshSurfaceOptimizer
             const label bpI,
             const bool transformIntoPlane = true
         ) const;
-        
+
         //- smooth the node using surface optimizer
         void nodeDisplacementSurfaceOptimizer
         (
             const label bpI,
             const scalar tol = 1.0
         );
-        
+
         //- smooth edge node
         void edgeNodeDisplacement(const label bpI) const;
-        
+
     // Functions needed for parallel runs
-    
+
         //- transfer data between processors
         void exchangeData
         (
             const labelListPMG& nodesToSmooth,
             std::map<label, DynList<parTriFace> >& m
         ) const;
-    
+
         //- laplacian smoothing of points at processor boundaries
         //- new position is the average of the neighbouring vertices
         void nodeDisplacementLaplacianParallel
@@ -194,20 +194,20 @@ class meshSurfaceOptimizer
             const labelListPMG& nodesToSmooth,
             const bool transformIntoPlane = true
         );
-        
+
         //- laplacian smoothing of points at processor boundaries
         void nodeDisplacementLaplacianFCParallel
         (
             const labelListPMG& nodesToSmooth,
             const bool transformIntoPlane = true
         );
-        
+
         //- smooth nodes at processor boundaries using surface optimizer
         void nodeDisplacementSurfaceOptimizerParallel
         (
             const labelListPMG& nodesToSmooth
         );
-        
+
         //- smooth edge nodes at processor boundaries
         void edgeNodeDisplacementParallel
         (
@@ -219,7 +219,7 @@ class meshSurfaceOptimizer
 
         //- Disallow default bitwise assignment
         void operator=(const meshSurfaceOptimizer&);
-        
+
     // Enumerators
         enum vertexTypes
         {
@@ -247,7 +247,7 @@ public:
 
     // Member Functions
         //- optimisation before the edges and corners are created
-        void preOptimizeSurface();
+        bool preOptimizeSurface();
 
         //- optimize boundary nodes after boundary regions are created
         void optimizeSurface(const label nIterations = 5);

meshLibrary/utilities/smoothers/geometry/meshSurfaceOptimizer/meshSurfaceOptimizerOptimizeSurface.C

@@ -146,10 +146,12 @@ label meshSurfaceOptimizer::findInvertedVertices
     return nInvertedTria;
 }