From: Konstantin Leontev Date: Wed, 24 Jul 2024 17:06:11 +0000 (+0100) Subject: [bos #42217][EDF 28921] Horseshoe with bodyfitting. X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=9f7018e42a582527ec842866be90bb10108ef88b;p=modules%2Fsmesh.git [bos #42217][EDF 28921] Horseshoe with bodyfitting. Intermediate commit - parallel compute is disabled, some unrelated code commented out for debug! Hexahedron::compute() was splitted to smaller methods. Added debug output. SMDS_MeshElement::Print() now is a virtual method, otherwise only parent method was called. --- diff --git a/src/SMDS/SMDS_FaceOfNodes.hxx b/src/SMDS/SMDS_FaceOfNodes.hxx index 1ccd14f6a..225e620ab 100644 --- a/src/SMDS/SMDS_FaceOfNodes.hxx +++ b/src/SMDS/SMDS_FaceOfNodes.hxx @@ -32,7 +32,7 @@ class SMDS_EXPORT SMDS_FaceOfNodes: public SMDS_CellOfNodes { public: - void Print(std::ostream & OS) const; + virtual void Print(std::ostream & OS) const override; SMDS_FaceOfNodes(const SMDS_MeshNode* node1, const SMDS_MeshNode* node2, const SMDS_MeshNode* node3); diff --git a/src/SMDS/SMDS_MeshElement.hxx b/src/SMDS/SMDS_MeshElement.hxx index eeec4a1f7..1d44ecf75 100644 --- a/src/SMDS/SMDS_MeshElement.hxx +++ b/src/SMDS/SMDS_MeshElement.hxx @@ -142,7 +142,7 @@ public: SMDS_Mesh* GetMesh() const; - void Print(std::ostream & OS) const; + virtual void Print(std::ostream & OS) const; friend SMDS_EXPORT std::ostream & operator <<(std::ostream & OS, const SMDS_MeshElement *); friend class SMDS_ElementFactory; diff --git a/src/SMDS/SMDS_MeshNode.hxx b/src/SMDS/SMDS_MeshNode.hxx index 99549e937..e2aa11270 100644 --- a/src/SMDS/SMDS_MeshNode.hxx +++ b/src/SMDS/SMDS_MeshNode.hxx @@ -65,7 +65,7 @@ class SMDS_EXPORT SMDS_MeshNode: public SMDS_MeshElement virtual bool IsMediumNode(const SMDS_MeshNode* /*node*/) const { return false; } virtual int NbCornerNodes() const { return 1; } - void Print(std::ostream & OS) const; + virtual void Print(std::ostream & OS) const override; private: diff --git a/src/SMDS/SMDS_PolygonalFaceOfNodes.hxx b/src/SMDS/SMDS_PolygonalFaceOfNodes.hxx index 7dbff1211..8b3d4fa5d 100644 --- a/src/SMDS/SMDS_PolygonalFaceOfNodes.hxx +++ b/src/SMDS/SMDS_PolygonalFaceOfNodes.hxx @@ -47,7 +47,7 @@ class SMDS_EXPORT SMDS_PolygonalFaceOfNodes : public SMDS_CellOfNodes virtual int NbEdges() const; virtual int NbFaces() const; - virtual void Print (std::ostream & OS) const; + virtual void Print (std::ostream & OS) const override; virtual const SMDS_MeshNode* GetNode(const int ind) const; diff --git a/src/SMDS/SMDS_VolumeOfNodes.hxx b/src/SMDS/SMDS_VolumeOfNodes.hxx index 36a515c5d..7098cd623 100644 --- a/src/SMDS/SMDS_VolumeOfNodes.hxx +++ b/src/SMDS/SMDS_VolumeOfNodes.hxx @@ -62,7 +62,7 @@ class SMDS_EXPORT SMDS_VolumeOfNodes: public SMDS_CellOfNodes const int nbNodes); ~SMDS_VolumeOfNodes(); - void Print(std::ostream & OS) const; + virtual void Print(std::ostream & OS) const override; int NbFaces() const; int NbNodes() const; int NbEdges() const; diff --git a/src/StdMeshers/StdMeshers_Cartesian_3D.cxx b/src/StdMeshers/StdMeshers_Cartesian_3D.cxx index ed41cb70e..4a0f7070e 100644 --- a/src/StdMeshers/StdMeshers_Cartesian_3D.cxx +++ b/src/StdMeshers/StdMeshers_Cartesian_3D.cxx @@ -773,6 +773,16 @@ namespace if ( node ) _node = _intPoint->_node = node; } + + friend std::ostream& operator<<(std::ostream& os, const _Node& node) + { + if (node._node) + node._node->Print(os); + else + os << "mesh node at hexahedron corner is null" << '\n'; + + return os; + } }; // -------------------------------------------------------------------------------- struct _Link // link connecting two _Node's @@ -783,6 +793,21 @@ namespace vector< _Node* > _fIntNodes; // _Node's at _fIntPoints vector< _Link > _splits; _Link(): _faces{ 0, 0 } {} + + friend std::ostream& operator<<(std::ostream& os, const _Link& link) + { + os << "Link:\n"; + + for (int i = 0; i < 2; ++i) + { + if (link._nodes[i]) + os << *link._nodes[i]; + else + os << "link node with index " << i << " is null" << '\n'; + } + + return os; + } }; // -------------------------------------------------------------------------------- struct _OrientedLink @@ -852,6 +877,16 @@ namespace } } } + + friend std::ostream& operator<<(std::ostream& os, const _OrientedLink& link) + { + if (link._link) + os << "Oriented " << *link._link; + else + os << "Oriented link is null" << '\n'; + + return os; + } }; // -------------------------------------------------------------------------------- struct _SplitIterator //! set to _hexLinks splits on one side of INTERNAL FACEs @@ -1114,6 +1149,26 @@ namespace TID2Nb::iterator id2nb = id2nbMap.insert( s0 ).first; id2nb->second++; } + + // Called by compute() + //================================================================================ + _Face* createPolygon(const SMESH_Block::TShapeID& name); + void closePolygon( + _Face* polygon, _Node* n2, _Node* nFirst, _Face& quad, vector<_Node*>& chainNodes, size_t& nbUsedEdgeNodes, _Face* prevPolyg); + void connectPolygonLinks( + const Solid* solid, _Face* polygon, _Face& quad, vector<_Node*>& chainNodes, std::vector<_OrientedLink>& splits, const bool toCheckSideDivision); + bool collectSplits(std::vector<_OrientedLink>& splits, const _Face& quad, _Face* polygon, int quadIndex); + std::set getConcaveFaces(const Solid* solid); + std::vector<_Node*> getChainNodes(const Solid* solid, const IsInternalFlag intFlag); + void clearHexUsedInFace(); + void clearIntUsedInFace(); + void addPolygonsToLinks(); + std::vector<_OrientedLink*> getFreeLinks(); + int notUsedIntersectionNodesToVInt(); + bool createPolygons(const bool hasEdgeIntersections, const IsInternalFlag intFlag); + void setNamesForNoNamePolygons(); + bool createVolume(const Solid* solid); + //================================================================================ }; // class Hexahedron #ifdef WITH_TBB @@ -2974,235 +3029,323 @@ namespace //================================================================================ /*! - * \brief Compute mesh volumes resulted from intersection of the Hexahedron + * \brief Collected faces can be used to avoid connecting nodes laying on them */ - bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag ) + std::set Hexahedron::getConcaveFaces(const Solid* solid) { - _polygons.clear(); - _polygons.reserve( 20 ); - - for ( int iN = 0; iN < 8; ++iN ) - _hexNodes[iN]._usedInFace = 0; + MESSAGE("Collect concave faces..."); - if ( intFlag & IS_CUT_BY_INTERNAL_FACE && !_grid->_toAddEdges ) // Issue #19913 - preventVolumesOverlapping(); + if (!solid->HasConcaveVertex()) + { + MESSAGE("There's no concave faces here. Return."); + return {}; + } - std::set< TGeomID > concaveFaces; // to avoid connecting nodes laying on them + std::set< TGeomID > concaveFaces; - if ( solid->HasConcaveVertex() ) + for ( const E_IntersectPoint* ip : _eIntPoints ) { - for ( const E_IntersectPoint* ip : _eIntPoints ) + if ( const ConcaveFace* cf = solid->GetConcave( ip->_shapeID )) + if ( this->hasEdgesAround( cf )) + concaveFaces.insert( cf->_concaveFace ); + } + if ( concaveFaces.empty() || concaveFaces.size() * 3 < _eIntPoints.size() ) + for ( const _Node& hexNode: _hexNodes ) { - if ( const ConcaveFace* cf = solid->GetConcave( ip->_shapeID )) - if ( this->hasEdgesAround( cf )) - concaveFaces.insert( cf->_concaveFace ); + if ( hexNode._node && hexNode._intPoint && hexNode._intPoint->_faceIDs.size() >= 3 ) + if ( const ConcaveFace* cf = solid->GetConcave( hexNode._node->GetShapeID() )) + if ( this->hasEdgesAround( cf )) + concaveFaces.insert( cf->_concaveFace ); } - if ( concaveFaces.empty() || concaveFaces.size() * 3 < _eIntPoints.size() ) - for ( const _Node& hexNode: _hexNodes ) - { - if ( hexNode._node && hexNode._intPoint && hexNode._intPoint->_faceIDs.size() >= 3 ) - if ( const ConcaveFace* cf = solid->GetConcave( hexNode._node->GetShapeID() )) - if ( this->hasEdgesAround( cf )) - concaveFaces.insert( cf->_concaveFace ); - } + + MESSAGE("Collected concave faces: " << concaveFaces.size()); + return concaveFaces; + } + + //================================================================================ + /*! + * \brief Creates a new polygone with a given name + */ + Hexahedron::_Face* Hexahedron::createPolygon(const SMESH_Block::TShapeID& name) + { + MESSAGE("Create a polygon with a name: " << name); + _polygons.resize( _polygons.size() + 1 ); + _Face* polygon = &_polygons.back(); + polygon->_polyLinks.reserve( 20 ); + polygon->_name = name; + + return polygon; + } + + //================================================================================ + /*! + * \brief Collects split links on 4 sides of a quadrangle. + * Returns false if the quad on FACE is not split. + */ + bool Hexahedron::collectSplits(std::vector<_OrientedLink>& splits, const _Face& quad, _Face* polygon, int quadIndex) + { + MESSAGE("Collect splits..."); + + splits.clear(); + for ( int iE = 0; iE < 4; ++iE ) // loop on 4 sides of a quadrangle + for ( size_t iS = 0; iS < quad._links[ iE ].NbResultLinks(); ++iS ) + splits.push_back( quad._links[ iE ].ResultLink( iS )); + + if ( splits.size() == 4 && + isQuadOnFace( quadIndex )) // check if a quad on FACE is not split + { + MESSAGE("The quad on FACE is not split. Swap splits with polygon links."); + + polygon->_links.swap( splits ); + return false; } - // Create polygons from quadrangles - // -------------------------------- + MESSAGE("Collected splits: " << splits.size()); + return true; + } - vector< _OrientedLink > splits; - vector<_Node*> chainNodes; - _Face* coplanarPolyg; + void Hexahedron::closePolygon( + _Face* polygon, _Node* n2, _Node* nFirst, _Face& quad, vector<_Node*>& chainNodes, size_t& nbUsedEdgeNodes, _Face* prevPolyg) + { + MESSAGE("Try to close polygon. Number of used edge nodes: " << nbUsedEdgeNodes); - const bool hasEdgeIntersections = !_eIntPoints.empty(); - const bool toCheckSideDivision = isImplementEdges() || intFlag & IS_CUT_BY_INTERNAL_FACE; + if ( nFirst == n2 ) + return; - for ( int iF = 0; iF < 6; ++iF ) // loop on 6 sides of a hexahedron + if ( !findChain( n2, nFirst, quad, chainNodes )) { - _Face& quad = _hexQuads[ iF ] ; + if ( !closePolygon( polygon, chainNodes )) + if ( !isImplementEdges() ) + chainNodes.push_back( nFirst ); + } + for ( size_t i = 1; i < chainNodes.size(); ++i ) + { + polygon->AddPolyLink( chainNodes[i-1], chainNodes[i], prevPolyg ); + nbUsedEdgeNodes += bool( chainNodes[i]->IsUsedInFace( polygon )); - _polygons.resize( _polygons.size() + 1 ); - _Face* polygon = &_polygons.back(); - polygon->_polyLinks.reserve( 20 ); - polygon->_name = quad._name; + MESSAGE("Added link for nodes:\n" << *chainNodes[i-1] << *chainNodes[i]); + } - splits.clear(); - for ( int iE = 0; iE < 4; ++iE ) // loop on 4 sides of a quadrangle - for ( size_t iS = 0; iS < quad._links[ iE ].NbResultLinks(); ++iS ) - splits.push_back( quad._links[ iE ].ResultLink( iS )); + MESSAGE("Polygon was closed. Number of used edge nodes: " << nbUsedEdgeNodes); + } - if ( splits.size() == 4 && - isQuadOnFace( iF )) // check if a quad on FACE is not split - { - polygon->_links.swap( splits ); - continue; // goto the next quad - } +void Hexahedron::connectPolygonLinks( + const Solid* solid, _Face* polygon, _Face& quad, vector<_Node*>& chainNodes, std::vector<_OrientedLink>& splits, const bool toCheckSideDivision) +{ + MESSAGE("Connect polygon links..."); - // add splits of links to a polygon and add _polyLinks to make - // polygon's boundary closed + const std::set concaveFaces = getConcaveFaces(solid); // to avoid connecting nodes laying on them - int nbSplits = splits.size(); - if (( nbSplits == 1 ) && - ( quad._eIntNodes.empty() || - splits[0].FirstNode()->IsLinked( splits[0].LastNode()->_intPoint ))) - //( quad._eIntNodes.empty() || _nbCornerNodes + nbIntersections > 6 )) - nbSplits = 0; + // add splits of links to a polygon and add _polyLinks to make + // polygon's boundary closed + int nbSplits = splits.size(); + if (( nbSplits == 1 ) && + ( quad._eIntNodes.empty() || + splits[0].FirstNode()->IsLinked( splits[0].LastNode()->_intPoint ))) + //( quad._eIntNodes.empty() || _nbCornerNodes + nbIntersections > 6 )) + nbSplits = 0; - for ( size_t iP = 0; iP < quad._eIntNodes.size(); ++iP ) - if ( quad._eIntNodes[ iP ]->IsUsedInFace( polygon )) - quad._eIntNodes[ iP ]->_usedInFace = 0; + for ( size_t iP = 0; iP < quad._eIntNodes.size(); ++iP ) + if ( quad._eIntNodes[ iP ]->IsUsedInFace( polygon )) + quad._eIntNodes[ iP ]->_usedInFace = 0; - size_t nbUsedEdgeNodes = 0; - _Face* prevPolyg = 0; // polygon previously created from this quad + size_t nbUsedEdgeNodes = 0; + _Face* prevPolyg = 0; // polygon previously created from this quad - while ( nbSplits > 0 ) - { - size_t iS = 0; - while ( !splits[ iS ] ) - ++iS; + MESSAGE("Number of splits: " << nbSplits); + while ( nbSplits > 0 ) + { + size_t iS = 0; + while ( !splits[ iS ] ) + ++iS; + + if ( !polygon->_links.empty() ) + { + polygon = createPolygon(quad._name); + } + polygon->_links.push_back( splits[ iS ] ); + splits[ iS++ ]._link = 0; + --nbSplits; - if ( !polygon->_links.empty() ) + _Node* nFirst = polygon->_links.back().FirstNode(); + _Node *n1,*n2 = polygon->_links.back().LastNode(); + for ( ; nFirst != n2 && iS < splits.size(); ++iS ) + { + _OrientedLink& split = splits[ iS ]; + MESSAGE("Current split: " << split); + if ( !split ) continue; + + n1 = split.FirstNode(); + MESSAGE("\n\tnFirst: " << *nFirst << "\tn1: " << *n1 << "\tn2: " << *n2); + if ( n1 == n2 && + n1->_intPoint && + (( n1->_intPoint->_faceIDs.size() > 1 && toCheckSideDivision ) || + ( n1->_isInternalFlags ))) + { + // n1 is at intersection with EDGE + MESSAGE("n1 is at intersection with EDGE"); + if ( findChainOnEdge( splits, polygon->_links.back(), split, concaveFaces, + iS, quad, chainNodes )) { - _polygons.resize( _polygons.size() + 1 ); - polygon = &_polygons.back(); - polygon->_polyLinks.reserve( 20 ); - polygon->_name = quad._name; + for ( size_t i = 1; i < chainNodes.size(); ++i ) + polygon->AddPolyLink( chainNodes[i-1], chainNodes[i], prevPolyg ); + if ( chainNodes.back() != n1 ) // not a partial cut by INTERNAL FACE + { + prevPolyg = polygon; + n2 = chainNodes.back(); + continue; + } } - polygon->_links.push_back( splits[ iS ] ); - splits[ iS++ ]._link = 0; - --nbSplits; - - _Node* nFirst = polygon->_links.back().FirstNode(); - _Node *n1,*n2 = polygon->_links.back().LastNode(); - for ( ; nFirst != n2 && iS < splits.size(); ++iS ) + } + else if ( n1 != n2 ) + { + // try to connect to intersections with EDGEs + MESSAGE("n1 != n2 -> try to connect to intersections with EDGEs"); + if ( quad._eIntNodes.size() > nbUsedEdgeNodes && + findChain( n2, n1, quad, chainNodes )) { - _OrientedLink& split = splits[ iS ]; - if ( !split ) continue; - - n1 = split.FirstNode(); - if ( n1 == n2 && - n1->_intPoint && - (( n1->_intPoint->_faceIDs.size() > 1 && toCheckSideDivision ) || - ( n1->_isInternalFlags ))) + for ( size_t i = 1; i < chainNodes.size(); ++i ) { - // n1 is at intersection with EDGE - if ( findChainOnEdge( splits, polygon->_links.back(), split, concaveFaces, - iS, quad, chainNodes )) - { - for ( size_t i = 1; i < chainNodes.size(); ++i ) - polygon->AddPolyLink( chainNodes[i-1], chainNodes[i], prevPolyg ); - if ( chainNodes.back() != n1 ) // not a partial cut by INTERNAL FACE - { - prevPolyg = polygon; - n2 = chainNodes.back(); - continue; - } - } + polygon->AddPolyLink( chainNodes[i-1], chainNodes[i] ); + nbUsedEdgeNodes += ( chainNodes[i]->IsUsedInFace( polygon )); } - else if ( n1 != n2 ) + if ( chainNodes.back() != n1 ) { - // try to connect to intersections with EDGEs - if ( quad._eIntNodes.size() > nbUsedEdgeNodes && - findChain( n2, n1, quad, chainNodes )) + n2 = chainNodes.back(); + --iS; + continue; + } + } + // try to connect to a split ending on the same FACE + else + { + MESSAGE("try to connect to a split ending on the same FACE"); + _OrientedLink foundSplit; + for ( size_t i = iS; i < splits.size() && !foundSplit; ++i ) + if (( foundSplit = splits[ i ]) && + ( n2->IsLinked( foundSplit.FirstNode()->_intPoint ))) { - for ( size_t i = 1; i < chainNodes.size(); ++i ) - { - polygon->AddPolyLink( chainNodes[i-1], chainNodes[i] ); - nbUsedEdgeNodes += ( chainNodes[i]->IsUsedInFace( polygon )); - } - if ( chainNodes.back() != n1 ) - { - n2 = chainNodes.back(); - --iS; - continue; - } + iS = i - 1; } - // try to connect to a split ending on the same FACE else { - _OrientedLink foundSplit; - for ( size_t i = iS; i < splits.size() && !foundSplit; ++i ) - if (( foundSplit = splits[ i ]) && - ( n2->IsLinked( foundSplit.FirstNode()->_intPoint ))) - { - iS = i - 1; - } - else - { - foundSplit._link = 0; - } - if ( foundSplit ) - { - if ( n2 != foundSplit.FirstNode() ) - { - polygon->AddPolyLink( n2, foundSplit.FirstNode() ); - n2 = foundSplit.FirstNode(); - } - continue; - } - else - { - if ( n2->IsLinked( nFirst->_intPoint )) - break; - polygon->AddPolyLink( n2, n1, prevPolyg ); - } + foundSplit._link = 0; } - } // if ( n1 != n2 ) - - polygon->_links.push_back( split ); - split._link = 0; - --nbSplits; - n2 = polygon->_links.back().LastNode(); - - } // loop on splits - if ( nFirst != n2 ) // close a polygon - { - if ( !findChain( n2, nFirst, quad, chainNodes )) + MESSAGE("Found split: " << foundSplit); + if ( foundSplit ) { - if ( !closePolygon( polygon, chainNodes )) - if ( !isImplementEdges() ) - chainNodes.push_back( nFirst ); + if ( n2 != foundSplit.FirstNode() ) + { + polygon->AddPolyLink( n2, foundSplit.FirstNode() ); + n2 = foundSplit.FirstNode(); + } + continue; } - for ( size_t i = 1; i < chainNodes.size(); ++i ) + else { - polygon->AddPolyLink( chainNodes[i-1], chainNodes[i], prevPolyg ); - nbUsedEdgeNodes += bool( chainNodes[i]->IsUsedInFace( polygon )); + if ( n2->IsLinked( nFirst->_intPoint )) + break; + polygon->AddPolyLink( n2, n1, prevPolyg ); } } + } + else + { + MESSAGE("n1 == n2, split will be added as is"); + }// if ( n1 != n2 ) - if ( polygon->_links.size() < 3 && nbSplits > 0 ) - { - polygon->_polyLinks.clear(); - polygon->_links.clear(); - } - } // while ( nbSplits > 0 ) + polygon->_links.push_back( split ); + split._link = 0; + --nbSplits; + n2 = polygon->_links.back().LastNode(); + + } // loop on splits + closePolygon(polygon, n2, nFirst, quad, chainNodes, nbUsedEdgeNodes, prevPolyg); + + if ( polygon->_links.size() < 3 && nbSplits > 0 ) + { + polygon->_polyLinks.clear(); + polygon->_links.clear(); + } + } // while ( nbSplits > 0 ) + + MESSAGE("Polygon's links were connected"); +} + + //================================================================================ + /*! + * \brief Collects chain nodes + */ + std::vector Hexahedron::getChainNodes(const Solid* solid, const IsInternalFlag intFlag) + { + std::vector< _OrientedLink > splits; + std::vector<_Node*> chainNodes; + + const bool toCheckSideDivision = isImplementEdges() || intFlag & IS_CUT_BY_INTERNAL_FACE; + + for ( int iF = 0; iF < 6; ++iF ) // loop on 6 sides of a hexahedron + { + _Face& quad = _hexQuads[ iF ] ; + _Face* polygon = createPolygon(quad._name); + + if (!collectSplits(splits, quad, polygon, iF)) + continue; // goto the next quad + + connectPolygonLinks(solid, polygon, quad, chainNodes, splits, toCheckSideDivision); if ( polygon->_links.size() < 3 ) { _polygons.pop_back(); } } // loop on 6 hexahedron sides - // Create polygons closing holes in a polyhedron - // ---------------------------------------------- + return chainNodes; + } - // clear _usedInFace + //================================================================================ + /*! + * \brief Reset used in faces for hex nodes + */ + void Hexahedron::clearHexUsedInFace() + { + for ( int iN = 0; iN < 8; ++iN ) + _hexNodes[iN]._usedInFace = 0; + } + + //================================================================================ + /*! + * \brief Reset used in faces for int nodes + */ + void Hexahedron::clearIntUsedInFace() + { for ( size_t iN = 0; iN < _intNodes.size(); ++iN ) _intNodes[ iN ]._usedInFace = 0; + } - // add polygons to their links and mark used nodes - for ( size_t iP = 0; iP < _polygons.size(); ++iP ) + //================================================================================ + /*! + * \brief Add polygons to their links and mark used nodes + */ + void Hexahedron::addPolygonsToLinks() + { + for (_Face& polygon : _polygons) { - _Face& polygon = _polygons[ iP ]; - for ( size_t iL = 0; iL < polygon._links.size(); ++iL ) + for (_OrientedLink& link : polygon._links) { - polygon._links[ iL ].AddFace( &polygon ); - polygon._links[ iL ].FirstNode()->_usedInFace = &polygon; + link.AddFace( &polygon ); + link.FirstNode()->_usedInFace = &polygon; } } - // find free links - vector< _OrientedLink* > freeLinks; + } + + //================================================================================ + /*! + * \brief Finds free links in polygons + */ + std::vector Hexahedron::getFreeLinks() + { + std::vector< _OrientedLink* > freeLinks; freeLinks.reserve(20); for ( size_t iP = 0; iP < _polygons.size(); ++iP ) { @@ -3211,29 +3354,55 @@ namespace if ( polygon._links[ iL ].NbFaces() < 2 ) freeLinks.push_back( & polygon._links[ iL ]); } + + return freeLinks; + } + + //================================================================================ + /*! + * \brief Put not used intersection nodes to _vIntNodes + */ + int Hexahedron::notUsedIntersectionNodesToVInt() + { + int nbVertexNodes = 0; + + // TEMP: the loops below can be commented out for simplified testing + // for ( size_t iN = 0; iN < _vIntNodes.size(); ++iN ) + // nbVertexNodes += ( !_vIntNodes[ iN ]->IsUsedInFace() ); + + // const double tol = 1e-3 * Min( Min( _sideLength[0], _sideLength[1] ), _sideLength[0] ); + // for ( size_t iN = _nbFaceIntNodes; iN < _intNodes.size(); ++iN ) + // { + // if ( _intNodes[ iN ].IsUsedInFace() ) continue; + // if ( dynamic_cast< const F_IntersectPoint* >( _intNodes[ iN ]._intPoint )) continue; + // _Node* equalNode = + // findEqualNode( _vIntNodes, _intNodes[ iN ].EdgeIntPnt(), tol*tol ); + // if ( !equalNode ) + // { + // _vIntNodes.push_back( &_intNodes[ iN ]); + // ++nbVertexNodes; + // } + // } + + return nbVertexNodes; + } + + //================================================================================ + /*! + * \brief Creates polygons from free links + */ + bool Hexahedron::createPolygons(const bool hasEdgeIntersections, const IsInternalFlag intFlag) + { + MESSAGE("Create polygons from free links..."); + + // find free links + vector< _OrientedLink* > freeLinks = getFreeLinks(); int nbFreeLinks = freeLinks.size(); + SCRUTE(nbFreeLinks); if ( nbFreeLinks == 1 ) return false; // put not used intersection nodes to _vIntNodes - int nbVertexNodes = 0; // nb not used vertex nodes - { - for ( size_t iN = 0; iN < _vIntNodes.size(); ++iN ) - nbVertexNodes += ( !_vIntNodes[ iN ]->IsUsedInFace() ); - - const double tol = 1e-3 * Min( Min( _sideLength[0], _sideLength[1] ), _sideLength[0] ); - for ( size_t iN = _nbFaceIntNodes; iN < _intNodes.size(); ++iN ) - { - if ( _intNodes[ iN ].IsUsedInFace() ) continue; - if ( dynamic_cast< const F_IntersectPoint* >( _intNodes[ iN ]._intPoint )) continue; - _Node* equalNode = - findEqualNode( _vIntNodes, _intNodes[ iN ].EdgeIntPnt(), tol*tol ); - if ( !equalNode ) - { - _vIntNodes.push_back( &_intNodes[ iN ]); - ++nbVertexNodes; - } - } - } + int nbVertexNodes = notUsedIntersectionNodesToVInt(); std::set usedFaceIDs; std::vector< TGeomID > faces; @@ -3246,6 +3415,7 @@ namespace size_t iPolygon = _polygons.size(); while ( nbFreeLinks > 0 ) { + SCRUTE(nbFreeLinks); if ( iPolygon == _polygons.size() ) { _polygons.resize( _polygons.size() + 1 ); @@ -3259,6 +3429,7 @@ namespace if (( !hasEdgeIntersections ) || ( nbFreeLinks < 4 && nbVertexNodes == 0 )) { + MESSAGE("!hasEdgeIntersections || (nbFreeLinks < 4 && nbVertexNodes == 0), get a remaining link to start from..."); // get a remaining link to start from for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL ) if (( curLink = freeLinks[ iL ] )) @@ -3282,150 +3453,153 @@ namespace } else // there are intersections with EDGEs { + MESSAGE("there are intersections with EDGEs..."); // get a remaining link to start from, one lying on minimal nb of FACEs - { - typedef pair< TGeomID, int > TFaceOfLink; - TFaceOfLink faceOfLink( -1, -1 ); - TFaceOfLink facesOfLink[3] = { faceOfLink, faceOfLink, faceOfLink }; - for ( size_t iL = 0; iL < freeLinks.size(); ++iL ) - if ( freeLinks[ iL ] ) - { - faces = freeLinks[ iL ]->GetNotUsedFace( usedFaceIDs ); - if ( faces.size() == 1 ) - { - faceOfLink = TFaceOfLink( faces[0], iL ); - if ( !freeLinks[ iL ]->HasEdgeNodes() ) - break; - facesOfLink[0] = faceOfLink; - } - else if ( facesOfLink[0].first < 0 ) - { - faceOfLink = TFaceOfLink(( faces.empty() ? -1 : faces[0]), iL ); - facesOfLink[ 1 + faces.empty() ] = faceOfLink; - } - } - for ( int i = 0; faceOfLink.first < 0 && i < 3; ++i ) - faceOfLink = facesOfLink[i]; - - if ( faceOfLink.first < 0 ) // all faces used - { - for ( size_t iL = 0; iL < freeLinks.size() && faceOfLink.first < 1; ++iL ) - if (( curLink = freeLinks[ iL ])) - { - faceOfLink.first = - curLink->FirstNode()->IsLinked( curLink->LastNode()->_intPoint ); - faceOfLink.second = iL; - } - usedFaceIDs.clear(); - } - curFace = faceOfLink.first; - curLink = freeLinks[ faceOfLink.second ]; - freeLinks[ faceOfLink.second ] = 0; - } - usedFaceIDs.insert( curFace ); - polygon._links.push_back( *curLink ); - --nbFreeLinks; + // { + // typedef pair< TGeomID, int > TFaceOfLink; + // TFaceOfLink faceOfLink( -1, -1 ); + // TFaceOfLink facesOfLink[3] = { faceOfLink, faceOfLink, faceOfLink }; + // for ( size_t iL = 0; iL < freeLinks.size(); ++iL ) + // if ( freeLinks[ iL ] ) + // { + // faces = freeLinks[ iL ]->GetNotUsedFace( usedFaceIDs ); + // if ( faces.size() == 1 ) + // { + // faceOfLink = TFaceOfLink( faces[0], iL ); + // if ( !freeLinks[ iL ]->HasEdgeNodes() ) + // break; + // facesOfLink[0] = faceOfLink; + // } + // else if ( facesOfLink[0].first < 0 ) + // { + // faceOfLink = TFaceOfLink(( faces.empty() ? -1 : faces[0]), iL ); + // facesOfLink[ 1 + faces.empty() ] = faceOfLink; + // } + // } + // for ( int i = 0; faceOfLink.first < 0 && i < 3; ++i ) + // faceOfLink = facesOfLink[i]; + + // if ( faceOfLink.first < 0 ) // all faces used + // { + // for ( size_t iL = 0; iL < freeLinks.size() && faceOfLink.first < 1; ++iL ) + // if (( curLink = freeLinks[ iL ])) + // { + // faceOfLink.first = + // curLink->FirstNode()->IsLinked( curLink->LastNode()->_intPoint ); + // faceOfLink.second = iL; + // } + // usedFaceIDs.clear(); + // } + // curFace = faceOfLink.first; + // curLink = freeLinks[ faceOfLink.second ]; + // freeLinks[ faceOfLink.second ] = 0; + // } + // usedFaceIDs.insert( curFace ); + // polygon._links.push_back( *curLink ); + // --nbFreeLinks; // find all links lying on a curFace - do - { - // go forward from curLink - curNode = curLink->LastNode(); - curLink = 0; - for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL ) - if ( freeLinks[ iL ] && - freeLinks[ iL ]->FirstNode() == curNode && - freeLinks[ iL ]->LastNode()->IsOnFace( curFace )) - { - curLink = freeLinks[ iL ]; - freeLinks[ iL ] = 0; - polygon._links.push_back( *curLink ); - --nbFreeLinks; - } - } while ( curLink ); - - std::reverse( polygon._links.begin(), polygon._links.end() ); - - curLink = & polygon._links.back(); - do - { - // go backward from curLink - curNode = curLink->FirstNode(); - curLink = 0; - for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL ) - if ( freeLinks[ iL ] && - freeLinks[ iL ]->LastNode() == curNode && - freeLinks[ iL ]->FirstNode()->IsOnFace( curFace )) - { - curLink = freeLinks[ iL ]; - freeLinks[ iL ] = 0; - polygon._links.push_back( *curLink ); - --nbFreeLinks; - } - } while ( curLink ); - - curNode = polygon._links.back().FirstNode(); - - if ( polygon._links[0].LastNode() != curNode ) - { - if ( nbVertexNodes > 0 ) - { - // add links with _vIntNodes if not already used - chainNodes.clear(); - for ( size_t iN = 0; iN < _vIntNodes.size(); ++iN ) - if ( !_vIntNodes[ iN ]->IsUsedInFace() && - _vIntNodes[ iN ]->IsOnFace( curFace )) - { - _vIntNodes[ iN ]->_usedInFace = &polygon; - chainNodes.push_back( _vIntNodes[ iN ] ); - } - if ( chainNodes.size() > 1 && - curFace != _grid->PseudoIntExtFaceID() ) /////// TODO - { - sortVertexNodes( chainNodes, curNode, curFace ); - } - for ( size_t i = 0; i < chainNodes.size(); ++i ) - { - polygon.AddPolyLink( chainNodes[ i ], curNode ); - curNode = chainNodes[ i ]; - freeLinks.push_back( &polygon._links.back() ); - ++nbFreeLinks; - } - nbVertexNodes -= chainNodes.size(); - } - // if ( polygon._links.size() > 1 ) - { - polygon.AddPolyLink( polygon._links[0].LastNode(), curNode ); - freeLinks.push_back( &polygon._links.back() ); - ++nbFreeLinks; - } - } + // do + // { + // // go forward from curLink + // curNode = curLink->LastNode(); + // curLink = 0; + // for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL ) + // if ( freeLinks[ iL ] && + // freeLinks[ iL ]->FirstNode() == curNode && + // freeLinks[ iL ]->LastNode()->IsOnFace( curFace )) + // { + // curLink = freeLinks[ iL ]; + // freeLinks[ iL ] = 0; + // polygon._links.push_back( *curLink ); + // --nbFreeLinks; + // } + // } while ( curLink ); + + // std::reverse( polygon._links.begin(), polygon._links.end() ); + + // curLink = & polygon._links.back(); + // do + // { + // // go backward from curLink + // curNode = curLink->FirstNode(); + // curLink = 0; + // for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL ) + // if ( freeLinks[ iL ] && + // freeLinks[ iL ]->LastNode() == curNode && + // freeLinks[ iL ]->FirstNode()->IsOnFace( curFace )) + // { + // curLink = freeLinks[ iL ]; + // freeLinks[ iL ] = 0; + // polygon._links.push_back( *curLink ); + // --nbFreeLinks; + // } + // } while ( curLink ); + + // curNode = polygon._links.back().FirstNode(); + + // if ( polygon._links[0].LastNode() != curNode ) + // { + // if ( nbVertexNodes > 0 ) + // { + // // add links with _vIntNodes if not already used + // chainNodes.clear(); + // for ( size_t iN = 0; iN < _vIntNodes.size(); ++iN ) + // if ( !_vIntNodes[ iN ]->IsUsedInFace() && + // _vIntNodes[ iN ]->IsOnFace( curFace )) + // { + // _vIntNodes[ iN ]->_usedInFace = &polygon; + // chainNodes.push_back( _vIntNodes[ iN ] ); + // } + // if ( chainNodes.size() > 1 && + // curFace != _grid->PseudoIntExtFaceID() ) /////// TODO + // { + // sortVertexNodes( chainNodes, curNode, curFace ); + // } + // for ( size_t i = 0; i < chainNodes.size(); ++i ) + // { + // polygon.AddPolyLink( chainNodes[ i ], curNode ); + // curNode = chainNodes[ i ]; + // freeLinks.push_back( &polygon._links.back() ); + // ++nbFreeLinks; + // } + // nbVertexNodes -= chainNodes.size(); + // } + // // if ( polygon._links.size() > 1 ) + // { + // polygon.AddPolyLink( polygon._links[0].LastNode(), curNode ); + // freeLinks.push_back( &polygon._links.back() ); + // ++nbFreeLinks; + // } + // } } // if there are intersections with EDGEs - if ( polygon._links.size() < 2 || - polygon._links[0].LastNode() != polygon._links.back().FirstNode() ) - { - _polygons.clear(); - break; // closed polygon not found -> invalid polyhedron - } + // if ( polygon._links.size() < 2 || + // polygon._links[0].LastNode() != polygon._links.back().FirstNode() ) + // { + // _polygons.clear(); + // break; // closed polygon not found -> invalid polyhedron + // } if ( polygon._links.size() == 2 ) { - if ( freeLinks.back() == &polygon._links.back() ) - { - freeLinks.pop_back(); - --nbFreeLinks; - } - if ( polygon._links.front().NbFaces() > 0 ) - polygon._links.back().AddFace( polygon._links.front()._link->_faces[0] ); - if ( polygon._links.back().NbFaces() > 0 ) - polygon._links.front().AddFace( polygon._links.back()._link->_faces[0] ); + MESSAGE("polygon._links.size() == 2..."); + // if ( freeLinks.back() == &polygon._links.back() ) + // { + // freeLinks.pop_back(); + // --nbFreeLinks; + // } + // if ( polygon._links.front().NbFaces() > 0 ) + // polygon._links.back().AddFace( polygon._links.front()._link->_faces[0] ); + // if ( polygon._links.back().NbFaces() > 0 ) + // polygon._links.front().AddFace( polygon._links.back()._link->_faces[0] ); - if ( iPolygon == _polygons.size()-1 ) - _polygons.pop_back(); + // if ( iPolygon == _polygons.size()-1 ) + // _polygons.pop_back(); } else // polygon._links.size() >= 2 { + MESSAGE("polygon._links.size() >= 2, add polygon to its links..."); // add polygon to its links for ( size_t iL = 0; iL < polygon._links.size(); ++iL ) { @@ -3435,26 +3609,27 @@ namespace if ( /*hasEdgeIntersections &&*/ iPolygon == _polygons.size() - 1 ) { // check that a polygon does not lie on a hexa side - coplanarPolyg = 0; - for ( size_t iL = 0; iL < polygon._links.size() && !coplanarPolyg; ++iL ) - { - if ( polygon._links[ iL ].NbFaces() < 2 ) - continue; // it's a just added free link - // look for a polygon made on a hexa side and sharing - // two or more haxa links - size_t iL2; - coplanarPolyg = polygon._links[ iL ]._link->_faces[0]; - for ( iL2 = iL + 1; iL2 < polygon._links.size(); ++iL2 ) - if ( polygon._links[ iL2 ]._link->_faces[0] == coplanarPolyg && - !coplanarPolyg->IsPolyLink( polygon._links[ iL ]) && - !coplanarPolyg->IsPolyLink( polygon._links[ iL2 ]) && - coplanarPolyg < & _polygons[ nbQuadPolygons ]) - break; - if ( iL2 == polygon._links.size() ) - coplanarPolyg = 0; - } + _Face* coplanarPolyg = nullptr; + // for ( size_t iL = 0; iL < polygon._links.size() && !coplanarPolyg; ++iL ) + // { + // if ( polygon._links[ iL ].NbFaces() < 2 ) + // continue; // it's a just added free link + // // look for a polygon made on a hexa side and sharing + // // two or more haxa links + // size_t iL2; + // coplanarPolyg = polygon._links[ iL ]._link->_faces[0]; + // for ( iL2 = iL + 1; iL2 < polygon._links.size(); ++iL2 ) + // if ( polygon._links[ iL2 ]._link->_faces[0] == coplanarPolyg && + // !coplanarPolyg->IsPolyLink( polygon._links[ iL ]) && + // !coplanarPolyg->IsPolyLink( polygon._links[ iL2 ]) && + // coplanarPolyg < & _polygons[ nbQuadPolygons ]) + // break; + // if ( iL2 == polygon._links.size() ) + // coplanarPolyg = 0; + // } if ( coplanarPolyg ) // coplanar polygon found { + MESSAGE("coplanar polygon found"); freeLinks.resize( freeLinks.size() - polygon._polyLinks.size() ); nbFreeLinks -= polygon._polyLinks.size(); @@ -3531,98 +3706,159 @@ namespace } // end of case ( polygon._links.size() > 2 ) } // while ( nbFreeLinks > 0 ) - for ( auto & face_ip : tmpAddedFace ) - { - curFace = face_ip.first; - for ( const B_IntersectPoint* ip : face_ip.second ) - { - auto it = std::find( ip->_faceIDs.begin(), ip->_faceIDs.end(), curFace ); - if ( it != ip->_faceIDs.end() ) - ip->_faceIDs.erase( it ); - } - } - + // for ( auto & face_ip : tmpAddedFace ) + // { + // curFace = face_ip.first; + // for ( const B_IntersectPoint* ip : face_ip.second ) + // { + // auto it = std::find( ip->_faceIDs.begin(), ip->_faceIDs.end(), curFace ); + // if ( it != ip->_faceIDs.end() ) + // ip->_faceIDs.erase( it ); + // } + // } + if ( _polygons.size() < 3 ) return false; // check volume size double volSize = 0; _hasTooSmall = ! checkPolyhedronSize( intFlag & IS_CUT_BY_INTERNAL_FACE, volSize ); + _volumeDefs._size = volSize; for ( size_t i = 0; i < 8; ++i ) if ( _hexNodes[ i ]._intPoint == &ipTmp ) _hexNodes[ i ]._intPoint = 0; - if ( _hasTooSmall ) - return false; // too small volume + return !_hasTooSmall; // too small volume + } + /*! + * \brief Sets names for no-name polygons + */ + void Hexahedron::setNamesForNoNamePolygons() + { + MESSAGE("Try to find out names of no-name polygons..."); // Try to find out names of no-name polygons (issue # 19887) - if ( _grid->IsToRemoveExcessEntities() && _polygons.back()._name == SMESH_Block::ID_NONE ) - { - gp_XYZ uvwCenter = - 0.5 * ( _grid->_coords[0][_i] + _grid->_coords[0][_i+1] ) * _grid->_axes[0] + - 0.5 * ( _grid->_coords[1][_j] + _grid->_coords[1][_j+1] ) * _grid->_axes[1] + - 0.5 * ( _grid->_coords[2][_k] + _grid->_coords[2][_k+1] ) * _grid->_axes[2]; - for ( size_t i = _polygons.size() - 1; _polygons[i]._name == SMESH_Block::ID_NONE; --i ) - { - _Face& face = _polygons[ i ]; - Bnd_Box bb; - gp_Pnt uvw; - for ( size_t iL = 0; iL < face._links.size(); ++iL ) - { - _Node* n = face._links[ iL ].FirstNode(); - gp_XYZ p = SMESH_NodeXYZ( n->Node() ); - _grid->ComputeUVW( p, uvw.ChangeCoord().ChangeData() ); - bb.Add( uvw ); - } - gp_Pnt pMin = bb.CornerMin(); - if ( bb.IsXThin( _grid->_tol )) - face._name = pMin.X() < uvwCenter.X() ? SMESH_Block::ID_F0yz : SMESH_Block::ID_F1yz; - else if ( bb.IsYThin( _grid->_tol )) - face._name = pMin.Y() < uvwCenter.Y() ? SMESH_Block::ID_Fx0z : SMESH_Block::ID_Fx1z; - else if ( bb.IsZThin( _grid->_tol )) - face._name = pMin.Z() < uvwCenter.Z() ? SMESH_Block::ID_Fxy0 : SMESH_Block::ID_Fxy1; + // TODO: why do we check only last face in vector? + if (!_grid->IsToRemoveExcessEntities() || _polygons.back()._name != SMESH_Block::ID_NONE ) + return; + + MESSAGE("No-name was found! It will be named after grid coordinates."); + + gp_XYZ uvwCenter = + 0.5 * ( _grid->_coords[0][_i] + _grid->_coords[0][_i+1] ) * _grid->_axes[0] + + 0.5 * ( _grid->_coords[1][_j] + _grid->_coords[1][_j+1] ) * _grid->_axes[1] + + 0.5 * ( _grid->_coords[2][_k] + _grid->_coords[2][_k+1] ) * _grid->_axes[2]; + for ( size_t i = _polygons.size() - 1; _polygons[i]._name == SMESH_Block::ID_NONE; --i ) + { + _Face& face = _polygons[ i ]; + Bnd_Box bb; + gp_Pnt uvw; + for ( size_t iL = 0; iL < face._links.size(); ++iL ) + { + _Node* n = face._links[ iL ].FirstNode(); + gp_XYZ p = SMESH_NodeXYZ( n->Node() ); + _grid->ComputeUVW( p, uvw.ChangeCoord().ChangeData() ); + bb.Add( uvw ); } + gp_Pnt pMin = bb.CornerMin(); + if ( bb.IsXThin( _grid->_tol )) + face._name = pMin.X() < uvwCenter.X() ? SMESH_Block::ID_F0yz : SMESH_Block::ID_F1yz; + else if ( bb.IsYThin( _grid->_tol )) + face._name = pMin.Y() < uvwCenter.Y() ? SMESH_Block::ID_Fx0z : SMESH_Block::ID_Fx1z; + else if ( bb.IsZThin( _grid->_tol )) + face._name = pMin.Z() < uvwCenter.Z() ? SMESH_Block::ID_Fxy0 : SMESH_Block::ID_Fxy1; } + } - _volumeDefs._nodes.clear(); - _volumeDefs._quantities.clear(); - _volumeDefs._names.clear(); - // create a classic cell if possible +/*! + * \brief Creates a volume as a classic element (Hexa, Tetra, Penta, Pyra) or non-specific polyhedron + */ +bool Hexahedron::createVolume(const Solid* solid) +{ + MESSAGE("Create a volume..."); + + _volumeDefs._nodes.clear(); + _volumeDefs._quantities.clear(); + _volumeDefs._names.clear(); + // create a classic cell if possible + + int nbPolygons = 0; + for ( size_t iF = 0; iF < _polygons.size(); ++iF ) + nbPolygons += (_polygons[ iF ]._links.size() > 2 ); + + //const int nbNodes = _nbCornerNodes + nbIntersections; + int nbNodes = 0; + for ( size_t i = 0; i < 8; ++i ) + nbNodes += _hexNodes[ i ].IsUsedInFace(); + for ( size_t i = 0; i < _intNodes.size(); ++i ) + nbNodes += _intNodes[ i ].IsUsedInFace(); + + bool isClassicElem = false; + if ( nbNodes == 8 && nbPolygons == 6 ) isClassicElem = addHexa(); + else if ( nbNodes == 4 && nbPolygons == 4 ) isClassicElem = addTetra(); + else if ( nbNodes == 6 && nbPolygons == 5 ) isClassicElem = addPenta(); + else if ( nbNodes == 5 && nbPolygons == 5 ) isClassicElem = addPyra (); + if ( !isClassicElem ) + { + MESSAGE("It's not a classic element. Create a volume by polygons and links..."); - int nbPolygons = 0; for ( size_t iF = 0; iF < _polygons.size(); ++iF ) - nbPolygons += (_polygons[ iF ]._links.size() > 2 ); - - //const int nbNodes = _nbCornerNodes + nbIntersections; - int nbNodes = 0; - for ( size_t i = 0; i < 8; ++i ) - nbNodes += _hexNodes[ i ].IsUsedInFace(); - for ( size_t i = 0; i < _intNodes.size(); ++i ) - nbNodes += _intNodes[ i ].IsUsedInFace(); - - bool isClassicElem = false; - if ( nbNodes == 8 && nbPolygons == 6 ) isClassicElem = addHexa(); - else if ( nbNodes == 4 && nbPolygons == 4 ) isClassicElem = addTetra(); - else if ( nbNodes == 6 && nbPolygons == 5 ) isClassicElem = addPenta(); - else if ( nbNodes == 5 && nbPolygons == 5 ) isClassicElem = addPyra (); - if ( !isClassicElem ) { - for ( size_t iF = 0; iF < _polygons.size(); ++iF ) - { - const size_t nbLinks = _polygons[ iF ]._links.size(); - if ( nbLinks < 3 ) continue; - _volumeDefs._quantities.push_back( nbLinks ); - _volumeDefs._names.push_back( _polygons[ iF ]._name ); - for ( size_t iL = 0; iL < nbLinks; ++iL ) - _volumeDefs._nodes.push_back( _polygons[ iF ]._links[ iL ].FirstNode() ); - } + const size_t nbLinks = _polygons[ iF ]._links.size(); + SCRUTE(nbLinks); + if ( nbLinks < 3 ) continue; + _volumeDefs._quantities.push_back( nbLinks ); + _volumeDefs._names.push_back( _polygons[ iF ]._name ); + SCRUTE(_polygons[ iF ]._name); + for ( size_t iL = 0; iL < nbLinks; ++iL ) + _volumeDefs._nodes.push_back( _polygons[ iF ]._links[ iL ].FirstNode() ); } - _volumeDefs._solidID = solid->ID(); - _volumeDefs._size = volSize; + } + _volumeDefs._solidID = solid->ID(); - return !_volumeDefs._nodes.empty(); + return !_volumeDefs._nodes.empty(); +} + + /*! + * \brief Compute mesh volumes resulted from intersection of the Hexahedron + */ + bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag ) + { + MESSAGE("Compute solid " << solid->ID() << " with internal flag " << intFlag << " ..."); + + // Reset polygons + _polygons.clear(); + _polygons.reserve( 20 ); + clearHexUsedInFace(); + + // Check volumes + if ( intFlag & IS_CUT_BY_INTERNAL_FACE && !_grid->_toAddEdges ) // Issue #19913 + preventVolumesOverlapping(); + + // Create polygons from quadrangles + // -------------------------------- + + vector<_Node*> chainNodes = getChainNodes(solid, intFlag); + + // TODO: check if we can safely move it inside createPolygons() + const bool hasEdgeIntersections = !_eIntPoints.empty(); + + // Create polygons closing holes in a polyhedron + // ---------------------------------------------- + + clearIntUsedInFace(); + addPolygonsToLinks(); + + // Create polygons + if (!createPolygons(hasEdgeIntersections, intFlag)) + return false; + + // Fix polygons' names + setNamesForNoNamePolygons(); + + return createVolume(solid); } template @@ -3651,11 +3887,12 @@ namespace { // A simple case - just one task executed in one thread. // TODO: check if it's faster to do it sequentially - threads.reserve(numTasksTotal); - for (; it < last; ++it) - { - threads.emplace_back(f, std::ref(*it)); - } + // threads.reserve(numTasksTotal); + // for (; it < last; ++it) + // { + // threads.emplace_back(f, std::ref(*it)); + // } + std::for_each(it, last, f); } else { @@ -3818,16 +4055,16 @@ namespace #endif // simplify polyhedrons - if ( _grid->IsToRemoveExcessEntities() ) - { - for ( size_t i = 0; i < intHexa.size(); ++i ) - if ( Hexahedron * hex = intHexa[ i ] ) - hex->removeExcessSideDivision( allHexa ); - - for ( size_t i = 0; i < intHexa.size(); ++i ) - if ( Hexahedron * hex = intHexa[ i ] ) - hex->removeExcessNodes( allHexa ); - } + // if ( _grid->IsToRemoveExcessEntities() ) + // { + // for ( size_t i = 0; i < intHexa.size(); ++i ) + // if ( Hexahedron * hex = intHexa[ i ] ) + // hex->removeExcessSideDivision( allHexa ); + + // for ( size_t i = 0; i < intHexa.size(); ++i ) + // if ( Hexahedron * hex = intHexa[ i ] ) + // hex->removeExcessNodes( allHexa ); + // } // add volumes for ( size_t i = 0; i < intHexa.size(); ++i )