From: eap Date: Tue, 13 Feb 2018 12:16:32 +0000 (+0300) Subject: GPUSPHGUI: Offset transformation X-Git-Tag: V8_5_0a1~8 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=ead27b045b448d98442ad5817fd47947178cd7e8;p=modules%2Fsmesh.git GPUSPHGUI: Offset transformation add missing files --- diff --git a/src/SMESHGUI/SMESHGUI_OffsetDlg.h b/src/SMESHGUI/SMESHGUI_OffsetDlg.h new file mode 100644 index 000000000..e8f4b9c60 --- /dev/null +++ b/src/SMESHGUI/SMESHGUI_OffsetDlg.h @@ -0,0 +1,134 @@ +// Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE +// +// 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 +// + +// SMESH SMESHGUI : GUI for SMESH component +// File : SMESHGUI_OffsetDlg.h +// +#ifndef SMESHGUI_OffsetDLG_H +#define SMESHGUI_OffsetDLG_H + +// SMESH includes +#include "SMESH_SMESHGUI.hxx" +#include "SMESHGUI_PreviewDlg.h" + +// IDL includes +#include +#include CORBA_SERVER_HEADER(SMESH_Mesh) + +class QButtonGroup; +class QGroupBox; +class QLabel; +class QLineEdit; +class QPushButton; +class QRadioButton; +class QCheckBox; +class SMESHGUI; +class SMESHGUI_IdValidator; +class SMESHGUI_SpinBox; +class SMESHGUI_FilterDlg; +class SMESH_Actor; +class SVTK_Selector; +class LightApp_SelectionMgr; +class SMESH_LogicalFilter; + +//================================================================================= +// class : SMESHGUI_OffsetDlg +// purpose : +//================================================================================= +class SMESHGUI_EXPORT SMESHGUI_OffsetDlg : public SMESHGUI_MultiPreviewDlg +{ + Q_OBJECT + + public: + SMESHGUI_OffsetDlg( SMESHGUI* ); + ~SMESHGUI_OffsetDlg(); + + private: + void Init( bool = true ); + void enterEvent( QEvent* ); /* mouse enter the QWidget */ + void keyPressEvent( QKeyEvent* ); + int GetConstructorId(); + void setNewMeshName(); + + bool isValid(); + void getOffset( SMESH::PointStruct& thePoint, + SMESH::double_array_var& theOffsetFact); + + SMESHGUI_IdValidator* myIdValidator; + LightApp_SelectionMgr* mySelectionMgr; /* User shape selection */ + QString myElementsId; + int myNbOkElements; /* to check when elements are defined */ + + SVTK_Selector* mySelector; + + QWidget* myEditCurrentArgument; + + bool myBusy; + SMESH_Actor* myActor; + SMESH_LogicalFilter* myMeshOrSubMeshOrGroupFilter; + + QList myObjects; + QList myObjectsNames; + QList myMeshes; + + + QGroupBox* ConstructorsBox; + QGroupBox* GroupButtons; + QPushButton* buttonOk; + QPushButton* buttonCancel; + QPushButton* buttonApply; + QPushButton* buttonHelp; + QGroupBox* GroupArguments; + QLabel* TextLabelElements; + //QPushButton* SelectElementsButton; + QLineEdit* LineEditElements; + QCheckBox* CheckBoxMesh; + SMESHGUI_SpinBox* SpinBox; + QGroupBox* ActionBox; + QButtonGroup* ActionGroup; + QCheckBox* MakeGroupsCheck; + QLineEdit* LineEditNewMesh; + + QString myHelpFileName; + + QPushButton* myFilterBtn; + SMESHGUI_FilterDlg* myFilterDlg; + + +protected slots: + virtual void onDisplaySimulation( bool ); + virtual void reject(); + void onFilterAccepted(); + +private slots: + void ClickOnOk(); + bool ClickOnApply(); + void ClickOnHelp(); + void SelectionIntoArgument(); + void DeactivateActiveDialog(); + void ActivateThisDialog(); + void onTextChange( const QString& ); + void onSelectMesh( bool ); + void onActionClicked( int ); + void onOpenView(); + void onCloseView(); + void setFilters(); +}; + +#endif // SMESHGUI_OffsetDLG_H diff --git a/src/SMESHUtils/SMESH_MeshAlgos.cxx b/src/SMESHUtils/SMESH_MeshAlgos.cxx index 1de6ba82a..fa65474f8 100644 --- a/src/SMESHUtils/SMESH_MeshAlgos.cxx +++ b/src/SMESHUtils/SMESH_MeshAlgos.cxx @@ -258,7 +258,7 @@ namespace // Utils used in SMESH_ElementSearcherImpl::FindElementsByPoint() //!< allocator of ElementBox's and SMESH_TreeLimit struct LimitAndPool : public SMESH_TreeLimit { - TElementBoxPool _elBoPool; + TElementBoxPool _elBoPool; LimitAndPool():SMESH_TreeLimit( MaxLevel, /*minSize=*/0. ) {} }; LimitAndPool* getLimitAndPool() const diff --git a/src/SMESHUtils/SMESH_Offset.cxx b/src/SMESHUtils/SMESH_Offset.cxx new file mode 100644 index 000000000..f571535f2 --- /dev/null +++ b/src/SMESHUtils/SMESH_Offset.cxx @@ -0,0 +1,2912 @@ +// 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_Offset.cxx +// Created : Mon Dec 25 15:52:38 2017 +// Author : Edward AGAPOV (eap) + +#include "SMESH_MeshAlgos.hxx" + +#include +#include "SMDS_Mesh.hxx" + +#include + +#include +#include +#include +#include + +#include +#include + +namespace +{ + const size_t theMaxNbFaces = 256; // max number of faces sharing a node + + typedef NCollection_DataMap< Standard_Address, const SMDS_MeshNode* > TNNMap; + typedef NCollection_Map< SMESH_Link, SMESH_Link > TLinkMap; + + //-------------------------------------------------------------------------------- + /*! + * \brief Intersected face side storing a node created at this intersection + * and a intersected face + */ + struct CutLink + { + bool myReverse; + const SMDS_MeshNode* myNode[2]; // side nodes + mutable SMESH_NodeXYZ myIntNode; // intersection node + const SMDS_MeshElement* myFace; // cutter face + int myIndex; // index of a node on the same link + + CutLink(const SMDS_MeshNode* node1=0, + const SMDS_MeshNode* node2=0, + const SMDS_MeshElement* face=0, + const int index=0) { Set ( node1, node2, face, index ); } + + void Set( const SMDS_MeshNode* node1, + const SMDS_MeshNode* node2, + const SMDS_MeshElement* face, + const int index=0) + { + myNode[0] = node1; myNode[1] = node2; myFace = face; myIndex = index; myReverse = false; + if ( myNode[0] && ( myReverse = ( myNode[0]->GetID() > myNode[1]->GetID() ))) + std::swap( myNode[0], myNode[1] ); + } + const SMDS_MeshNode* IntNode() const { return myIntNode.Node(); } + const SMDS_MeshNode* Node1() const { return myNode[ myReverse ]; } + const SMDS_MeshNode* Node2() const { return myNode[ !myReverse ]; } + + static Standard_Integer HashCode(const CutLink& link, + const Standard_Integer upper) + { + Standard_Integer n = ( link.myNode[0]->GetID() + + link.myNode[1]->GetID() + + link.myIndex ); + return ::HashCode( n, upper ); + } + static Standard_Boolean IsEqual(const CutLink& link1, const CutLink& link2 ) + { + return ( link1.myNode[0] == link2.myNode[0] && + link1.myNode[1] == link2.myNode[1] && + link1.myIndex == link2.myIndex ); + } + }; + + typedef NCollection_Map< CutLink, CutLink > TCutLinkMap; + + //-------------------------------------------------------------------------------- + /*! + * \brief Part of a divided face edge + */ + struct EdgePart + { + const SMDS_MeshNode* myNode1; + const SMDS_MeshNode* myNode2; + int myIndex; // positive -> side index, negative -> State + const SMDS_MeshElement* myFace; + + enum State { _INTERNAL = -1, _COPLANAR = -2 }; + + void Set( const SMDS_MeshNode* Node1, + const SMDS_MeshNode* Node2, + const SMDS_MeshElement* Face = 0, + int EdgeIndex = _INTERNAL ) + { myNode1 = Node1; myNode2 = Node2; myIndex = EdgeIndex; myFace = Face; } + + // bool HasSameNode( const EdgePart& other ) { return ( myNode1 == other.myNode1 || + // myNode1 == other.myNode2 || + // myNode2 == other.myNode1 || + // myNode2 == other.myNode2 ); + // } + bool IsInternal() const { return myIndex < 0; } + bool IsTwin( const EdgePart& e ) const { return myNode1 == e.myNode2 && myNode2 == e.myNode1; } + bool IsSame( const EdgePart& e ) const { + return (( myNode1 == e.myNode2 && myNode2 == e.myNode1 ) || + ( myNode1 == e.myNode1 && myNode2 == e.myNode2 )); } + bool ReplaceCoplanar( const EdgePart& e ); + operator SMESH_Link() const { return SMESH_Link( myNode1, myNode2 ); } + operator gp_Vec() const { return SMESH_NodeXYZ( myNode2 ) - SMESH_NodeXYZ( myNode1 ); } + }; + + //-------------------------------------------------------------------------------- + /*! + * \brief Loop of EdgePart's forming a new face which is a part of CutFace + */ + struct EdgeLoop : public SMDS_PolygonalFaceOfNodes + { + std::vector< const EdgePart* > myLinks; + bool myIsBndConnected; //!< is there a path to CutFace side edges + bool myHasPending; //!< an edge encounters twice + + EdgeLoop() : SMDS_PolygonalFaceOfNodes( std::vector() ) {} + void Clear() { myLinks.clear(); myIsBndConnected = false; myHasPending = false; } + bool SetConnected() { bool was = myIsBndConnected; myIsBndConnected = true; return !was; } + bool Contains( const SMDS_MeshNode* n ) const + { + for ( size_t i = 0; i < myLinks.size(); ++i ) + if ( myLinks[i]->myNode1 == n ) return true; + return false; + } + virtual int NbNodes() const { return myLinks.size(); } + virtual SMDS_ElemIteratorPtr nodesIterator() const + { + return setNodes(), SMDS_PolygonalFaceOfNodes::nodesIterator(); + } + virtual SMDS_NodeIteratorPtr nodeIterator() const + { + return setNodes(), SMDS_PolygonalFaceOfNodes::nodeIterator(); + } + void setNodes() const //!< set nodes to SMDS_PolygonalFaceOfNodes + { + EdgeLoop* me = const_cast( this ); + me->myNodes.resize( NbNodes() ); + size_t iMin = 0; + for ( size_t i = 1; i < myNodes.size(); ++i ) { + if ( myLinks[ i ]->myNode1->GetID() < myLinks[ iMin ]->myNode1->GetID() ) + iMin = i; + } + for ( size_t i = 0; i < myNodes.size(); ++i ) + me->myNodes[ i ] = myLinks[ ( iMin + i ) % myNodes.size() ]->myNode1; + } + }; + + //-------------------------------------------------------------------------------- + /*! + * \brief Set of EdgeLoop's constructed from a CutFace + */ + struct EdgeLoopSet + { + std::vector< EdgeLoop > myLoops; //!< buffer of EdgeLoop's + size_t myNbLoops; //!< number of constructed loops + + const EdgePart* myEdge0; //!< & CutFace.myLinks[0] + size_t myNbUsedEdges; //!< nb of EdgePart's added to myLoops + boost::dynamic_bitset<> myIsUsedEdge; //!< is i-th EdgePart of CutFace is in any EdgeLoop + std::vector< EdgeLoop* > myLoopOfEdge; //!< EdgeLoop of CutFace.myLinks[i] + std::vector< EdgePart* > myCandidates; //!< EdgePart's starting at the same node + + EdgeLoopSet(): myLoops(100) {} + + void Init( const std::vector< EdgePart >& edges ) + { + size_t nb = edges.size(); + myEdge0 = & edges[0]; + myNbLoops = 0; + myNbUsedEdges = 0; + myIsUsedEdge.reset(); + myIsUsedEdge.resize( nb, false ); + myLoopOfEdge.clear(); + myLoopOfEdge.resize( nb, (EdgeLoop*) 0 ); + } + EdgeLoop& AddNewLoop() + { + if ( ++myNbLoops >= myLoops.size() ) + myLoops.resize( myNbLoops + 10 ); + myLoops[ myNbLoops-1 ].Clear(); + return myLoops[ myNbLoops-1 ]; + } + bool AllEdgesUsed() const { return myNbUsedEdges == myLoopOfEdge.size(); } + + bool AddEdge( EdgePart& edge ) + { + size_t i = Index( edge ); + if ( myIsUsedEdge[ i ]) + return false; + myLoops[ myNbLoops-1 ].myLinks.push_back( &edge ); + myLoopOfEdge[ i ] = & myLoops[ myNbLoops-1 ]; + myIsUsedEdge[ i ] = true; + ++myNbUsedEdges; + return true; + } + void Erase( EdgeLoop* loop ) + { + for ( size_t iE = 0; iE < loop->myLinks.size(); ++iE ) + myLoopOfEdge[ Index( *loop->myLinks[ iE ] )] = 0; + loop->Clear(); + } + size_t Index( const EdgePart& edge ) const { return &edge - myEdge0; } + EdgeLoop* GetLoopOf( const EdgePart* edge ) { return myLoopOfEdge[ Index( *edge )]; } + }; + + //-------------------------------------------------------------------------------- + /*! + * \brief Intersections of a face + */ + struct CutFace + { + mutable std::vector< EdgePart > myLinks; + const SMDS_MeshElement* myInitFace; + + CutFace( const SMDS_MeshElement* face ): myInitFace( face ) {} + void AddEdge( const CutLink& p1, + const CutLink& p2, + const SMDS_MeshElement* cutter, + const int nbOnPlane, + const int iNotOnPlane = -1) const; + void AddPoint( const CutLink& p1, const CutLink& p2, double tol ) const; + bool ReplaceNodes( const TNNMap& theRm2KeepMap ) const; + bool IsCut() const; + int NbInternalEdges() const; + void MakeLoops( EdgeLoopSet& loops, const gp_XYZ& theFaceNorm ) const; + bool RemoveInternalLoops( EdgeLoopSet& theLoops ) const; + void CutOffLoops( EdgeLoopSet& theLoops, + const double theSign, + const std::vector< gp_XYZ >& theNormals, + std::vector< EdgePart >& theCutOffLinks, + TLinkMap& theCutOffCoplanarLinks) const; + void InitLinks() const; + bool IsCoplanar( const EdgePart* edge ) const; + + static Standard_Integer HashCode(const CutFace& f, const Standard_Integer upper) + { + return ::HashCode( f.myInitFace->GetID(), upper ); + } + static Standard_Boolean IsEqual(const CutFace& f1, const CutFace& f2 ) + { + return f1.myInitFace == f2.myInitFace; + } + void Dump() const; + + private: + + EdgePart* getTwin( const EdgePart* edge ) const; + }; + + typedef NCollection_Map< CutFace, CutFace > TCutFaceMap; + + //-------------------------------------------------------------------------------- + /*! + * \brief Intersection point of two edges of co-planar triangles + */ + struct IntPoint2D + { + size_t myEdgeInd[2]; //!< edge indices of triangles + double myU [2]; //!< parameter [0,1] on edges of triangles + SMESH_NodeXYZ myNode; //!< intersection node + bool myIsCollinear;//!< edges are collinear + + IntPoint2D() : myIsCollinear( false ) {} + + void InitLink( CutLink& link, int iFace, const std::vector< SMESH_NodeXYZ >& nodes ) const + { + link.Set( nodes[ myEdgeInd[ iFace ] ].Node(), + nodes[( myEdgeInd[ iFace ] + 1 ) % nodes.size() ].Node(), + link.myFace ); + link.myIntNode = myNode; + } + const SMDS_MeshNode* Node() const { return myNode.Node(); } + }; + struct IntPoint2DCompare + { + int myI; + IntPoint2DCompare( int iFace=0 ): myI( iFace ) {} + bool operator() ( const IntPoint2D* ip1, const IntPoint2D* ip2 ) const + { + return ip1->myU[ myI ] < ip2->myU[ myI ]; + } + bool operator() ( const IntPoint2D& ip1, const IntPoint2D& ip2 ) const + { + return ip1.myU[ myI ] < ip2.myU[ myI ]; + } + }; + typedef boost::container::flat_set< IntPoint2D, IntPoint2DCompare > TIntPointSet; + typedef boost::container::flat_set< IntPoint2D*, IntPoint2DCompare > TIntPointPtrSet; + + //-------------------------------------------------------------------------------- + /*! + * \brief Face used to find translated position of the node + */ + struct Face + { + const SMDS_MeshElement* myFace; + SMESH_TNodeXYZ myNode1; //!< nodes neighboring another node of myFace + SMESH_TNodeXYZ myNode2; + const gp_XYZ* myNorm; + bool myNodeRightOrder; + void operator=(const SMDS_MeshElement* f) { myFace = f; } + const SMDS_MeshElement* operator->() { return myFace; } + void SetNodes( int i0, int i1 ) //!< set myNode's + { + myNode1.Set( myFace->GetNode( i1 )); + int i2 = ( i0 - 1 + myFace->NbCornerNodes() ) % myFace->NbCornerNodes(); + if ( i2 == i1 ) + i2 = ( i0 + 1 ) % myFace->NbCornerNodes(); + myNode2.Set( myFace->GetNode( i2 )); + myNodeRightOrder = ( Abs( i2-i1 ) == 1 ) ? i2 > i1 : i2 < i1; + } + void SetOldNodes( const SMDS_Mesh& theSrcMesh ) + { + myNode1.Set( theSrcMesh.FindNode( myNode1->GetID() )); + myNode2.Set( theSrcMesh.FindNode( myNode2->GetID() )); + } + bool SetNormal( const std::vector< gp_XYZ >& faceNormals ) + { + myNorm = & faceNormals[ myFace->GetID() ]; + return ( myNorm->SquareModulus() > gp::Resolution() * gp::Resolution() ); + } + const gp_XYZ& Norm() const { return *myNorm; } + }; + + //-------------------------------------------------------------------------------- + /*! + * \brief Offset plane used to find translated position of the node + */ + struct OffsetPlane + { + gp_XYZ myNode; + Face* myFace; + gp_Pln myPln; + gp_Lin myLines[2]; //!< line of intersection with neighbor OffsetPlane's + bool myIsLineOk[2]; + double myWeight[2]; + + void Init( const gp_XYZ& node, Face& tria, double offset ) + { + myNode = node; + myFace = & tria; + myPln = gp_Pln( node + tria.Norm() * offset, tria.Norm() ); + myIsLineOk[0] = myIsLineOk[1] = false; + myWeight[0] = myWeight[1] = 0; + } + bool ComputeIntersectionLine( OffsetPlane& pln ); + void SetSkewLine( const gp_Lin& line ); + gp_XYZ GetCommonPoint( int & nbOkPoints, double& sumWeight ); + gp_XYZ ProjectNodeOnLine( int & nbOkPoints ); + double Weight() const { return myWeight[0] + myWeight[1]; } + }; + + //================================================================================ + /*! + * \brief Set the second line + */ + //================================================================================ + + void OffsetPlane::SetSkewLine( const gp_Lin& line ) + { + myLines[1] = line; + gp_XYZ n = myLines[0].Direction().XYZ() ^ myLines[1].Direction().XYZ(); + if (( myIsLineOk[1] = n.SquareModulus() > gp::Resolution() )) + myPln = gp_Pln( myPln.Location(), n ); + } + + //================================================================================ + /*! + * \brief Project myNode on myLine[0] + */ + //================================================================================ + + gp_XYZ OffsetPlane::ProjectNodeOnLine( int & nbOkPoints ) + { + gp_XYZ p = gp::Origin().XYZ(); + if ( myIsLineOk[0] ) + { + gp_Vec l2n( myLines[0].Location(), myNode ); + double u = l2n * myLines[0].Direction(); + p = myLines[0].Location().XYZ() + u * myLines[0].Direction().XYZ(); + ++nbOkPoints; + } + return p; + } + + //================================================================================ + /*! + * \brief Computes intersection point of myLines + */ + //================================================================================ + + gp_XYZ OffsetPlane::GetCommonPoint( int & nbOkPoints, double& sumWeight ) + { + if ( !myIsLineOk[0] || !myIsLineOk[1] ) + { + // sumWeight += myWeight[0]; + // return ProjectNodeOnLine( nbOkPoints ) * myWeight[0]; + return gp::Origin().XYZ(); + } + + gp_XYZ p; + + gp_Vec lPerp0 = myLines[0].Direction().XYZ() ^ myPln.Axis().Direction().XYZ(); + double dot01 = lPerp0 * myLines[1].Direction().XYZ(); + if ( Abs( dot01 ) > 0.05 ) + { + gp_Vec l0l1 = myLines[1].Location().XYZ() - myLines[0].Location().XYZ(); + double u1 = - ( lPerp0 * l0l1 ) / dot01; + p = ( myLines[1].Location().XYZ() + myLines[1].Direction().XYZ() * u1 ); + } + else + { + gp_Vec lv0( myLines[0].Location(), myNode), lv1(myLines[1].Location(), myNode ); + double dot0( lv0 * myLines[0].Direction() ), dot1( lv1 * myLines[1].Direction() ); + p = 0.5 * ( myLines[0].Location().XYZ() + myLines[0].Direction().XYZ() * dot0 ); + p += 0.5 * ( myLines[1].Location().XYZ() + myLines[1].Direction().XYZ() * dot1 ); + } + + sumWeight += Weight(); + ++nbOkPoints; + + return p * Weight(); + } + + //================================================================================ + /*! + * \brief Compute line of intersection of 2 planes + */ + //================================================================================ + + bool OffsetPlane::ComputeIntersectionLine( OffsetPlane& theNextPln ) + { + const gp_XYZ& n1 = myFace->Norm(); + const gp_XYZ& n2 = theNextPln.myFace->Norm(); + + gp_XYZ lineDir = n1 ^ n2; + gp_Pnt linePos; + + double x = Abs( lineDir.X() ); + double y = Abs( lineDir.Y() ); + double z = Abs( lineDir.Z() ); + + int cooMax; // max coordinate + if (x > y) { + if (x > z) cooMax = 1; + else cooMax = 3; + } + else { + if (y > z) cooMax = 2; + else cooMax = 3; + } + + bool ok = true; + if ( Abs( lineDir.Coord( cooMax )) < 0.05 ) + { + // parallel planes - intersection is an offset of the common edge + linePos = 0.5 * ( myPln.Location().XYZ() + theNextPln.myPln.Location().XYZ() ); + lineDir = myNode - myFace->myNode2; + ok = false; + myWeight[0] = 0; + } + else + { + // the constants in the 2 plane equations + double d1 = - ( n1 * myPln.Location().XYZ() ); + double d2 = - ( n2 * theNextPln.myPln.Location().XYZ() ); + + switch ( cooMax ) { + case 1: + linePos.SetX( 0 ); + linePos.SetY(( d2*n1.Z() - d1*n2.Z()) / lineDir.X() ); + linePos.SetZ(( d1*n2.Y() - d2*n1.Y()) / lineDir.X() ); + break; + case 2: + linePos.SetX(( d1*n2.Z() - d2*n1.Z()) / lineDir.Y() ); + linePos.SetY( 0 ); + linePos.SetZ(( d2*n1.X() - d1*n2.X()) / lineDir.Y() ); + break; + case 3: + linePos.SetX(( d2*n1.Y() - d1*n2.Y()) / lineDir.Z() ); + linePos.SetY(( d1*n2.X() - d2*n1.X()) / lineDir.Z() ); + linePos.SetZ( 0 ); + } + myWeight[0] = lineDir.SquareModulus(); + if ( n1 * n2 < 0 ) + myWeight[0] = 2. - myWeight[0]; + } + myLines [ 0 ].SetDirection( lineDir ); + myLines [ 0 ].SetLocation ( linePos ); + myIsLineOk[ 0 ] = ok; + + theNextPln.myLines [ 1 ] = myLines[ 0 ]; + theNextPln.myIsLineOk[ 1 ] = ok; + theNextPln.myWeight [ 1 ] = myWeight[ 0 ]; + + return ok; + } + + //================================================================================ + /*! + * \brief Return a translated position of a node + * \param [in] new2OldNodes - new and old nodes + * \param [in] faceNormals - normals to input faces + * \param [in] theSrcMesh - initial mesh + * \param [in] theNewPos - a computed normal + * \return bool - true if theNewPos is computed + */ + //================================================================================ + + bool getTranslatedPosition( const SMDS_MeshNode* theNewNode, + const double theOffset, + const double theTol, + const double theSign, + const std::vector< gp_XYZ >& theFaceNormals, + SMDS_Mesh& theSrcMesh, + gp_XYZ& theNewPos) + { + bool useOneNormal = true; + + // check if theNewNode needs an average position, i.e. theNewNode is convex + // SMDS_ElemIteratorPtr faceIt = theNewNode->GetInverseElementIterator(); + // const SMDS_MeshElement* f0 = faceIt->next(); + // const gp_XYZ& norm0 = theFaceNormals[ f0->GetID() ]; + // const SMESH_NodeXYZ nodePos = theNewNode; + // while ( faceIt->more() ) + // { + // const SMDS_MeshElement* f = faceIt->next(); + // const int nodeInd = f->GetNodeIndex( theNewNode ); + // SMESH_NodeXYZ nodePos2 = f->GetWrapNode( nodeInd + 1 ); + // try { + // const gp_XYZ nnDir = ( nodePos2 - nodePos ).Normalized(); + // } + // catch { + // continue; + // } + // const double dot = norm0 * nnDir; + // bool isConvex = + + + + // get faces surrounding theNewNode and sort them + Face faces[ theMaxNbFaces ]; + SMDS_ElemIteratorPtr faceIt = theNewNode->GetInverseElementIterator(); + faces[0] = faceIt->next(); + while ( !faces[0].SetNormal( theFaceNormals ) && faceIt->more() ) + faces[0] = faceIt->next(); + int i0 = faces[0]->GetNodeIndex( theNewNode ); + int i1 = ( i0 + 1 ) % faces[0]->NbCornerNodes(); + faces[0].SetNodes( i0, i1 ); + TIDSortedElemSet elemSet, avoidSet; + int iFace = 0; + const SMDS_MeshElement* f; + for ( ; faceIt->more(); faceIt->next() ) + { + avoidSet.insert( faces[ iFace ].myFace ); + f = SMESH_MeshAlgos::FindFaceInSet( theNewNode, faces[ iFace ].myNode2.Node(), + elemSet, avoidSet, &i0, &i1 ); + if ( !f ) + { + std::reverse( &faces[0], &faces[0] + iFace + 1 ); + for ( int i = 0; i <= iFace; ++i ) + { + std::swap( faces[i].myNode1, faces[i].myNode2 ); + faces[i].myNodeRightOrder = !faces[i].myNodeRightOrder; + } + f = SMESH_MeshAlgos::FindFaceInSet( theNewNode, faces[ iFace ].myNode2.Node(), + elemSet, avoidSet, &i0, &i1 ); + if ( !f ) + break; + } + faces[ ++iFace ] = f; + faces[ iFace ].SetNodes( i0, i1 ); + faces[ iFace ].SetNormal( theFaceNormals ); + } + int nbFaces = Min( iFace + 1, (int)theMaxNbFaces ); + + theNewPos.SetCoord( 0, 0, 0 ); + gp_XYZ oldXYZ = SMESH_NodeXYZ( theNewNode ); + + // check if all faces are co-planar + bool isPlanar = true; + const double tol = 1e-2; + for ( int i = 1; i < nbFaces && isPlanar; ++i ) + isPlanar = ( faces[i].Norm() - faces[i-1].Norm() ).SquareModulus() < tol*tol; + + if ( isPlanar ) + { + theNewPos = oldXYZ + faces[0].Norm() * theOffset; + return useOneNormal; + } + + // prepare OffsetPlane's + OffsetPlane pln[ theMaxNbFaces ]; + for ( int i = 0; i < nbFaces; ++i ) + { + faces[i].SetOldNodes( theSrcMesh ); + pln[i].Init( oldXYZ, faces[i], theOffset ); + } + // intersect neighboring OffsetPlane's + int nbOkPoints = 0; + for ( int i = 1; i < nbFaces; ++i ) + nbOkPoints += pln[ i-1 ].ComputeIntersectionLine( pln[ i ]); + nbOkPoints += pln[ nbFaces-1 ].ComputeIntersectionLine( pln[ 0 ]); + + // move intersection lines to over parallel planes + if ( nbOkPoints > 1 ) + for ( int i = 0; i < nbFaces; ++i ) + if ( pln[i].myIsLineOk[0] && !pln[i].myIsLineOk[1] ) + for ( int j = 1; j < nbFaces && !pln[i].myIsLineOk[1]; ++j ) + { + int i2 = ( i + j ) % nbFaces; + if ( pln[i2].myIsLineOk[0] ) + pln[i].SetSkewLine( pln[i2].myLines[0] ); + } + + // get the translated position + nbOkPoints = 0; + double sumWegith = 0; + const double minWeight = Sin( 30 * M_PI / 180. ) * Sin( 30 * M_PI / 180. ); + for ( int i = 0; i < nbFaces; ++i ) + if ( pln[ i ].Weight() > minWeight ) + theNewPos += pln[ i ].GetCommonPoint( nbOkPoints, sumWegith ); + + if ( nbOkPoints == 0 ) + { + // there is only one feature edge; + // find the theNewPos by projecting oldXYZ to any intersection line + for ( int i = 0; i < nbFaces; ++i ) + theNewPos += pln[ i ].ProjectNodeOnLine( nbOkPoints ); + + if ( nbOkPoints == 0 ) + { + theNewPos = oldXYZ + faces[0].Norm() * theOffset; + return useOneNormal; + } + sumWegith = nbOkPoints; + } + theNewPos /= sumWegith; + + + // mark theNewNode if it is concave + useOneNormal = false; + gp_Vec moveVec( oldXYZ, theNewPos ); + for ( int i = 0, iPrev = nbFaces-1; i < nbFaces; iPrev = i++ ) + { + gp_Vec nodeVec( oldXYZ, faces[ i ].myNode1 ); + double u = ( moveVec * nodeVec ) / nodeVec.SquareMagnitude(); + if ( u > 0.5 ) // param [0,1] on nodeVec + { + theNewNode->setIsMarked( true ); + } + if ( !useOneNormal ) + { + gp_XYZ inFaceVec = faces[ i ].Norm() ^ nodeVec.XYZ(); + double dot = inFaceVec * faces[ iPrev ].Norm(); + if ( !faces[ i ].myNodeRightOrder ) + dot *= -1; + if ( dot * theSign < 0 ) + { + gp_XYZ p1 = oldXYZ + faces[ i ].Norm() * theOffset; + gp_XYZ p2 = oldXYZ + faces[ iPrev ].Norm() * theOffset; + useOneNormal = ( p1 - p2 ).SquareModulus() > theTol * theTol; + } + } + if ( useOneNormal && theNewNode->isMarked() ) + break; + } + + return useOneNormal; + } + + //-------------------------------------------------------------------------------- + /*! + * \brief Intersect faces of a mesh + */ + struct Intersector + { + SMDS_Mesh* myMesh; + double myTol, myEps; + const std::vector< gp_XYZ >& myNormals; + TCutLinkMap myCutLinks; //!< assure sharing of new nodes + TCutFaceMap myCutFaces; + TNNMap myRemove2KeepNodes; //!< node merge map + + // data to intersect 2 faces + const SMDS_MeshElement* myFace1; + const SMDS_MeshElement* myFace2; + std::vector< SMESH_NodeXYZ > myNodes1, myNodes2; + std::vector< double > myDist1, myDist2; + int myInd1, myInd2; // coordinate indices on an axis-aligned plane + int myNbOnPlane1, myNbOnPlane2; + TIntPointSet myIntPointSet; + + Intersector( SMDS_Mesh* mesh, double tol, const std::vector< gp_XYZ >& normals ) + : myMesh( mesh ), + myTol( tol ), + myEps( 1e-100 ), + //myEps( Sqrt( std::numeric_limits::min() )), + //myEps( gp::Resolution() ), + myNormals( normals ) + {} + void Cut( const SMDS_MeshElement* face1, + const SMDS_MeshElement* face2, + const int nbCommonNodes ); + void MakeNewFaces( SMESH_MeshAlgos::TEPairVec& theNew2OldFaces, + SMESH_MeshAlgos::TNPairVec& theNew2OldNodes, + const double theSign ); + + private: + + bool isPlaneIntersected( const gp_XYZ& n2, + const double d2, + const std::vector< SMESH_NodeXYZ >& nodes1, + std::vector< double > & dist1, + int & nbOnPlane1, + int & iNotOnPlane1); + void computeIntervals( const std::vector< SMESH_NodeXYZ >& nodes, + const std::vector< double >& dist, + const int nbOnPln, + const int iMaxCoo, + double * u, + int* iE); + void cutCoplanar(); + void addLink ( CutLink& link ); + bool findLink( CutLink& link ); + bool coincide( const gp_XYZ& p1, const gp_XYZ& p2, const double tol ) const + { + return ( p1 - p2 ).SquareModulus() < tol * tol; + } + gp_XY p2D( const gp_XYZ& p ) const { return gp_XY( p.Coord( myInd1 ), p.Coord( myInd2 )); } + + void intersectLink( const std::vector< SMESH_NodeXYZ >& nodes1, + const std::vector< double > & dist1, + const int iEdge1, + const SMDS_MeshElement* face2, + CutLink& link1); + void findIntPointOnPlane( const std::vector< SMESH_NodeXYZ >& nodes, + const std::vector< double > & dist, + CutLink& link ); + void replaceIntNode( const SMDS_MeshNode* nToKeep, const SMDS_MeshNode* nToRemove ); + void computeIntPoint( const double u1, + const double u2, + const int iE1, + const int iE2, + CutLink & link, + const SMDS_MeshNode* & node1, + const SMDS_MeshNode* & node2); + void cutCollinearLink( const int iNotOnPlane1, + const std::vector< SMESH_NodeXYZ >& nodes1, + const SMDS_MeshElement* face2, + const CutLink& link1, + const CutLink& link2); + void setPlaneIndices( const gp_XYZ& planeNorm ); + bool intersectEdgeEdge( const gp_XY s1p0, const gp_XY s1p1, + const gp_XY s2p0, const gp_XY s2p1, + double & t1, double & t2, + bool & isCollinear ); + bool intersectEdgeEdge( int iE1, int iE2, IntPoint2D& intPoint ); + bool isPointInTriangle( const gp_XYZ& p, const std::vector< SMESH_NodeXYZ >& nodes ); + void intersectNewEdges( const CutFace& theCFace ); + const SMDS_MeshNode* createNode( const gp_XYZ& p ); + }; + + //================================================================================ + /*! + * \brief Return coordinate index with maximal abs value + */ + //================================================================================ + + int MaxIndex( const gp_XYZ& x ) + { + int iMaxCoo = ( Abs( x.X()) < Abs( x.Y() )) + 1; + if ( Abs( x.Coord( iMaxCoo )) < Abs( x.Z() )) + iMaxCoo = 3; + return iMaxCoo; + } + //================================================================================ + /*! + * \brief Store a CutLink + */ + //================================================================================ + + const SMDS_MeshNode* Intersector::createNode( const gp_XYZ& p ) + { + const SMDS_MeshNode* n = myMesh->AddNode( p.X(), p.Y(), p.Z() ); + n->setIsMarked( true ); // cut nodes are marked + return n; + } + + //================================================================================ + /*! + * \brief Store a CutLink + */ + //================================================================================ + + void Intersector::addLink( CutLink& link ) + { + link.myIndex = 0; + const CutLink* added = & myCutLinks.Added( link ); + while ( added->myIntNode.Node() != link.myIntNode.Node() ) + { + if ( !added->myIntNode ) + { + added->myIntNode = link.myIntNode; + break; + } + else + { + link.myIndex++; + added = & myCutLinks.Added( link ); + } + } + link.myIndex = 0; + } + + //================================================================================ + /*! + * \brief Find a CutLink with an intersection point coincident with that of a given link + */ + //================================================================================ + + bool Intersector::findLink( CutLink& link ) + { + link.myIndex = 0; + while ( myCutLinks.Contains( link )) + { + const CutLink* added = & myCutLinks.Added( link ); + if ( !!added->myIntNode && coincide( added->myIntNode, link.myIntNode, myTol )) + { + link.myIntNode = added->myIntNode; + return true; + } + link.myIndex++; + } + return false; + } + + //================================================================================ + /*! + * \brief Check if a triangle intersects the plane of another triangle + * \param [in] nodes1 - nodes of triangle 1 + * \param [in] n2 - normal of triangle 2 + * \param [in] d2 - a constant of the plane equation 2 + * \param [out] dist1 - distance of nodes1 from the plane 2 + * \param [out] nbOnPlane - number of nodes1 lying on the plane 2 + * \return bool - true if the triangle intersects the plane 2 + */ + //================================================================================ + + bool Intersector::isPlaneIntersected( const gp_XYZ& n2, + const double d2, + const std::vector< SMESH_NodeXYZ >& nodes1, + std::vector< double > & dist1, + int & nbOnPlane1, + int & iNotOnPlane1) + { + iNotOnPlane1 = nbOnPlane1 = 0; + dist1.resize( nodes1.size() ); + for ( size_t i = 0; i < nodes1.size(); ++i ) + { + dist1[i] = n2 * nodes1[i] + d2; + if ( Abs( dist1[i] ) < myTol ) + { + ++nbOnPlane1; + dist1[i] = 0.; + } + else + { + iNotOnPlane1 = i; + } + } + if ( nbOnPlane1 == 0 ) + for ( size_t i = 0; i < nodes1.size(); ++i ) + if ( dist1[iNotOnPlane1] * dist1[i] < 0 ) + return true; + + return nbOnPlane1; + } + + //================================================================================ + /*! + * \brief Compute parameters on the plane intersection line of intersections + * of edges of a triangle + * \param [in] nodes - triangle nodes + * \param [in] dist - distance of triangle nodes from the plane of another triangle + * \param [in] nbOnPln - number of nodes lying on the plane of another triangle + * \param [in] iMaxCoo - index of coordinate of max component of the plane intersection line + * \param [out] u - two computed parameters on the plane intersection line + * \param [out] iE - indices of intersected edges + */ + //================================================================================ + + void Intersector::computeIntervals( const std::vector< SMESH_NodeXYZ >& nodes, + const std::vector< double >& dist, + const int nbOnPln, + const int iMaxCoo, + double * u, + int* iE) + { + if ( nbOnPln == 3 ) + { + u[0] = u[1] = 1e+100; + return; + } + int nb = 0; + int i1 = 2, i2 = 0; + if ( nbOnPln == 1 && ( dist[i1] == 0. || dist[i2] == 0 )) + { + int i = dist[i1] == 0 ? i1 : i2; + u [ 1 ] = nodes[ i ].Coord( iMaxCoo ); + iE[ 1 ] = i; + i1 = i2++; + } + for ( ; i2 < 3 && nb < 2; i1 = i2++ ) + { + double dd = dist[i1] - dist[i2]; + if ( dd != 0. && dist[i2] * dist[i1] <= 0. ) + { + double x1 = nodes[i1].Coord( iMaxCoo ); + double x2 = nodes[i2].Coord( iMaxCoo ); + u [ nb ] = x1 + ( x2 - x1 ) * dist[i1] / dd; + iE[ nb ] = i1; + ++nb; + } + } + if ( u[0] > u[1] ) + { + std::swap( u [0], u [1] ); + std::swap( iE[0], iE[1] ); + } + } + + //================================================================================ + /*! + * \brief Try to find an intersection node on a link collinear with the plane intersection line + */ + //================================================================================ + + void Intersector::findIntPointOnPlane( const std::vector< SMESH_NodeXYZ >& nodes, + const std::vector< double > & dist, + CutLink& link ) + { + int i1 = ( dist[0] == 0 ? 0 : 1 ), i2 = ( dist[2] == 0 ? 2 : 1 ); + CutLink link2 = link; + link2.Set( nodes[i1].Node(), nodes[i2].Node(), 0 ); + if ( findLink( link2 )) + link.myIntNode = link2.myIntNode; + } + + //================================================================================ + /*! + * \brief Compute intersection point of a link1 with a face2 + */ + //================================================================================ + + void Intersector::intersectLink( const std::vector< SMESH_NodeXYZ >& nodes1, + const std::vector< double > & dist1, + const int iEdge1, + const SMDS_MeshElement* face2, + CutLink& link1) + { + const int iEdge2 = ( iEdge1 + 1 ) % nodes1.size(); + const SMESH_NodeXYZ& p1 = nodes1[ iEdge1 ]; + const SMESH_NodeXYZ& p2 = nodes1[ iEdge2 ]; + + link1.Set( p1.Node(), p2.Node(), face2 ); + const CutLink* link = & myCutLinks.Added( link1 ); + if ( !link->IntNode() ) + { + if ( dist1[ iEdge1 ] == 0. ) link1.myIntNode = p1; + else if ( dist1[ iEdge2 ] == 0. ) link1.myIntNode = p2; + else + { + gp_XYZ p = p1 + ( p2 - p1 ) * dist1[ iEdge1 ] / ( dist1[ iEdge1 ] - dist1[ iEdge2 ]); + (gp_XYZ&)link1.myIntNode = p; + } + } + else + { + gp_XYZ p = p1 + ( p2 - p1 ) * dist1[ iEdge1 ] / ( dist1[ iEdge1 ] - dist1[ iEdge2 ]); + while ( link->IntNode() ) + { + if ( coincide( p, link->myIntNode, myTol )) + { + link1.myIntNode = link->myIntNode; + break; + } + link1.myIndex++; + link = & myCutLinks.Added( link1 ); + } + if ( !link1.IntNode() ) + { + if ( dist1[ iEdge1 ] == 0. ) link1.myIntNode = p1; + else if ( dist1[ iEdge2 ] == 0. ) link1.myIntNode = p2; + else (gp_XYZ&)link1.myIntNode = p; + } + } + } + + //================================================================================ + /*! + * \brief Store node replacement in myCutFaces + */ + //================================================================================ + + void Intersector::replaceIntNode( const SMDS_MeshNode* nToKeep, + const SMDS_MeshNode* nToRemove ) + { + if ( nToKeep == nToRemove ) + return; + if ( nToRemove->GetID() < nToKeep->GetID() ) // keep node with lower ID + myRemove2KeepNodes.Bind((void*) nToKeep, nToRemove ); + else + myRemove2KeepNodes.Bind((void*) nToRemove, nToKeep ); + } + + //================================================================================ + /*! + * \brief Compute intersection point on a link of either of faces by choosing + * a link whose parameter on the intersection line in maximal + * \param [in] u1 - parameter on the intersection line of link iE1 of myFace1 + * \param [in] u2 - parameter on the intersection line of link iE2 of myFace2 + * \param [in] iE1 - index of a link myFace1 + * \param [in] iE2 - index of a link myFace2 + * \param [out] link - CutLink storing the intersection point + * \param [out] node1 - a node of the 2nd link if two links intersect + * \param [out] node2 - a node of the 2nd link if two links intersect + */ + //================================================================================ + + void Intersector::computeIntPoint( const double u1, + const double u2, + const int iE1, + const int iE2, + CutLink & link, + const SMDS_MeshNode* & node1, + const SMDS_MeshNode* & node2) + { + if ( u1 > u2 + myTol ) + { + intersectLink( myNodes1, myDist1, iE1, myFace2, link ); + node1 = node2 = 0; + if ( myNbOnPlane2 == 2 ) + findIntPointOnPlane( myNodes2, myDist2, link ); + } + else if ( u2 > u1 + myTol ) + { + intersectLink( myNodes2, myDist2, iE2, myFace1, link ); + node1 = node2 = 0; + if ( myNbOnPlane1 == 2 ) + findIntPointOnPlane( myNodes1, myDist1, link ); + } + else // edges of two faces intersect the line at the same point + { + CutLink link2; + intersectLink( myNodes1, myDist1, iE1, myFace2, link ); + intersectLink( myNodes2, myDist2, iE2, myFace1, link2 ); + node1 = link2.Node1(); + node2 = link2.Node2(); + + if ( !link.IntNode() && link2.IntNode() ) + link.myIntNode = link2.myIntNode; + + else if ( !link.IntNode() && !link2.IntNode() ) + (gp_XYZ&)link.myIntNode = 0.5 * ( link.myIntNode + link2.myIntNode ); + + else if ( link.IntNode() && link2.IntNode() ) + replaceIntNode( link.IntNode(), link2.IntNode() ); + } + } + + //================================================================================ + /*! + * \brief Add intersections to a link collinear with the intersection line + */ + //================================================================================ + + void Intersector::cutCollinearLink( const int iNotOnPlane1, + const std::vector< SMESH_NodeXYZ >& nodes1, + const SMDS_MeshElement* face2, + const CutLink& link1, + const CutLink& link2) + + { + int iN1 = ( iNotOnPlane1 + 1 ) % 3; + int iN2 = ( iNotOnPlane1 + 2 ) % 3; + CutLink link( nodes1[ iN1 ].Node(), nodes1[ iN2 ].Node(), face2 ); + if ( link1.myFace != face2 ) + { + link.myIntNode = link1.myIntNode; + addLink( link ); + } + if ( link2.myFace != face2 ) + { + link.myIntNode = link2.myIntNode; + addLink( link ); + } + } + + //================================================================================ + /*! + * \brief Choose indices on an axis-aligned plane + */ + //================================================================================ + + void Intersector::setPlaneIndices( const gp_XYZ& planeNorm ) + { + switch ( MaxIndex( planeNorm )) { + case 1: myInd1 = 2; myInd2 = 3; break; + case 2: myInd1 = 3; myInd2 = 1; break; + case 3: myInd1 = 1; myInd2 = 2; break; + } + } + + //================================================================================ + /*! + * \brief Intersect two faces + */ + //================================================================================ + + void Intersector::Cut( const SMDS_MeshElement* face1, + const SMDS_MeshElement* face2, + const int nbCommonNodes) + { + myFace1 = face1; + myFace2 = face2; + myNodes1.assign( face1->begin_nodes(), face1->end_nodes() ); + myNodes2.assign( face2->begin_nodes(), face2->end_nodes() ); + + const gp_XYZ& n1 = myNormals[ face1->GetID() ]; + const gp_XYZ& n2 = myNormals[ face2->GetID() ]; + + // check if triangles intersect + int iNotOnPlane1, iNotOnPlane2; + const double d2 = -( n2 * myNodes2[0]); + if ( !isPlaneIntersected( n2, d2, myNodes1, myDist1, myNbOnPlane1, iNotOnPlane1 )) + return; + const double d1 = -( n1 * myNodes1[0]); + if ( !isPlaneIntersected( n1, d1, myNodes2, myDist2, myNbOnPlane2, iNotOnPlane2 )) + return; + + if ( myNbOnPlane1 == 3 || myNbOnPlane2 == 3 )// triangles are co-planar + { + setPlaneIndices( myNbOnPlane1 == 3 ? n2 : n1 ); // choose indices on an axis-aligned plane + cutCoplanar(); + } + else if ( nbCommonNodes < 2 ) // triangle planes intersect + { + gp_XYZ lineDir = n1 ^ n2; // intersection line + + // check if intervals of intersections of triangles with lineDir overlap + + double u1[2], u2 [2]; // parameters on lineDir of edge intersection points { minU, maxU } + int iE1[2], iE2[2]; // indices of edges + int iMaxCoo = MaxIndex( lineDir ); + computeIntervals( myNodes1, myDist1, myNbOnPlane1, iMaxCoo, u1, iE1 ); + computeIntervals( myNodes2, myDist2, myNbOnPlane2, iMaxCoo, u2, iE2 ); + if ( u1[1] < u2[0] - myTol || u2[1] < u1[0] - myTol ) + return; // intervals do not overlap + + // make intersection nodes + + const SMDS_MeshNode *l1n1, *l1n2, *l2n1, *l2n2; + CutLink link1; // intersection with smaller u on lineDir + computeIntPoint( u1[0], u2[0], iE1[0], iE2[0], link1, l1n1, l1n2 ); + CutLink link2; // intersection with larger u on lineDir + computeIntPoint( -u1[1], -u2[1], iE1[1], iE2[1], link2, l2n1, l2n2 ); + + const CutFace& cf1 = myCutFaces.Added( CutFace( face1 )); + const CutFace& cf2 = myCutFaces.Added( CutFace( face2 )); + + if ( coincide( link1.myIntNode, link2.myIntNode, myTol )) + { + // intersection is a point + if ( link1.IntNode() && link2.IntNode() ) + replaceIntNode( link1.IntNode(), link2.IntNode() ); + + CutLink* link = link2.IntNode() ? &link2 : &link1; + if ( !link->IntNode() ) + { + gp_XYZ p = 0.5 * ( link1.myIntNode + link2.myIntNode ); + link->myIntNode.Set( createNode( p )); + } + if ( !link1.IntNode() ) link1.myIntNode = link2.myIntNode; + if ( !link2.IntNode() ) link2.myIntNode = link1.myIntNode; + + cf1.AddPoint( link1, link2, myTol ); + cf2.AddPoint( link1, link2, myTol ); + } + else + { + // intersection is a line segment + if ( !link1.IntNode() ) + link1.myIntNode.Set( createNode( link1.myIntNode )); + if ( !link2.IntNode() ) + link2.myIntNode.Set( createNode( link2.myIntNode )); + + cf1.AddEdge( link1, link2, face2, myNbOnPlane1, iNotOnPlane1 ); + if ( l1n1 ) link1.Set( l1n1, l1n2, face2 ); + if ( l2n1 ) link2.Set( l2n1, l2n2, face2 ); + cf2.AddEdge( link1, link2, face1, myNbOnPlane2, iNotOnPlane2 ); + + // add intersections to a link collinear with the intersection line + if ( myNbOnPlane1 == 2 && ( link1.myFace != face2 || link2.myFace != face2 )) + cutCollinearLink( iNotOnPlane1, myNodes1, face2, link1, link2 ); + + if ( myNbOnPlane2 == 2 && ( link1.myFace != face1 || link2.myFace != face1 )) + cutCollinearLink( iNotOnPlane2, myNodes2, face1, link1, link2 ); + } + + addLink( link1 ); + addLink( link2 ); + + } // non co-planar case + + return; + } + + //================================================================================ + /*! + * \brief Intersect two 2D line segments + */ + //================================================================================ + + bool Intersector::intersectEdgeEdge( const gp_XY s1p0, const gp_XY s1p1, + const gp_XY s2p0, const gp_XY s2p1, + double & t1, double & t2, + bool & isCollinear ) + { + gp_XY u = s1p1 - s1p0; + gp_XY v = s2p1 - s2p0; + gp_XY w = s1p0 - s2p0; + double perpDotUV = u * gp_XY( -v.Y(), v.X() ); + double perpDotVW = v * gp_XY( -w.Y(), w.X() ); + double perpDotUW = u * gp_XY( -w.Y(), w.X() ); + double u2 = u.SquareModulus(); + double v2 = v.SquareModulus(); + if ( u2 < myEps * myEps || v2 < myEps * myEps ) + return false; + if ( perpDotUV * perpDotUV / u2 / v2 < 1e-6 ) // cos ^ 2 + { + if ( !isCollinear ) + return false; // no need in collinear solution + if ( perpDotUW * perpDotUW / u2 > myTol * myTol ) + return false; // parallel + + // collinear + gp_XY w2 = s1p1 - s2p0; + if ( Abs( v.X()) + Abs( u.X()) > Abs( v.Y()) + Abs( u.Y())) { + t1 = w.X() / v.X(); // params on segment 2 + t2 = w2.X() / v.X(); + } + else { + t1 = w.Y() / v.Y(); + t2 = w2.Y() / v.Y(); + } + if ( Max( t1,t2 ) <= 0 || Min( t1,t2 ) >= 1 ) + return false; // no overlap + return true; + } + isCollinear = false; + + t1 = perpDotVW / perpDotUV; // param on segment 1 + if ( t1 < 0. || t1 > 1. ) + return false; // intersection not within the segment + + t2 = perpDotUW / perpDotUV; // param on segment 2 + if ( t2 < 0. || t2 > 1. ) + return false; // intersection not within the segment + + return true; + } + + //================================================================================ + /*! + * \brief Intersect two edges of co-planar triangles + * \param [inout] iE1 - edge index of triangle 1 + * \param [inout] iE2 - edge index of triangle 2 + * \param [inout] intPoints - intersection points + * \param [inout] nbIntPoints - nb of found intersection points + */ + //================================================================================ + + bool Intersector::intersectEdgeEdge( int iE1, int iE2, IntPoint2D& intPoint ) + { + int i01 = iE1, i11 = ( iE1 + 1 ) % 3; + int i02 = iE2, i12 = ( iE2 + 1 ) % 3; + if (( !intPoint.myIsCollinear ) && + ( myNodes1[ i01 ] == myNodes2[ i02 ] || + myNodes1[ i01 ] == myNodes2[ i12 ] || + myNodes1[ i11 ] == myNodes2[ i02 ] || + myNodes1[ i11 ] == myNodes2[ i12 ] )) + return false; + + // segment 1 + gp_XY s1p0 = p2D( myNodes1[ i01 ]); + gp_XY s1p1 = p2D( myNodes1[ i11 ]); + + // segment 2 + gp_XY s2p0 = p2D( myNodes2[ i02 ]); + gp_XY s2p1 = p2D( myNodes2[ i12 ]); + + double t1, t2; + if ( !intersectEdgeEdge( s1p0,s1p1, s2p0,s2p1, t1, t2, intPoint.myIsCollinear )) + return false; + + intPoint.myEdgeInd[0] = iE1; + intPoint.myEdgeInd[1] = iE2; + intPoint.myU[0] = t1; + intPoint.myU[1] = t2; + (gp_XYZ&)intPoint.myNode = myNodes1[i01] * ( 1 - t1 ) + myNodes1[i11] * t1; + + if ( intPoint.myIsCollinear ) + return true; + + // try to find existing node at intPoint.myNode + + if ( myNodes1[ i01 ] == myNodes2[ i02 ] || + myNodes1[ i01 ] == myNodes2[ i12 ] || + myNodes1[ i11 ] == myNodes2[ i02 ] || + myNodes1[ i11 ] == myNodes2[ i12 ] ) + return false; + + const double coincTol = myTol * 1e-3; + + CutLink link1( myNodes1[i01].Node(), myNodes1[i11].Node(), myFace2 ); + CutLink link2( myNodes2[i02].Node(), myNodes2[i12].Node(), myFace1 ); + + SMESH_NodeXYZ& n1 = myNodes1[ t1 < 0.5 ? i01 : i11 ]; + bool same1 = coincide( n1, intPoint.myNode, coincTol ); + if ( same1 ) + { + link2.myIntNode = intPoint.myNode = n1; + addLink( link2 ); + } + SMESH_NodeXYZ& n2 = myNodes2[ t2 < 0.5 ? i02 : i12 ]; + bool same2 = coincide( n2, intPoint.myNode, coincTol ); + if ( same2 ) + { + link1.myIntNode = intPoint.myNode = n2; + addLink( link1 ); + if ( same1 ) + { + replaceIntNode( n1.Node(), n2.Node() ); + return false; + } + return true; + } + if ( same1 ) + return true; + + link1.myIntNode = intPoint.myNode; + if ( findLink( link1 )) + { + intPoint.myNode = link2.myIntNode = link1.myIntNode; + addLink( link2 ); + return true; + } + + link2.myIntNode = intPoint.myNode; + if ( findLink( link2 )) + { + intPoint.myNode = link1.myIntNode = link2.myIntNode; + addLink( link1 ); + return true; + } + + for ( int is2nd = 0; is2nd < 2; ++is2nd ) + { + const SMDS_MeshElement* f = is2nd ? myFace1 : myFace2; + const CutFace& cf = myCutFaces.Added( CutFace( is2nd ? myFace2 : myFace1 )); + for ( size_t i = 0; i < cf.myLinks.size(); ++i ) + if ( cf.myLinks[i].myFace == f && + //cf.myLinks[i].myIndex != EdgePart::_COPLANAR && + coincide( intPoint.myNode, SMESH_NodeXYZ( cf.myLinks[i].myNode1 ), coincTol )) + { + intPoint.myNode.Set( cf.myLinks[i].myNode1 ); + return true; + } + } + + // make a new node + + intPoint.myNode._node = createNode( intPoint.myNode ); + link1.myIntNode = link2.myIntNode = intPoint.myNode; + addLink( link1 ); + addLink( link2 ); + + return true; + } + + + //================================================================================ + /*! + * \brief Check if a point is contained in a triangle + */ + //================================================================================ + + bool Intersector::isPointInTriangle( const gp_XYZ& p, const std::vector< SMESH_NodeXYZ >& nodes ) + { + double bc1, bc2; + SMESH_MeshAlgos::GetBarycentricCoords( p2D( p ), + p2D( nodes[0] ), p2D( nodes[1] ), p2D( nodes[2] ), + bc1, bc2 ); + return ( 0. < bc1 && 0. < bc2 && bc1 + bc2 < 1. ); + } + + //================================================================================ + /*! + * \brief Intersect two co-planar faces + */ + //================================================================================ + + void Intersector::cutCoplanar() + { + // find intersections of edges + + IntPoint2D intPoints[ 6 ]; + int nbIntPoints = 0; + for ( int iE1 = 0; iE1 < 3; ++iE1 ) + { + int maxNbIntPoints = nbIntPoints + 2; + for ( int iE2 = 0; iE2 < 3 && nbIntPoints < maxNbIntPoints; ++iE2 ) + nbIntPoints += intersectEdgeEdge( iE1, iE2, intPoints[ nbIntPoints ]); + } + const int minNbOnPlane = Min( myNbOnPlane1, myNbOnPlane2 ); + + if ( nbIntPoints == 0 ) // no intersections of edges + { + bool is1in2; + if ( isPointInTriangle( myNodes1[0], myNodes2 )) // face2 includes face1 + is1in2 = true; + else if ( isPointInTriangle( myNodes2[0], myNodes1 )) // face1 includes face2 + is1in2 = false; + else + return; + + // add edges of an inner triangle to an outer one + + const std::vector< SMESH_NodeXYZ >& nodesIn = is1in2 ? myNodes1 : myNodes2; + const SMDS_MeshElement* faceOut = is1in2 ? myFace2 : myFace1; + const SMDS_MeshElement* faceIn = is1in2 ? myFace1 : myFace2; + + const CutFace& outFace = myCutFaces.Added( CutFace( faceOut )); + CutLink link1( nodesIn.back().Node(), nodesIn.back().Node(), faceOut ); + CutLink link2( nodesIn.back().Node(), nodesIn.back().Node(), faceOut ); + + link1.myIntNode = nodesIn.back(); + for ( size_t i = 0; i < nodesIn.size(); ++i ) + { + link2.myIntNode = nodesIn[ i ]; + outFace.AddEdge( link1, link2, faceIn, minNbOnPlane ); + link1.myIntNode = link2.myIntNode; + } + } + else + { + // add parts of edges to a triangle including them + + CutLink link1, link2; + IntPoint2D ip0, ip1; + ip0.myU[0] = ip0.myU[1] = 0.; + ip1.myU[0] = ip1.myU[1] = 1.; + ip0.myEdgeInd[0] = ip0.myEdgeInd[1] = ip1.myEdgeInd[0] = ip1.myEdgeInd[1] = 0; + + for ( int isFromFace1 = 0; isFromFace1 < 2; ++isFromFace1 ) + { + const SMDS_MeshElement* faceTo = isFromFace1 ? myFace2 : myFace1; + const SMDS_MeshElement* faceFrom = isFromFace1 ? myFace1 : myFace2; + const std::vector< SMESH_NodeXYZ >& nodesTo = isFromFace1 ? myNodes2 : myNodes1; + const std::vector< SMESH_NodeXYZ >& nodesFrom = isFromFace1 ? myNodes1 : myNodes2; + const int iTo = isFromFace1 ? 1 : 0; + const int iFrom = isFromFace1 ? 0 : 1; + //const int nbOnPlaneFrom = isFromFace1 ? myNbOnPlane1 : myNbOnPlane2; + + const CutFace* cutFaceTo = & myCutFaces.Added( CutFace( faceTo )); + // const CutFace* cutFaceFrom = 0; + // if ( nbOnPlaneFrom > minNbOnPlane ) + // cutFaceFrom = & myCutFaces.Added( CutFace( faceTo )); + + link1.myFace = link2.myFace = faceTo; + + IntPoint2DCompare ipCompare( iFrom ); + TIntPointPtrSet pointsOnEdge( ipCompare ); // IntPoint2D sorted by parameter on edge + + for ( size_t iE = 0; iE < nodesFrom.size(); ++iE ) + { + // get parts of an edge iE + + ip0.myEdgeInd[ iTo ] = iE; + ip1.myEdgeInd[ iTo ] = ( iE + 1 ) % nodesFrom.size(); + ip0.myNode = nodesFrom[ ip0.myEdgeInd[ iTo ]]; + ip1.myNode = nodesFrom[ ip1.myEdgeInd[ iTo ]]; + + pointsOnEdge.clear(); + + for ( int iP = 0; iP < nbIntPoints; ++iP ) + if ( intPoints[ iP ].myEdgeInd[ iFrom ] == iE ) + pointsOnEdge.insert( & intPoints[ iP ] ); + + pointsOnEdge.insert( pointsOnEdge.begin(), & ip0 ); + pointsOnEdge.insert( pointsOnEdge.end(), & ip1 ); + + // add edge parts to faceTo + + TIntPointPtrSet::iterator ipIt = pointsOnEdge.begin() + 1; + for ( ; ipIt != pointsOnEdge.end(); ++ipIt ) + { + const IntPoint2D* p1 = *(ipIt-1); + const IntPoint2D* p2 = *ipIt; + gp_XYZ middle = 0.5 * ( p1->myNode + p2->myNode ); + if ( isPointInTriangle( middle, nodesTo )) + { + p1->InitLink( link1, iTo, ( p1 != & ip0 ) ? nodesTo : nodesFrom ); + p2->InitLink( link2, iTo, ( p2 != & ip1 ) ? nodesTo : nodesFrom ); + cutFaceTo->AddEdge( link1, link2, faceFrom, minNbOnPlane ); + + // if ( cutFaceFrom ) + // { + // p1->InitLink( link1, iFrom, nodesFrom ); + // p2->InitLink( link2, iFrom, nodesFrom ); + // cutFaceTo->AddEdge( link1, link2, faceTo, minNbOnPlane ); + // } + } + } + } + } + } + return; + + } // Intersector::cutCoplanar() + + //================================================================================ + /*! + * \brief Intersect edges added to myCutFaces + */ + //================================================================================ + + void Intersector::intersectNewEdges( const CutFace& cf ) + { + IntPoint2D intPoint; + + if ( cf.NbInternalEdges() < 2 ) + return; + + const gp_XYZ& faceNorm = myNormals[ cf.myInitFace->GetID() ]; + setPlaneIndices( faceNorm ); // choose indices on an axis-aligned plane + + size_t limit = cf.myLinks.size() * cf.myLinks.size() * 2; + + for ( size_t i1 = 3; i1 < cf.myLinks.size(); ++i1 ) + { + if ( !cf.myLinks[i1].IsInternal() ) + continue; + + myIntPointSet.clear(); + for ( size_t i2 = i1 + 2; i2 < cf.myLinks.size(); ++i2 ) + { + if ( !cf.myLinks[i2].IsInternal() ) + continue; + + // prepare to intersection + myFace1 = cf.myLinks[i1].myFace; + myNodes1[0] = cf.myLinks[i1].myNode1; + myNodes1[1] = cf.myLinks[i1].myNode2; + myFace2 = cf.myLinks[i2].myFace; + myNodes2[0] = cf.myLinks[i2].myNode1; + myNodes2[1] = cf.myLinks[i2].myNode2; + + // intersect + intPoint.myIsCollinear = true; // to find collinear solutions + if ( intersectEdgeEdge( 0, 0, intPoint )) + { + if ( cf.myLinks[i1].IsSame( cf.myLinks[i2] )) // remove i2 + { + cf.myLinks[i1].ReplaceCoplanar( cf.myLinks[i2] ); + cf.myLinks.erase( cf.myLinks.begin() + i2, cf.myLinks.begin() + i2 + 2 ); + --i2; + continue; + } + if ( !intPoint.myIsCollinear ) + { + intPoint.myEdgeInd[1] = i2; + myIntPointSet.insert( intPoint ); + } + else // if ( intPoint.myIsCollinear ) // overlapping edges + { + myIntPointSet.clear(); // to recompute + + if ( intPoint.myU[0] > intPoint.myU[1] ) // orient in same direction + { + std::swap( intPoint.myU[0], intPoint.myU[1] ); + std::swap( myNodes1[0], myNodes1[1] ); + } + // replace _COPLANAR by _INTERNAL + cf.myLinks[i1].ReplaceCoplanar( cf.myLinks[i1+1] ); + cf.myLinks[i2].ReplaceCoplanar( cf.myLinks[i2+1] ); + + if ( coincide( myNodes1[0], myNodes2[0], myTol ) && + coincide( myNodes1[1], myNodes2[1], myTol )) + { + cf.myLinks.erase( cf.myLinks.begin() + i2, cf.myLinks.begin() + i2 + 2 ); + --i2; + continue; + } + + EdgePart common = cf.myLinks[i1]; + common.ReplaceCoplanar( cf.myLinks[i2] ); + + const SMDS_MeshNode* n1 = myNodes1[0].Node(); // end nodes of an overlapping part + const SMDS_MeshNode* n2 = myNodes1[1].Node(); + size_t i3 = cf.myLinks.size(); + + if ( myNodes1[0] != myNodes2[0] ) // a part before the overlapping one + { + if ( intPoint.myU[0] < 0 ) + cf.myLinks[i1].Set( myNodes1[0].Node(), myNodes2[0].Node(), + cf.myLinks[i1].myFace, cf.myLinks[i1].myIndex ); + else + cf.myLinks[i1].Set( myNodes2[0].Node(), myNodes1[0].Node(), + cf.myLinks[i2].myFace, cf.myLinks[i2].myIndex ); + + cf.myLinks[i1+1].Set( cf.myLinks[i1].myNode2, + cf.myLinks[i1].myNode1, + cf.myLinks[i1].myFace, + cf.myLinks[i1].myIndex); + n1 = cf.myLinks[i1].myNode2; + } + else + i3 = i1; + + if ( myNodes1[1] != myNodes2[1] ) // a part after the overlapping one + { + if ( intPoint.myU[1] < 1 ) + cf.myLinks[i2].Set( myNodes1[1].Node(), myNodes2[1].Node(), + cf.myLinks[i2].myFace, cf.myLinks[i2].myIndex ); + else + cf.myLinks[i2].Set( myNodes2[1].Node(), myNodes1[1].Node(), + cf.myLinks[i1].myFace, cf.myLinks[i1].myIndex ); + + cf.myLinks[i2+1].Set( cf.myLinks[i2].myNode2, + cf.myLinks[i2].myNode1, + cf.myLinks[i2].myFace, + cf.myLinks[i2].myIndex); + n2 = cf.myLinks[i2].myNode1; + } + else + i3 = i2; + + if ( i3 == cf.myLinks.size() ) + cf.myLinks.resize( i3 + 2 ); + + cf.myLinks[i3].Set ( n1, n2, common.myFace, common.myIndex ); + cf.myLinks[i3+1].Set( n2, n1, common.myFace, common.myIndex ); + + i2 = i1 + 1; // recheck modified i1 + continue; + } + //else + // { + // // remember a new node + // CutLink link1( myNodes1[0].Node(), myNodes1[1].Node(), cf.myInitFace ); + // CutLink link2( myNodes2[0].Node(), myNodes2[1].Node(), cf.myInitFace ); + // link2.myIntNode = link1.myIntNode = intPoint.myNode; + // addLink( link1 ); + // addLink( link2 ); + + // // split edges + // size_t i = cf.myLinks.size(); + // if ( intPoint.myNode != cf.myLinks[ i1 ].myNode1 && + // intPoint.myNode != cf.myLinks[ i1 ].myNode2 ) + // { + // cf.myLinks.push_back( cf.myLinks[ i1 ]); + // cf.myLinks.push_back( cf.myLinks[ i1 + 1 ]); + // cf.myLinks[ i1 ].myNode2 = cf.myLinks[ i1 + 1 ].myNode1 = intPoint.Node(); + // cf.myLinks[ i ].myNode1 = cf.myLinks[ i + 1 ].myNode2 = intPoint.Node(); + // } + // if ( intPoint.myNode != cf.myLinks[ i2 ].myNode1 && + // intPoint.myNode != cf.myLinks[ i2 ].myNode2 ) + // { + // i = cf.myLinks.size(); + // cf.myLinks.push_back( cf.myLinks[ i2 ]); + // cf.myLinks.push_back( cf.myLinks[ i2 + 1 ]); + // cf.myLinks[ i2 ].myNode2 = cf.myLinks[ i2 + 1 ].myNode1 = intPoint.Node(); + // cf.myLinks[ i ].myNode1 = cf.myLinks[ i + 1 ].myNode2 = intPoint.Node(); + // } + // } + + } // if ( intersectEdgeEdge( 0, 0, intPoint )) + + ++i2; + --limit; + } + + // split i1 edge and all edges it intersects + // don't do it inside intersection loop in order not to loose direction of i1 edge + if ( !myIntPointSet.empty() ) + { + cf.myLinks.reserve( cf.myLinks.size() + myIntPointSet.size() * 2 + 2 ); + + EdgePart* edge1 = &cf.myLinks[ i1 ]; + EdgePart* twin1 = &cf.myLinks[ i1 + 1 ]; + + TIntPointSet::iterator ipIt = myIntPointSet.begin(); + for ( ; ipIt != myIntPointSet.end(); ++ipIt ) // int points sorted on i1 edge + { + size_t i = cf.myLinks.size(); + if ( ipIt->myNode != edge1->myNode1 && + ipIt->myNode != edge1->myNode2 ) + { + cf.myLinks.push_back( *edge1 ); + cf.myLinks.push_back( *twin1 ); + edge1->myNode2 = twin1->myNode1 = ipIt->Node(); + cf.myLinks[ i ].myNode1 = cf.myLinks[ i + 1 ].myNode2 = ipIt->Node(); + edge1 = & cf.myLinks[ i ]; + twin1 = & cf.myLinks[ i + 1 ]; + } + size_t i2 = ipIt->myEdgeInd[1]; + if ( ipIt->myNode != cf.myLinks[ i2 ].myNode1 && + ipIt->myNode != cf.myLinks[ i2 ].myNode2 ) + { + i = cf.myLinks.size(); + cf.myLinks.push_back( cf.myLinks[ i2 ]); + cf.myLinks.push_back( cf.myLinks[ i2 + 1 ]); + cf.myLinks[ i2 ].myNode2 = cf.myLinks[ i2 + 1 ].myNode1 = ipIt->Node(); + cf.myLinks[ i ].myNode1 = cf.myLinks[ i + 1 ].myNode2 = ipIt->Node(); + } + } + if ( cf.myLinks.size() >= limit ) + throw SALOME_Exception( "Infinite loop in Intersector::intersectNewEdges()" ); + } + ++i1; // each internal edge encounters twice + } + return; + } + + //================================================================================ + /*! + * \brief Split intersected faces + */ + //================================================================================ + + void Intersector::MakeNewFaces( SMESH_MeshAlgos::TEPairVec& theNew2OldFaces, + SMESH_MeshAlgos::TNPairVec& theNew2OldNodes, + const double theSign) + { + // unmark all nodes except intersection ones + + for ( SMDS_NodeIteratorPtr nIt = myMesh->nodesIterator(); nIt->more(); ) + { + const SMDS_MeshNode* n = nIt->next(); + if ( n->isMarked() && n->GetID()-1 < (int) theNew2OldNodes.size() ) + n->setIsMarked( false ); + } + // SMESH_MeshAlgos::MarkElems( myMesh->nodesIterator(), false ); + + TCutLinkMap::const_iterator cutLinksIt = myCutLinks.cbegin(); + // for ( ; cutLinksIt != myCutLinks.cend(); ++cutLinksIt ) + // { + // const CutLink& link = *cutLinksIt; + // if ( link.IntNode() && link.IntNode()->GetID()-1 < (int) theNew2OldNodes.size() ) + // link.IntNode()->setIsMarked( true ); + // } + + // intersect edges added to myCutFaces + + TCutFaceMap::const_iterator cutFacesIt = myCutFaces.cbegin(); + for ( ; cutFacesIt != myCutFaces.cend(); ++cutFacesIt ) + { + const CutFace& cf = *cutFacesIt; + cf.ReplaceNodes( myRemove2KeepNodes ); + intersectNewEdges( cf ); + } + + // make new faces + + EdgeLoopSet loopSet; + SMESH_MeshAlgos::Triangulate triangulator; + std::vector< EdgePart > cutOffLinks; + TLinkMap cutOffCoplanarLinks; + std::vector< const CutFace* > touchedFaces; + SMESH_MeshAlgos::TEPairVec::value_type new2OldTria; + CutFace cutFace(0); + std::vector< const SMDS_MeshNode* > nodes; + std::vector faces; + + cutOffLinks.reserve( myCutFaces.Extent() * 2 ); + + for ( cutFacesIt = myCutFaces.cbegin(); cutFacesIt != myCutFaces.cend(); ++cutFacesIt ) + { + const CutFace& cf = *cutFacesIt; + if ( !cf.IsCut() ) + { + touchedFaces.push_back( & cf ); + continue; + } + + const gp_XYZ& normal = myNormals[ cf.myInitFace->GetID() ]; + + // form loops of new faces + cf.ReplaceNodes( myRemove2KeepNodes ); + cf.MakeLoops( loopSet, normal ); + + // avoid loops that are not connected to boundary edges of cf.myInitFace + if ( cf.RemoveInternalLoops( loopSet )) + { + intersectNewEdges( cf ); + cf.MakeLoops( loopSet, normal ); + } + // erase loops that are cut off by face intersections + cf.CutOffLoops( loopSet, theSign, myNormals, cutOffLinks, cutOffCoplanarLinks ); + + int index = cf.myInitFace->GetID(); // index in theNew2OldFaces + + const SMDS_MeshElement* tria; + for ( size_t iL = 0; iL < loopSet.myNbLoops; ++iL ) + { + EdgeLoop& loop = loopSet.myLoops[ iL ]; + if ( loop.myLinks.size() == 0 ) + continue; + + int nbTria = triangulator.GetTriangles( &loop, nodes ); + int nbNodes = 3 * nbTria; + for ( int i = 0; i < nbNodes; i += 3 ) + { + if ( nodes[i] == nodes[i+1] || nodes[i] == nodes[i+2] || nodes[i+1] == nodes[i+2] ) + { +#ifdef _DEBUG_ + std::cerr << "BAD tria" << std::endl; + cf.Dump(); +#endif + continue; + } + if (!( tria = myMesh->FindFace( nodes[i], nodes[i+1], nodes[i+2] ))) + tria = myMesh->AddFace( nodes[i], nodes[i+1], nodes[i+2] ); + tria->setIsMarked( true ); // not to remove it + + new2OldTria = std::make_pair( tria, theNew2OldFaces[ index ].second ); + if ( tria->GetID() < (int)theNew2OldFaces.size() ) + theNew2OldFaces[ tria->GetID() ] = new2OldTria; + else + theNew2OldFaces.push_back( new2OldTria ); + + if ( index == tria->GetID() ) + index = 0; // do not remove tria + } + } + theNew2OldFaces[ index ].first = 0; + } + + // remove split faces + for ( size_t id = 1; id < theNew2OldFaces.size(); ++id ) + { + if ( theNew2OldFaces[id].first ) + continue; + if ( const SMDS_MeshElement* f = myMesh->FindElement( id )) + myMesh->RemoveFreeElement( f ); + } + + // remove face connected to cut off parts of cf.myInitFace + + nodes.resize(2); + for ( size_t i = 0; i < cutOffLinks.size(); ++i ) + { + //break; + nodes[0] = cutOffLinks[i].myNode1; + nodes[1] = cutOffLinks[i].myNode2; + + if ( nodes[0] != nodes[1] && + myMesh->GetElementsByNodes( nodes, faces )) + { + if ( cutOffLinks[i].myFace && + cutOffLinks[i].myIndex != EdgePart::_COPLANAR && + faces.size() == 2 ) + continue; + for ( size_t iF = 0; iF < faces.size(); ++iF ) + { + int index = faces[iF]->GetID(); + // if ( //faces[iF]->isMarked() || // kept part of cutFace + // !theNew2OldFaces[ index ].first ) // already removed + // continue; + cutFace.myInitFace = faces[iF]; + // if ( myCutFaces.Contains( cutFace )) // keep cutting faces needed in CutOffLoops() + // { + // if ( !myCutFaces.Added( cutFace ).IsCut() ) + // theNew2OldFaces[ index ].first = 0; + // continue; + // } + cutFace.myLinks.clear(); + cutFace.InitLinks(); + for ( size_t iL = 0; iL < cutFace.myLinks.size(); ++iL ) + if ( !cutOffLinks[i].IsSame( cutFace.myLinks[ iL ])) + cutOffLinks.push_back( cutFace.myLinks[ iL ]); + + theNew2OldFaces[ index ].first = 0; + myMesh->RemoveFreeElement( faces[iF] ); + } + } + } + + // replace nodes in touched faces + + // treat touched faces + for ( size_t i = 0; i < touchedFaces.size(); ++i ) + { + const CutFace& cf = *touchedFaces[i]; + + int index = cf.myInitFace->GetID(); // index in theNew2OldFaces + if ( !theNew2OldFaces[ index ].first ) + continue; // already cut off + + if ( !cf.ReplaceNodes( myRemove2KeepNodes )) + continue; // just keep as is + + if ( cf.myLinks.size() == 3 ) + { + const SMDS_MeshElement* tria = myMesh->AddFace( cf.myLinks[0].myNode1, + cf.myLinks[1].myNode1, + cf.myLinks[2].myNode1 ); + new2OldTria = std::make_pair( tria, theNew2OldFaces[ index ].second ); + if ( tria->GetID() < (int)theNew2OldFaces.size() ) + theNew2OldFaces[ tria->GetID() ] = new2OldTria; + else + theNew2OldFaces.push_back( new2OldTria ); + } + theNew2OldFaces[ index ].first = 0; + } + + + // add used new nodes to theNew2OldNodes + SMESH_MeshAlgos::TNPairVec::value_type new2OldNode; + new2OldNode.second = NULL; + for ( cutLinksIt = myCutLinks.cbegin(); cutLinksIt != myCutLinks.cend(); ++cutLinksIt ) + { + const CutLink& link = *cutLinksIt; + if ( link.IntNode() ) // && link.IntNode()->NbInverseElements() > 0 ) + { + new2OldNode.first = link.IntNode(); + theNew2OldNodes.push_back( new2OldNode ); + } + } + + return; + } + + //================================================================================ + /*! + * \brief Debug + */ + //================================================================================ + + void CutFace::Dump() const + { + std::cout << std::endl << "INI F " << myInitFace->GetID() << std::endl; + for ( size_t i = 0; i < myLinks.size(); ++i ) + std::cout << "[" << i << "] (" + << char(( myLinks[i].IsInternal() ? 'j' : '0' ) + myLinks[i].myIndex ) << ") " + << myLinks[i].myNode1->GetID() << " - " << myLinks[i].myNode2->GetID() + << " " << ( myLinks[i].myFace ? 'F' : 'C' ) + << ( myLinks[i].myFace ? myLinks[i].myFace->GetID() : 0 ) << " " << std::endl; + } + + //================================================================================ + /*! + * \brief Add an edge cutting this face + * \param [in] p1 - start point of the edge + * \param [in] p2 - end point of the edge + * \param [in] cutter - a face producing the added cut edge. + * \param [in] nbOnPlane - nb of triangle nodes lying on the plane of the cutter face + */ + //================================================================================ + + void CutFace::AddEdge( const CutLink& p1, + const CutLink& p2, + const SMDS_MeshElement* cutterFace, + const int nbOnPlane, + const int iNotOnPlane) const + { + int iN[2] = { myInitFace->GetNodeIndex( p1.IntNode() ), + myInitFace->GetNodeIndex( p2.IntNode() ) }; + if ( iN[0] >= 0 && iN[1] >= 0 ) + { + // the cutting edge is a whole side + if (( cutterFace && nbOnPlane < 3 ) && + !( cutterFace->GetNodeIndex( p1.IntNode() ) >= 0 && + cutterFace->GetNodeIndex( p2.IntNode() ) >= 0 )) + { + InitLinks(); + myLinks[ Abs( iN[0] - iN[1] ) == 1 ? Min( iN[0], iN[1] ) : 2 ].myFace = cutterFace; + } + return; + } + + if ( p1.IntNode() == p2.IntNode() ) + { + AddPoint( p1, p2, 1e-10 ); + return; + } + + InitLinks(); + + // cut side edges by a new one + + int iEOnPlane = ( nbOnPlane == 2 ) ? ( iNotOnPlane + 1 ) % 3 : -1; + + double dist[2]; + for ( int is2nd = 0; is2nd < 2; ++is2nd ) + { + const CutLink& p = is2nd ? p2 : p1; + dist[ is2nd ] = 0; + if ( iN[ is2nd ] >= 0 ) + continue; + + int iE = Max( iEOnPlane, myInitFace->GetNodeIndex( p.Node1() )); + if ( iE < 0 ) + continue; // link of other face + + SMESH_NodeXYZ n0 = myLinks[iE].myNode1; + dist[ is2nd ] = ( n0 - p.myIntNode ).SquareModulus(); + + for ( size_t i = 0; i < myLinks.size(); ++i ) + if ( myLinks[i].myIndex == iE ) + { + double d1 = n0.SquareDistance( myLinks[i].myNode1 ); + if ( d1 < dist[ is2nd ] ) + { + double d2 = n0.SquareDistance( myLinks[i].myNode2 ); + if ( dist[ is2nd ] < d2 ) + { + myLinks.push_back( myLinks[i] ); + myLinks.back().myNode1 = myLinks[i].myNode2 = p.IntNode(); + break; + } + } + } + } + + int state = nbOnPlane == 3 ? EdgePart::_COPLANAR : EdgePart::_INTERNAL; + + // look for an existing equal edge + if ( nbOnPlane == 2 ) + { + SMESH_NodeXYZ n0 = myLinks[ iEOnPlane ].myNode1; + if ( iN[0] >= 0 ) dist[0] = ( n0 - p1.myIntNode ).SquareModulus(); + if ( iN[1] >= 0 ) dist[1] = ( n0 - p2.myIntNode ).SquareModulus(); + if ( dist[0] > dist[1] ) + std::swap( dist[0], dist[1] ); + for ( size_t i = 0; i < myLinks.size(); ++i ) + { + if ( myLinks[i].myIndex != iEOnPlane ) + continue; + gp_XYZ mid = 0.5 * ( SMESH_NodeXYZ( myLinks[i].myNode1 ) + + SMESH_NodeXYZ( myLinks[i].myNode2 )); + double d = ( n0 - mid ).SquareModulus(); + if ( dist[0] < d && d < dist[1] ) + myLinks[i].myFace = cutterFace; + } + return; + } + else + { + EdgePart newEdge; newEdge.Set( p1.IntNode(), p2.IntNode(), cutterFace, state ); + for ( size_t i = 0; i < myLinks.size(); ++i ) + { + if ( myLinks[i].IsSame( newEdge )) + { + // if ( !myLinks[i].IsInternal() ) + // myLinks[ i ].myFace = cutterFace; + // else + myLinks[ i ].ReplaceCoplanar( newEdge ); + myLinks[ i+1 ].ReplaceCoplanar( newEdge ); + return; + } + i += myLinks[i].IsInternal(); + } + } + + size_t i = myLinks.size(); + myLinks.resize( i + 2 ); + myLinks[ i ].Set( p1.IntNode(), p2.IntNode(), cutterFace, state ); + myLinks[ i+1 ].Set( p2.IntNode(), p1.IntNode(), cutterFace, state ); + } + + //================================================================================ + /*! + * \brief Add a point cutting this face + */ + //================================================================================ + + void CutFace::AddPoint( const CutLink& p1, const CutLink& p2, double tol ) const + { + if ( myInitFace->GetNodeIndex( p1.IntNode() ) >= 0 || + myInitFace->GetNodeIndex( p2.IntNode() ) >= 0 ) + return; + + InitLinks(); + + const CutLink* link = &p1; + int iE = myInitFace->GetNodeIndex( link->Node1() ); + if ( iE < 0 ) + { + link = &p2; + iE = myInitFace->GetNodeIndex( link->Node1() ); + } + if ( iE >= 0 ) + { + // cut an existing edge by the point + SMESH_NodeXYZ n0 = link->Node1(); + double d = ( n0 - link->myIntNode ).SquareModulus(); + + for ( size_t i = 0; i < myLinks.size(); ++i ) + if ( myLinks[i].myIndex == iE ) + { + double d1 = n0.SquareDistance( myLinks[i].myNode1 ); + if ( d1 < d ) + { + double d2 = n0.SquareDistance( myLinks[i].myNode2 ); + if ( d < d2 ) + { + myLinks.push_back( myLinks[i] ); + myLinks.back().myNode1 = myLinks[i].myNode2 = link->IntNode(); + return; + } + } + } + } + else // point is inside the triangle + { + // // check if a point already added + // for ( size_t i = 3; i < myLinks.size(); ++i ) + // if ( myLinks[i].myNode1 == p1.IntNode() ) + // return; + + // // create a link between the point and the closest corner node + // const SMDS_MeshNode* closeNode = myLinks[0].myNode1; + // double minDist = p1.myIntNode.SquareDistance( closeNode ); + // for ( int i = 1; i < 3; ++i ) + // { + // double dist = p1.myIntNode.SquareDistance( myLinks[i].myNode1 ); + // if ( dist < minDist ) + // { + // minDist = dist; + // closeNode = myLinks[i].myNode1; + // } + // } + // if ( minDist > tol * tol ) + // { + // size_t i = myLinks.size(); + // myLinks.resize( i + 2 ); + // myLinks[ i ].Set( p1.IntNode(), closeNode ); + // myLinks[ i+1 ].Set( closeNode, p1.IntNode() ); + // } + } + } + + //================================================================================ + /*! + * \brief Perform node replacement + */ + //================================================================================ + + bool CutFace::ReplaceNodes( const TNNMap& theRm2KeepMap ) const + { + bool replaced = false; + for ( size_t i = 0; i < myLinks.size(); ++i ) + { + while ( theRm2KeepMap.IsBound((Standard_Address) myLinks[i].myNode1 )) + replaced = ( myLinks[i].myNode1 = theRm2KeepMap((Standard_Address) myLinks[i].myNode1 )); + + while ( theRm2KeepMap.IsBound((Standard_Address) myLinks[i].myNode2 )) + replaced = ( myLinks[i].myNode2 = theRm2KeepMap((Standard_Address) myLinks[i].myNode2 )); + } + + //if ( replaced ) // remove equal links + { + for ( size_t i1 = 0; i1 < myLinks.size(); ++i1 ) + { + if ( myLinks[i1].myNode1 == myLinks[i1].myNode2 ) + { + myLinks.erase( myLinks.begin() + i1, + myLinks.begin() + i1 + 1 + myLinks[i1].IsInternal() ); + --i1; + continue; + } + size_t i2 = i1 + 1 + myLinks[i1].IsInternal(); + for ( ; i2 < myLinks.size(); ++i2 ) + { + if ( !myLinks[i2].IsInternal() ) + continue; + if ( myLinks[i1].IsSame( myLinks[i2] )) + { + myLinks[i1]. ReplaceCoplanar( myLinks[i2] ); + if ( myLinks[i1].IsInternal() ) + myLinks[i1+1].ReplaceCoplanar( myLinks[i2+1] ); + if ( !myLinks[i1].myFace && myLinks[i2].myFace ) + { + myLinks[i1]. myFace = myLinks[i2].myFace; + if ( myLinks[i1].IsInternal() ) + myLinks[i1+1].myFace = myLinks[i2+1].myFace; + } + myLinks.erase( myLinks.begin() + i2, + myLinks.begin() + i2 + 2 ); + --i2; + continue; + } + ++i2; + } + i1 += myLinks[i1].IsInternal(); + } + } + + return replaced; + } + + //================================================================================ + /*! + * \brief Initialize myLinks with edges of myInitFace + */ + //================================================================================ + + void CutFace::InitLinks() const + { + if ( !myLinks.empty() ) return; + + int nbNodes = myInitFace->NbNodes(); + myLinks.reserve( nbNodes * 2 ); + myLinks.resize( nbNodes ); + + for ( int i = 0; i < nbNodes; ++i ) + { + const SMDS_MeshNode* n1 = myInitFace->GetNode( i ); + const SMDS_MeshNode* n2 = myInitFace->GetNodeWrap( i + 1); + myLinks[i].Set( n1, n2, 0, i ); + } + } + + //================================================================================ + /*! + * \brief Return number of internal edges + */ + //================================================================================ + + int CutFace::NbInternalEdges() const + { + int nb = 0; + for ( size_t i = 3; i < myLinks.size(); ++i ) + nb += myLinks[i].IsInternal(); + + return nb / 2; // each internal edge encounters twice + } + + //================================================================================ + /*! + * \brief Remove loops that are not connected to boundary edges of myFace by + * adding edges connecting these loops to the boundary + */ + //================================================================================ + + bool CutFace::RemoveInternalLoops( EdgeLoopSet& theLoops ) const + { + size_t nbReachedLoops = 0; + + // count loops including boundary EdgeParts + for ( size_t iL = 0; iL < theLoops.myNbLoops; ++iL ) + { + EdgeLoop& loop = theLoops.myLoops[ iL ]; + + for ( size_t iE = 0; iE < loop.myLinks.size(); ++iE ) + if ( !loop.myLinks[ iE ]->IsInternal() ) + { + nbReachedLoops += loop.SetConnected(); + break; + } + } + if ( nbReachedLoops == theLoops.myNbLoops ) + return false; // no unreachable loops + + + // try to reach all loops by propagating via internal edges shared by loops + size_t prevNbReached; + do + { + prevNbReached = nbReachedLoops; + + for ( size_t iL = 0; iL < theLoops.myNbLoops; ++iL ) + { + EdgeLoop& loop = theLoops.myLoops[ iL ]; + if ( !loop.myIsBndConnected ) + continue; + + for ( size_t iE = 0; iE < loop.myLinks.size(); ++iE ) + if ( loop.myLinks[ iE ]->IsInternal() ) + { + const EdgePart* twinEdge = getTwin( loop.myLinks[ iE ]); + EdgeLoop* loop2 = theLoops.GetLoopOf( twinEdge ); + if ( loop2->SetConnected() && ++nbReachedLoops == theLoops.myNbLoops ) + return false; // no unreachable loops + } + } + } + while ( prevNbReached < nbReachedLoops ); + + + // add links connecting internal loops with the boundary ones + + for ( size_t iL = 0; iL < theLoops.myNbLoops; ++iL ) + { + EdgeLoop& loop = theLoops.myLoops[ iL ]; + if ( loop.myIsBndConnected ) + continue; + + // find a pair of closest nodes + const SMDS_MeshNode *closestNode1, *closestNode2; + double minDist = 1e100; + for ( size_t iE = 0; iE < loop.myLinks.size(); ++iE ) + { + SMESH_NodeXYZ n1 = loop.myLinks[ iE ]->myNode1; + + for ( size_t i = 0; i < myLinks.size(); ++i ) + { + if ( !loop.Contains( myLinks[i].myNode1 )) + { + double dist = n1.SquareDistance( myLinks[i].myNode1 ); + if ( dist < minDist ) + { + minDist = dist; + closestNode1 = loop.myLinks[ iE ]->myNode1; + closestNode2 = myLinks[i].myNode1; + } + } + if ( myLinks[i].IsInternal() ) + ++i; + } + } + + size_t i = myLinks.size(); + myLinks.resize( i + 2 ); + myLinks[ i ].Set( closestNode1, closestNode2 ); + myLinks[ i+1 ].Set( closestNode2, closestNode1 ); + } + + return true; + } + + //================================================================================ + /*! + * \brief Return equal reversed edge + */ + //================================================================================ + + EdgePart* CutFace::getTwin( const EdgePart* edge ) const + { + size_t i = edge - & myLinks[0]; + + if ( i > 2 && myLinks[ i-1 ].IsTwin( *edge )) + return & myLinks[ i-1 ]; + + if ( i+1 < myLinks.size() && + myLinks[ i+1 ].IsTwin( *edge )) + return & myLinks[ i+1 ]; + + return 0; + } + + //================================================================================ + /*! + * \brief Fill loops of edges + */ + //================================================================================ + + void CutFace::MakeLoops( EdgeLoopSet& theLoops, const gp_XYZ& theFaceNorm ) const + { + theLoops.Init( myLinks ); + + if ( myLinks.size() == 3 ) + { + theLoops.AddNewLoop(); + theLoops.AddEdge( myLinks[0] ); + theLoops.AddEdge( myLinks[1] ); + theLoops.AddEdge( myLinks[2] ); + return; + } + + while ( !theLoops.AllEdgesUsed() ) + { + theLoops.AddNewLoop(); + + // add 1st edge to a new loop + size_t i1; + for ( i1 = theLoops.myNbLoops - 1; i1 < myLinks.size(); ++i1 ) + if ( theLoops.AddEdge( myLinks[i1] )) + break; + + EdgePart* lastEdge = & myLinks[ i1 ]; + EdgePart* twinEdge = getTwin( lastEdge ); + const SMDS_MeshNode* firstNode = lastEdge->myNode1; + const SMDS_MeshNode* lastNode = lastEdge->myNode2; + + do // add the rest edges + { + theLoops.myCandidates.clear(); // edges starting at lastNode + int nbInternal = 0; + + // find candidate edges + for ( size_t i = i1 + 1; i < myLinks.size(); ++i ) + if ( myLinks[ i ].myNode1 == lastNode && + &myLinks[ i ] != twinEdge && + !theLoops.myIsUsedEdge[ i ]) + { + theLoops.myCandidates.push_back( & myLinks[ i ]); + nbInternal += myLinks[ i ].IsInternal(); + } + + // choose among candidates + if ( theLoops.myCandidates.size() == 0 ) + { + theLoops.GetLoopOf( lastEdge )->myHasPending = true; + lastEdge = twinEdge; + } + else if ( theLoops.myCandidates.size() == 1 ) + { + lastEdge = theLoops.myCandidates[0]; + } + else if ( nbInternal == 1 && !lastEdge->IsInternal() ) + { + lastEdge = theLoops.myCandidates[ !theLoops.myCandidates[0]->IsInternal() ]; + } + else + { + gp_Vec lastVec = *lastEdge; + double maxAngle = -2 * M_PI; + for ( size_t i = 0; i < theLoops.myCandidates.size(); ++i ) + { + double angle = lastVec.AngleWithRef( *theLoops.myCandidates[i], theFaceNorm ); + if ( angle > maxAngle ) + { + maxAngle = angle; + lastEdge = theLoops.myCandidates[i]; + } + } + } + theLoops.AddEdge( *lastEdge ); + lastNode = lastEdge->myNode2; + twinEdge = getTwin( lastEdge ); + } + while ( lastNode != firstNode ); + + } // while ( !theLoops.AllEdgesUsed() ) + + return; + } + + //================================================================================ + /*! + * \brief Erase loops that are cut off by face intersections + */ + //================================================================================ + + void CutFace::CutOffLoops( EdgeLoopSet& theLoops, + const double theSign, + const std::vector< gp_XYZ >& theNormals, + std::vector< EdgePart >& theCutOffLinks, + TLinkMap& theCutOffCoplanarLinks) const + { + EdgePart sideEdge; + for ( size_t i = 0; i < myLinks.size(); ++i ) + { + if ( !myLinks[i].myFace ) + continue; + + EdgeLoop* loop = theLoops.GetLoopOf( & myLinks[i] ); + if ( !loop || loop->myLinks.empty() || loop->myHasPending ) + continue; + + bool toErase, isCoplanar = ( myLinks[i].myIndex == EdgePart::_COPLANAR ); + + gp_Vec iniNorm = theNormals[ myInitFace->GetID() ]; + if ( isCoplanar ) + { + toErase = ( myLinks[i].myFace->GetID() > myInitFace->GetID() ); + + const EdgePart* twin = getTwin( & myLinks[i] ); + if ( !twin || twin->myFace == myLinks[i].myFace ) + { + // only one co-planar face includes myLinks[i] + gp_Vec inFaceDir = iniNorm ^ myLinks[i]; + gp_XYZ edgePnt = SMESH_NodeXYZ( myLinks[i].myNode1 ); + for ( int iN = 0; iN < 3; ++iN ) + { + gp_Vec inCutFaceDir = ( SMESH_NodeXYZ( myLinks[i].myFace->GetNode( iN )) - edgePnt ); + if ( inCutFaceDir * inFaceDir < 0 ) + { + toErase = false; + break; + } + } + } + } + else + { + gp_Vec cutNorm = theNormals[ myLinks[i].myFace->GetID() ]; + gp_Vec inFaceDir = iniNorm ^ myLinks[i]; + + toErase = inFaceDir * cutNorm * theSign < 0; + if ( !toErase ) + { + // erase a neighboring loop + loop = 0; + if ( const EdgePart* twin = getTwin( & myLinks[i] )) + loop = theLoops.GetLoopOf( twin ); + toErase = ( loop && !loop->myLinks.empty() ); + } + } + + if ( toErase ) + { + if ( !isCoplanar ) + { + // remember whole sides of myInitFace that are cut off + for ( size_t iE = 0; iE < loop->myLinks.size(); ++iE ) + { + if ( !loop->myLinks[ iE ]->myFace && + !loop->myLinks[ iE ]->IsInternal() )// && + // !loop->myLinks[ iE ]->myNode1->isMarked() && // cut nodes are marked + // !loop->myLinks[ iE ]->myNode2->isMarked() ) + { + int i = loop->myLinks[ iE ]->myIndex; + sideEdge.Set( myInitFace->GetNode ( i ), + myInitFace->GetNodeWrap( i+1 )); + theCutOffLinks.push_back( sideEdge ); + + if ( !sideEdge.IsSame( *loop->myLinks[ iE ] )) // nodes replaced + { + theCutOffLinks.push_back( *loop->myLinks[ iE ] ); + } + } + else if ( IsCoplanar( loop->myLinks[ iE ])) + { + // propagate erasure to a co-planar face + theCutOffLinks.push_back( *loop->myLinks[ iE ]); + } + else if ( loop->myLinks[ iE ]->myFace && + loop->myLinks[ iE ]->IsInternal() ) + theCutOffLinks.push_back( *loop->myLinks[ iE ]); + } + + // clear the loop + theLoops.Erase( loop ); + } + } + } + return; + } + + //================================================================================ + /*! + * \brief Check if the face has cut edges + */ + //================================================================================ + + bool CutFace::IsCut() const + { + if ( myLinks.size() > 3 ) + return true; + + if ( myLinks.size() == 3 ) + for ( size_t i = 0; i < 3; ++i ) + if ( myLinks[i].myFace ) + return true; + + return false; + } + + //================================================================================ + /*! + * \brief Check if an edge is produced by a co-planar cut + */ + //================================================================================ + + bool CutFace::IsCoplanar( const EdgePart* edge ) const + { + if ( edge->myIndex == EdgePart::_COPLANAR ) + { + const EdgePart* twin = getTwin( edge ); + return ( !twin || twin->myIndex == EdgePart::_COPLANAR ); + } + return false; + } + + //================================================================================ + /*! + * \brief Replace _COPLANAR cut edge by _INTERNAL oe vice versa + */ + //================================================================================ + + bool EdgePart::ReplaceCoplanar( const EdgePart& e ) + { + if ( myIndex + e.myIndex == _COPLANAR + _INTERNAL ) + { + //check if the faces are connected + int nbCommonNodes = SMESH_MeshAlgos::GetCommonNodes( e.myFace, myFace ).size(); + bool toReplace = (( myIndex == _INTERNAL && nbCommonNodes > 1 ) || + ( myIndex == _COPLANAR && nbCommonNodes < 2 )); + if ( toReplace ) + { + myIndex = e.myIndex; + myFace = e.myFace; + return true; + } + } + return false; + } + +} // namespace + +//================================================================================ +/*! + * \brief Create an offsetMesh of given faces + * \param [in] faceIt - the input faces + * \param [out] new2OldFaces - history of faces + * \param [out] new2OldNodes - history of nodes + * \return SMDS_Mesh* - the new offset mesh, a caller should delete + */ +//================================================================================ + +SMDS_Mesh* SMESH_MeshAlgos::MakeOffset( SMDS_ElemIteratorPtr theFaceIt, + SMDS_Mesh& theSrcMesh, + const double theOffset, + const bool theFixIntersections, + TEPairVec& theNew2OldFaces, + TNPairVec& theNew2OldNodes) +{ + SMDS_Mesh* newMesh = new SMDS_Mesh; + theNew2OldFaces.clear(); + theNew2OldNodes.clear(); + theNew2OldFaces.push_back + ( std::make_pair(( const SMDS_MeshElement*) 0, + ( const SMDS_MeshElement*) 0)); // to have index == face->GetID() + + if ( theSrcMesh.GetMeshInfo().NbFaces( ORDER_QUADRATIC ) > 0 ) + throw SALOME_Exception( "Offset of quadratic mesh not supported" ); + if ( theSrcMesh.GetMeshInfo().NbFaces() > theSrcMesh.GetMeshInfo().NbTriangles() ) + throw SALOME_Exception( "Offset of non-triangular mesh not supported" ); + + // copy input faces to the newMesh keeping IDs of nodes + + double minNodeDist = 1e100; + + std::vector< const SMDS_MeshNode* > nodes; + while ( theFaceIt->more() ) + { + const SMDS_MeshElement* face = theFaceIt->next(); + if ( face->GetType() != SMDSAbs_Face ) continue; + + // copy nodes + nodes.assign( face->begin_nodes(), face->end_nodes() ); + for ( size_t i = 0; i < nodes.size(); ++i ) + { + const SMDS_MeshNode* newNode = newMesh->FindNode( nodes[i]->GetID() ); + if ( !newNode ) + { + SMESH_NodeXYZ xyz( nodes[i] ); + newNode = newMesh->AddNodeWithID( xyz.X(), xyz.Y(), xyz.Z(), nodes[i]->GetID() ); + theNew2OldNodes.push_back( std::make_pair( newNode, nodes[i] )); + nodes[i] = newNode; + } + } + const SMDS_MeshElement* newFace = 0; + switch ( face->GetEntityType() ) + { + case SMDSEntity_Triangle: + newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2] ); + break; + case SMDSEntity_Quad_Triangle: + newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2], + nodes[3],nodes[4],nodes[5] ); + break; + case SMDSEntity_BiQuad_Triangle: + newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2], + nodes[3],nodes[4],nodes[5],nodes[6] ); + break; + case SMDSEntity_Quadrangle: + newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],nodes[3] ); + break; + case SMDSEntity_Quad_Quadrangle: + newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],nodes[3], + nodes[4],nodes[5],nodes[6],nodes[7] ); + break; + case SMDSEntity_BiQuad_Quadrangle: + newFace = newMesh->AddFace( nodes[0],nodes[1],nodes[2],nodes[3],nodes[4], + nodes[5],nodes[6],nodes[7],nodes[8] ); + break; + case SMDSEntity_Polygon: + newFace = newMesh->AddPolygonalFace( nodes ); + break; + case SMDSEntity_Quad_Polygon: + newFace = newMesh->AddQuadPolygonalFace( nodes ); + break; + default: + continue; + } + theNew2OldFaces.push_back( std::make_pair( newFace, face )); + + SMESH_NodeXYZ pPrev = nodes.back(), p; + for ( size_t i = 0; i < nodes.size(); ++i ) + { + p.Set( nodes[i] ); + double dist = ( pPrev - p ).SquareModulus(); + if ( dist > std::numeric_limits::min() ) + minNodeDist = dist; + pPrev = p; + } + } // while ( faceIt->more() ) + + + // compute normals to faces + std::vector< gp_XYZ > normals( theNew2OldFaces.size() ); + for ( size_t i = 1; i < normals.size(); ++i ) + { + if ( !SMESH_MeshAlgos::FaceNormal( theNew2OldFaces[i].second, normals[i] )) + normals[i].SetCoord( 0,0,0 ); // TODO find norm by neighbors + } + + const double tol = 1e-3 * Sqrt( minNodeDist ); + const double sign = ( theOffset < 0 ? -1 : +1 ); + + // translate new nodes by normal to input faces + gp_XYZ newXYZ; + std::vector< const SMDS_MeshNode* > multiNormalNodes; + for ( size_t i = 0; i < theNew2OldNodes.size(); ++i ) + { + const SMDS_MeshNode* newNode = theNew2OldNodes[i].first; + + if ( getTranslatedPosition( newNode, theOffset, tol*10., sign, normals, theSrcMesh, newXYZ )) + newMesh->MoveNode( newNode, newXYZ.X(), newXYZ.Y(), newXYZ.Z() ); + else + multiNormalNodes.push_back( newNode ); + } + // make multi-normal translation + std::vector< SMESH_NodeXYZ > multiPos(10); + for ( size_t i = 0; i < multiNormalNodes.size(); ++i ) + { + const SMDS_MeshNode* newNode = multiNormalNodes[i]; + newNode->setIsMarked( true ); + SMESH_NodeXYZ oldXYZ = newNode; + multiPos.clear(); + for ( SMDS_ElemIteratorPtr fIt = newNode->GetInverseElementIterator(); fIt->more(); ) + { + const SMDS_MeshElement* newFace = fIt->next(); + const int faceIndex = newFace->GetID(); + const gp_XYZ& oldNorm = normals[ faceIndex ]; + const gp_XYZ newXYZ = oldXYZ + oldNorm * theOffset; + if ( multiPos.empty() ) + { + newMesh->MoveNode( newNode, newXYZ.X(), newXYZ.Y(), newXYZ.Z() ); + multiPos.emplace_back( newNode ); + } + else + { + newNode = 0; + for ( size_t iP = 0; iP < multiPos.size() && !newNode; ++iP ) + if (( multiPos[iP] - newXYZ ).SquareModulus() < tol * tol ) + newNode = multiPos[iP].Node(); + if ( !newNode ) + { + newNode = newMesh->AddNode( newXYZ.X(), newXYZ.Y(), newXYZ.Z() ); + newNode->setIsMarked( true ); + theNew2OldNodes.push_back( std::make_pair( newNode, theNew2OldNodes[i].second )); + multiPos.emplace_back( newNode ); + } + } + if ( newNode != oldXYZ.Node() ) + { + nodes.assign( newFace->begin_nodes(), newFace->end_nodes() ); + nodes[ newFace->GetNodeIndex( oldXYZ.Node() )] = newNode; + newMesh->ChangeElementNodes( newFace, & nodes[0], nodes.size() ); + } + } + } + + if ( !theFixIntersections ) + return newMesh; + + + // remove new faces around concave nodes (they are marked) if the faces are inverted + gp_XYZ faceNorm; + for ( size_t i = 0; i < theNew2OldNodes.size(); ++i ) + { + const SMDS_MeshNode* newNode = theNew2OldNodes[i].first; + //const SMDS_MeshNode* oldNode = theNew2OldNodes[i].second; + if ( newNode->isMarked() ) + { + //gp_XYZ moveVec = sign * ( SMESH_NodeXYZ( newNode ) - SMESH_NodeXYZ( oldNode )); + + //bool haveInverseFace = false; + for ( SMDS_ElemIteratorPtr fIt = newNode->GetInverseElementIterator(); fIt->more(); ) + { + const SMDS_MeshElement* newFace = fIt->next(); + const int faceIndex = newFace->GetID(); + const gp_XYZ& oldNorm = normals[ faceIndex ]; + if ( !SMESH_MeshAlgos::FaceNormal( newFace, faceNorm, /*normalize=*/false ) || + //faceNorm * moveVec < 0 ) + faceNorm * oldNorm < 0 ) + { + //haveInverseFace = true; + theNew2OldFaces[ faceIndex ].first = 0; + newMesh->RemoveFreeElement( newFace ); + //break; + } + } + // if ( haveInverseFace ) + // { + // newMesh->MoveNode( newNode, oldNode->X(), oldNode->Y(), oldNode->Z() ); + + // for ( SMDS_ElemIteratorPtr fIt = newNode->GetInverseElementIterator(); fIt->more(); ) + // { + // const SMDS_MeshElement* newFace = fIt->next(); + // if ( !SMESH_MeshAlgos::FaceNormal( newFace, normals[ newFace->GetID() ] )) + // normals[i].SetCoord( 0,0,0 ); // TODO find norm by neighbors + // } + // } + } + // mark all new nodes located closer than theOffset from theSrcMesh + } + + // ================================================== + // find self-intersections of new faces and fix them + // ================================================== + + std::unique_ptr< SMESH_ElementSearcher > fSearcher + ( SMESH_MeshAlgos::GetElementSearcher( *newMesh, tol )); + + Intersector intersector( newMesh, tol, normals ); + + std::vector< const SMDS_MeshElement* > closeFaces; + std::vector< const SMDS_MeshNode* > faceNodes; + Bnd_B3d faceBox; + for ( size_t iF = 1; iF < theNew2OldFaces.size(); ++iF ) + { + const SMDS_MeshElement* newFace = theNew2OldFaces[iF].first; + if ( !newFace ) continue; + faceNodes.assign( newFace->begin_nodes(), newFace->end_nodes() ); + + bool isConcaveNode1 = false; + for ( size_t iN = 0; iN < faceNodes.size() && !isConcaveNode1; ++iN ) + isConcaveNode1 = faceNodes[iN]->isMarked(); + + // get faces close to a newFace + closeFaces.clear(); + faceBox.Clear(); + for ( size_t i = 0; i < faceNodes.size(); ++i ) + faceBox.Add( SMESH_NodeXYZ( faceNodes[i] )); + faceBox.Enlarge( tol ); + + fSearcher->GetElementsInBox( faceBox, SMDSAbs_Face, closeFaces ); + + // intersect the newFace with closeFaces + + for ( size_t i = 0; i < closeFaces.size(); ++i ) + { + const SMDS_MeshElement* closeFace = closeFaces[i]; + if ( closeFace->GetID() <= newFace->GetID() ) + continue; // this pair already treated + + // do not intersect connected faces if they have no concave nodes + int nbCommonNodes = 0; + for ( size_t iN = 0; iN < faceNodes.size(); ++iN ) + nbCommonNodes += ( closeFace->GetNodeIndex( faceNodes[iN] ) >= 0 ); + + if ( !isConcaveNode1 ) + { + bool isConcaveNode2 = false; + for ( SMDS_ElemIteratorPtr nIt = closeFace->nodesIterator(); nIt->more(); ) + if (( isConcaveNode2 = nIt->next()->isMarked() )) + break; + + if ( !isConcaveNode2 && nbCommonNodes > 0 ) + continue; + } + + intersector.Cut( newFace, closeFace, nbCommonNodes ); + } + } + intersector.MakeNewFaces( theNew2OldFaces, theNew2OldNodes, sign ); + + return newMesh; +} diff --git a/src/SMESHUtils/SMESH_Triangulate.cxx b/src/SMESHUtils/SMESH_Triangulate.cxx new file mode 100644 index 000000000..6652aeae1 --- /dev/null +++ b/src/SMESHUtils/SMESH_Triangulate.cxx @@ -0,0 +1,341 @@ +// 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_Triangulate.cxx +// Created : Thu Jan 18 18:00:13 2018 +// Author : Edward AGAPOV (eap) + +// Extracted from ../DriverSTL/DriverSTL_W_SMDS_Mesh.cxx + +#include "SMESH_MeshAlgos.hxx" + +#include +#include +#include + +using namespace SMESH_MeshAlgos; + +//================================================================================ +/*! + * \brief Initialization + */ +//================================================================================ + +void Triangulate::PolyVertex::SetNodeAndNext( const SMDS_MeshNode* n, + PolyVertex& v ) +{ + _nxyz.Set( n ); + _next = &v; + v._prev = this; +} +//================================================================================ +/*! + * \brief Remove self from a polygon + */ +//================================================================================ + +Triangulate::PolyVertex* Triangulate::PolyVertex::Delete() +{ + _prev->_next = _next; + _next->_prev = _prev; + return _next; +} + +//================================================================================ +/*! + * \brief Return nodes of a triangle + */ +//================================================================================ + +void Triangulate::PolyVertex::GetTriaNodes( const SMDS_MeshNode** nodes) const +{ + nodes[0] = _prev->_nxyz._node; + nodes[1] = this->_nxyz._node; + nodes[2] = _next->_nxyz._node; +} + +//================================================================================ +/*! + * \brief Compute triangle area + */ +//================================================================================ + +inline static double Area( const gp_XY& xy0, const gp_XY& xy1, const gp_XY& xy2 ) +{ + gp_XY vPrev = xy0 - xy1; + gp_XY vNext = xy2 - xy1; + return vNext ^ vPrev; +} + +//================================================================================ +/*! + * \brief Compute triangle area + */ +//================================================================================ + +double Triangulate::PolyVertex::TriaArea() const +{ + return Area( _prev->_xy, this->_xy, _next->_xy ); +} + +//================================================================================ +/*! + * \brief Check if a vertex is inside a triangle + */ +//================================================================================ + +bool Triangulate::PolyVertex::IsInsideTria( const PolyVertex* v ) +{ + if ( this ->_nxyz == v->_nxyz || + _prev->_nxyz == v->_nxyz || + _next->_nxyz == v->_nxyz ) + return false; + + gp_XY p = _prev->_xy - v->_xy; + gp_XY t = this->_xy - v->_xy; + gp_XY n = _next->_xy - v->_xy; + const double tol = -1e-12; + return (( p ^ t ) >= tol && + ( t ^ n ) >= tol && + ( n ^ p ) >= tol ); + // return ( Area( _prev, this, v ) > 0 && + // Area( this, _next, v ) > 0 && + // Area( _next, _prev, v ) > 0 ); +} + +//================================================================================ +/*! + * \brief Triangulate a polygon. Assure correct orientation for concave polygons + */ +//================================================================================ + +bool Triangulate::triangulate( std::vector< const SMDS_MeshNode*>& nodes, + const size_t nbNodes ) +{ + // connect nodes into a ring + _pv.resize( nbNodes ); + for ( size_t i = 1; i < nbNodes; ++i ) + _pv[i-1].SetNodeAndNext( nodes[i-1], _pv[i] ); + _pv[ nbNodes-1 ].SetNodeAndNext( nodes[ nbNodes-1 ], _pv[0] ); + + // get a polygon normal + gp_XYZ normal(0,0,0), p0,v01,v02; + p0 = _pv[0]._nxyz; + v01 = _pv[1]._nxyz - p0; + for ( size_t i = 2; i < nbNodes; ++i ) + { + v02 = _pv[i]._nxyz - p0; + normal += v01 ^ v02; + v01 = v02; + } + // project nodes to the found plane + gp_Ax2 axes; + try { + axes = gp_Ax2( p0, normal, v01 ); + } + catch ( Standard_Failure ) { + return false; + } + for ( size_t i = 0; i < nbNodes; ++i ) + { + gp_XYZ p = _pv[i]._nxyz - p0; + _pv[i]._xy.SetX( axes.XDirection().XYZ() * p ); + _pv[i]._xy.SetY( axes.YDirection().XYZ() * p ); + } + + // in a loop, find triangles with positive area and having no vertices inside + int iN = 0, nbTria = nbNodes - 2; + nodes.reserve( nbTria * 3 ); + const double minArea = 1e-6; + PolyVertex* v = &_pv[0], *vi; + int nbVertices = nbNodes, nbBadTria = 0, isGoodTria; + while ( nbBadTria < nbVertices ) + { + if (( isGoodTria = v->TriaArea() > minArea )) + { + for ( vi = v->_next->_next; + vi != v->_prev; + vi = vi->_next ) + { + if ( v->IsInsideTria( vi )) + break; + } + isGoodTria = ( vi == v->_prev ); + } + if ( isGoodTria ) + { + v->GetTriaNodes( &nodes[ iN ] ); + iN += 3; + v = v->Delete(); + if ( --nbVertices == 3 ) + { + // last triangle remains + v->GetTriaNodes( &nodes[ iN ] ); + return true; + } + nbBadTria = 0; + } + else + { + v = v->_next; + ++nbBadTria; + } + } + + // the polygon is invalid; add triangles with positive area + nbBadTria = 0; + while ( nbBadTria < nbVertices ) + { + isGoodTria = v->TriaArea() > minArea; + if ( isGoodTria ) + { + v->GetTriaNodes( &nodes[ iN ] ); + iN += 3; + v = v->Delete(); + if ( --nbVertices == 3 ) + { + // last triangle remains + v->GetTriaNodes( &nodes[ iN ] ); + return true; + } + nbBadTria = 0; + } + else + { + v = v->_next; + ++nbBadTria; + } + } + + // add all the rest triangles + while ( nbVertices >= 3 ) + { + v->GetTriaNodes( &nodes[ iN ] ); + iN += 3; + v = v->Delete(); + --nbVertices; + } + + return true; + +} // triangulate() + +//================================================================================ +/*! + * \brief Return nb triangles in a decomposed mesh face + * \retval int - number of triangles + */ +//================================================================================ + +int Triangulate::GetNbTriangles( const SMDS_MeshElement* face ) +{ + // WARNING: counting triangles must be coherent with GetTriangles() + switch ( face->GetEntityType() ) + { + case SMDSEntity_BiQuad_Triangle: + case SMDSEntity_BiQuad_Quadrangle: + return face->NbNodes() - 1; + // case SMDSEntity_Triangle: + // case SMDSEntity_Quad_Triangle: + // case SMDSEntity_Quadrangle: + // case SMDSEntity_Quad_Quadrangle: + // case SMDSEntity_Polygon: + // case SMDSEntity_Quad_Polygon: + default: + return face->NbNodes() - 2; + } + return 0; +} + +//================================================================================ +/*! + * \brief Decompose a mesh face into triangles + * \retval int - number of triangles + */ +//================================================================================ + +int Triangulate::GetTriangles( const SMDS_MeshElement* face, + std::vector< const SMDS_MeshNode*>& nodes) +{ + if ( face->GetType() != SMDSAbs_Face ) + return 0; + + // WARNING: decomposing into triangles must be coherent with getNbTriangles() + int nbTria, i = 0, nbNodes = face->NbNodes(); + SMDS_NodeIteratorPtr nIt = face->interlacedNodesIterator(); + nodes.resize( nbNodes * 3 ); + nodes[ i++ ] = nIt->next(); + nodes[ i++ ] = nIt->next(); + + const SMDSAbs_EntityType type = face->GetEntityType(); + switch ( type ) + { + case SMDSEntity_BiQuad_Triangle: + case SMDSEntity_BiQuad_Quadrangle: + + nbTria = ( type == SMDSEntity_BiQuad_Triangle ) ? 6 : 8; + nodes[ i++ ] = face->GetNode( nbTria ); + for ( i = 3; i < 3*(nbTria-1); i += 3 ) + { + nodes[ i+0 ] = nodes[ i-2 ]; + nodes[ i+1 ] = nIt->next(); + nodes[ i+2 ] = nodes[ 2 ]; + } + nodes[ i+0 ] = nodes[ i-2 ]; + nodes[ i+1 ] = nodes[ 0 ]; + nodes[ i+2 ] = nodes[ 2 ]; + break; + + case SMDSEntity_Triangle: + + nbTria = 1; + nodes[ i++ ] = nIt->next(); + break; + + default: + + // case SMDSEntity_Quad_Triangle: + // case SMDSEntity_Quadrangle: + // case SMDSEntity_Quad_Quadrangle: + // case SMDSEntity_Polygon: + // case SMDSEntity_Quad_Polygon: + + nbTria = nbNodes - 2; + while ( nIt->more() ) + nodes[ i++ ] = nIt->next(); + + if ( nbTria > 1 && !triangulate( nodes, nbNodes )) + { + nIt = face->interlacedNodesIterator(); + nodes[ 0 ] = nIt->next(); + nodes[ 1 ] = nIt->next(); + nodes[ 2 ] = nIt->next(); + for ( i = 3; i < 3*nbTria; i += 3 ) + { + nodes[ i+0 ] = nodes[ 0 ]; + nodes[ i+1 ] = nodes[ i-1 ]; + nodes[ i+2 ] = nIt->next(); + } + } + } + + return nbTria; +} diff --git a/src/SMESHUtils/SMESH_Triangulate.hxx b/src/SMESHUtils/SMESH_Triangulate.hxx new file mode 100644 index 000000000..8734c4b7a --- /dev/null +++ b/src/SMESHUtils/SMESH_Triangulate.hxx @@ -0,0 +1,51 @@ +// 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_Triangulate.hxx +// Created : Thu Jan 18 17:51:34 2018 +// Author : Edward AGAPOV (eap) + + +#ifndef __SMESH_Triangulate_HXX__ +#define __SMESH_Triangulate_HXX__ + +/*! + * \brief Divide a mesh face into triangles + */ +class SMESHUtils_EXPORT SMESH_Triangulate +{ + public: + + static int GetNbTriangles( const SMDS_MeshElement* face ); + + int GetTriangles( const SMDS_MeshElement* face, + std::vector< const SMDS_MeshNode*>& nodes); + + private: + + bool triangulate( std::vector< const SMDS_MeshNode*>& nodes, const size_t nbNodes ); + + struct PolyVertex; + std::vector< PolyVertex > _pv; +}; + + +#endif