+ //================================================================================
+ /*!
+ * \brief Remove edges and nodes dividing a hexa side in the case if an adjacent
+ * volume also sharing the dividing edge is missing due to its small side.
+ * Issue #19887.
+ */
+ //================================================================================
+
+ void Hexahedron::removeExcessSideDivision(const vector< Hexahedron* >& allHexa)
+ {
+ if ( ! _volumeDefs.IsPolyhedron() )
+ return; // not a polyhedron
+
+ // look for a divided side adjacent to a small hexahedron
+
+ int di[6] = { 0, 0, 0, 0,-1, 1 };
+ int dj[6] = { 0, 0,-1, 1, 0, 0 };
+ int dk[6] = {-1, 1, 0, 0, 0, 0 };
+
+ for ( int iF = 0; iF < 6; ++iF ) // loop on 6 sides of a hexahedron
+ {
+ size_t neighborIndex = _grid->CellIndex( _i + di[iF],
+ _j + dj[iF],
+ _k + dk[iF] );
+ if ( neighborIndex >= allHexa.size() ||
+ !allHexa[ neighborIndex ] ||
+ !allHexa[ neighborIndex ]->_hasTooSmall )
+ continue;
+
+ // check if a side is divided into several polygons
+ for ( _volumeDef* volDef = &_volumeDefs; volDef; volDef = volDef->_next )
+ {
+ int nbPolygons = 0, nbNodes = 0;
+ for ( size_t i = 0; i < volDef->_names.size(); ++i )
+ if ( volDef->_names[ i ] == _hexQuads[ iF ]._name )
+ {
+ ++nbPolygons;
+ nbNodes += volDef->_quantities[ i ];
+ }
+ if ( nbPolygons < 2 )
+ continue;
+
+ // construct loops from polygons
+ typedef _volumeDef::_linkDef TLinkDef;
+ std::vector< TLinkDef* > loops;
+ std::vector< TLinkDef > links( nbNodes );
+ for ( size_t i = 0, iN = 0, iLoop = 0; iLoop < volDef->_quantities.size(); ++iLoop )
+ {
+ size_t nbLinks = volDef->_quantities[ iLoop ];
+ if ( volDef->_names[ iLoop ] != _hexQuads[ iF ]._name )
+ {
+ iN += nbLinks;
+ continue;
+ }
+ loops.push_back( & links[i] );
+ for ( size_t n = 0; n < nbLinks-1; ++n, ++i, ++iN )
+ {
+ links[i].init( volDef->_nodes[iN], volDef->_nodes[iN+1], iLoop );
+ links[i].setNext( &links[i+1] );
+ }
+ links[i].init( volDef->_nodes[iN], volDef->_nodes[iN-nbLinks+1], iLoop );
+ links[i].setNext( &links[i-nbLinks+1] );
+ ++i; ++iN;
+ }
+
+ // look for equal links in different loops and join such loops
+ bool loopsJoined = false;
+ std::set< TLinkDef > linkSet;
+ for ( size_t iLoop = 0; iLoop < loops.size(); ++iLoop )
+ {
+ TLinkDef* beg = 0;
+ for ( TLinkDef* l = loops[ iLoop ]; l != beg; l = l->_next ) // walk around the iLoop
+ {
+ std::pair< std::set< TLinkDef >::iterator, bool > it2new = linkSet.insert( *l );
+ if ( !it2new.second ) // equal found, join loops
+ {
+ const TLinkDef* equal = &(*it2new.first);
+ if ( equal->_loopIndex == l->_loopIndex )
+ continue; // error?
+
+ loopsJoined = true;
+
+ for ( size_t i = iLoop - 1; i < loops.size(); --i )
+ if ( loops[ i ] && loops[ i ]->_loopIndex == equal->_loopIndex )
+ loops[ i ] = 0;
+
+ // exclude l and equal and join two loops
+ if ( l->_prev != equal )
+ l->_prev->setNext( equal->_next );
+ if ( equal->_prev != l )
+ equal->_prev->setNext( l->_next );
+
+ if ( volDef->_quantities[ l->_loopIndex ] > 0 )
+ volDef->_quantities[ l->_loopIndex ] *= -1;
+ if ( volDef->_quantities[ equal->_loopIndex ] > 0 )
+ volDef->_quantities[ equal->_loopIndex ] *= -1;
+
+ if ( loops[ iLoop ] == l )
+ loops[ iLoop ] = l->_prev->_next;
+ }
+ beg = loops[ iLoop ];
+ }
+ }
+ // update volDef
+ if ( loopsJoined )
+ {
+ // set unchanged polygons
+ std::vector< int > newQuantities;
+ std::vector< _volumeDef::_nodeDef > newNodes;
+ vector< SMESH_Block::TShapeID > newNames;
+ newQuantities.reserve( volDef->_quantities.size() );
+ newNodes.reserve ( volDef->_nodes.size() );
+ newNames.reserve ( volDef->_names.size() );
+ for ( size_t i = 0, iLoop = 0; iLoop < volDef->_quantities.size(); ++iLoop )
+ {
+ if ( volDef->_quantities[ iLoop ] < 0 )
+ {
+ i -= volDef->_quantities[ iLoop ];
+ continue;
+ }
+ newQuantities.push_back( volDef->_quantities[ iLoop ]);
+ newNodes.insert( newNodes.end(),
+ volDef->_nodes.begin() + i,
+ volDef->_nodes.begin() + i + newQuantities.back() );
+ newNames.push_back( volDef->_names[ iLoop ]);
+ i += volDef->_quantities[ iLoop ];
+ }
+
+ // set joined loops
+ for ( size_t iLoop = 0; iLoop < loops.size(); ++iLoop )
+ {
+ if ( !loops[ iLoop ] )
+ continue;
+ newQuantities.push_back( 0 );
+ TLinkDef* beg = 0;
+ for ( TLinkDef* l = loops[ iLoop ]; l != beg; l = l->_next, ++newQuantities.back() )
+ {
+ newNodes.push_back( l->_node1 );
+ beg = loops[ iLoop ];
+ }
+ newNames.push_back( _hexQuads[ iF ]._name );
+ }
+ volDef->_quantities.swap( newQuantities );
+ volDef->_nodes.swap( newNodes );
+ volDef->_names.swap( newNames );
+ }
+ } // loop on volDef's
+ } // loop on hex sides
+
+ return;
+ } // removeExcessSideDivision()
+
+
+ //================================================================================
+ /*!
+ * \brief Remove nodes splitting Cartesian cell edges in the case if a node
+ * is used in every cells only by two polygons sharing the edge
+ * Issue #19887.
+ */
+ //================================================================================
+
+ void Hexahedron::removeExcessNodes(vector< Hexahedron* >& allHexa)
+ {
+ if ( ! _volumeDefs.IsPolyhedron() )
+ return; // not a polyhedron
+
+ typedef vector< _volumeDef::_nodeDef >::iterator TNodeIt;
+ vector< int > nodesInPoly[ 4 ]; // node index in _volumeDefs._nodes
+ vector< int > volDefInd [ 4 ]; // index of a _volumeDefs
+ Hexahedron* hexa [ 4 ];
+ int i,j,k, cellIndex, iLink = 0, iCellLink;
+ for ( int iDir = 0; iDir < 3; ++iDir )
+ {
+ CellsAroundLink fourCells( _grid, iDir );
+ for ( int iL = 0; iL < 4; ++iL, ++iLink ) // 4 links in a direction
+ {
+ _Link& link = _hexLinks[ iLink ];
+ fourCells.Init( _i, _j, _k, iLink );
+
+ for ( size_t iP = 0; iP < link._fIntPoints.size(); ++iP ) // loop on nodes on the link
+ {
+ bool nodeRemoved = true;
+ _volumeDef::_nodeDef node; node._intPoint = link._fIntPoints[iP];
+
+ for ( size_t i = 0, nb = _volumeDefs.size(); i < nb && nodeRemoved; ++i )
+ if ( _volumeDef* vol = _volumeDefs.at( i ))
+ nodeRemoved =
+ ( std::find( vol->_nodes.begin(), vol->_nodes.end(), node ) == vol->_nodes.end() );
+ if ( nodeRemoved )
+ continue; // node already removed
+
+ // check if a node encounters zero or two times in 4 cells sharing iLink
+ // if so, the node can be removed from the cells
+ bool nodeIsOnEdge = true;
+ int nbPolyhedraWithNode = 0;
+ for ( int iC = 0; iC < 4; ++iC ) // loop on 4 cells sharing a link
+ {
+ nodesInPoly[ iC ].clear();
+ volDefInd [ iC ].clear();
+ hexa [ iC ] = 0;
+ if ( !fourCells.GetCell( iC, i,j,k, cellIndex, iCellLink ))
+ continue;
+ hexa[ iC ] = allHexa[ cellIndex ];
+ if ( !hexa[ iC ])
+ continue;
+ for ( size_t i = 0, nb = hexa[ iC ]->_volumeDefs.size(); i < nb; ++i )
+ if ( _volumeDef* vol = hexa[ iC ]->_volumeDefs.at( i ))
+ {
+ for ( TNodeIt nIt = vol->_nodes.begin(); nIt != vol->_nodes.end(); ++nIt )
+ {
+ nIt = std::find( nIt, vol->_nodes.end(), node );
+ if ( nIt != vol->_nodes.end() )
+ {
+ nodesInPoly[ iC ].push_back( std::distance( vol->_nodes.begin(), nIt ));
+ volDefInd [ iC ].push_back( i );
+ }
+ else
+ break;
+ }
+ nbPolyhedraWithNode += ( !nodesInPoly[ iC ].empty() );
+ }
+ if ( nodesInPoly[ iC ].size() != 0 &&
+ nodesInPoly[ iC ].size() != 2 )
+ {
+ nodeIsOnEdge = false;
+ break;
+ }
+ } // loop on 4 cells
+
+ // remove nodes from polyhedra
+ if ( nbPolyhedraWithNode > 0 && nodeIsOnEdge )
+ {
+ for ( int iC = 0; iC < 4; ++iC ) // loop on 4 cells sharing the link
+ {
+ if ( nodesInPoly[ iC ].empty() )
+ continue;
+ for ( int i = volDefInd[ iC ].size() - 1; i >= 0; --i )
+ {
+ _volumeDef* vol = hexa[ iC ]->_volumeDefs.at( volDefInd[ iC ][ i ]);
+ int nIndex = nodesInPoly[ iC ][ i ];
+ // decrement _quantities
+ for ( size_t iQ = 0; iQ < vol->_quantities.size(); ++iQ )
+ if ( nIndex < vol->_quantities[ iQ ])
+ {
+ vol->_quantities[ iQ ]--;
+ break;
+ }
+ else
+ {
+ nIndex -= vol->_quantities[ iQ ];
+ }
+ vol->_nodes.erase( vol->_nodes.begin() + nodesInPoly[ iC ][ i ]);
+
+ if ( i == 0 &&
+ vol->_nodes.size() == 6 * 4 &&
+ vol->_quantities.size() == 6 ) // polyhedron becomes hexahedron?
+ {
+ bool allQuads = true;
+ for ( size_t iQ = 0; iQ < vol->_quantities.size() && allQuads; ++iQ )
+ allQuads = ( vol->_quantities[ iQ ] == 4 );
+ if ( allQuads )
+ {
+ // set side nodes as this: bottom, top, top, ...
+ int iTop = 0, iBot = 0; // side indices
+ for ( int iS = 0; iS < 6; ++iS )
+ {
+ if ( vol->_names[ iS ] == SMESH_Block::ID_Fxy0 )
+ iBot = iS;
+ if ( vol->_names[ iS ] == SMESH_Block::ID_Fxy1 )
+ iTop = iS;
+ }
+ if ( iBot != 0 )
+ {
+ if ( iTop == 0 )
+ {
+ std::copy( vol->_nodes.begin(),
+ vol->_nodes.begin() + 4,
+ vol->_nodes.begin() + 4 );
+ iTop = 1;
+ }
+ std::copy( vol->_nodes.begin() + 4 * iBot,
+ vol->_nodes.begin() + 4 * ( iBot + 1),
+ vol->_nodes.begin() );
+ }
+ if ( iTop != 1 )
+ std::copy( vol->_nodes.begin() + 4 * iTop,
+ vol->_nodes.begin() + 4 * ( iTop + 1),
+ vol->_nodes.begin() + 4 );
+
+ std::copy( vol->_nodes.begin() + 4,
+ vol->_nodes.begin() + 8,
+ vol->_nodes.begin() + 8 );
+ // set up top facet nodes by comparing their uvw with bottom nodes
+ E_IntersectPoint ip[8];
+ for ( int iN = 0; iN < 8; ++iN )
+ {
+ SMESH_NodeXYZ p = vol->_nodes[ iN ].Node();
+ _grid->ComputeUVW( p, ip[ iN ]._uvw );
+ }
+ const double tol2 = _grid->_tol * _grid->_tol;
+ for ( int iN = 0; iN < 4; ++iN )
+ {
+ gp_Pnt2d pBot( ip[ iN ]._uvw[0], ip[ iN ]._uvw[1] );
+ for ( int iT = 4; iT < 8; ++iT )
+ {
+ gp_Pnt2d pTop( ip[ iT ]._uvw[0], ip[ iT ]._uvw[1] );
+ if ( pBot.SquareDistance( pTop ) < tol2 )
+ {
+ // vol->_nodes[ iN + 4 ]._node = ip[ iT ]._node;
+ // vol->_nodes[ iN + 4 ]._intPoint = 0;
+ vol->_nodes[ iN + 4 ] = vol->_nodes[ iT + 4 ];
+ break;
+ }
+ }
+ }
+ vol->_nodes.resize( 8 );
+ vol->_quantities.clear();
+ //vol->_names.clear();
+ }
+ }
+ } // loop on _volumeDefs
+ } // loop on 4 cell abound a link
+ } // if ( nodeIsOnEdge )
+ } // loop on intersection points of a link
+ } // loop on 4 links of a direction
+ } // loop on 3 directions
+
+ return;
+
+ } // removeExcessNodes()
+
+ //================================================================================
+ /*!
+ * \brief [Issue #19913] Modify _hexLinks._splits to prevent creating overlapping volumes
+ */
+ //================================================================================
+
+ void Hexahedron::preventVolumesOverlapping()
+ {
+ // Cut off a quadrangle corner if two links sharing the corner
+ // are shared by same two solids, in this case each of solids gets
+ // a triangle for it-self.
+ std::vector< TGeomID > soIDs[4];
+ for ( int iF = 0; iF < 6; ++iF ) // loop on 6 sides of a hexahedron
+ {
+ _Face& quad = _hexQuads[ iF ] ;
+
+ int iFOpposite = iF + ( iF % 2 ? -1 : 1 );
+ _Face& quadOpp = _hexQuads[ iFOpposite ] ;
+
+ int nbSides = 0, nbSidesOpp = 0;
+ for ( int iE = 0; iE < 4; ++iE ) // loop on 4 sides of a quadrangle
+ {
+ nbSides += ( quad._links [ iE ].NbResultLinks() > 0 );
+ nbSidesOpp += ( quadOpp._links[ iE ].NbResultLinks() > 0 );
+ }
+ if ( nbSides < 4 || nbSidesOpp != 2 )
+ continue;
+
+ for ( int iE = 0; iE < 4; ++iE )
+ {
+ soIDs[ iE ].clear();
+ _Node* n = quad._links[ iE ].FirstNode();
+ if ( n->_intPoint && n->_intPoint->_faceIDs.size() )
+ soIDs[ iE ] = _grid->GetSolidIDs( n->_intPoint->_faceIDs[0] );
+ }
+ if ((( soIDs[0].size() >= 2 ) +
+ ( soIDs[1].size() >= 2 ) +
+ ( soIDs[2].size() >= 2 ) +
+ ( soIDs[3].size() >= 2 ) ) < 3 )
+ continue;
+
+ bool done = false;
+ for ( int i = 0; i < 4; ++i )
+ {
+ int i1 = _grid->_helper->WrapIndex( i + 1, 4 );
+ int i2 = _grid->_helper->WrapIndex( i + 2, 4 );
+ int i3 = _grid->_helper->WrapIndex( i + 3, 4 );
+ if ( soIDs[i1].size() == 2 && soIDs[i ] != soIDs[i1] &&
+ soIDs[i2].size() == 2 && soIDs[i1] == soIDs[i2] &&
+ soIDs[i3].size() == 2 && soIDs[i2] == soIDs[i3] )
+ {
+ quad._links[ i1 ]._link->_splits.clear();
+ quad._links[ i2 ]._link->_splits.clear();
+ done = true;
+ break;
+ }
+ }
+ if ( done )
+ break;
+ }
+ return;
+ } // preventVolumesOverlapping()
+