#include "SMESH_OctreeNode.hxx"
#include "SMESH_subMesh.hxx"
-#include <Basics_OCCTVersion.hxx>
-
#include "utilities.h"
#include "chrono.hxx"
#include <sstream>
#include <boost/tuple/tuple.hpp>
+#include <boost/container/flat_set.hpp>
#include <Standard_Failure.hxx>
#include <Standard_ErrorHandler.hxx>
+#include <OSD_Parallel.hxx>
+
+#include "SMESH_TryCatch.hxx" // include after OCCT headers!
#define cast2Node(elem) static_cast<const SMDS_MeshNode*>( elem )
}
if ( otherFace && otherFace != theFace)
{
- // link must be reverse in otherFace if orientation ot otherFace
+ // link must be reverse in otherFace if orientation to otherFace
// is same as that of theFace
if ( abs(nodeInd2-nodeInd1) == 1 ? nodeInd2 > nodeInd1 : nodeInd1 > nodeInd2 )
{
//=======================================================================
/*!
* \brief Split each of given quadrangles into 4 triangles.
- * \param theElems - The faces to be splitted. If empty all faces are split.
+ * \param theElems - The faces to be split. If empty all faces are split.
*/
//=======================================================================
break;
case SMDSEntity_Penta:
case SMDSEntity_Quad_Penta:
+ case SMDSEntity_BiQuad_Penta:
connVariants = thePentaTo3; nbTet = 3; nbVariants = 6;
break;
default:
std::swap( itNN[0], itNN[1] );
std::swap( prevNod[0], prevNod[1] );
std::swap( nextNod[0], nextNod[1] );
-#if defined(__APPLE__)
std::swap( isSingleNode[0], isSingleNode[1] );
-#else
- isSingleNode.swap( isSingleNode[0], isSingleNode[1] );
-#endif
if ( nbSame > 0 )
sames[0] = 1 - sames[0];
iNotSameNode = 1 - iNotSameNode;
makeEdgePathPoints(aPrms, aTrackEdge,(aN1->GetID()==startNid), LPP);
LLPPs.push_back(LPP);
UsedNums.Add(k);
- // update startN for search following egde
+ // update startN for search following edge
if( aN1->GetID() == startNid ) startNid = aN2->GetID();
else startNid = aN1->GetID();
break;
makeEdgePathPoints(aPrms, aTrackEdge, aN1isOK, LPP);
LLPPs.push_back(LPP);
UsedNums.Add(k);
- // update startN for search following egde
+ // update startN for search following edge
if ( aN1isOK ) aVprev = aV2;
else aVprev = aV1;
break;
typedef map< SortableElement, int > TMapOfNodeSet;
typedef list<int> TGroupOfElems;
- if ( theElements.empty() )
- { // get all elements in the mesh
- SMDS_ElemIteratorPtr eIt = GetMeshDS()->elementsIterator();
- while ( eIt->more() )
- theElements.insert( theElements.end(), eIt->next() );
- }
+ SMDS_ElemIteratorPtr elemIt;
+ if ( theElements.empty() ) elemIt = GetMeshDS()->elementsIterator();
+ else elemIt = elemSetIterator( theElements );
vector< TGroupOfElems > arrayOfGroups;
TGroupOfElems groupOfElems;
TMapOfNodeSet mapOfNodeSet;
- TIDSortedElemSet::iterator elemIt = theElements.begin();
- for ( int i = 0; elemIt != theElements.end(); ++elemIt )
+ for ( int iGroup = 0; elemIt->more(); )
{
- const SMDS_MeshElement* curElem = *elemIt;
+ const SMDS_MeshElement* curElem = elemIt->next();
SortableElement SE(curElem);
// check uniqueness
- pair< TMapOfNodeSet::iterator, bool> pp = mapOfNodeSet.insert(make_pair(SE, i));
+ pair< TMapOfNodeSet::iterator, bool> pp = mapOfNodeSet.insert(make_pair(SE, iGroup));
if ( !pp.second ) { // one more coincident elem
TMapOfNodeSet::iterator& itSE = pp.first;
- int ind = (*itSE).second;
- arrayOfGroups[ind].push_back( curElem->GetID() );
+ int iG = itSE->second;
+ arrayOfGroups[ iG ].push_back( curElem->GetID() );
}
else {
arrayOfGroups.push_back( groupOfElems );
arrayOfGroups.back().push_back( curElem->GetID() );
- i++;
+ iGroup++;
}
}
SMESH_MesherHelper& theHelper,
const bool theForce3d)
{
+ //MESSAGE("convertElemToQuadratic");
int nbElem = 0;
if( !theSm ) return nbElem;
case SMDSEntity_Quad_Triangle:
case SMDSEntity_Quad_Quadrangle:
case SMDSEntity_Quad_Hexa:
+ case SMDSEntity_Quad_Penta:
alreadyOK = !theHelper.GetIsBiQuadratic(); break;
case SMDSEntity_BiQuad_Triangle:
case SMDSEntity_BiQuad_Quadrangle:
case SMDSEntity_TriQuad_Hexa:
+ case SMDSEntity_BiQuad_Penta:
alreadyOK = theHelper.GetIsBiQuadratic();
hasCentralNodes = true;
break;
default:
alreadyOK = true;
}
- // take into account already present modium nodes
+ // take into account already present medium nodes
switch ( aType ) {
case SMDSAbs_Volume:
theHelper.AddTLinks( static_cast< const SMDS_MeshVolume* >( elem )); break;
NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], id, theForce3d);
break;
case SMDSEntity_Penta:
+ case SMDSEntity_Quad_Penta:
+ case SMDSEntity_BiQuad_Penta:
NewElem = theHelper.AddVolume(nodes[0], nodes[1], nodes[2], nodes[3], nodes[4], nodes[5], id, theForce3d);
break;
case SMDSEntity_Hexa:
void SMESH_MeshEditor::ConvertToQuadratic(const bool theForce3d, const bool theToBiQuad)
{
+ //MESSAGE("ConvertToQuadratic "<< theForce3d << " " << theToBiQuad);
SMESHDS_Mesh* meshDS = GetMeshDS();
SMESH_MesherHelper aHelper(*myMesh);
{
case SMDSEntity_Quad_Hexa: alreadyOK = !theToBiQuad; break;
case SMDSEntity_TriQuad_Hexa: alreadyOK = theToBiQuad; break;
+ case SMDSEntity_Quad_Penta: alreadyOK = !theToBiQuad; break;
+ case SMDSEntity_BiQuad_Penta: alreadyOK = theToBiQuad; break;
default: alreadyOK = true;
}
if ( alreadyOK )
nodes[3], nodes[4], id, theForce3d);
break;
case SMDSEntity_Penta:
+ case SMDSEntity_Quad_Penta:
+ case SMDSEntity_BiQuad_Penta:
NewVolume = aHelper.AddVolume(nodes[0], nodes[1], nodes[2],
nodes[3], nodes[4], nodes[5], id, theForce3d);
+ for ( size_t i = 15; i < nodes.size(); ++i ) // rm central nodes
+ if ( nodes[i]->NbInverseElements() == 0 )
+ GetMeshDS()->RemoveFreeNode( nodes[i], /*sm=*/0, /*fromGroups=*/true );
break;
case SMDSEntity_Hexagonal_Prism:
default:
return SEW_OK;
}
+namespace // automatically find theAffectedElems for DoubleNodes()
+{
+ bool isOut( const SMDS_MeshNode* n, const gp_XYZ& norm, const SMDS_MeshElement* elem );
+
+ //--------------------------------------------------------------------------------
+ // Nodes shared by adjacent FissureBorder's.
+ // 1 node if FissureBorder separates faces
+ // 2 nodes if FissureBorder separates volumes
+ struct SubBorder
+ {
+ const SMDS_MeshNode* _nodes[2];
+ int _nbNodes;
+
+ SubBorder( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2 = 0 )
+ {
+ _nodes[0] = n1;
+ _nodes[1] = n2;
+ _nbNodes = bool( n1 ) + bool( n2 );
+ if ( _nbNodes == 2 && n1 > n2 )
+ std::swap( _nodes[0], _nodes[1] );
+ }
+ bool operator<( const SubBorder& other ) const
+ {
+ for ( int i = 0; i < _nbNodes; ++i )
+ {
+ if ( _nodes[i] < other._nodes[i] ) return true;
+ if ( _nodes[i] > other._nodes[i] ) return false;
+ }
+ return false;
+ }
+ };
+
+ //--------------------------------------------------------------------------------
+ // Map a SubBorder to all FissureBorder it bounds
+ struct FissureBorder;
+ typedef std::map< SubBorder, std::vector< FissureBorder* > > TBorderLinks;
+ typedef TBorderLinks::iterator TMappedSub;
+
+ //--------------------------------------------------------------------------------
+ /*!
+ * \brief Element border (volume facet or face edge) at a fissure
+ */
+ struct FissureBorder
+ {
+ std::vector< const SMDS_MeshNode* > _nodes; // border nodes
+ const SMDS_MeshElement* _elems[2]; // volume or face adjacent to fissure
+
+ std::vector< TMappedSub > _mappedSubs; // Sub() in TBorderLinks map
+ std::vector< const SMDS_MeshNode* > _sortedNodes; // to compare FissureBorder's
+
+ FissureBorder( FissureBorder && from ) // move constructor
+ {
+ std::swap( _nodes, from._nodes );
+ std::swap( _sortedNodes, from._sortedNodes );
+ _elems[0] = from._elems[0];
+ _elems[1] = from._elems[1];
+ }
+
+ FissureBorder( const SMDS_MeshElement* elemToDuplicate,
+ std::vector< const SMDS_MeshElement* > & adjElems)
+ : _nodes( elemToDuplicate->NbCornerNodes() )
+ {
+ for ( size_t i = 0; i < _nodes.size(); ++i )
+ _nodes[i] = elemToDuplicate->GetNode( i );
+
+ SMDSAbs_ElementType type = SMDSAbs_ElementType( elemToDuplicate->GetType() + 1 );
+ findAdjacent( type, adjElems );
+ }
+
+ FissureBorder( const SMDS_MeshNode** nodes,
+ const size_t nbNodes,
+ const SMDSAbs_ElementType adjElemsType,
+ std::vector< const SMDS_MeshElement* > & adjElems)
+ : _nodes( nodes, nodes + nbNodes )
+ {
+ findAdjacent( adjElemsType, adjElems );
+ }
+
+ void findAdjacent( const SMDSAbs_ElementType adjElemsType,
+ std::vector< const SMDS_MeshElement* > & adjElems)
+ {
+ _elems[0] = _elems[1] = 0;
+ adjElems.clear();
+ if ( SMDS_Mesh::GetElementsByNodes( _nodes, adjElems, adjElemsType ))
+ for ( size_t i = 0; i < adjElems.size() && i < 2; ++i )
+ _elems[i] = adjElems[i];
+ }
+
+ bool operator<( const FissureBorder& other ) const
+ {
+ return GetSortedNodes() < other.GetSortedNodes();
+ }
+
+ const std::vector< const SMDS_MeshNode* >& GetSortedNodes() const
+ {
+ if ( _sortedNodes.empty() && !_nodes.empty() )
+ {
+ FissureBorder* me = const_cast<FissureBorder*>( this );
+ me->_sortedNodes = me->_nodes;
+ std::sort( me->_sortedNodes.begin(), me->_sortedNodes.end() );
+ }
+ return _sortedNodes;
+ }
+
+ size_t NbSub() const
+ {
+ return _nodes.size();
+ }
+
+ SubBorder Sub(size_t i) const
+ {
+ return SubBorder( _nodes[i], NbSub() > 2 ? _nodes[ (i+1)%NbSub() ] : 0 );
+ }
+
+ void AddSelfTo( TBorderLinks& borderLinks )
+ {
+ _mappedSubs.resize( NbSub() );
+ for ( size_t i = 0; i < NbSub(); ++i )
+ {
+ TBorderLinks::iterator s2b =
+ borderLinks.insert( std::make_pair( Sub(i), TBorderLinks::mapped_type() )).first;
+ s2b->second.push_back( this );
+ _mappedSubs[ i ] = s2b;
+ }
+ }
+
+ void Clear()
+ {
+ _nodes.clear();
+ }
+
+ const SMDS_MeshElement* GetMarkedElem() const
+ {
+ if ( _nodes.empty() ) return 0; // cleared
+ if ( _elems[0] && _elems[0]->isMarked() ) return _elems[0];
+ if ( _elems[1] && _elems[1]->isMarked() ) return _elems[1];
+ return 0;
+ }
+
+ gp_XYZ GetNorm() const // normal to the border
+ {
+ gp_XYZ norm;
+ if ( _nodes.size() == 2 )
+ {
+ gp_XYZ avgNorm( 0,0,0 ); // sum of normals of adjacent faces
+ if ( SMESH_MeshAlgos::FaceNormal( _elems[0], norm ))
+ avgNorm += norm;
+ if ( SMESH_MeshAlgos::FaceNormal( _elems[1], norm ))
+ avgNorm += norm;
+
+ gp_XYZ bordDir( SMESH_NodeXYZ( _nodes[0] ) - SMESH_NodeXYZ( _nodes[1] ));
+ norm = bordDir ^ avgNorm;
+ }
+ else
+ {
+ SMESH_NodeXYZ p0( _nodes[0] );
+ SMESH_NodeXYZ p1( _nodes[1] );
+ SMESH_NodeXYZ p2( _nodes[2] );
+ norm = ( p0 - p1 ) ^ ( p2 - p1 );
+ }
+ if ( isOut( _nodes[0], norm, GetMarkedElem() ))
+ norm.Reverse();
+
+ return norm;
+ }
+
+ void ChooseSide() // mark an _elem located at positive side of fissure
+ {
+ _elems[0]->setIsMarked( true );
+ gp_XYZ norm = GetNorm();
+ double maxX = norm.Coord(1);
+ if ( Abs( maxX ) < Abs( norm.Coord(2)) ) maxX = norm.Coord(2);
+ if ( Abs( maxX ) < Abs( norm.Coord(3)) ) maxX = norm.Coord(3);
+ if ( maxX < 0 )
+ {
+ _elems[0]->setIsMarked( false );
+ _elems[1]->setIsMarked( true );
+ }
+ }
+
+ }; // struct FissureBorder
+
+ //--------------------------------------------------------------------------------
+ /*!
+ * \brief Classifier of elements at fissure edge
+ */
+ class FissureNormal
+ {
+ std::vector< gp_XYZ > _normals;
+ bool _bothIn;
+
+ public:
+ void Add( const SMDS_MeshNode* n, const FissureBorder& bord )
+ {
+ _bothIn = false;
+ _normals.reserve(2);
+ _normals.push_back( bord.GetNorm() );
+ if ( _normals.size() == 2 )
+ _bothIn = !isOut( n, _normals[0], bord.GetMarkedElem() );
+ }
+
+ bool IsIn( const SMDS_MeshNode* n, const SMDS_MeshElement* elem ) const
+ {
+ bool isIn = false;
+ switch ( _normals.size() ) {
+ case 1:
+ {
+ isIn = !isOut( n, _normals[0], elem );
+ break;
+ }
+ case 2:
+ {
+ bool in1 = !isOut( n, _normals[0], elem );
+ bool in2 = !isOut( n, _normals[1], elem );
+ isIn = _bothIn ? ( in1 && in2 ) : ( in1 || in2 );
+ }
+ }
+ return isIn;
+ }
+ };
+
+ //================================================================================
+ /*!
+ * \brief Classify an element by a plane passing through a node
+ */
+ //================================================================================
+
+ bool isOut( const SMDS_MeshNode* n, const gp_XYZ& norm, const SMDS_MeshElement* elem )
+ {
+ SMESH_NodeXYZ p = n;
+ double sumDot = 0;
+ for ( int i = 0, nb = elem->NbCornerNodes(); i < nb; ++i )
+ {
+ SMESH_NodeXYZ pi = elem->GetNode( i );
+ sumDot += norm * ( pi - p );
+ }
+ return sumDot < -1e-100;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Find FissureBorder's by nodes to duplicate
+ */
+ //================================================================================
+
+ void findFissureBorders( const TIDSortedElemSet& theNodes,
+ std::vector< FissureBorder > & theFissureBorders )
+ {
+ TIDSortedElemSet::const_iterator nIt = theNodes.begin();
+ const SMDS_MeshNode* n = dynamic_cast< const SMDS_MeshNode*>( *nIt );
+ if ( !n ) return;
+ SMDSAbs_ElementType elemType = SMDSAbs_Volume;
+ if ( n->NbInverseElements( elemType ) == 0 )
+ {
+ elemType = SMDSAbs_Face;
+ if ( n->NbInverseElements( elemType ) == 0 )
+ return;
+ }
+ // unmark elements touching the fissure
+ for ( ; nIt != theNodes.end(); ++nIt )
+ SMESH_MeshAlgos::MarkElems( cast2Node(*nIt)->GetInverseElementIterator(), false );
+
+ // loop on elements touching the fissure to get their borders belonging to the fissure
+ std::set< FissureBorder > fissureBorders;
+ std::vector< const SMDS_MeshElement* > adjElems;
+ std::vector< const SMDS_MeshNode* > nodes;
+ SMDS_VolumeTool volTool;
+ for ( nIt = theNodes.begin(); nIt != theNodes.end(); ++nIt )
+ {
+ SMDS_ElemIteratorPtr invIt = cast2Node(*nIt)->GetInverseElementIterator( elemType );
+ while ( invIt->more() )
+ {
+ const SMDS_MeshElement* eInv = invIt->next();
+ if ( eInv->isMarked() ) continue;
+ eInv->setIsMarked( true );
+
+ if ( elemType == SMDSAbs_Volume )
+ {
+ volTool.Set( eInv );
+ int iQuad = eInv->IsQuadratic() ? 2 : 1;
+ for ( int iF = 0, nbF = volTool.NbFaces(); iF < nbF; ++iF )
+ {
+ const SMDS_MeshNode** nn = volTool.GetFaceNodes( iF );
+ int nbN = volTool.NbFaceNodes( iF ) / iQuad;
+ nodes.clear();
+ bool allOnFissure = true;
+ for ( int iN = 0; iN < nbN && allOnFissure; iN += iQuad )
+ if (( allOnFissure = theNodes.count( nn[ iN ])))
+ nodes.push_back( nn[ iN ]);
+ if ( allOnFissure )
+ fissureBorders.insert( std::move( FissureBorder( &nodes[0], nodes.size(),
+ elemType, adjElems )));
+ }
+ }
+ else // elemType == SMDSAbs_Face
+ {
+ const SMDS_MeshNode* nn[2] = { eInv->GetNode( eInv->NbCornerNodes()-1 ), 0 };
+ bool onFissure0 = theNodes.count( nn[0] ), onFissure1;
+ for ( int iN = 0, nbN = eInv->NbCornerNodes(); iN < nbN; ++iN )
+ {
+ nn[1] = eInv->GetNode( iN );
+ onFissure1 = theNodes.count( nn[1] );
+ if ( onFissure0 && onFissure1 )
+ fissureBorders.insert( std::move( FissureBorder( nn, 2, elemType, adjElems )));
+ nn[0] = nn[1];
+ onFissure0 = onFissure1;
+ }
+ }
+ }
+ }
+
+ theFissureBorders.reserve( theFissureBorders.size() + fissureBorders.size());
+ std::set< FissureBorder >::iterator bord = fissureBorders.begin();
+ for ( ; bord != fissureBorders.end(); ++bord )
+ {
+ theFissureBorders.push_back( std::move( const_cast<FissureBorder&>( *bord ) ));
+ }
+ return;
+ } // findFissureBorders()
+
+ //================================================================================
+ /*!
+ * \brief Find elements on one side of a fissure defined by elements or nodes to duplicate
+ * \param [in] theElemsOrNodes - elements or nodes to duplicate
+ * \param [in] theNodesNot - nodes not to duplicate
+ * \param [out] theAffectedElems - the found elements
+ */
+ //================================================================================
+
+ void findAffectedElems( const TIDSortedElemSet& theElemsOrNodes,
+ TIDSortedElemSet& theAffectedElems)
+ {
+ if ( theElemsOrNodes.empty() ) return;
+
+ // find FissureBorder's
+
+ std::vector< FissureBorder > fissure;
+ std::vector< const SMDS_MeshElement* > elemsByFacet;
+
+ TIDSortedElemSet::const_iterator elIt = theElemsOrNodes.begin();
+ if ( (*elIt)->GetType() == SMDSAbs_Node )
+ {
+ findFissureBorders( theElemsOrNodes, fissure );
+ }
+ else
+ {
+ fissure.reserve( theElemsOrNodes.size() );
+ for ( ; elIt != theElemsOrNodes.end(); ++elIt )
+ fissure.push_back( std::move( FissureBorder( *elIt, elemsByFacet )));
+ }
+ if ( fissure.empty() )
+ return;
+
+ // fill borderLinks
+
+ TBorderLinks borderLinks;
+
+ for ( size_t i = 0; i < fissure.size(); ++i )
+ {
+ fissure[i].AddSelfTo( borderLinks );
+ }
+
+ // get theAffectedElems
+
+ // unmark elements having nodes on the fissure, theAffectedElems elements will be marked
+ for ( size_t i = 0; i < fissure.size(); ++i )
+ for ( size_t j = 0; j < fissure[i]._nodes.size(); ++j )
+ {
+ SMESH_MeshAlgos::MarkElemNodes( fissure[i]._nodes[j]->GetInverseElementIterator(),
+ false, /*markElem=*/true );
+ }
+
+ std::vector<const SMDS_MeshNode *> facetNodes;
+ std::map< const SMDS_MeshNode*, FissureNormal > fissEdgeNodes2Norm;
+ boost::container::flat_set< const SMDS_MeshNode* > fissureNodes;
+
+ // choose a side of fissure
+ fissure[0].ChooseSide();
+ theAffectedElems.insert( fissure[0].GetMarkedElem() );
+
+ size_t nbCheckedBorders = 0;
+ while ( nbCheckedBorders < fissure.size() )
+ {
+ // find a FissureBorder to treat
+ FissureBorder* bord = 0;
+ for ( size_t i = 0; i < fissure.size() && !bord; ++i )
+ if ( fissure[i].GetMarkedElem() )
+ bord = & fissure[i];
+ for ( size_t i = 0; i < fissure.size() && !bord; ++i )
+ if ( fissure[i].NbSub() > 0 && fissure[i]._elems[0] )
+ {
+ bord = & fissure[i];
+ bord->ChooseSide();
+ theAffectedElems.insert( bord->GetMarkedElem() );
+ }
+ if ( !bord ) return;
+ ++nbCheckedBorders;
+
+ // treat FissureBorder's linked to bord
+ fissureNodes.clear();
+ fissureNodes.insert( bord->_nodes.begin(), bord->_nodes.end() );
+ for ( size_t i = 0; i < bord->NbSub(); ++i )
+ {
+ TBorderLinks::iterator l2b = bord->_mappedSubs[ i ];
+ if ( l2b == borderLinks.end() || l2b->second.empty() ) continue;
+ std::vector< FissureBorder* >& linkedBorders = l2b->second;
+ const SubBorder& sb = l2b->first;
+ const SMDS_MeshElement* bordElem = bord->GetMarkedElem();
+
+ if ( linkedBorders.size() == 1 ) // fissure edge reached, fill fissEdgeNodes2Norm
+ {
+ for ( int j = 0; j < sb._nbNodes; ++j )
+ fissEdgeNodes2Norm[ sb._nodes[j] ].Add( sb._nodes[j], *bord );
+ continue;
+ }
+
+ // add to theAffectedElems elems sharing nodes of a SubBorder and a node of bordElem
+ // until an elem adjacent to a neighbour FissureBorder is found
+ facetNodes.clear();
+ facetNodes.insert( facetNodes.end(), sb._nodes, sb._nodes + sb._nbNodes );
+ facetNodes.resize( sb._nbNodes + 1 );
+
+ while ( bordElem )
+ {
+ // check if bordElem is adjacent to a neighbour FissureBorder
+ for ( size_t j = 0; j < linkedBorders.size(); ++j )
+ {
+ FissureBorder* bord2 = linkedBorders[j];
+ if ( bord2 == bord ) continue;
+ if ( bordElem == bord2->_elems[0] || bordElem == bord2->_elems[1] )
+ bordElem = 0;
+ else
+ fissureNodes.insert( bord2->_nodes.begin(), bord2->_nodes.end() );
+ }
+ if ( !bordElem )
+ break;
+
+ // find the next bordElem
+ const SMDS_MeshElement* nextBordElem = 0;
+ for ( int iN = 0, nbN = bordElem->NbCornerNodes(); iN < nbN && !nextBordElem; ++iN )
+ {
+ const SMDS_MeshNode* n = bordElem->GetNode( iN );
+ if ( fissureNodes.count( n )) continue;
+
+ facetNodes[ sb._nbNodes ] = n;
+ elemsByFacet.clear();
+ if ( SMDS_Mesh::GetElementsByNodes( facetNodes, elemsByFacet ) > 1 )
+ {
+ for ( size_t iE = 0; iE < elemsByFacet.size(); ++iE )
+ if ( elemsByFacet[ iE ] != bordElem &&
+ !elemsByFacet[ iE ]->isMarked() )
+ {
+ theAffectedElems.insert( elemsByFacet[ iE ]);
+ elemsByFacet[ iE ]->setIsMarked( true );
+ if ( elemsByFacet[ iE ]->GetType() == bordElem->GetType() )
+ nextBordElem = elemsByFacet[ iE ];
+ }
+ }
+ }
+ bordElem = nextBordElem;
+
+ } // while ( bordElem )
+
+ linkedBorders.clear(); // not to treat this link any more
+
+ } // loop on SubBorder's of a FissureBorder
+
+ bord->Clear();
+
+ } // loop on FissureBorder's
+
+
+ // add elements sharing only one node of the fissure, except those sharing fissure edge nodes
+
+ // mark nodes of theAffectedElems
+ SMESH_MeshAlgos::MarkElemNodes( theAffectedElems.begin(), theAffectedElems.end(), true );
+
+ // unmark nodes of the fissure
+ elIt = theElemsOrNodes.begin();
+ if ( (*elIt)->GetType() == SMDSAbs_Node )
+ SMESH_MeshAlgos::MarkElems( elIt, theElemsOrNodes.end(), false );
+ else
+ SMESH_MeshAlgos::MarkElemNodes( elIt, theElemsOrNodes.end(), false );
+
+ std::vector< gp_XYZ > normVec;
+
+ // loop on nodes of the fissure, add elements having marked nodes
+ for ( elIt = theElemsOrNodes.begin(); elIt != theElemsOrNodes.end(); ++elIt )
+ {
+ const SMDS_MeshElement* e = (*elIt);
+ if ( e->GetType() != SMDSAbs_Node )
+ e->setIsMarked( true ); // avoid adding a fissure element
+
+ for ( int iN = 0, nbN = e->NbCornerNodes(); iN < nbN; ++iN )
+ {
+ const SMDS_MeshNode* n = e->GetNode( iN );
+ if ( fissEdgeNodes2Norm.count( n ))
+ continue;
+
+ SMDS_ElemIteratorPtr invIt = n->GetInverseElementIterator();
+ while ( invIt->more() )
+ {
+ const SMDS_MeshElement* eInv = invIt->next();
+ if ( eInv->isMarked() ) continue;
+ eInv->setIsMarked( true );
+
+ SMDS_ElemIteratorPtr nIt = eInv->nodesIterator();
+ while( nIt->more() )
+ if ( nIt->next()->isMarked())
+ {
+ theAffectedElems.insert( eInv );
+ SMESH_MeshAlgos::MarkElems( eInv->nodesIterator(), true );
+ n->setIsMarked( false );
+ break;
+ }
+ }
+ }
+ }
+
+ // add elements on the fissure edge
+ std::map< const SMDS_MeshNode*, FissureNormal >::iterator n2N;
+ for ( n2N = fissEdgeNodes2Norm.begin(); n2N != fissEdgeNodes2Norm.end(); ++n2N )
+ {
+ const SMDS_MeshNode* edgeNode = n2N->first;
+ const FissureNormal & normals = n2N->second;
+
+ SMDS_ElemIteratorPtr invIt = edgeNode->GetInverseElementIterator();
+ while ( invIt->more() )
+ {
+ const SMDS_MeshElement* eInv = invIt->next();
+ if ( eInv->isMarked() ) continue;
+ eInv->setIsMarked( true );
+
+ // classify eInv using normals
+ bool toAdd = normals.IsIn( edgeNode, eInv );
+ if ( toAdd ) // check if all nodes lie on the fissure edge
+ {
+ bool notOnEdge = false;
+ for ( int iN = 0, nbN = eInv->NbCornerNodes(); iN < nbN && !notOnEdge; ++iN )
+ notOnEdge = !fissEdgeNodes2Norm.count( eInv->GetNode( iN ));
+ toAdd = notOnEdge;
+ }
+ if ( toAdd )
+ {
+ theAffectedElems.insert( eInv );
+ }
+ }
+ }
+
+ return;
+ } // findAffectedElems()
+} // namespace
+
//================================================================================
/*!
* \brief Create elements equal (on same nodes) to given ones
SMDSAbs_ElementType type = SMDSAbs_All;
SMDS_ElemIteratorPtr elemIt;
- vector< const SMDS_MeshElement* > allElems;
if ( theElements.empty() )
{
if ( mesh->NbNodes() == 0 )
type = types[i];
break;
}
- // put all elements in the vector <allElems>
- allElems.reserve( mesh->GetMeshInfo().NbElements( type ));
elemIt = mesh->elementsIterator( type );
- while ( elemIt->more() )
- allElems.push_back( elemIt->next());
- elemIt = elemSetIterator( allElems );
}
else
{
elemIt = elemSetIterator( theElements );
}
+ // un-mark all elements to avoid duplicating just created elements
+ SMESH_MeshAlgos::MarkElems( mesh->elementsIterator( type ), false );
+
// duplicate elements
ElemFeatures elemType;
while ( elemIt->more() )
{
const SMDS_MeshElement* elem = elemIt->next();
- if ( elem->GetType() != type )
+ if ( elem->GetType() != type || elem->isMarked() )
continue;
elemType.Init( elem, /*basicOnly=*/false );
nodes.assign( elem->begin_nodes(), elem->end_nodes() );
- AddElement( nodes, elemType );
+ if ( const SMDS_MeshElement* newElem = AddElement( nodes, elemType ))
+ newElem->setIsMarked( true );
}
}
const TIDSortedElemSet& theNodesNot,
const TIDSortedElemSet& theAffectedElems )
{
- myLastCreatedElems.Clear();
- myLastCreatedNodes.Clear();
+ ClearLastCreated();
if ( theElems.size() == 0 )
return false;
for ( ; elemItr != theElems.end(); ++elemItr )
{
const SMDS_MeshElement* anElem = *elemItr;
- if (!anElem)
- continue;
+ // if (!anElem)
+ // continue;
// duplicate nodes to duplicate element
bool isDuplicate = false;
bool SMESH_MeshEditor::DoubleNodes( const std::list< int >& theListOfNodes,
const std::list< int >& theListOfModifiedElems )
{
- myLastCreatedElems.Clear();
- myLastCreatedNodes.Clear();
+ ClearLastCreated();
if ( theListOfNodes.size() == 0 )
return false;
std::list< int >::const_iterator aNodeIter;
for ( aNodeIter = theListOfNodes.begin(); aNodeIter != theListOfNodes.end(); ++aNodeIter )
{
- int aCurr = *aNodeIter;
- SMDS_MeshNode* aNode = (SMDS_MeshNode*)aMeshDS->FindNode( aCurr );
+ const SMDS_MeshNode* aNode = aMeshDS->FindNode( *aNodeIter );
if ( !aNode )
continue;
}
}
- // Create map of new nodes for modified elements
+ // Change nodes of elements
- std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> > anElemToNodes;
+ std::vector<const SMDS_MeshNode*> aNodeArr;
std::list< int >::const_iterator anElemIter;
- for ( anElemIter = theListOfModifiedElems.begin();
- anElemIter != theListOfModifiedElems.end(); ++anElemIter )
+ for ( anElemIter = theListOfModifiedElems.begin();
+ anElemIter != theListOfModifiedElems.end();
+ anElemIter++ )
{
- int aCurr = *anElemIter;
- SMDS_MeshElement* anElem = (SMDS_MeshElement*)aMeshDS->FindElement( aCurr );
+ const SMDS_MeshElement* anElem = aMeshDS->FindElement( *anElemIter );
if ( !anElem )
continue;
- vector<const SMDS_MeshNode*> aNodeArr( anElem->NbNodes() );
-
- SMDS_ElemIteratorPtr anIter = anElem->nodesIterator();
- int ind = 0;
- while ( anIter->more() )
+ aNodeArr.assign( anElem->begin_nodes(), anElem->end_nodes() );
+ for( size_t i = 0; i < aNodeArr.size(); ++i )
{
- SMDS_MeshNode* aCurrNode = (SMDS_MeshNode*)anIter->next();
- if ( aCurr && anOldNodeToNewNode.find( aCurrNode ) != anOldNodeToNewNode.end() )
- {
- const SMDS_MeshNode* aNewNode = anOldNodeToNewNode[ aCurrNode ];
- aNodeArr[ ind++ ] = aNewNode;
- }
- else
- aNodeArr[ ind++ ] = aCurrNode;
- }
- anElemToNodes[ anElem ] = aNodeArr;
- }
-
- // Change nodes of elements
-
- std::map< SMDS_MeshElement*, vector<const SMDS_MeshNode*> >::iterator
- anElemToNodesIter = anElemToNodes.begin();
- for ( ; anElemToNodesIter != anElemToNodes.end(); ++anElemToNodesIter )
- {
- const SMDS_MeshElement* anElem = anElemToNodesIter->first;
- vector<const SMDS_MeshNode*> aNodeArr = anElemToNodesIter->second;
- if ( anElem )
- {
- aMeshDS->ChangeElementNodes( anElem, &aNodeArr[ 0 ], anElem->NbNodes() );
+ std::map< const SMDS_MeshNode*, const SMDS_MeshNode* >::iterator n2n =
+ anOldNodeToNewNode.find( aNodeArr[ i ]);
+ if ( n2n != anOldNodeToNewNode.end() )
+ aNodeArr[ i ] = n2n->second;
}
+ aMeshDS->ChangeElementNodes( anElem, &aNodeArr[ 0 ], aNodeArr.size() );
}
return true;
{
gp_XYZ centerXYZ (0, 0, 0);
SMDS_ElemIteratorPtr aNodeItr = theElem->nodesIterator();
- while (aNodeItr->more())
- centerXYZ += SMESH_TNodeXYZ(cast2Node( aNodeItr->next()));
+ while ( aNodeItr->more() )
+ centerXYZ += SMESH_NodeXYZ( aNodeItr->next() );
gp_Pnt aPnt = centerXYZ / theElem->NbNodes();
theClassifier.Perform(aPnt, theTol);
(select elements with a gravity center on the side given by faces normals).
This mode (null shape) is faster, but works only when theElems are faces, with coherents orientations.
The replicated nodes should be associated to affected elements.
- \return groups of affected elements
+ \return true
\sa DoubleNodeElemGroupsInRegion()
- */
+*/
//================================================================================
bool SMESH_MeshEditor::AffectedElemGroupsInRegion( const TIDSortedElemSet& theElems,
{
if ( theShape.IsNull() )
{
- std::set<const SMDS_MeshNode*> alreadyCheckedNodes;
- std::set<const SMDS_MeshElement*> alreadyCheckedElems;
- std::set<const SMDS_MeshElement*> edgesToCheck;
- alreadyCheckedNodes.clear();
- alreadyCheckedElems.clear();
- edgesToCheck.clear();
-
- // --- iterates on elements to be replicated and get elements by back references from their nodes
-
- TIDSortedElemSet::const_iterator elemItr = theElems.begin();
- for ( ; elemItr != theElems.end(); ++elemItr )
- {
- SMDS_MeshElement* anElem = (SMDS_MeshElement*)*elemItr;
- if (!anElem || (anElem->GetType() != SMDSAbs_Face))
- continue;
- gp_XYZ normal;
- SMESH_MeshAlgos::FaceNormal( anElem, normal, /*normalized=*/true );
- std::set<const SMDS_MeshNode*> nodesElem;
- nodesElem.clear();
- SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator();
- while ( nodeItr->more() )
- {
- const SMDS_MeshNode* aNode = cast2Node(nodeItr->next());
- nodesElem.insert(aNode);
- }
- std::set<const SMDS_MeshNode*>::iterator nodit = nodesElem.begin();
- for (; nodit != nodesElem.end(); nodit++)
- {
- const SMDS_MeshNode* aNode = *nodit;
- if ( !aNode || theNodesNot.find(aNode) != theNodesNot.end() )
- continue;
- if (alreadyCheckedNodes.find(aNode) != alreadyCheckedNodes.end())
- continue;
- alreadyCheckedNodes.insert(aNode);
- SMDS_ElemIteratorPtr backElemItr = aNode->GetInverseElementIterator();
- while ( backElemItr->more() )
- {
- const SMDS_MeshElement* curElem = backElemItr->next();
- if (alreadyCheckedElems.find(curElem) != alreadyCheckedElems.end())
- continue;
- if (theElems.find(curElem) != theElems.end())
- continue;
- alreadyCheckedElems.insert(curElem);
- double x=0, y=0, z=0;
- int nb = 0;
- SMDS_ElemIteratorPtr nodeItr2 = curElem->nodesIterator();
- while ( nodeItr2->more() )
- {
- const SMDS_MeshNode* anotherNode = cast2Node(nodeItr2->next());
- x += anotherNode->X();
- y += anotherNode->Y();
- z += anotherNode->Z();
- nb++;
- }
- gp_XYZ p;
- p.SetCoord( x/nb -aNode->X(),
- y/nb -aNode->Y(),
- z/nb -aNode->Z() );
- if (normal*p > 0)
- {
- theAffectedElems.insert( curElem );
- }
- else if (curElem->GetType() == SMDSAbs_Edge)
- edgesToCheck.insert(curElem);
- }
- }
- }
- // --- add also edges lying on the set of faces (all nodes in alreadyCheckedNodes)
- std::set<const SMDS_MeshElement*>::iterator eit = edgesToCheck.begin();
- for( ; eit != edgesToCheck.end(); eit++)
- {
- bool onside = true;
- const SMDS_MeshElement* anEdge = *eit;
- SMDS_ElemIteratorPtr nodeItr = anEdge->nodesIterator();
- while ( nodeItr->more() )
- {
- const SMDS_MeshNode* aNode = cast2Node(nodeItr->next());
- if (alreadyCheckedNodes.find(aNode) == alreadyCheckedNodes.end())
- {
- onside = false;
- break;
- }
- }
- if (onside)
- {
- theAffectedElems.insert(anEdge);
- }
- }
+ findAffectedElems( theElems, theAffectedElems );
}
else
{
TIDSortedElemSet::const_iterator elemItr = theElems.begin();
for ( ; elemItr != theElems.end(); ++elemItr )
{
- SMDS_MeshElement* anElem = (SMDS_MeshElement*)*elemItr;
- if (!anElem)
- continue;
+ SMDS_MeshElement* anElem = (SMDS_MeshElement*)*elemItr;
SMDS_ElemIteratorPtr nodeItr = anElem->nodesIterator();
while ( nodeItr->more() )
{
{
const SMDS_MeshElement* curElem = backElemItr->next();
if ( curElem && theElems.find(curElem) == theElems.end() &&
- ( bsc3d.get() ?
- isInside( curElem, *bsc3d, aTol ) :
- isInside( curElem, *aFaceClassifier, aTol )))
+ ( bsc3d.get() ?
+ isInside( curElem, *bsc3d, aTol ) :
+ isInside( curElem, *aFaceClassifier, aTol )))
theAffectedElems.insert( curElem );
}
}
default:;
}
}
+
+namespace // utils for MakePolyLine
+{
+ //================================================================================
+ /*!
+ * \brief Sequence of found points and a current point data
+ */
+ struct Path
+ {
+ std::vector< gp_XYZ > myPoints;
+ double myLength;
+
+ int mySrcPntInd; //!< start point index
+ const SMDS_MeshElement* myFace;
+ SMESH_NodeXYZ myNode1;
+ SMESH_NodeXYZ myNode2;
+ int myNodeInd1;
+ int myNodeInd2;
+ double myDot1;
+ double myDot2;
+ TIDSortedElemSet myElemSet, myAvoidSet;
+
+ Path(): myLength(0.0), myFace(0) {}
+
+ bool SetCutAtCorner( const SMESH_NodeXYZ& cornerNode,
+ const SMDS_MeshElement* face,
+ const gp_XYZ& plnNorm,
+ const gp_XYZ& plnOrig );
+
+ void AddPoint( const gp_XYZ& p );
+
+ bool Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig );
+
+ bool ReachSamePoint( const Path& other );
+
+ static void Remove( std::vector< Path > & paths, size_t& i );
+ };
+
+ //================================================================================
+ /*!
+ * \brief Return true if this Path meats another
+ */
+ //================================================================================
+
+ bool Path::ReachSamePoint( const Path& other )
+ {
+ return ( mySrcPntInd != other.mySrcPntInd &&
+ myFace == other.myFace );
+ }
+
+ //================================================================================
+ /*!
+ * \brief Remove a path from a vector
+ */
+ //================================================================================
+
+ void Path::Remove( std::vector< Path > & paths, size_t& i )
+ {
+ if ( paths.size() > 1 )
+ {
+ size_t j = paths.size() - 1; // last item to be removed
+ if ( i < j )
+ {
+ paths[ i ].myPoints.swap( paths[ j ].myPoints );
+ paths[ i ].myLength = paths[ j ].myLength;
+ paths[ i ].mySrcPntInd = paths[ j ].mySrcPntInd;
+ paths[ i ].myFace = paths[ j ].myFace;
+ paths[ i ].myNode1 = paths[ j ].myNode1;
+ paths[ i ].myNode2 = paths[ j ].myNode2;
+ paths[ i ].myNodeInd1 = paths[ j ].myNodeInd1;
+ paths[ i ].myNodeInd2 = paths[ j ].myNodeInd2;
+ paths[ i ].myDot1 = paths[ j ].myDot1;
+ paths[ i ].myDot2 = paths[ j ].myDot2;
+ }
+ }
+ paths.pop_back();
+ if ( i > 0 )
+ --i;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Store a point that is at a node of a face if the face is intersected by plane.
+ * Return false if the node is a sole intersection point of the face and the plane
+ */
+ //================================================================================
+
+ bool Path::SetCutAtCorner( const SMESH_NodeXYZ& cornerNode,
+ const SMDS_MeshElement* face,
+ const gp_XYZ& plnNorm,
+ const gp_XYZ& plnOrig )
+ {
+ if ( face == myFace )
+ return false;
+ myNodeInd1 = face->GetNodeIndex( cornerNode._node );
+ myNodeInd2 = ( myNodeInd1 + 1 ) % face->NbCornerNodes();
+ int ind3 = ( myNodeInd1 + 2 ) % face->NbCornerNodes();
+ myNode1.Set( face->GetNode( ind3 ));
+ myNode2.Set( face->GetNode( myNodeInd2 ));
+
+ myDot1 = plnNorm * ( myNode1 - plnOrig );
+ myDot2 = plnNorm * ( myNode2 - plnOrig );
+
+ bool ok = ( myDot1 * myDot2 < 0 );
+ if ( !ok && myDot1 * myDot2 == 0 )
+ {
+ ok = ( myDot1 != myDot2 );
+ if ( ok && myFace )
+ ok = ( myFace->GetNodeIndex(( myDot1 == 0 ? myNode1 : myNode2 )._node ) < 0 );
+ }
+ if ( ok )
+ {
+ myFace = face;
+ myDot1 = 0;
+ AddPoint( cornerNode );
+ }
+ return ok;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Store a point and update myLength
+ */
+ //================================================================================
+
+ void Path::AddPoint( const gp_XYZ& p )
+ {
+ if ( !myPoints.empty() )
+ myLength += ( p - myPoints.back() ).Modulus();
+ else
+ myLength = 0;
+ myPoints.push_back( p );
+ }
+
+ //================================================================================
+ /*!
+ * \brief Try to find the next point
+ * \param [in] plnNorm - cutting plane normal
+ * \param [in] plnOrig - cutting plane origin
+ */
+ //================================================================================
+
+ bool Path::Extend( const gp_XYZ& plnNorm, const gp_XYZ& plnOrig )
+ {
+ int nodeInd3 = ( myNodeInd1 + 1 ) % myFace->NbCornerNodes();
+ if ( myNodeInd2 == nodeInd3 )
+ nodeInd3 = ( myNodeInd1 + 2 ) % myFace->NbCornerNodes();
+
+ SMESH_NodeXYZ node3 = myFace->GetNode( nodeInd3 );
+ double dot3 = plnNorm * ( node3 - plnOrig );
+
+ if ( dot3 * myDot1 < 0. )
+ {
+ myNode2 = node3;
+ myNodeInd2 = nodeInd3;
+ myDot2 = dot3;
+ }
+ else if ( dot3 * myDot2 < 0. )
+ {
+ myNode1 = node3;
+ myNodeInd1 = nodeInd3;
+ myDot1 = dot3;
+ }
+ else if ( dot3 == 0. )
+ {
+ SMDS_ElemIteratorPtr fIt = node3._node->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() )
+ if ( SetCutAtCorner( node3, fIt->next(), plnNorm, plnOrig ))
+ return true;
+ return false;
+ }
+ else if ( myDot2 == 0. )
+ {
+ SMESH_NodeXYZ node2 = myNode2; // copy as myNode2 changes in SetCutAtCorner()
+ SMDS_ElemIteratorPtr fIt = node2._node->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() )
+ if ( SetCutAtCorner( node2, fIt->next(), plnNorm, plnOrig ))
+ return true;
+ return false;
+ }
+
+ double r = Abs( myDot1 / ( myDot2 - myDot1 ));
+ AddPoint( myNode1 * ( 1 - r ) + myNode2 * r );
+
+ myAvoidSet.clear();
+ myAvoidSet.insert( myFace );
+ myFace = SMESH_MeshAlgos::FindFaceInSet( myNode1._node, myNode2._node,
+ myElemSet, myAvoidSet,
+ &myNodeInd1, &myNodeInd2 );
+ return myFace;
+ }
+
+ //================================================================================
+ /*!
+ * \brief Compute a path between two points of PolySegment
+ */
+ struct PolyPathCompute
+ {
+ SMESH_MeshEditor::TListOfPolySegments& mySegments; //!< inout PolySegment's
+ std::vector< Path >& myPaths; //!< path of each of segments to compute
+ SMESH_Mesh* myMesh;
+ mutable std::vector< std::string > myErrors;
+
+ PolyPathCompute( SMESH_MeshEditor::TListOfPolySegments& theSegments,
+ std::vector< Path >& thePaths,
+ SMESH_Mesh* theMesh):
+ mySegments( theSegments ),
+ myPaths( thePaths ),
+ myMesh( theMesh ),
+ myErrors( theSegments.size() )
+ {
+ }
+#undef SMESH_CAUGHT
+#define SMESH_CAUGHT myErrors[i] =
+ void operator() ( const int i ) const
+ {
+ SMESH_TRY;
+ const_cast< PolyPathCompute* >( this )->Compute( i );
+ SMESH_CATCH( SMESH::returnError );
+ }
+#undef SMESH_CAUGHT
+ //================================================================================
+ /*!
+ * \brief Compute a path of a given segment
+ */
+ //================================================================================
+
+ void Compute( const int iSeg )
+ {
+ SMESH_MeshEditor::PolySegment& polySeg = mySegments[ iSeg ];
+
+ // get a cutting plane
+
+ gp_XYZ p1 = SMESH_NodeXYZ( polySeg.myNode1[0] );
+ gp_XYZ p2 = SMESH_NodeXYZ( polySeg.myNode1[1] );
+ if ( polySeg.myNode2[0] ) p1 = 0.5 * ( p1 + SMESH_NodeXYZ( polySeg.myNode2[0] ));
+ if ( polySeg.myNode2[1] ) p2 = 0.5 * ( p2 + SMESH_NodeXYZ( polySeg.myNode2[1] ));
+
+ gp_XYZ plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
+ gp_XYZ plnOrig = p2;
+
+ // find paths connecting the 2 end points of polySeg
+
+ std::vector< Path > paths; paths.reserve(10);
+
+ // initialize paths
+
+ for ( int iP = 0; iP < 2; ++iP ) // loop on the polySeg end points
+ {
+ Path path;
+ path.mySrcPntInd = iP;
+ size_t nbPaths = paths.size();
+
+ if ( polySeg.myNode2[ iP ] && polySeg.myNode2[ iP ] != polySeg.myNode1[ iP ] )
+ {
+ while (( path.myFace = SMESH_MeshAlgos::FindFaceInSet( polySeg.myNode1[ iP ],
+ polySeg.myNode2[ iP ],
+ path.myElemSet,
+ path.myAvoidSet,
+ &path.myNodeInd1,
+ &path.myNodeInd2 )))
+ {
+ path.myNode1.Set( polySeg.myNode1[ iP ]);
+ path.myNode2.Set( polySeg.myNode2[ iP ]);
+ path.myDot1 = plnNorm * ( path.myNode1 - plnOrig );
+ path.myDot2 = plnNorm * ( path.myNode2 - plnOrig );
+ path.myPoints.clear();
+ path.AddPoint( 0.5 * ( path.myNode1 + path.myNode2 ));
+ path.myAvoidSet.insert( path.myFace );
+ paths.push_back( path );
+ }
+ if ( nbPaths == paths.size() )
+ throw SALOME_Exception ( SMESH_Comment("No face edge found by point ") << iP+1
+ << " in a PolySegment " << iSeg );
+ }
+ else // an end point is at node
+ {
+ std::set<const SMDS_MeshNode* > nodes;
+ SMDS_ElemIteratorPtr fIt = polySeg.myNode1[ iP ]->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() )
+ {
+ path.myPoints.clear();
+ if ( path.SetCutAtCorner( polySeg.myNode1[ iP ], fIt->next(), plnNorm, plnOrig ))
+ {
+ if (( path.myDot1 * path.myDot2 != 0 ) ||
+ ( nodes.insert( path.myDot1 == 0 ? path.myNode1._node : path.myNode2._node ).second ))
+ paths.push_back( path );
+ }
+ }
+ }
+
+ // look for a one-segment path
+ for ( size_t i = 0; i < nbPaths; ++i )
+ for ( size_t j = nbPaths; j < paths.size(); ++j )
+ if ( paths[i].myFace == paths[j].myFace )
+ {
+ myPaths[ iSeg ].myPoints.push_back( paths[i].myPoints[0] );
+ myPaths[ iSeg ].myPoints.push_back( paths[j].myPoints[0] );
+ paths.clear();
+ }
+ }
+
+ // extend paths
+
+ myPaths[ iSeg ].myLength = 1e100;
+
+ while ( paths.size() >= 2 )
+ {
+ for ( size_t i = 0; i < paths.size(); ++i )
+ {
+ Path& path = paths[ i ];
+ if ( !path.Extend( plnNorm, plnOrig ) || // path reached a mesh boundary
+ path.myLength > myPaths[ iSeg ].myLength ) // path is longer than others
+ {
+ Path::Remove( paths, i );
+ continue;
+ }
+
+ // join paths that reach same point
+ for ( size_t j = 0; j < paths.size(); ++j )
+ {
+ if ( i != j && paths[i].ReachSamePoint( paths[j] ))
+ {
+ double distLast = ( paths[i].myPoints.back() - paths[j].myPoints.back() ).Modulus();
+ double fullLength = ( paths[i].myLength + paths[j].myLength + distLast );
+ if ( fullLength < myPaths[ iSeg ].myLength )
+ {
+ myPaths[ iSeg ].myLength = fullLength;
+ std::vector< gp_XYZ > & allPoints = myPaths[ iSeg ].myPoints;
+ allPoints.swap( paths[i].myPoints );
+ allPoints.insert( allPoints.end(),
+ paths[j].myPoints.rbegin(),
+ paths[j].myPoints.rend() );
+ }
+ Path::Remove( paths, i );
+ Path::Remove( paths, j );
+ }
+ }
+ }
+ if ( !paths.empty() && (int) paths[0].myPoints.size() > myMesh->NbFaces() )
+ throw SALOME_Exception(LOCALIZED( "Infinite loop in MakePolyLine()"));
+ }
+
+ if ( myPaths[ iSeg ].myPoints.empty() )
+ throw SALOME_Exception( SMESH_Comment("Can't find a full path for PolySegment #") << iSeg );
+
+ } // PolyPathCompute::Compute()
+
+ }; // struct PolyPathCompute
+
+} // namespace
+
+//=======================================================================
+//function : MakePolyLine
+//purpose : Create a polyline consisting of 1D mesh elements each lying on a 2D element of
+// the initial mesh
+//=======================================================================
+
+void SMESH_MeshEditor::MakePolyLine( TListOfPolySegments& theSegments,
+ SMESHDS_Group* theGroup,
+ SMESH_ElementSearcher* theSearcher)
+{
+ std::vector< Path > segPaths( theSegments.size() ); // path of each of segments
+
+ SMESH_ElementSearcher* searcher = theSearcher;
+ SMESHUtils::Deleter<SMESH_ElementSearcher> delSearcher;
+ if ( !searcher )
+ {
+ searcher = SMESH_MeshAlgos::GetElementSearcher( *GetMeshDS() );
+ delSearcher._obj = searcher;
+ }
+
+ // get cutting planes
+
+ std::vector< bool > isVectorOK( theSegments.size(), true );
+ const double planarCoef = 0.333; // plane height in planar case
+
+ for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
+ {
+ PolySegment& polySeg = theSegments[ iSeg ];
+
+ gp_XYZ p1 = SMESH_NodeXYZ( polySeg.myNode1[0] );
+ gp_XYZ p2 = SMESH_NodeXYZ( polySeg.myNode1[1] );
+ if ( polySeg.myNode2[0] ) p1 = 0.5 * ( p1 + SMESH_NodeXYZ( polySeg.myNode2[0] ));
+ if ( polySeg.myNode2[1] ) p2 = 0.5 * ( p2 + SMESH_NodeXYZ( polySeg.myNode2[1] ));
+
+ gp_XYZ plnNorm = ( p1 - p2 ) ^ polySeg.myVector.XYZ();
+
+ isVectorOK[ iSeg ] = ( plnNorm.Modulus() > std::numeric_limits<double>::min() );
+ if ( !isVectorOK[ iSeg ])
+ {
+ gp_XYZ pMid = 0.5 * ( p1 + p2 );
+ const SMDS_MeshElement* face;
+ polySeg.myMidProjPoint = searcher->Project( pMid, SMDSAbs_Face, &face );
+ polySeg.myVector = polySeg.myMidProjPoint.XYZ() - pMid;
+
+ gp_XYZ faceNorm;
+ SMESH_MeshAlgos::FaceNormal( face, faceNorm );
+
+ if ( polySeg.myVector.Magnitude() < Precision::Confusion() ||
+ polySeg.myVector * faceNorm < Precision::Confusion() )
+ {
+ polySeg.myVector = faceNorm;
+ polySeg.myMidProjPoint = pMid + faceNorm * ( p1 - p2 ).Modulus() * planarCoef;
+ }
+ }
+ else
+ {
+ polySeg.myVector = plnNorm ^ ( p1 - p2 );
+ }
+ }
+
+ // assure that inverse elements are constructed, avoid their concurrent building in threads
+ GetMeshDS()->nodesIterator()->next()->NbInverseElements();
+
+ // find paths
+
+ PolyPathCompute algo( theSegments, segPaths, myMesh );
+ OSD_Parallel::For( 0, theSegments.size(), algo, theSegments.size() == 1 );
+
+ for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
+ if ( !algo.myErrors[ iSeg ].empty() )
+ throw SALOME_Exception( algo.myErrors[ iSeg ].c_str() );
+
+ // create an 1D mesh
+
+ const SMDS_MeshNode *n, *nPrev = 0;
+ SMESHDS_Mesh* mesh = GetMeshDS();
+
+ for ( size_t iSeg = 0; iSeg < theSegments.size(); ++iSeg )
+ {
+ const Path& path = segPaths[iSeg];
+ if ( path.myPoints.size() < 2 )
+ continue;
+
+ double tol = path.myLength / path.myPoints.size() / 1000.;
+ if ( !nPrev || ( SMESH_NodeXYZ( nPrev ) - path.myPoints[0] ).SquareModulus() > tol*tol )
+ {
+ nPrev = mesh->AddNode( path.myPoints[0].X(), path.myPoints[0].Y(), path.myPoints[0].Z() );
+ myLastCreatedNodes.Append( nPrev );
+ }
+ for ( size_t iP = 1; iP < path.myPoints.size(); ++iP )
+ {
+ n = mesh->AddNode( path.myPoints[iP].X(), path.myPoints[iP].Y(), path.myPoints[iP].Z() );
+ myLastCreatedNodes.Append( n );
+
+ const SMDS_MeshElement* elem = mesh->AddEdge( nPrev, n );
+ myLastCreatedElems.Append( elem );
+ if ( theGroup )
+ theGroup->Add( elem );
+
+ nPrev = n;
+ }
+
+ // return a vector
+
+ gp_XYZ pMid = 0.5 * ( path.myPoints[0] + path.myPoints.back() );
+ if ( isVectorOK[ iSeg ])
+ {
+ // find the most distance point of a path
+ double maxDist = 0;
+ for ( size_t iP = 1; iP < path.myPoints.size(); ++iP )
+ {
+ double dist = Abs( theSegments[iSeg].myVector * ( path.myPoints[iP] - path.myPoints[0] ));
+ if ( dist > maxDist )
+ {
+ maxDist = dist;
+ theSegments[iSeg].myMidProjPoint = path.myPoints[iP];
+ }
+ }
+ if ( maxDist < Precision::Confusion() ) // planar case
+ theSegments[iSeg].myMidProjPoint =
+ pMid + theSegments[iSeg].myVector.XYZ().Normalized() * path.myLength * planarCoef;
+ }
+ theSegments[iSeg].myVector = gp_Vec( pMid, theSegments[iSeg].myMidProjPoint );
+ }
+
+ return;
+}
