Salome HOME
54122: Bad quality prismatic mesh
authoreap <eap@opencascade.com>
Thu, 20 Apr 2017 11:16:23 +0000 (14:16 +0300)
committereap <eap@opencascade.com>
Thu, 20 Apr 2017 11:16:23 +0000 (14:16 +0300)
   Extract class Delaunay from class Morph

+ Print warning when creating an algorithm with unexpected arguments
  (expected ones are algo type and geometry) -- smeshBuilder.py

src/SMESHUtils/CMakeLists.txt
src/SMESHUtils/SMESH_Delaunay.cxx [new file with mode: 0644]
src/SMESHUtils/SMESH_Delaunay.hxx [new file with mode: 0644]
src/SMESH_SWIG/smeshBuilder.py
src/StdMeshers/StdMeshers_Prism_3D.cxx
src/StdMeshers/StdMeshers_Prism_3D.hxx
src/StdMeshers/StdMeshers_ProjectionUtils.cxx
src/StdMeshers/StdMeshers_ProjectionUtils.hxx

index 3641d91..c84715e 100644 (file)
@@ -67,6 +67,7 @@ SET(SMESHUtils_HEADERS
   SMESH_MeshAlgos.hxx
   SMESH_MAT2d.hxx
   SMESH_ControlPnt.hxx
+  SMESH_Delaunay.hxx
 )
 
 # --- sources ---
@@ -84,6 +85,7 @@ SET(SMESHUtils_SOURCES
   SMESH_FreeBorders.cxx
   SMESH_ControlPnt.cxx
   SMESH_DeMerge.cxx
+  SMESH_Delaunay.cxx
   )
 
 # --- rules ---
diff --git a/src/SMESHUtils/SMESH_Delaunay.cxx b/src/SMESHUtils/SMESH_Delaunay.cxx
new file mode 100644 (file)
index 0000000..d5bc2f6
--- /dev/null
@@ -0,0 +1,296 @@
+// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library 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
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File      : SMESH_Delaunay.cxx
+// Created   : Wed Apr 19 15:41:15 2017
+// Author    : Edward AGAPOV (eap)
+
+#include "SMESH_Delaunay.hxx"
+#include "SMESH_MeshAlgos.hxx"
+
+#include <BRepAdaptor_Surface.hxx>
+#include <BRepMesh_Delaun.hxx>
+
+//================================================================================
+/*!
+ * \brief Construct a Delaunay triangulation of given boundary nodes
+ *  \param [in] boundaryNodes - vector of nodes of a wire
+ *  \param [in] face - the face
+ *  \param [in] faceID - the face ID
+ *  \param [in] nbNodesToVisit - nb of non-marked nodes on the face
+ */
+//================================================================================
+
+SMESH_Delaunay::SMESH_Delaunay(const std::vector< const UVPtStructVec* > & boundaryNodes,
+                               const TopoDS_Face&                          face,
+                               const int                                   faceID)
+  : _face( face ), _faceID( faceID ), _scale( 1., 1. )
+{
+  // compute _scale
+  BRepAdaptor_Surface surf( face );
+  if ( surf.GetType() != GeomAbs_Plane )
+  {
+    const int nbDiv = 100;
+    const double uRange = surf.LastUParameter() - surf.FirstUParameter();
+    const double vRange = surf.LastVParameter() - surf.FirstVParameter();
+    const double dU = uRange / nbDiv;
+    const double dV = vRange / nbDiv;
+    double u = surf.FirstUParameter(), v = surf.FirstVParameter();
+    gp_Pnt p0U = surf.Value( u, v ), p0V = p0U;
+    double lenU = 0, lenV = 0;
+    for ( ; u < surf.LastUParameter(); u += dU, v += dV )
+    {
+      gp_Pnt p1U = surf.Value( u, surf.FirstVParameter() );
+      lenU += p1U.Distance( p0U );
+      p0U = p1U;
+      gp_Pnt p1V = surf.Value( surf.FirstUParameter(), v );
+      lenV += p1V.Distance( p0V );
+      p0V = p1V;
+    }
+    _scale.SetCoord( lenU / uRange, lenV / vRange );
+  }
+
+  // count boundary points
+  int iP = 1, nbP = 0;
+  for ( size_t iW = 0; iW < boundaryNodes.size(); ++iW ) // loop on wires
+  {
+    nbP += boundaryNodes[iW]->size();
+    if ( boundaryNodes[iW]->front().node == boundaryNodes[iW]->back().node )
+      --nbP; // 1st and last points coincide
+  }
+  _bndNodes.resize( nbP );
+
+  // fill boundary points
+  BRepMesh::Array1OfVertexOfDelaun bndVert( 1, 1 + nbP );
+  BRepMesh_Vertex v( 0, 0, BRepMesh_Frontier );
+  for ( size_t iW = 0; iW < boundaryNodes.size(); ++iW )
+  {
+    const UVPtStructVec& bndPnt = *boundaryNodes[iW];
+    int i = 0, nb = bndPnt.size();
+    if ( bndPnt[0].node == bndPnt.back().node )
+      --nb;
+    for ( ; i < nb;  ++i, ++iP )
+    {
+      _bndNodes[ iP-1 ] = bndPnt[i].node;
+      bndPnt[i].node->setIsMarked( true );
+
+      v.ChangeCoord() = bndPnt[i].UV().Multiplied( _scale );
+      bndVert( iP )   = v;
+    }
+  }
+
+  // triangulate the srcFace in 2D
+  BRepMesh_Delaun Delaunay( bndVert );
+  _triaDS = Delaunay.Result();
+}
+
+//================================================================================
+/*!
+ * \brief Prepare to the exploration of nodes
+ */
+//================================================================================
+
+void SMESH_Delaunay::InitTraversal(const int nbNodesToVisit)
+{
+  _nbNodesToVisit = (size_t) nbNodesToVisit;
+  _nbVisitedNodes = _iBndNode = 0;
+  _noTriQueue.clear();
+}
+
+//================================================================================
+/*!
+ * \brief Return a node with its Barycentric Coordinates within the triangle
+ *        defined by its node indices (zero based)
+ *  \param [out] bc - Barycentric Coordinates of the returned node
+ *  \param [out] triaNodes - indices of triangle nodes
+ *  \return const SMDS_MeshNode* - the next node or NULL
+ */
+//================================================================================
+
+const SMDS_MeshNode* SMESH_Delaunay::NextNode( double bc[3], int triaNodes[3] )
+{
+  while ( _nbVisitedNodes < _nbNodesToVisit )
+  {
+    while ( !_noTriQueue.empty() )
+    {
+      const SMDS_MeshNode*     node = _noTriQueue.front().first;
+      const BRepMesh_Triangle* tria = _noTriQueue.front().second;
+      _noTriQueue.pop_front();
+      if ( node->isMarked() )
+        continue;
+      ++_nbVisitedNodes;
+      node->setIsMarked( true );
+
+      // find a Delaunay triangle containing the src node
+      gp_XY uv = getNodeUV( _face, node );
+      tria = FindTriangle( uv, tria, bc, triaNodes );
+      if ( tria )
+      {
+        addCloseNodes( node, tria, _faceID, _noTriQueue );
+        return node;
+      }
+    }
+    for ( ; _iBndNode < _bndNodes.size() &&  _noTriQueue.empty();  ++_iBndNode )
+    {
+      if ( const BRepMesh_Triangle* tria = GetTriangleNear( _iBndNode ))
+        addCloseNodes( _bndNodes[ _iBndNode ], tria, _faceID, _noTriQueue );
+    }
+    if ( _noTriQueue.empty() )
+      break;
+  }
+
+  // if ( _nbVisitedNodes < _nbNodesToVisit )
+  //   _nbVisitedNodes = std::numeric_limits<int>::max();
+  return NULL;
+}
+
+//================================================================================
+/*!
+ * \brief Find a Delaunay triangle containing a given 2D point and return
+ *        barycentric coordinates within the found triangle
+ */
+//================================================================================
+
+const BRepMesh_Triangle* SMESH_Delaunay::FindTriangle( const gp_XY&             UV,
+                                                       const BRepMesh_Triangle* tria,
+                                                       double                   bc[3],
+                                                       int                      triaNodes[3] )
+{
+  int   nodeIDs[3];
+  gp_XY nodeUVs[3];
+  int   linkIDs[3];
+  Standard_Boolean ori[3];
+
+  gp_XY uv = UV.Multiplied( _scale );
+
+  while ( tria )
+  {
+    // check if the uv is in tria
+
+    _triaDS->ElementNodes( *tria, nodeIDs );
+    nodeUVs[0] = _triaDS->GetNode( nodeIDs[0] ).Coord();
+    nodeUVs[1] = _triaDS->GetNode( nodeIDs[1] ).Coord();
+    nodeUVs[2] = _triaDS->GetNode( nodeIDs[2] ).Coord();
+
+    SMESH_MeshAlgos::GetBarycentricCoords( uv,
+                                           nodeUVs[0], nodeUVs[1], nodeUVs[2],
+                                           bc[0], bc[1] );
+    if ( bc[0] >= 0 && bc[1] >= 0 && bc[0] + bc[1] <= 1 )
+    {
+      bc[2] = 1 - bc[0] - bc[1];
+      triaNodes[0] = nodeIDs[0] - 1;
+      triaNodes[1] = nodeIDs[1] - 1;
+      triaNodes[2] = nodeIDs[2] - 1;
+      return tria;
+    }
+
+    // look for a neighbor triangle, which is adjacent to a link intersected
+    // by a segment( triangle center -> uv )
+
+    gp_XY gc = ( nodeUVs[0] + nodeUVs[1] + nodeUVs[2] ) / 3.;
+    gp_XY seg = uv - gc;
+
+    tria->Edges( linkIDs, ori );
+    int triaID = _triaDS->IndexOf( *tria );
+    tria = 0;
+
+    for ( int i = 0; i < 3; ++i )
+    {
+      const BRepMesh_PairOfIndex & triIDs = _triaDS->ElementsConnectedTo( linkIDs[i] );
+      if ( triIDs.Extent() < 2 )
+        continue; // no neighbor triangle
+
+      // check if a link intersects gc2uv
+      const BRepMesh_Edge & link = _triaDS->GetLink( linkIDs[i] );
+      const BRepMesh_Vertex & n1 = _triaDS->GetNode( link.FirstNode() );
+      const BRepMesh_Vertex & n2 = _triaDS->GetNode( link.LastNode() );
+      gp_XY uv1 = n1.Coord();
+      gp_XY lin = n2.Coord() - uv1; // link direction
+
+      double crossSegLin = seg ^ lin;
+      if ( Abs( crossSegLin ) < std::numeric_limits<double>::min() )
+        continue; // parallel
+
+      double uSeg = ( uv1 - gc ) ^ lin / crossSegLin;
+      if ( 0. <= uSeg && uSeg <= 1. )
+      {
+        tria = & _triaDS->GetElement( triIDs.Index( 1 + ( triIDs.Index(1) == triaID )));
+        break;
+      }
+    }
+  }
+  return tria;
+}
+
+//================================================================================
+/*!
+ * \brief Return a triangle sharing a given boundary node
+ *  \param [in] iBndNode - index of the boundary node
+ *  \return const BRepMesh_Triangle* - a found triangle
+ */
+//================================================================================
+
+const BRepMesh_Triangle* SMESH_Delaunay::GetTriangleNear( int iBndNode )
+{
+  const BRepMesh::ListOfInteger & linkIds = _triaDS->LinksConnectedTo( iBndNode + 1 );
+  const BRepMesh_PairOfIndex &    triaIds = _triaDS->ElementsConnectedTo( linkIds.First() );
+  const BRepMesh_Triangle&           tria = _triaDS->GetElement( triaIds.Index(1) );
+  return &tria;
+}
+
+//================================================================================
+/*!
+ * \brief Return UV of the i-th source boundary node (zero based)
+ */
+//================================================================================
+
+gp_XY SMESH_Delaunay::GetBndUV(const int iNode) const
+{
+  return _triaDS->GetNode( iNode+1 ).Coord();
+}
+
+//================================================================================
+/*!
+ * \brief Add non-marked nodes surrounding a given one to a queue
+ */
+//================================================================================
+
+void SMESH_Delaunay::addCloseNodes( const SMDS_MeshNode*     node,
+                                    const BRepMesh_Triangle* tria,
+                                    const int                faceID,
+                                    TNodeTriaList &          _noTriQueue )
+{
+  // find in-FACE nodes
+  SMDS_ElemIteratorPtr elems = node->GetInverseElementIterator(SMDSAbs_Face);
+  while ( elems->more() )
+  {
+    const SMDS_MeshElement* elem = elems->next();
+    if ( elem->getshapeId() == faceID )
+    {
+      for ( int i = 0, nb = elem->NbNodes(); i < nb; ++i )
+      {
+        const SMDS_MeshNode* n = elem->GetNode( i );
+        if ( !n->isMarked() /*&& n->getshapeId() == faceID*/ )
+          _noTriQueue.push_back( std::make_pair( n, tria ));
+      }
+    }
+  }
+}
diff --git a/src/SMESHUtils/SMESH_Delaunay.hxx b/src/SMESHUtils/SMESH_Delaunay.hxx
new file mode 100644 (file)
index 0000000..8d04f50
--- /dev/null
@@ -0,0 +1,112 @@
+// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
+//
+// Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+//
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+//
+// This library 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
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+//
+// File      : SMESH_Delaunay.hxx
+// Created   : Wed Apr 19 15:42:54 2017
+// Author    : Edward AGAPOV (eap)
+
+
+#ifndef __SMESH_Delaunay_HXX__
+#define __SMESH_Delaunay_HXX__
+
+#include "SMESH_TypeDefs.hxx"
+
+#include <BRepMesh_DataStructureOfDelaun.hxx>
+
+/*!
+ * \brief Create a Delaunay triangulation of nodes on a face boundary
+ *        and provide exploration of nodes shared by elements lying on
+ *        the face. For a returned node, also return a Delaunay triangle
+ *        the node lies in and its Barycentric Coordinates within the triangle.
+ *        Only non-marked nodes are visited. Boundary nodes given at the construction
+ *        are not returned.
+ *
+ *        For usage, this class needs to be subclassed to implement getNodeUV();
+ */
+class SMESHUtils_EXPORT SMESH_Delaunay
+{
+ public:
+
+  // construct a Delaunay triangulation of given boundary nodes
+  SMESH_Delaunay(const std::vector< const UVPtStructVec* > & boundaryNodes,
+                 const TopoDS_Face&                          face,
+                 const int                                   faceID);
+
+  virtual ~SMESH_Delaunay() {}
+
+  // prepare to the exploration of nodes
+  void InitTraversal(const int nbNodesToVisit = -1);
+
+  // return a node with its Barycentric Coordinates within the triangle
+  // defined by its node indices (zero based)
+  const SMDS_MeshNode* NextNode( double bc[3], int triaNodes[3] );
+
+  // return nb of nodes returned by NextNode()
+  int NbVisitedNodes() const { return _nbVisitedNodes; }
+
+
+  // find a triangle containing an UV, starting from a given triangle;
+  // return barycentric coordinates of the UV and the found triangle (indices are zero based).
+  const BRepMesh_Triangle* FindTriangle( const gp_XY&             uv,
+                                         const BRepMesh_Triangle* bmTria,
+                                         double                   bc[3],
+                                         int                      triaNodes[3]);
+
+  // return any Delaunay triangle neighboring a given boundary node (zero based)
+  const BRepMesh_Triangle* GetTriangleNear( int iBndNode );
+
+  // return source boundary nodes
+  const std::vector< const SMDS_MeshNode* >& GetBndNodes() const { return _bndNodes; }
+
+  // return UV of the i-th source boundary node (zero based)
+  gp_XY GetBndUV(const int iNode) const;
+
+  // return scale factor to convert real UV to/from UV used for Delauney meshing:
+  // delauney_UV = real_UV * scale
+  const gp_XY& GetScale() const { return _scale; }
+
+
+ protected:
+
+  // container of a node and a triangle serving as a start while searching a
+  // triangle including the node UV
+  typedef std::list< std::pair< const SMDS_MeshNode*, const BRepMesh_Triangle* > > TNodeTriaList;
+
+  // return UV of a node on the face
+  virtual gp_XY getNodeUV( const TopoDS_Face& face, const SMDS_MeshNode* node ) const = 0;
+
+  // add non-marked nodes surrounding a given one to a queue
+  static void addCloseNodes( const SMDS_MeshNode*     node,
+                             const BRepMesh_Triangle* bmTria,
+                             const int                faceID,
+                             TNodeTriaList &          noTriQueue );
+
+  const TopoDS_Face&                     _face;
+  int                                    _faceID;
+  std::vector< const SMDS_MeshNode* >    _bndNodes;
+  gp_XY                                  _scale;
+  Handle(BRepMesh_DataStructureOfDelaun) _triaDS;
+  size_t                                 _nbNodesToVisit, _nbVisitedNodes, _iBndNode;
+  TNodeTriaList                          _noTriQueue;
+
+};
+
+#endif
index 900aecf..88d2b0b 100644 (file)
@@ -5182,10 +5182,11 @@ omniORB.registerObjref(SMESH._objref_SMESH_Pattern._NP_RepositoryId, Pattern)
 ## Private class used to bind methods creating algorithms to the class Mesh
 #
 class algoCreator:
-    def __init__(self):
+    def __init__(self, method):
         self.mesh = None
         self.defaultAlgoType = ""
         self.algoTypeToClass = {}
+        self.method = method
 
     # Store a python class of algorithm
     def add(self, algoClass):
@@ -5199,25 +5200,48 @@ class algoCreator:
 
     # Create a copy of self and assign mesh to the copy
     def copy(self, mesh):
-        other = algoCreator()
+        other = algoCreator( self.method )
         other.defaultAlgoType = self.defaultAlgoType
-        other.algoTypeToClass  = self.algoTypeToClass
+        other.algoTypeToClass = self.algoTypeToClass
         other.mesh = mesh
         return other
 
     # Create an instance of algorithm
     def __call__(self,algo="",geom=0,*args):
-        algoType = self.defaultAlgoType
-        for arg in args + (algo,geom):
-            if isinstance( arg, geomBuilder.GEOM._objref_GEOM_Object ):
-                geom = arg
-            if isinstance( arg, str ) and arg:
+        algoType = ""
+        shape = 0
+        if isinstance( algo, str ):
+            algoType = algo
+        elif ( isinstance( algo, geomBuilder.GEOM._objref_GEOM_Object ) and \
+               not isinstance( geom, geomBuilder.GEOM._objref_GEOM_Object )):
+            shape = algo
+        elif algo:
+            args += (algo,)
+
+        if isinstance( geom, geomBuilder.GEOM._objref_GEOM_Object ):
+            shape = geom
+        elif not algoType and isinstance( geom, str ):
+            algoType = geom
+        elif geom:
+            args += (geom,)
+        for arg in args:
+            if isinstance( arg, geomBuilder.GEOM._objref_GEOM_Object ) and not shape:
+                shape = arg
+            elif isinstance( arg, str ) and not algoType:
                 algoType = arg
+            else:
+                import traceback, sys
+                msg = "Warning. Unexpected argument in mesh.%s() --->  %s" % ( self.method, arg )
+                sys.stderr.write( msg + '\n' )
+                tb = traceback.extract_stack(None,2)
+                traceback.print_list( [tb[0]] )
+        if not algoType:
+            algoType = self.defaultAlgoType
         if not algoType and self.algoTypeToClass:
             algoType = self.algoTypeToClass.keys()[0]
         if self.algoTypeToClass.has_key( algoType ):
             #print "Create algo",algoType
-            return self.algoTypeToClass[ algoType ]( self.mesh, geom )
+            return self.algoTypeToClass[ algoType ]( self.mesh, shape )
         raise RuntimeError, "No class found for algo type %s" % algoType
         return None
 
@@ -5299,7 +5323,7 @@ for pluginName in os.environ[ "SMESH_MeshersList" ].split( ":" ):
         if type( algo ).__name__ == 'classobj' and hasattr( algo, "meshMethod" ):
             #print "                     meshMethod:" , str(algo.meshMethod)
             if not hasattr( Mesh, algo.meshMethod ):
-                setattr( Mesh, algo.meshMethod, algoCreator()
+                setattr( Mesh, algo.meshMethod, algoCreator( algo.meshMethod ))
                 pass
             getattr( Mesh, algo.meshMethod ).add( algo )
             pass
index 5ae0a4b..0ee4041 100644 (file)
@@ -52,7 +52,6 @@
 #include <GeomLib_IsPlanarSurface.hxx>
 #include <Geom_Curve.hxx>
 #include <Standard_ErrorHandler.hxx>
-#include <TColStd_DataMapOfIntegerInteger.hxx>
 #include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
 #include <TopTools_ListIteratorOfListOfShape.hxx>
@@ -1177,6 +1176,9 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
   {
     // use transformation (issue 0020680, IPAL0052499) or a "straight line" approach
     StdMeshers_Sweeper sweeper;
+    sweeper.myHelper  = myHelper;
+    sweeper.myBotFace = thePrism.myBottom;
+    sweeper.myTopFace = thePrism.myTop;
 
     // load boundary nodes into sweeper
     bool dummy;
@@ -1213,15 +1215,15 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
     {
       double tol = getSweepTolerance( thePrism );
       bool allowHighBndError = !isSimpleBottom( thePrism );
-      myUseBlock = !sweeper.ComputeNodes( *myHelper, tol, allowHighBndError );
+      myUseBlock = !sweeper.ComputeNodesByTrsf( tol, allowHighBndError );
     }
     else if ( sweeper.CheckSameZ() )
     {
-      myUseBlock = !sweeper.ComputeNodesOnStraightSameZ( *myHelper );
+      myUseBlock = !sweeper.ComputeNodesOnStraightSameZ();
     }
     else
     {
-      myUseBlock = !sweeper.ComputeNodesOnStraight( *myHelper, thePrism.myBottom, thePrism.myTop );
+      myUseBlock = !sweeper.ComputeNodesOnStraight();
     }
     myHelper->SetElementsOnShape( false );
   }
@@ -4956,13 +4958,13 @@ void StdMeshers_Sweeper::applyBoundaryError(const vector< gp_XYZ >& bndPoints,
 
 //================================================================================
 /*!
- * \brief Create internal nodes of the prism
+ * \brief Create internal nodes of the prism by computing an affine transformation
+ *        from layer to layer
  */
 //================================================================================
 
-bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper,
-                                       const double        tol,
-                                       const bool          allowHighBndError)
+bool StdMeshers_Sweeper::ComputeNodesByTrsf( const double tol,
+                                             const bool   allowHighBndError)
 {
   const size_t zSize = myBndColumns[0]->size();
   const size_t zSrc = 0, zTgt = zSize-1;
@@ -5229,7 +5231,7 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper,
     for ( size_t z = zSrc + 1; z < zTgt; ++z ) // vertical loop on layers
     {
       const gp_XYZ & xyz = intPntsOfLayer[ z ][ iP ];
-      if ( !( nodeCol[ z ] = helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() )))
+      if ( !( nodeCol[ z ] = myHelper->AddNode( xyz.X(), xyz.Y(), xyz.Z() )))
         return false;
     }
   }
@@ -5301,7 +5303,7 @@ bool StdMeshers_Sweeper::CheckSameZ()
  */
 //================================================================================
 
-bool StdMeshers_Sweeper::ComputeNodesOnStraightSameZ( SMESH_MesherHelper& helper )
+bool StdMeshers_Sweeper::ComputeNodesOnStraightSameZ()
 {
   TZColumn& z = myZColumns[0];
 
@@ -5313,7 +5315,7 @@ bool StdMeshers_Sweeper::ComputeNodesOnStraightSameZ( SMESH_MesherHelper& helper
     for ( size_t iZ = 0; iZ < z.size(); ++iZ )
     {
       gp_XYZ p = n0 * ( 1 - z[iZ] ) + n1 * z[iZ];
-      nodes[ iZ+1 ] = helper.AddNode( p.X(), p.Y(), p.Z() );
+      nodes[ iZ+1 ] = myHelper->AddNode( p.X(), p.Y(), p.Z() );
     }
   }
 
@@ -5327,144 +5329,51 @@ bool StdMeshers_Sweeper::ComputeNodesOnStraightSameZ( SMESH_MesherHelper& helper
  */
 //================================================================================
 
-bool StdMeshers_Sweeper::ComputeNodesOnStraight( SMESH_MesherHelper& helper,
-                                                 const TopoDS_Face&  botFace,
-                                                 const TopoDS_Face&  topFace )
+bool StdMeshers_Sweeper::ComputeNodesOnStraight()
 {
-  // get data to create a Morph
-  UVPtStructVec botUV( myBndColumns.size() + 1 );
-  UVPtStructVec topUV( myBndColumns.size() + 1 );
-  for ( size_t i = 0; i < myBndColumns.size(); ++i )
-  {
-    TNodeColumn& nodes = *myBndColumns[i];
-    botUV[i].node = nodes[0];
-    botUV[i].SetUV( helper.GetNodeUV( botFace, nodes[0] ));
-    topUV[i].node = nodes.back();
-    topUV[i].SetUV( helper.GetNodeUV( topFace, nodes.back() ));
-    botUV[i].node->setIsMarked( true );
-  }
-  botUV.back() = botUV[0];
-  topUV.back() = topUV[0];
-
-  TopoDS_Edge dummyE;
-  TSideVector botWires( 1, StdMeshers_FaceSide::New( botUV, botFace, dummyE, helper.GetMesh() ));
-  TSideVector topWires( 1, StdMeshers_FaceSide::New( topUV, topFace, dummyE, helper.GetMesh() ));
-
-  // use Morph to make delauney mesh on the FACEs. Locating of a node within a
-  // delauney triangle will be used to get a weighted Z.
-  NSProjUtils::Morph botDelauney( botWires );
-  NSProjUtils::Morph topDelauney( topWires );
-
-  if ( helper.GetIsQuadratic() )
-  {
-    // mark all medium nodes of faces on botFace to avoid their treating
-    SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( botFace );
-    SMDS_ElemIteratorPtr eIt = smDS->GetElements();
-    while ( eIt->more() )
-    {
-      const SMDS_MeshElement* e = eIt->next();
-      for ( int i = e->NbCornerNodes(), nb = e->NbNodes(); i < nb; ++i )
-        e->GetNode( i )->setIsMarked( true );
-    }
-  }
+  prepareTopBotDelaunay();
 
-  // map to get a node column by a bottom node
-  TColStd_DataMapOfIntegerInteger iNode2iCol( myIntColumns.size() );
-
-  // un-mark nodes to treat (internal bottom nodes); later we will mark treated nodes
-  for ( size_t i = 0; i < myIntColumns.size(); ++i )
-  {
-    const SMDS_MeshNode* botNode = myIntColumns[i]->front();
-    botNode->setIsMarked( false );
-    iNode2iCol.Bind( botNode->GetID(), i );
-  }
-
-  const int botFaceID = helper.GetMesh()->GetSubMesh( botFace )->GetId();
   const SMDS_MeshNode     *botNode, *topNode;
-  const BRepMesh_Triangle *botTria, *topTria;
+  const BRepMesh_Triangle *topTria;
   double botBC[3], topBC[3]; // barycentric coordinates
   int    botTriaNodes[3], topTriaNodes[3];
   bool   checkUV = true;
 
-  // a queue of bottom nodes with starting delauney triangles
-  NSProjUtils::Morph::TNodeTriaList botNoTriQueue;
+  int nbInternalNodes = myIntColumns.size();
+  myBotDelaunay->InitTraversal( nbInternalNodes );
 
-  size_t iBndN = 1; // index of a bottom boundary node
-  int nbNodesToProcess = myIntColumns.size();
-  while ( nbNodesToProcess > 0 )
+  while (( botNode = myBotDelaunay->NextNode( botBC, botTriaNodes )))
   {
-    while ( !botNoTriQueue.empty() ) // treat all nodes in the queue
-    {
-      botNode = botNoTriQueue.front().first;
-      botTria = botNoTriQueue.front().second;
-      botNoTriQueue.pop_front();
-      if ( botNode->isMarked() )
-        continue;
-      --nbNodesToProcess;
-      botNode->setIsMarked( true );
-
-      TNodeColumn* column = myIntColumns[ iNode2iCol( botNode->GetID() )];
-
-      // find a delauney triangle containing the botNode
-      gp_XY botUV = helper.GetNodeUV( botFace, botNode, NULL, &checkUV );
-      botUV  *= botDelauney.GetScale();
-      botTria = botDelauney.FindTriangle( botUV, botTria, botBC, botTriaNodes );
-      if ( !botTria )
-        return false;
-
-      // find a delauney triangle containing the topNode
-      topNode = column->back();
-      gp_XY topUV = helper.GetNodeUV( topFace, topNode, NULL, &checkUV );
-      topUV *= topDelauney.GetScale();
-      // get a starting triangle basing on that top and bot boundary nodes have same index
-      topTria = topDelauney.GetTriangleNear( botTriaNodes[0] );
-      topTria = topDelauney.FindTriangle( topUV, topTria, topBC, topTriaNodes );
-      if ( !topTria )
-        return false;
-
-      // create nodes along a line
-      SMESH_NodeXYZ botP( botNode ), topP( topNode);
-      for ( size_t iZ = 0; iZ < myZColumns[0].size(); ++iZ )
-      {
-        // use barycentric coordinates as weight of Z of boundary columns
-        double botZ = 0, topZ = 0;
-        for ( int i = 0; i < 3; ++i )
-        {
-          botZ += botBC[i] * myZColumns[ botTriaNodes[i]-1 ][ iZ ];
-          topZ += topBC[i] * myZColumns[ topTriaNodes[i]-1 ][ iZ ];
-        }
-        double rZ = double( iZ + 1 ) / ( myZColumns[0].size() + 1 );
-        double z = botZ * ( 1 - rZ ) + topZ * rZ;
-        gp_XYZ p = botP * ( 1 - z  ) + topP * z;
-        (*column)[ iZ+1 ] = helper.AddNode( p.X(), p.Y(), p.Z() );
-      }
-
-      // add neighbor nodes to the queue
-      botDelauney.AddCloseNodes( botNode, botTria, botFaceID, botNoTriQueue );
-    }
+    TNodeColumn* column = myIntColumns[ myNodeID2ColID( botNode->GetID() )];
+
+    // find a Delaunay triangle containing the topNode
+    topNode = column->back();
+    gp_XY topUV = myHelper->GetNodeUV( myTopFace, topNode, NULL, &checkUV );
+    // get a starting triangle basing on that top and bot boundary nodes have same index
+    topTria = myTopDelaunay->GetTriangleNear( botTriaNodes[0] );
+    topTria = myTopDelaunay->FindTriangle( topUV, topTria, topBC, topTriaNodes );
+    if ( !topTria )
+      return false;
 
-    if ( nbNodesToProcess > 0 ) // fill the queue
+    // create nodes along a line
+    SMESH_NodeXYZ botP( botNode ), topP( topNode);
+    for ( size_t iZ = 0; iZ < myZColumns[0].size(); ++iZ )
     {
-      // assure that all bot nodes are visited
-      for ( ; iBndN-1 < myBndColumns.size() &&  botNoTriQueue.empty();  ++iBndN )
-      {
-        botTria = botDelauney.GetTriangleNear( iBndN );
-        const SMDS_MeshNode*  bndNode = botDelauney.GetBndNodes()[ iBndN ];
-        botDelauney.AddCloseNodes( bndNode, botTria, botFaceID, botNoTriQueue );
-      }
-      if ( botNoTriQueue.empty() )
+      // use barycentric coordinates as weight of Z of boundary columns
+      double botZ = 0, topZ = 0;
+      for ( int i = 0; i < 3; ++i )
       {
-        for ( size_t i = 0; i < myIntColumns.size(); ++i )
-        {
-          botNode = myIntColumns[i]->front();
-          if ( !botNode->isMarked() )
-            botNoTriQueue.push_back( make_pair( botNode, botTria ));
-        }
+        botZ += botBC[i] * myZColumns[ botTriaNodes[i] ][ iZ ];
+        topZ += topBC[i] * myZColumns[ topTriaNodes[i] ][ iZ ];
       }
+      double rZ = double( iZ + 1 ) / ( myZColumns[0].size() + 1 );
+      double z = botZ * ( 1 - rZ ) + topZ * rZ;
+      gp_XYZ p = botP * ( 1 - z  ) + topP * z;
+      (*column)[ iZ+1 ] = myHelper->AddNode( p.X(), p.Y(), p.Z() );
     }
   }
 
-  return true;
+  return myBotDelaunay->NbVisitedNodes() == nbInternalNodes;
 }
 
 //================================================================================
@@ -5490,3 +5399,58 @@ void StdMeshers_Sweeper::fillZColumn( TZColumn&    zColumn,
     zColumn[i] = ( line * vec ) / len2; // param [0,1] on the line
   }
 }
+
+//================================================================================
+/*!
+ * \brief Initialize *Delaunay members
+ */
+//================================================================================
+
+void StdMeshers_Sweeper::prepareTopBotDelaunay()
+{
+  UVPtStructVec botUV( myBndColumns.size() );
+  UVPtStructVec topUV( myBndColumns.size() );
+  for ( size_t i = 0; i < myBndColumns.size(); ++i )
+  {
+    TNodeColumn& nodes = *myBndColumns[i];
+    botUV[i].node = nodes[0];
+    botUV[i].SetUV( myHelper->GetNodeUV( myBotFace, nodes[0] ));
+    topUV[i].node = nodes.back();
+    topUV[i].SetUV( myHelper->GetNodeUV( myTopFace, nodes.back() ));
+    botUV[i].node->setIsMarked( true );
+  }
+  TopoDS_Edge dummyE;
+  SMESH_Mesh* mesh = myHelper->GetMesh();
+  TSideVector botWires( 1, StdMeshers_FaceSide::New( botUV, myBotFace, dummyE, mesh ));
+  TSideVector topWires( 1, StdMeshers_FaceSide::New( topUV, myTopFace, dummyE, mesh ));
+
+  // Delaunay mesh on the FACEs.
+  bool checkUV = false;
+  myBotDelaunay.reset( new NSProjUtils::Delaunay( botWires, checkUV ));
+  myTopDelaunay.reset( new NSProjUtils::Delaunay( topWires, checkUV ));
+
+  if ( myHelper->GetIsQuadratic() )
+  {
+    // mark all medium nodes of faces on botFace to avoid their treating
+    SMESHDS_SubMesh* smDS = myHelper->GetMeshDS()->MeshElements( myBotFace );
+    SMDS_ElemIteratorPtr eIt = smDS->GetElements();
+    while ( eIt->more() )
+    {
+      const SMDS_MeshElement* e = eIt->next();
+      for ( int i = e->NbCornerNodes(), nb = e->NbNodes(); i < nb; ++i )
+        e->GetNode( i )->setIsMarked( true );
+    }
+  }
+
+  // map to get a node column by a bottom node
+  myNodeID2ColID.Clear(/*doReleaseMemory=*/false);
+  myNodeID2ColID.ReSize( myIntColumns.size() );
+
+  // un-mark nodes to treat (internal bottom nodes) to be returned by myBotDelaunay
+  for ( size_t i = 0; i < myIntColumns.size(); ++i )
+  {
+    const SMDS_MeshNode* botNode = myIntColumns[i]->front();
+    botNode->setIsMarked( false );
+    myNodeID2ColID.Bind( botNode->GetID(), i );
+  }
+}
index ccdf545..51a9c4e 100644 (file)
 
 #include "SMESHDS_Mesh.hxx"
 #include "SMESH_Block.hxx"
+#include "SMESH_Comment.hxx"
 #include "SMESH_Mesh.hxx"
 #include "SMESH_MesherHelper.hxx"
 #include "SMESH_TypeDefs.hxx"
-#include "SMESH_Comment.hxx"
 #include "SMESH_subMesh.hxx"
+#include "StdMeshers_ProjectionUtils.hxx"
 
 #include <Adaptor2d_Curve2d.hxx>
 #include <Adaptor3d_Curve.hxx>
 #include <Adaptor3d_Surface.hxx>
 #include <BRepAdaptor_Surface.hxx>
+#include <TColStd_DataMapOfIntegerInteger.hxx>
 #include <TopTools_IndexedMapOfOrientedShape.hxx>
 #include <TopoDS_Face.hxx>
 #include <gp_Trsf.hxx>
@@ -53,10 +55,6 @@ namespace Prism_3D
   struct TNode;
   struct TPrismTopo;
 }
-namespace StdMeshers_ProjectionUtils
-{
-  class TrsfFinder3D;
-}
 class SMESHDS_SubMesh;
 class TopoDS_Edge;
 
@@ -417,23 +415,29 @@ private:
  */
 struct StdMeshers_Sweeper
 {
-  std::vector< TNodeColumn* > myBndColumns; // boundary nodes
-  std::vector< TNodeColumn* > myIntColumns; // internal nodes
+  SMESH_MesherHelper*                     myHelper;
+  TopoDS_Face                             myBotFace;
+  TopoDS_Face                             myTopFace;
+
+  std::vector< TNodeColumn* >             myBndColumns; // boundary nodes
+  std::vector< TNodeColumn* >             myIntColumns; // internal nodes
 
   typedef std::vector< double > TZColumn;
-  std::vector< TZColumn > myZColumns; // Z distribution of boundary nodes
+  std::vector< TZColumn >                 myZColumns; // Z distribution of boundary nodes
+
+  StdMeshers_ProjectionUtils::DelaunayPtr myTopDelaunay;
+  StdMeshers_ProjectionUtils::DelaunayPtr myBotDelaunay;
+  TColStd_DataMapOfIntegerInteger         myNodeID2ColID;
 
-  bool ComputeNodes( SMESH_MesherHelper& helper,
-                     const double        tol,
-                     const bool          allowHighBndError );
+
+  bool ComputeNodesByTrsf( const double tol,
+                           const bool   allowHighBndError );
 
   bool CheckSameZ();
 
-  bool ComputeNodesOnStraightSameZ( SMESH_MesherHelper& helper );
+  bool ComputeNodesOnStraightSameZ();
 
-  bool ComputeNodesOnStraight( SMESH_MesherHelper& helper,
-                               const TopoDS_Face&  bottom,
-                               const TopoDS_Face&  top);
+  bool ComputeNodesOnStraight();
 
 private:
 
@@ -459,6 +463,8 @@ private:
 
   static void fillZColumn( TZColumn&    zColumn,
                            TNodeColumn& nodes );
+
+  void prepareTopBotDelaunay();
 };
 
 // ===============================================
index cda4346..11f2b5e 100644 (file)
@@ -478,6 +478,26 @@ namespace {
     return !allBndEdges.empty();
   }
 
+  /*!
+   * \brief Convertor used in Delaunay constructor
+   */
+  struct SideVector2UVPtStructVec
+  {
+    std::vector< const UVPtStructVec* > _uvVecs;
+
+    SideVector2UVPtStructVec( const TSideVector& wires )
+    {
+      _uvVecs.resize( wires.size() );
+      for ( size_t i = 0; i < wires.size(); ++i )
+        _uvVecs[ i ] = & wires[i]->GetUVPtStruct();
+    }
+
+    operator const std::vector< const UVPtStructVec* > & () const
+    {
+      return _uvVecs;
+    }
+  };
+
 } // namespace
 
 //=======================================================================
@@ -2797,183 +2817,15 @@ namespace StdMeshers_ProjectionUtils
 
   //================================================================================
   /*!
-   * \brief Add in-FACE nodes surrounding a given node to a queue
-   */
-  //================================================================================
-
-  void Morph::AddCloseNodes( const SMDS_MeshNode*     srcNode,
-                             const BRepMesh_Triangle* bmTria,
-                             const int                srcFaceID,
-                             TNodeTriaList &          noTriQueue )
-  {
-    // find in-FACE nodes
-    SMDS_ElemIteratorPtr elems = srcNode->GetInverseElementIterator(SMDSAbs_Face);
-    while ( elems->more() )
-    {
-      const SMDS_MeshElement* elem = elems->next();
-      if ( elem->getshapeId() == srcFaceID )
-      {
-        for ( int i = 0, nb = elem->NbNodes(); i < nb; ++i )
-        {
-          const SMDS_MeshNode* n = elem->GetNode( i );
-          if ( !n->isMarked() /*&& n->getshapeId() == srcFaceID*/ )
-            noTriQueue.push_back( make_pair( n, bmTria ));
-        }
-      }
-    }
-  }
-
-  //================================================================================
-  /*!
-   * \brief Find a delauney triangle containing a given 2D point and return
-   *        barycentric coordinates within the found triangle
-   */
-  //================================================================================
-
-  const BRepMesh_Triangle* Morph::FindTriangle( const gp_XY&             uv,
-                                                const BRepMesh_Triangle* bmTria,
-                                                double                   bc[3],
-                                                int                      triaNodes[3] )
-  {
-    int   nodeIDs[3];
-    gp_XY nodeUVs[3];
-    int   linkIDs[3];
-    Standard_Boolean ori[3];
-
-    while ( bmTria )
-    {
-      // check bmTria
-
-      _triaDS->ElementNodes( *bmTria, nodeIDs );
-      nodeUVs[0] = _triaDS->GetNode( nodeIDs[0] ).Coord();
-      nodeUVs[1] = _triaDS->GetNode( nodeIDs[1] ).Coord();
-      nodeUVs[2] = _triaDS->GetNode( nodeIDs[2] ).Coord();
-
-      SMESH_MeshAlgos::GetBarycentricCoords( uv,
-                                             nodeUVs[0], nodeUVs[1], nodeUVs[2],
-                                             bc[0], bc[1] );
-      if ( bc[0] >= 0 && bc[1] >= 0 && bc[0] + bc[1] <= 1 )
-      {
-        bc[2] = 1 - bc[0] - bc[1];
-        triaNodes[0] = nodeIDs[0];
-        triaNodes[1] = nodeIDs[1];
-        triaNodes[2] = nodeIDs[2];
-        return bmTria;
-      }
-
-      // look for a neighbor triangle, which is adjacent to a link intersected
-      // by a segment( triangle center -> uv )
-
-      gp_XY gc = ( nodeUVs[0] + nodeUVs[1] + nodeUVs[2] ) / 3.;
-      gp_XY seg = uv - gc;
-
-      bmTria->Edges( linkIDs, ori );
-      int triaID = _triaDS->IndexOf( *bmTria );
-      bmTria = 0;
-
-      for ( int i = 0; i < 3; ++i )
-      {
-        const BRepMesh_PairOfIndex & triIDs = _triaDS->ElementsConnectedTo( linkIDs[i] );
-        if ( triIDs.Extent() < 2 )
-          continue; // no neighbor triangle
-
-        // check if a link intersects gc2uv
-        const BRepMesh_Edge & link = _triaDS->GetLink( linkIDs[i] );
-        const BRepMesh_Vertex & n1 = _triaDS->GetNode( link.FirstNode() );
-        const BRepMesh_Vertex & n2 = _triaDS->GetNode( link.LastNode() );
-        gp_XY uv1 = n1.Coord();
-        gp_XY lin = n2.Coord() - uv1; // link direction
-
-        double crossSegLin = seg ^ lin;
-        if ( Abs( crossSegLin ) < std::numeric_limits<double>::min() )
-          continue; // parallel
-
-        double uSeg = ( uv1 - gc ) ^ lin / crossSegLin;
-        if ( 0. <= uSeg && uSeg <= 1. )
-        {
-          bmTria = & _triaDS->GetElement( triIDs.Index( 1 + ( triIDs.Index(1) == triaID )));
-          break;
-        }
-      }
-    }
-    return bmTria;
-  }
-
-  //================================================================================
-  /*!
-   * \brief Return a triangle sharing a given boundary node
-   *  \param [in] iBndNode - index of the boundary node
-   *  \return const BRepMesh_Triangle* - a found triangle
-   */
-  //================================================================================
-
-  const BRepMesh_Triangle* Morph::GetTriangleNear( int iBndNode )
-  {
-    const BRepMesh::ListOfInteger & linkIds = _triaDS->LinksConnectedTo( iBndNode );
-    const BRepMesh_PairOfIndex &    triaIds = _triaDS->ElementsConnectedTo( linkIds.First() );
-    const BRepMesh_Triangle&           tria = _triaDS->GetElement( triaIds.Index(1) );
-    return &tria;
-  }
-
-  //================================================================================
-  /*!
    * \brief triangulate the srcFace in 2D
    *  \param [in] srcWires - boundary of the src FACE
    */
   //================================================================================
 
-  Morph::Morph(const TSideVector& srcWires)
+  Morph::Morph(const TSideVector& srcWires):
+    _delaunay( srcWires, /*checkUV=*/true )
   {
     _srcSubMesh = srcWires[0]->GetMesh()->GetSubMesh( srcWires[0]->Face() );
-
-    // compute _scale
-    {
-      BRepAdaptor_Surface surf( srcWires[0]->Face() );
-      const int nbDiv = 100;
-      const double uRange = surf.LastUParameter() - surf.FirstUParameter();
-      const double vRange = surf.LastVParameter() - surf.FirstVParameter();
-      const double dU = uRange / nbDiv;
-      const double dV = vRange / nbDiv;
-      double u = surf.FirstUParameter(), v = surf.FirstVParameter();
-      gp_Pnt p0U = surf.Value( u, v ), p0V = p0U;
-      double lenU = 0, lenV = 0;
-      for ( ; u < surf.LastUParameter(); u += dU, v += dV )
-      {
-        gp_Pnt p1U = surf.Value( u, surf.FirstVParameter() );
-        lenU += p1U.Distance( p0U );
-        p0U = p1U;
-        gp_Pnt p1V = surf.Value( surf.FirstUParameter(), v );
-        lenV += p1V.Distance( p0V );
-        p0V = p1V;
-      }
-      _scale.SetCoord( lenU / uRange, lenV / vRange );
-    }
-
-    // count boundary points
-    int iP = 1, nbP = 0;
-    for ( size_t iW = 0; iW < srcWires.size(); ++iW )
-      nbP += srcWires[iW]->NbPoints() - 1; // 1st and last points coincide
-
-    _bndSrcNodes.resize( nbP + 1 ); _bndSrcNodes[0] = 0;
-
-    // fill boundary points
-    BRepMesh::Array1OfVertexOfDelaun srcVert( 1, 1 + nbP );
-    BRepMesh_Vertex v( 0, 0, BRepMesh_Frontier );
-    for ( size_t iW = 0; iW < srcWires.size(); ++iW )
-    {
-      const UVPtStructVec& srcPnt = srcWires[iW]->GetUVPtStruct();
-      for ( int i = 0, nb = srcPnt.size() - 1;  i < nb;  ++i, ++iP )
-      {
-        _bndSrcNodes[ iP ]  = srcPnt[i].node;
-        srcPnt[i].node->setIsMarked( true );
-
-        v.ChangeCoord() = srcPnt[i].UV().Multiplied( _scale );
-        srcVert( iP )   = v;
-      }
-    }
-    // triangulate the srcFace in 2D
-    BRepMesh_Delaun delauney( srcVert );
-    _triaDS = delauney.Result();
   }
 
   //================================================================================
@@ -2994,32 +2846,27 @@ namespace StdMeshers_ProjectionUtils
                       const TNodeNodeMap&           src2tgtNodes,
                       const bool                    moveAll)
   {
-    // get tgt boundary points corresponding to _bndSrcNodes
+    // get tgt boundary points corresponding to src boundary nodes
     size_t nbP = 0;
     for ( size_t iW = 0; iW < tgtWires.size(); ++iW )
       nbP += tgtWires[iW]->NbPoints() - 1; // 1st and last points coincide
-    if ( nbP != _bndSrcNodes.size() - 1 )
+    if ( nbP != _delaunay.GetBndNodes().size() )
       return false;
 
-    BRepMesh::Array1OfVertexOfDelaun tgtVert( 1, 1 + nbP );
-    BRepMesh_Vertex v( 0, 0, BRepMesh_Frontier );
-    for ( size_t iW = 0, iP = 1; iW < tgtWires.size(); ++iW )
+    std::vector< gp_XY > tgtUV( nbP );
+    for ( size_t iW = 0, iP = 0; iW < tgtWires.size(); ++iW )
     {
       const UVPtStructVec& tgtPnt = tgtWires[iW]->GetUVPtStruct();
       for ( int i = 0, nb = tgtPnt.size() - 1;  i < nb;  ++i, ++iP )
       {
-        v.ChangeCoord() = tgtPnt[i].UV().Multiplied( _scale );
-        tgtVert( iP )   = v;
+        tgtUV[ iP ] = tgtPnt[i].UV();
       }
     }
 
-    const TopoDS_Face& srcFace = TopoDS::Face( _srcSubMesh->GetSubShape() );
-    const int        srcFaceID = _srcSubMesh->GetId();
     SMESHDS_Mesh* tgtMesh = tgtHelper.GetMeshDS();
     const SMDS_MeshNode *srcNode, *tgtNode;
-    const BRepMesh_Triangle *bmTria;
 
-    // un-mark internal src nodes; later we will mark moved nodes
+    // un-mark internal src nodes in order iterate them using _delaunay
     int nbSrcNodes = 0;
     SMDS_NodeIteratorPtr nIt = _srcSubMesh->GetSubMeshDS()->GetNodes();
     if ( !nIt || !nIt->more() ) return true;
@@ -3038,82 +2885,75 @@ namespace StdMeshers_ProjectionUtils
     // Move tgt nodes
 
     double bc[3]; // barycentric coordinates
-    int    nodeIDs[3];
-    bool   checkUV = true;
+    int    nodeIDs[3]; // nodes of a delaunay triangle
     const SMDS_FacePosition* pos;
 
-    // a queue of nodes with starting triangles
-    TNodeTriaList noTriQueue;
-    size_t iBndSrcN = 1;
+    _delaunay.InitTraversal( nbSrcNodes );
 
-    while ( nbSrcNodes > 0 )
+    while (( srcNode = _delaunay.NextNode( bc, nodeIDs )))
     {
-      while ( !noTriQueue.empty() )
-      {
-        srcNode = noTriQueue.front().first;
-        bmTria  = noTriQueue.front().second;
-        noTriQueue.pop_front();
-        if ( srcNode->isMarked() )
-          continue;
-        --nbSrcNodes;
-        srcNode->setIsMarked( true );
-
-        // find a delauney triangle containing the src node
-        gp_XY uv = tgtHelper.GetNodeUV( srcFace, srcNode, NULL, &checkUV );
-        uv *= _scale;
-        bmTria = FindTriangle( uv, bmTria, bc, nodeIDs );
-        if ( !bmTria )
-          continue;
-
-        // compute new coordinates for a corresponding tgt node
-        gp_XY uvNew( 0., 0. ), nodeUV;
-        for ( int i = 0; i < 3; ++i )
-          uvNew += bc[i] * tgtVert( nodeIDs[i]).Coord();
-        uvNew.SetCoord( uvNew.X() / _scale.X(),
-                        uvNew.Y() / _scale.Y() );
-        gp_Pnt xyz = tgtSurface->Value( uvNew );
-
-        // find and move tgt node
-        TNodeNodeMap::const_iterator n2n = src2tgtNodes.find( srcNode );
-        if ( n2n == src2tgtNodes.end() ) continue;
-        tgtNode = n2n->second;
-        tgtMesh->MoveNode( tgtNode, xyz.X(), xyz.Y(), xyz.Z() );
-
-        if (( pos = dynamic_cast< const SMDS_FacePosition* >( tgtNode->GetPosition() )))
-          const_cast<SMDS_FacePosition*>( pos )->SetParameters( uvNew.X(), uvNew.Y() );
-
-        AddCloseNodes( srcNode, bmTria, srcFaceID, noTriQueue );
-      }
+      // compute new coordinates for a corresponding tgt node
+      gp_XY uvNew( 0., 0. ), nodeUV;
+      for ( int i = 0; i < 3; ++i )
+        uvNew += bc[i] * tgtUV[ nodeIDs[i]];
+      gp_Pnt xyz = tgtSurface->Value( uvNew );
 
-      if ( nbSrcNodes > 0 )
-      {
-        // assure that all src nodes are visited
-        for ( ; iBndSrcN < _bndSrcNodes.size() &&  noTriQueue.empty();  ++iBndSrcN )
-        {
-          const BRepMesh_Triangle* tria = GetTriangleNear( iBndSrcN );
-          AddCloseNodes( _bndSrcNodes[ iBndSrcN ], tria, srcFaceID, noTriQueue );
-        }
-        if ( noTriQueue.empty() )
-        {
-          SMDS_NodeIteratorPtr nIt = _srcSubMesh->GetSubMeshDS()->GetNodes();
-          while ( nIt->more() )
-          {
-            srcNode = nIt->next();
-            if ( !srcNode->isMarked() )
-              noTriQueue.push_back( make_pair( srcNode, bmTria ));
-          }
-        }
-      }
+      // find and move tgt node
+      TNodeNodeMap::const_iterator n2n = src2tgtNodes.find( srcNode );
+      if ( n2n == src2tgtNodes.end() ) continue;
+      tgtNode = n2n->second;
+      tgtMesh->MoveNode( tgtNode, xyz.X(), xyz.Y(), xyz.Z() );
+
+      if (( pos = dynamic_cast< const SMDS_FacePosition* >( tgtNode->GetPosition() )))
+        const_cast<SMDS_FacePosition*>( pos )->SetParameters( uvNew.X(), uvNew.Y() );
+
+      --nbSrcNodes;
     }
 
-    return true;
+    return nbSrcNodes == 0;
 
   } // Morph::Perform
 
-  gp_XY Morph::GetBndUV(const int iNode) const
+  //=======================================================================
+  //function : Delaunay
+  //purpose  : construct from face sides
+  //=======================================================================
+
+  Delaunay::Delaunay( const TSideVector& wires, bool checkUV ):
+    SMESH_Delaunay( SideVector2UVPtStructVec( wires ),
+                    TopoDS::Face( wires[0]->FaceHelper()->GetSubShape() ),
+                    wires[0]->FaceHelper()->GetSubShapeID() )
   {
-    return _triaDS->GetNode( iNode ).Coord();
+    _wire = wires[0]; // keep a wire to assure _helper to keep alive
+    _helper = _wire->FaceHelper();
+    _checkUVPtr = checkUV ? & _checkUV : 0;
   }
 
+  //=======================================================================
+  //function : Delaunay
+  //purpose  : construct from UVPtStructVec's
+  //=======================================================================
+
+  Delaunay::Delaunay( const std::vector< const UVPtStructVec* > & boundaryNodes,
+                      SMESH_MesherHelper&                         faceHelper,
+                      bool                                        checkUV):
+    SMESH_Delaunay( boundaryNodes,
+                    TopoDS::Face( faceHelper.GetSubShape() ),
+                    faceHelper.GetSubShapeID() )
+  {
+    _helper = & faceHelper;
+    _checkUVPtr = checkUV ? & _checkUV : 0;
+  }
+
+  //=======================================================================
+  //function : getNodeUV
+  //purpose  : 
+  //=======================================================================
+
+  gp_XY Delaunay::getNodeUV( const TopoDS_Face& face, const SMDS_MeshNode* node ) const
+  {
+    return _helper->GetNodeUV( face, node, 0, _checkUVPtr );
+  }
+  
 
 } // namespace StdMeshers_ProjectionUtils
index 6c94e45..aa38284 100644 (file)
 
 #include "SMESH_StdMeshers.hxx"
 
-#include "StdMeshers_FaceSide.hxx"
 #include "SMDS_MeshElement.hxx"
+#include "SMESH_Delaunay.hxx"
+#include "StdMeshers_FaceSide.hxx"
 
-#include <BRepMesh_DataStructureOfDelaun.hxx>
 #include <ShapeAnalysis_Surface.hxx>
 #include <TopTools_DataMapOfShapeShape.hxx>
 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
@@ -54,6 +54,7 @@ class SMESH_Mesh;
 class SMESH_subMesh;
 class TopoDS_Shape;
 
+//-----------------------------------------------------------------------------------------
 /*!
  * \brief Struct used instead of a sole TopTools_DataMapOfShapeShape to avoid
  *        problems with bidirectional bindings
@@ -94,6 +95,7 @@ namespace StdMeshers_ProjectionUtils
                    TIDCompare>                                 TNodeNodeMap;
 
 
+  //-----------------------------------------------------------------------------------------
   /*!
    * \brief Finds transformation between two sets of 2D points using
    *        a least square approximation
@@ -114,6 +116,7 @@ namespace StdMeshers_ProjectionUtils
 
     bool IsIdentity() const { return ( _trsf.Form() == gp_Identity ); }
   };
+  //-----------------------------------------------------------------------------------------
   /*!
    * \brief Finds transformation between two sets of 3D points using
    *        a least square approximation
@@ -139,16 +142,44 @@ namespace StdMeshers_ProjectionUtils
     bool Invert();
   };
 
+  //-----------------------------------------------------------------------------------------
+  /*!
+   * \brief Create a Delaunay triangulation of nodes on a face boundary
+   *        and provide exploration of nodes shared by elements lying on
+   *        the face. For a returned node, also return a Delaunay triangle
+   *        the node lies in and its Barycentric Coordinates within the triangle.
+   *        Only non-marked nodes are visited.
+   *
+   * The main methods are defined in ../SMESHUtils/SMESH_Delaunay.hxx
+   */
+  class Delaunay : public SMESH_Delaunay
+  {
+  public:
+
+    Delaunay( const TSideVector& wires, bool checkUV = false );
+
+    Delaunay( const std::vector< const UVPtStructVec* > & boundaryNodes,
+              SMESH_MesherHelper&                         faceHelper,
+              bool                                        checkUV = false);
+
+  protected:
+    virtual gp_XY getNodeUV( const TopoDS_Face& face, const SMDS_MeshNode* node ) const;
+
+  private:
+    SMESH_MesherHelper*    _helper;
+    StdMeshers_FaceSidePtr _wire;
+    bool                  *_checkUVPtr, _checkUV;
+  };
+  typedef boost::shared_ptr< Delaunay > DelaunayPtr;
+
+  //-----------------------------------------------------------------------------------------
   /*!
    * \brief Morph mesh on the target FACE to lie within FACE boundary w/o distortion
    */
   class Morph
   {
-    std::vector< const SMDS_MeshNode* >    _bndSrcNodes;
-    Handle(BRepMesh_DataStructureOfDelaun) _triaDS;
-    SMESH_subMesh*                         _srcSubMesh;
-    bool                                   _moveAll;
-    gp_XY                                  _scale;
+    Delaunay       _delaunay;
+    SMESH_subMesh* _srcSubMesh;
   public:
 
     Morph(const TSideVector& srcWires);
@@ -158,37 +189,9 @@ namespace StdMeshers_ProjectionUtils
                  Handle(ShapeAnalysis_Surface) tgtSurface,
                  const TNodeNodeMap&           src2tgtNodes,
                  const bool                    moveAll);
-
-    // find a triangle containing an UV starting from a given triangle;
-    // return barycentric coordinates of the UV in the found triangle
-    const BRepMesh_Triangle* FindTriangle( const gp_XY&             uv,
-                                           const BRepMesh_Triangle* bmTria,
-                                           double                   bc[3],
-                                           int                      triaNodes[3]);
-
-    // return any delauney triangle neighboring a given boundary node
-    const BRepMesh_Triangle* GetTriangleNear( int iBndNode );
-
-    // return source boundary nodes. 0-th node is NULL so that indices of
-    // boundary nodes correspond to indices used by Delauney mesh
-    const std::vector< const SMDS_MeshNode* >& GetBndNodes() const { return _bndSrcNodes; }
-
-    // return UV of the i-th source boundary node
-    gp_XY GetBndUV(const int iNode) const;
-
-    // return scale factor to convert real UV to/from UV used for Delauney meshing:
-    // delauney_UV = real_UV * scale
-    const gp_XY& GetScale() const { return _scale; }
-
-    typedef std::list< std::pair< const SMDS_MeshNode*, const BRepMesh_Triangle* > > TNodeTriaList;
-
-    // add non-marked nodes surrounding a given one to a list
-    static void AddCloseNodes( const SMDS_MeshNode*     node,
-                               const BRepMesh_Triangle* bmTria,
-                               const int                faceID,
-                               TNodeTriaList &          noTriQueue );
   };
 
+  //-----------------------------------------------------------------------------------------
   /*!
    * \brief Looks for association of all sub-shapes of two shapes
    * \param theShape1 - shape 1