[bos #42217][EDF 28921] Horseshoe with bodyfitting.

Nodes linked to a null nodes without any possible splits now being cleared befor splits initialization.

src/StdMeshers/StdMeshers_Cartesian_3D.cxx

@@ -203,6 +203,8 @@ namespace
 
   const TGeomID theUndefID = 1e+9;
 
+  const std::string debugSepLine = "\n===========================================\n";
+
   //=============================================================================
   // Definitions of internal utils
   // --------------------------------------------------------------------------
@@ -774,12 +776,23 @@ namespace
           _node = _intPoint->_node = node;
       }
 
+      void clear()
+      {
+        _node = nullptr;
+        _boundaryCornerNode = nullptr;
+        _intPoint = nullptr;
+        _usedInFace = nullptr;
+        _isInternalFlags = IS_NOT_INTERNAL;
+      }
+
       friend std::ostream& operator<<(std::ostream& os, const _Node& node)
       {
         if (node._node)
-          node._node->Print(os);
+          node._node->Print(os); // mesh node at hexahedron corner
+        else if (node._intPoint && node._intPoint->_node)
+          node._intPoint->_node->Print(os); // intersection point
         else
-          os << "mesh node at hexahedron corner is null" << '\n';
+          os << "mesh node is null" << '\n';
 
         return os;
       }
@@ -787,18 +800,33 @@ namespace
     // --------------------------------------------------------------------------------
     struct _Link // link connecting two _Node's
     {
-      _Node* _nodes[2];
+      static const std::size_t arrSize = 2;
-      _Face* _faces[2]; // polygons sharing a link
+
+      _Node* _nodes[arrSize];
+      _Face* _faces[arrSize]; // polygons sharing a link
       vector< const F_IntersectPoint* > _fIntPoints; // GridLine intersections with FACEs
       vector< _Node* >                  _fIntNodes;   // _Node's at _fIntPoints
       vector< _Link >                   _splits;
       _Link(): _faces{ 0, 0 } {}
 
+      void clear()
+      {
+        for (std::size_t i = 0; i < arrSize; ++i)
+        {
+          if (_nodes[i])
+            _nodes[i]->clear();
+        }
+
+        _fIntPoints.clear();
+        _fIntNodes.clear();
+        _splits.clear();
+      }
+
       friend std::ostream& operator<<(std::ostream& os, const _Link& link)
       {
         os << "Link:\n";
 
-        for (int i = 0; i < 2; ++i)
+        for (std::size_t i = 0; i < arrSize; ++i)
         {
           if (link._nodes[i])
             os << *link._nodes[i];
@@ -806,6 +834,10 @@ namespace
             os << "link node with index " << i << " is null" << '\n';
         }
 
+        os << "_fIntPoints: " << link._fIntPoints.size() << '\n';
+        os << "_fIntNodes: " << link._fIntNodes.size() << '\n';
+        os << "_splits: " << link._splits.size() << '\n';
+
         return os;
       }
     };
@@ -958,6 +990,34 @@ namespace
         _polyLinks.push_back( l );
         _links.push_back( _OrientedLink( &_polyLinks.back() ));
       }
+
+      friend std::ostream& operator<<(std::ostream& os, const _Face& face)
+      {
+        os << "Face " << face._name << '\n';
+
+        os << "Links on GridLines: \n";
+        for (const auto& link : face._links)
+        {
+          os << link;
+        }
+
+        os << "Links added to close a polygonal face: \n";
+        for (const auto& link : face._polyLinks)
+        {
+          os << link;
+        }
+
+        os << "Nodes at intersection with EDGEs: \n";
+        for (const auto node : face._eIntNodes)
+        {
+          if (node)
+          {
+            os << *node;
+          }
+        }
+
+        return os;
+      }
     };
     // --------------------------------------------------------------------------------
     struct _volumeDef // holder of nodes of a volume mesh element
@@ -1037,9 +1097,12 @@ namespace
     };
 
     // topology of a hexahedron
-    _Node _hexNodes [8];
+    static const std::size_t HEX_NODES_NUM = 8;
-    _Link _hexLinks [12];
+    static const std::size_t HEX_LINKS_NUM = 12;
-    _Face _hexQuads [6];
+    static const std::size_t HEX_QUADS_NUM = 6;
+    _Node _hexNodes [HEX_NODES_NUM];
+    _Link _hexLinks [HEX_LINKS_NUM];
+    _Face _hexQuads [HEX_QUADS_NUM];
 
     // faces resulted from hexahedron intersection
     vector< _Face > _polygons;
@@ -1074,6 +1137,8 @@ namespace
     Hexahedron(const Hexahedron& other, size_t i, size_t j, size_t k, int cellID );
     void init( size_t i, size_t j, size_t k, const Solid* solid=0 );
     void init( size_t i );
+    void clearNodesLinkedToNull(const Solid* solid, SMESH_MesherHelper& helper);
+    bool isSplittedLink(const Solid* solid, SMESH_MesherHelper& helper, const Hexahedron::_Link& linkIn) const;
     void setIJK( size_t i );
     bool compute( const Solid* solid, const IsInternalFlag intFlag );
     size_t getSolids( TGeomID ids[] );
@@ -2726,12 +2791,108 @@ namespace
     init( _i, _j, _k );
   }
 
+  //================================================================================
+  /*!
+   * \brief Clear nodes linked to null by not splitted links
+   */
+  void Hexahedron::clearNodesLinkedToNull(const Solid* solid, SMESH_MesherHelper& helper)
+  {
+    for (std::size_t i = 0; i < HEX_LINKS_NUM; ++i)
+    {
+      _Link& link = _hexLinks[i];
+
+      if (isSplittedLink(solid, helper, link))
+        continue;
+
+      // Links where both nodes are valid should have itself as a split.
+      // Check if we have at least one null node here for a chance if we have
+      // some unexpected edge case.
+      _Node* n1 = link._nodes[0];
+      _Node* n2 = link._nodes[1];
+
+      assert(!n1->_node || !n2->_node);
+      if (!n1->_node || !n2->_node)
+      {
+        link.clear();
+      }
+    }
+  }
+
+  //================================================================================
+  /*!
+   * \brief Returns true if the link will be splitted on the Hexaedron initialization
+   */
+  bool Hexahedron::isSplittedLink(const Solid* solid, SMESH_MesherHelper& helper, const Hexahedron::_Link& linkIn) const
+  {
+    _Link split;
+    std::vector<_Node> intNodes;
+
+    _Link link = linkIn;
+    link._fIntNodes.clear();
+    link._fIntNodes.reserve( link._fIntPoints.size() );
+
+    for ( size_t i = 0; i < link._fIntPoints.size(); ++i )
+      if ( solid->ContainsAny( link._fIntPoints[i]->_faceIDs ))
+      {
+        intNodes.push_back( _Node( 0, link._fIntPoints[i] ));
+        link._fIntNodes.push_back( & intNodes.back() );
+      }
+
+    link._splits.clear();
+    split._nodes[ 0 ] = link._nodes[0];
+    bool isOut = ( ! link._nodes[0]->Node() );
+    bool checkTransition;
+    for ( size_t i = 0; i < link._fIntNodes.size(); ++i )
+    {
+      const bool isGridNode = ( ! link._fIntNodes[i]->Node() );
+      if ( !isGridNode ) // intersection non-coincident with a grid node
+      {
+        if ( split._nodes[ 0 ]->Node() && !isOut )
+        {
+          return true;
+        }
+        split._nodes[ 0 ] = link._fIntNodes[i];
+        checkTransition = true;
+      }
+      else // FACE intersection coincident with a grid node (at link ends)
+      {
+        checkTransition = ( i == 0 && link._nodes[0]->Node() );
+      }
+      if ( checkTransition )
+      {
+        const vector< TGeomID >& faceIDs = link._fIntNodes[i]->_intPoint->_faceIDs;
+        if ( _grid->IsInternal( faceIDs.back() ))
+          isOut = false;
+        else if ( faceIDs.size() > 1 || _eIntPoints.size() > 0 )
+          isOut = isOutPoint( link, i, helper, solid );
+        else
+        {
+          bool okTransi = _grid->IsCorrectTransition( faceIDs[0], solid );
+          switch ( link._fIntNodes[i]->FaceIntPnt()->_transition ) {
+          case Trans_OUT: isOut = okTransi;  break;
+          case Trans_IN : isOut = !okTransi; break;
+          default:
+            isOut = isOutPoint( link, i, helper, solid );
+          }
+        }
+      }
+    }
+    if ( link._nodes[ 1 ]->Node() && split._nodes[ 0 ]->Node() && !isOut )
+    {
+      return true;
+    }
+
+    return false;
+  }
+
   //================================================================================
   /*!
    * \brief Initializes its data by given grid cell nodes and intersections
    */
   void Hexahedron::init( size_t i, size_t j, size_t k, const Solid* solid )
   {
+    MESSAGE(debugSepLine << "START Hexahedron::init()" << debugSepLine);
+
     _i = i; _j = j; _k = k;
 
     bool isCompute = solid;
@@ -2780,6 +2941,8 @@ namespace
       // this method can be called in parallel, so use own helper
       SMESH_MesherHelper helper( *_grid->_helper->GetMesh() );
 
+      clearNodesLinkedToNull(solid, helper);
+
       // Create sub-links (_Link::_splits) by splitting links with _Link::_fIntPoints
       // ---------------------------------------------------------------
       _Link split;
@@ -2968,8 +3131,8 @@ namespace
         }
       }
     }
-    return;
 
+    MESSAGE(debugSepLine << "END Hexahedron::init()" << debugSepLine);
   } // init( _i, _j, _k )
 
   //================================================================================
@@ -3084,23 +3247,26 @@ namespace
    */
   bool Hexahedron::collectSplits(std::vector<_OrientedLink>& splits, const _Face& quad, _Face* polygon, int quadIndex)
   {
-    MESSAGE("Collect splits...");
+    MESSAGE(debugSepLine << "START Collect splits..." << debugSepLine);
 
     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 ));
+        // SCRUTE(splits.back());
+      }
 
     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.");
+      MESSAGE(debugSepLine << "END. The quad on FACE is not split. Swap splits with polygon links." << debugSepLine);
 
       polygon->_links.swap( splits );
       return false;
     }
 
-    MESSAGE("Collected splits: " << splits.size());
+    MESSAGE(debugSepLine << "Collected splits: " << splits.size() << debugSepLine);
     return true;
   }
 
@@ -3132,7 +3298,7 @@ namespace
 void Hexahedron::connectPolygonLinks(
   const Solid* solid, _Face* polygon, _Face& quad, vector<_Node*>& chainNodes, std::vector<_OrientedLink>& splits, const bool toCheckSideDivision)
 {
-  MESSAGE("Connect polygon links...");
+  MESSAGE(debugSepLine << "START connectPolygonLinks()..." << debugSepLine);
 
   const std::set<TGeomID> concaveFaces = getConcaveFaces(solid); // to avoid connecting nodes laying on them
 
@@ -3164,6 +3330,7 @@ void Hexahedron::connectPolygonLinks(
       polygon = createPolygon(quad._name);
     }
     polygon->_links.push_back( splits[ iS ] );
+    MESSAGE("Split link was added to a polygon: " << splits[iS]);
     splits[ iS++ ]._link = 0;
     --nbSplits;
 
@@ -3200,10 +3367,11 @@ void Hexahedron::connectPolygonLinks(
       else if ( n1 != n2 )
       {
         // try to connect to intersections with EDGEs
-        MESSAGE("n1 != n2 -> try to connect to intersections with EDGEs");
+        MESSAGE("n1 != n2...");
         if ( quad._eIntNodes.size() > nbUsedEdgeNodes  &&
               findChain( n2, n1, quad, chainNodes ))
         {
+          MESSAGE("try to connect to intersections with EDGEs");
           for ( size_t i = 1; i < chainNodes.size(); ++i )
           {
             polygon->AddPolyLink( chainNodes[i-1], chainNodes[i] );
@@ -3237,6 +3405,7 @@ void Hexahedron::connectPolygonLinks(
           {
             if ( n2 != foundSplit.FirstNode() )
             {
+              MESSAGE("Add link n2 -> foundSplit.FirstNode()");
               polygon->AddPolyLink( n2, foundSplit.FirstNode() );
               n2 = foundSplit.FirstNode();
             }
@@ -3271,7 +3440,7 @@ void Hexahedron::connectPolygonLinks(
     }
   } // while ( nbSplits > 0 )
 
-    MESSAGE("Polygon's links were connected");
+    MESSAGE(debugSepLine << "END connectPolygonLinks(): " << *polygon << debugSepLine);
 }
 
   //================================================================================
@@ -3280,6 +3449,8 @@ void Hexahedron::connectPolygonLinks(
    */
   std::vector<Hexahedron::_Node*> Hexahedron::getChainNodes(const Solid* solid, const IsInternalFlag intFlag)
   {
+    MESSAGE(debugSepLine << "START Get chain nodes..." << debugSepLine);
+
     std::vector< _OrientedLink > splits;
     std::vector<_Node*>          chainNodes;
 
@@ -3289,6 +3460,7 @@ void Hexahedron::connectPolygonLinks(
     {
       _Face& quad = _hexQuads[ iF ] ;
       _Face* polygon = createPolygon(quad._name);
+      SCRUTE(quad);
 
       if (!collectSplits(splits, quad, polygon, iF))
         continue; // goto the next quad
@@ -3300,6 +3472,7 @@ void Hexahedron::connectPolygonLinks(
       }
     }  // loop on 6 hexahedron sides
 
+    MESSAGE(debugSepLine << "END Collected " << chainNodes.size() << " chain nodes" << debugSepLine);
     return chainNodes;
   }
 
@@ -3329,14 +3502,17 @@ void Hexahedron::connectPolygonLinks(
    */
   void Hexahedron::addPolygonsToLinks()
   {
+    MESSAGE(debugSepLine << "START addPolygonsToLinks()..." << debugSepLine);
     for (_Face& polygon : _polygons)
     {
+      // MESSAGE("START add polygon to its links : " << polygon << debugSepLine);
       for (_OrientedLink& link : polygon._links)
       {
         link.AddFace( &polygon );
         link.FirstNode()->_usedInFace = &polygon;
       }
     }
+    MESSAGE(debugSepLine << "END addPolygonsToLinks()" << debugSepLine);
   }
 
   //================================================================================
@@ -3393,7 +3569,7 @@ void Hexahedron::connectPolygonLinks(
    */
   bool Hexahedron::createPolygons(const bool hasEdgeIntersections, const IsInternalFlag intFlag)
   {
-    MESSAGE("Create polygons from free links...");
+    MESSAGE(debugSepLine << "START createPolygons()..." << debugSepLine);
 
     // find free links
     vector< _OrientedLink* > freeLinks = getFreeLinks();
@@ -3408,6 +3584,7 @@ void Hexahedron::connectPolygonLinks(
     std::vector< TGeomID > faces;
     TGeomID curFace = 0;
     const size_t nbQuadPolygons = _polygons.size();
+    SCRUTE(nbQuadPolygons);
     E_IntersectPoint ipTmp;
     std::map< TGeomID, std::vector< const B_IntersectPoint* > > tmpAddedFace; // face added to _intPoint
 
@@ -3423,6 +3600,7 @@ void Hexahedron::connectPolygonLinks(
         _polygons[ iPolygon ]._links.reserve( 20 );
       }
       _Face& polygon = _polygons[ iPolygon ];
+      SCRUTE(polygon);
 
       _OrientedLink* curLink = 0;
       _Node*         curNode;
@@ -3450,6 +3628,7 @@ void Hexahedron::connectPolygonLinks(
               polygon._links.push_back( *curLink );
             }
         } while ( curLink );
+        SCRUTE(polygon);
       }
       else // there are intersections with EDGEs
       {
@@ -3605,6 +3784,7 @@ void Hexahedron::connectPolygonLinks(
         {
           polygon._links[ iL ].AddFace( &polygon );
           polygon._links[ iL ].Reverse();
+          SCRUTE(polygon);
         }
         if ( /*hasEdgeIntersections &&*/ iPolygon == _polygons.size() - 1 )
         {
@@ -3729,6 +3909,7 @@ void Hexahedron::connectPolygonLinks(
       if ( _hexNodes[ i ]._intPoint == &ipTmp )
         _hexNodes[ i ]._intPoint = 0;
 
+    MESSAGE(debugSepLine << "END createPolygons()" << debugSepLine);
     return !_hasTooSmall; // too small volume
   }
 
@@ -3777,7 +3958,7 @@ void Hexahedron::connectPolygonLinks(
   */
 bool Hexahedron::createVolume(const Solid* solid)
 {
-  MESSAGE("Create a volume...");
+  MESSAGE(debugSepLine << "START createVolume()" << debugSepLine);
 
   _volumeDefs._nodes.clear();
   _volumeDefs._quantities.clear();
@@ -3785,9 +3966,12 @@ bool Hexahedron::createVolume(const Solid* solid)
   // create a classic cell if possible
 
   int nbPolygons = 0;
+  SCRUTE(_polygons.size());
   for ( size_t iF = 0; iF < _polygons.size(); ++iF )
     nbPolygons += (_polygons[ iF ]._links.size() > 2 );
 
+  SCRUTE(nbPolygons);
+
   //const int nbNodes = _nbCornerNodes + nbIntersections;
   int nbNodes = 0;
   for ( size_t i = 0; i < 8; ++i )
@@ -3795,6 +3979,7 @@ bool Hexahedron::createVolume(const Solid* solid)
   for ( size_t i = 0; i < _intNodes.size(); ++i )
     nbNodes += _intNodes[ i ].IsUsedInFace();
 
+  SCRUTE(nbNodes);
   bool isClassicElem = false;
   if (      nbNodes == 8 && nbPolygons == 6 ) isClassicElem = addHexa();
   else if ( nbNodes == 4 && nbPolygons == 4 ) isClassicElem = addTetra();
@@ -3807,17 +3992,18 @@ bool Hexahedron::createVolume(const Solid* solid)
     for ( size_t iF = 0; iF < _polygons.size(); ++iF )
     {
       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() );
+
+      // SCRUTE(_polygons[ iF ]);
     }
   }
   _volumeDefs._solidID = solid->ID();
 
+  MESSAGE(debugSepLine << "END createVolume()" << debugSepLine);
   return !_volumeDefs._nodes.empty();
 }
 
@@ -3826,7 +4012,7 @@ bool Hexahedron::createVolume(const Solid* solid)
    */
   bool Hexahedron::compute( const Solid* solid, const IsInternalFlag intFlag )
   {
-    MESSAGE("Compute solid " << solid->ID() << " with internal flag " << intFlag << " ...");
+    MESSAGE(debugSepLine << "START Compute solid " << solid->ID() << " ..." << debugSepLine);
 
     // Reset polygons
     _polygons.clear();
@@ -3858,7 +4044,11 @@ bool Hexahedron::createVolume(const Solid* solid)
     // Fix polygons' names
     setNamesForNoNamePolygons();
 
-    return createVolume(solid);
+    const bool isVolumeCreated = createVolume(solid);
+    SCRUTE(isVolumeCreated);
+    MESSAGE(debugSepLine << "END Compute solid " << solid->ID() << debugSepLine);
+
+    return isVolumeCreated;
   }
 
   template<typename Type>