Salome HOME
23422: EDF 14312 - Strange behavior of Viscous Layer Surface offset + smooth
[modules/smesh.git] / src / StdMeshers / StdMeshers_ProjectionUtils.cxx
index a0164ef01e073439895d4a43bc8d8c62b6ab092f..e7051db9b8efa3751ed70d85c13137066fb2c113 100644 (file)
@@ -28,6 +28,7 @@
 #include "StdMeshers_ProjectionUtils.hxx"
 
 #include "SMDS_EdgePosition.hxx"
+#include "SMDS_FacePosition.hxx"
 #include "SMESHDS_Mesh.hxx"
 #include "SMESH_Algo.hxx"
 #include "SMESH_Block.hxx"
@@ -46,6 +47,7 @@
 #include "utilities.h"
 
 #include <BRepAdaptor_Surface.hxx>
+#include <BRepMesh_Delaun.hxx>
 #include <BRepTools.hxx>
 #include <BRepTools_WireExplorer.hxx>
 #include <BRep_Builder.hxx>
@@ -2787,4 +2789,320 @@ namespace StdMeshers_ProjectionUtils
     }
     return true;
   }
-}
+
+  //================================================================================
+  /*!
+   * \brief Add in-FACE nodes surrounding a given node to a queue
+   */
+  //================================================================================
+
+  typedef list< pair< const SMDS_MeshNode*, const BRepMesh_Triangle* > > TNodeTriaList;
+
+  void 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() )
+            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* findTriangle( const gp_XY&                            uv,
+                                         const BRepMesh_Triangle*                bmTria,
+                                         Handle(BRepMesh_DataStructureOfDelaun)& triaDS,
+                                         double                                  bc[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];
+        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 triangulate the srcFace in 2D
+   *  \param [in] srcWires - boundary of the src FACE
+   */
+  //================================================================================
+
+  Morph::Morph(const TSideVector& srcWires)
+  {
+    _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();
+  }
+
+  //================================================================================
+  /*!
+   * \brief Move non-marked target nodes
+   *  \param [in,out] tgtHelper - helper
+   *  \param [in] tgtWires - boundary nodes of the target FACE; must be in the
+   *         same order as the nodes in srcWires given in the constructor
+   *  \param [in] src2tgtNodes - map of src -> tgt nodes
+   *  \param [in] moveAll - to move all nodes; if \c false, move only non-marked nodes
+   *  \return bool - Ok or not
+   */
+  //================================================================================
+
+  bool Morph::Perform(SMESH_MesherHelper&           tgtHelper,
+                      const TSideVector&            tgtWires,
+                      Handle(ShapeAnalysis_Surface) tgtSurface,
+                      const TNodeNodeMap&           src2tgtNodes,
+                      const bool                    moveAll)
+  {
+    // get tgt boundary points corresponding to _bndSrcNodes
+    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 )
+      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 )
+    {
+      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;
+      }
+    }
+
+    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;
+
+    // initialize a queue of nodes with starting triangles
+    TNodeTriaList noTriQueue;
+    size_t iBndSrcN = 1;
+    for ( ; iBndSrcN < _bndSrcNodes.size() &&  noTriQueue.empty();  ++iBndSrcN )
+    {
+      // get a triangle
+      const BRepMesh::ListOfInteger & linkIds = _triaDS->LinksConnectedTo( iBndSrcN );
+      const BRepMesh_PairOfIndex &    triaIds = _triaDS->ElementsConnectedTo( linkIds.First() );
+      const BRepMesh_Triangle&           tria = _triaDS->GetElement( triaIds.Index(1) );
+
+      addCloseNodes( _bndSrcNodes[ iBndSrcN ], &tria, srcFaceID, noTriQueue );
+    }
+
+    // un-mark internal src nodes; later we will mark moved nodes
+    int nbSrcNodes = 0;
+    SMDS_NodeIteratorPtr nIt = _srcSubMesh->GetSubMeshDS()->GetNodes();
+    if ( !nIt || !nIt->more() ) return true;
+    if ( moveAll )
+    {
+      nbSrcNodes = _srcSubMesh->GetSubMeshDS()->NbNodes();
+      while ( nIt->more() )
+        nIt->next()->setIsMarked( false );
+    }
+    else
+    {
+      while ( nIt->more() )
+        nbSrcNodes += int( !nIt->next()->isMarked() );
+    }
+
+    // Move tgt nodes
+
+    double bc[3]; // barycentric coordinates
+    int    nodeIDs[3];
+    bool   checkUV = true;
+    const SMDS_FacePosition* pos;
+
+    while ( nbSrcNodes > 0 )
+    {
+      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, _triaDS, bc );
+        if ( !bmTria )
+          continue;
+
+        // compute new coordinates for a corresponding tgt node
+        gp_XY uvNew( 0., 0. ), nodeUV;
+        _triaDS->ElementNodes( *bmTria, nodeIDs );
+        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 );
+      }
+
+      if ( nbSrcNodes > 0 )
+      {
+        // assure that all src nodes are visited
+        for ( ; iBndSrcN < _bndSrcNodes.size() &&  noTriQueue.empty();  ++iBndSrcN )
+        {
+          const BRepMesh::ListOfInteger & linkIds = _triaDS->LinksConnectedTo( iBndSrcN );
+          const BRepMesh_PairOfIndex &    triaIds = _triaDS->ElementsConnectedTo( linkIds.First() );
+          const BRepMesh_Triangle&           tria = _triaDS->GetElement( triaIds.Index(1) );
+          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 ));
+          }
+        }
+      }
+    }
+
+    return true;
+
+  } // Morph::Perform
+
+gp_XY Morph::GetBndUV(const int iNode) const
+{
+    return _triaDS->GetNode( iNode ).Coord();
+  }
+
+
+} // namespace StdMeshers_ProjectionUtils