#19913 [CEA] Crash when compute mesh with body fitting and shared faces

fix elements on top of cylinder

src/StdMeshers/StdMeshers_Cartesian_3D.cxx

@@ -731,11 +731,11 @@ namespace
     {
       struct _Split // data of a link split
       {
-        int    _linkID;      // hex link ID
+        int    _linkID;          // hex link ID
         _Node* _nodes[2];
         int    _iCheckIteration; // iteration where split is tried as Hexahedron split
         _Link* _checkedSplit;    // split set to hex links
-        bool   _isUsed;    // used in a volume
+        bool   _isUsed;          // used in a volume
 
         _Split( _Link & split, int iLink ):
           _linkID( iLink ), _nodes{ split._nodes[0], split._nodes[1] },
@@ -944,6 +944,7 @@ namespace
     void getBoundaryElems( vector< const SMDS_MeshElement* > & boundaryVolumes );
     void removeExcessSideDivision(const vector< Hexahedron* >& allHexa);
     void removeExcessNodes(vector< Hexahedron* >& allHexa);
+    void preventVolumesOverlapping();
     TGeomID getAnyFace() const;
     void cutByExtendedInternal( std::vector< Hexahedron* >& hexes,
                                 const TColStd_MapOfInteger& intEdgeIDs );
@@ -2688,6 +2689,9 @@ namespace
     for ( int iN = 0; iN < 8; ++iN )
       _hexNodes[iN]._usedInFace = 0;
 
+    if ( intFlag & IS_CUT_BY_INTERNAL_FACE && !_grid->_toAddEdges ) // Issue #19913
+      preventVolumesOverlapping();
+
     // Create polygons from quadrangles
     // --------------------------------
 
@@ -2695,7 +2699,8 @@ namespace
     vector<_Node*>          chainNodes;
     _Face*                  coplanarPolyg;
 
-    bool hasEdgeIntersections = !_eIntPoints.empty();
+    const bool hasEdgeIntersections = !_eIntPoints.empty();
+    const bool toCheckSideDivision = isImplementEdges() || intFlag & IS_CUT_BY_INTERNAL_FACE;
 
     for ( int iF = 0; iF < 6; ++iF ) // loop on 6 sides of a hexahedron
     {
@@ -2755,7 +2760,7 @@ namespace
           n1 = split.FirstNode();
           if ( n1 == n2 &&
                n1->_intPoint &&
-               (( n1->_intPoint->_faceIDs.size() > 1 && isImplementEdges() ) ||
+               (( n1->_intPoint->_faceIDs.size() > 1 && toCheckSideDivision ) ||
                 ( n1->_isInternalFlags )))
           {
             // n1 is at intersection with EDGE
@@ -5472,6 +5477,69 @@ namespace
 
   } // 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()
+
   //================================================================================
   /*!
    * \brief Set to _hexLinks a next portion of splits located on one side of INTERNAL FACEs