Salome HOME
IPAL22856 2D quadrangle mesher of reduced type works wrong
[modules/smesh.git] / src / StdMeshers / StdMeshers_Cartesian_3D.cxx
index 5956468089c6e80d1c16d72ba28ad738dbd63763..b22843958ce647a57ac52b964290dbd3a426d25b 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2011  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
 //
 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
@@ -85,6 +85,14 @@ using namespace std;
 
 //#define _MY_DEBUG_
 
+#define ELLIPSOLID_WORKAROUND // remove it as soon as http://tracker.dev.opencascade.org/view.php?id=22809 is solved
+
+#ifdef ELLIPSOLID_WORKAROUND
+#include <BRepIntCurveSurface_Inter.hxx>
+#include <BRepTopAdaptor_TopolTool.hxx>
+#include <BRepAdaptor_HSurface.hxx>
+#endif
+
 //=============================================================================
 /*!
  * Constructor
@@ -246,6 +254,41 @@ namespace
                         const TopoDS_Shape&   shape );
     void ComputeNodes(SMESH_MesherHelper& helper);
   };
+#ifdef ELLIPSOLID_WORKAROUND
+  // --------------------------------------------------------------------------
+  /*!
+   * \brief struct temporary replacing IntCurvesFace_Intersector until
+   *        OCCT bug 0022809 is fixed
+   *        http://tracker.dev.opencascade.org/view.php?id=22809
+   */
+  struct TMP_IntCurvesFace_Intersector
+  {
+    BRepAdaptor_Surface                       _surf;
+    double                                    _tol;
+    BRepIntCurveSurface_Inter                 _intcs;
+    vector<IntCurveSurface_IntersectionPoint> _points;
+    BRepTopAdaptor_TopolTool                  _clsf;
+
+    TMP_IntCurvesFace_Intersector(const TopoDS_Face& face, const double tol)
+      :_surf( face ), _tol( tol ), _clsf( new BRepAdaptor_HSurface(_surf) ) {}
+    Bnd_Box Bounding() const { Bnd_Box b; BRepBndLib::Add (_surf.Face(), b); return b; }
+    void Perform( const gp_Lin& line, const double w0, const double w1 )
+    {
+      _points.clear();
+      for ( _intcs.Init( _surf.Face(), line, _tol ); _intcs.More(); _intcs.Next() )
+        if ( w0 <= _intcs.W() && _intcs.W() <= w1 )
+          _points.push_back( _intcs.Point() );
+    }
+    bool IsDone() const { return true; }
+    int  NbPnt()  const { return _points.size(); }
+    IntCurveSurface_TransitionOnCurve Transition( const int i ) const { return _points[ i-1 ].Transition(); }
+    double       WParameter( const int i ) const { return _points[ i-1 ].W(); }
+    TopAbs_State ClassifyUVPoint(const gp_Pnt2d& p) { return _clsf.Classify( p, _tol ); }
+  };
+#define __IntCurvesFace_Intersector TMP_IntCurvesFace_Intersector
+#else
+#define __IntCurvesFace_Intersector IntCurvesFace_Intersector
+#endif
   // --------------------------------------------------------------------------
   /*!
    * \brief Intersector of TopoDS_Face with all GridLine's
@@ -255,7 +298,7 @@ namespace
     TopoDS_Face _face;
     Grid*       _grid;
     Bnd_Box     _bndBox;
-    IntCurvesFace_Intersector* _surfaceInt;
+    __IntCurvesFace_Intersector* _surfaceInt;
     vector< std::pair< GridLine*, IntersectionPoint > > _intersections;
 
     FaceGridIntersector(): _grid(0), _surfaceInt(0) {}
@@ -272,18 +315,18 @@ namespace
       GetCurveFaceIntersector();
       return _bndBox;
     }
-    IntCurvesFace_Intersector* GetCurveFaceIntersector()
+    __IntCurvesFace_Intersector* GetCurveFaceIntersector()
     {
       if ( !_surfaceInt )
       {
-        _surfaceInt = new IntCurvesFace_Intersector( _face, Precision::PConfusion() );
+        _surfaceInt = new __IntCurvesFace_Intersector( _face, Precision::PConfusion() );
         _bndBox     = _surfaceInt->Bounding();
         if ( _bndBox.IsVoid() )
           BRepBndLib::Add (_face, _bndBox);
       }
       return _surfaceInt;
     }
-    bool IsThreadSafe() const;
+    bool IsThreadSafe(set< const Standard_Transient* >& noSafeTShapes) const;
   };
   // --------------------------------------------------------------------------
   /*!
@@ -301,7 +344,7 @@ namespace
     gp_Cone     _cone;
     gp_Sphere   _sphere;
     gp_Torus    _torus;
-    IntCurvesFace_Intersector* _surfaceInt;
+    __IntCurvesFace_Intersector* _surfaceInt;
 
     vector< IntersectionPoint > _intPoints;
 
@@ -312,6 +355,7 @@ namespace
     void IntersectWithTorus   (const GridLine& gridLine);
     void IntersectWithSurface (const GridLine& gridLine);
 
+    bool UVIsOnFace() const;
     void addIntPoint(const bool toClassify=true);
     bool isParamOnLineOK( const double linLength )
     {
@@ -347,8 +391,6 @@ namespace
       vector< _Node>  _intNodes; // _Node's at GridLine intersections
       vector< _Link > _splits;
       vector< _Face*> _faces;
-      const GridLine* _line;
-      _Link(): _line(0) {}
     };
     // --------------------------------------------------------------------------------
     struct _OrientedLink
@@ -881,12 +923,20 @@ namespace
   }
   //================================================================================
   /*
-   * Store an intersection if it is In or ON the face
+   * Return true if (_u,_v) is on the face
+   */
+  bool FaceLineIntersector::UVIsOnFace() const
+  {
+    TopAbs_State state = _surfaceInt->ClassifyUVPoint(gp_Pnt2d( _u,_v ));
+    return ( state == TopAbs_IN || state == TopAbs_ON );
+  }
+  //================================================================================
+  /*
+   * Store an intersection if it is IN or ON the face
    */
   void FaceLineIntersector::addIntPoint(const bool toClassify)
   {
-    TopAbs_State state = toClassify ? _surfaceInt->ClassifyUVPoint(gp_Pnt2d( _u, _v )) : TopAbs_IN;
-    if ( state == TopAbs_IN || state == TopAbs_ON )
+    if ( !toClassify || UVIsOnFace() )
     {
       IntersectionPoint p;
       p._paramOnLine = _w;
@@ -954,8 +1004,7 @@ namespace
       _w = linCone.ParamOnConic( i );
       if ( !isParamOnLineOK( gridLine._length )) continue;
       ElSLib::Parameters(_cone, linCone.Point(i) ,_u,_v);
-      TopAbs_State state = _surfaceInt->ClassifyUVPoint(gp_Pnt2d( _u, _v ));
-      if ( state == TopAbs_IN || state == TopAbs_ON )
+      if ( UVIsOnFace() )
       {
         ElSLib::D1( _u, _v, _cone, P, du, dv );
         norm = du ^ dv;
@@ -1025,8 +1074,7 @@ namespace
       _w = linTorus.ParamOnLine( i );
       if ( !isParamOnLineOK( gridLine._length )) continue;
       linTorus.ParamOnTorus( i, _u,_v );
-      TopAbs_State state = _surfaceInt->ClassifyUVPoint(gp_Pnt2d( _u, _v ));
-      if ( state == TopAbs_IN || state == TopAbs_ON )
+      if ( UVIsOnFace() )
       {
         ElSLib::D1( _u, _v, _torus, P, du, dv );
         norm = du ^ dv;
@@ -1062,8 +1110,10 @@ namespace
   /*
    * check if its face can be safely intersected in a thread
    */
-  bool FaceGridIntersector::IsThreadSafe() const
+  bool FaceGridIntersector::IsThreadSafe(set< const Standard_Transient* >& noSafeTShapes) const
   {
+    bool isSafe = true;
+
     // check surface
     TopLoc_Location loc;
     Handle(Geom_Surface) surf = BRep_Tool::Surface( _face, loc );
@@ -1075,12 +1125,14 @@ namespace
     }
     if ( surf->IsKind( STANDARD_TYPE(Geom_BSplineSurface )) ||
          surf->IsKind( STANDARD_TYPE(Geom_BezierSurface )))
-      return false;
+      if ( !noSafeTShapes.insert((const Standard_Transient*) _face.TShape() ).second )
+        isSafe = false;
 
     double f, l;
     TopExp_Explorer exp( _face, TopAbs_EDGE );
     for ( ; exp.More(); exp.Next() )
     {
+      bool edgeIsSafe = true;
       const TopoDS_Edge& e = TopoDS::Edge( exp.Current() );
       // check 3d curve
       {
@@ -1094,10 +1146,11 @@ namespace
           }
           if ( c->IsKind( STANDARD_TYPE(Geom_BSplineCurve )) ||
                c->IsKind( STANDARD_TYPE(Geom_BezierCurve )))
-            return false;
+            edgeIsSafe = false;
         }
       }
       // check 2d curve
+      if ( edgeIsSafe )
       {
         Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( e, surf, loc, f, l);
         if ( !c2.IsNull() )
@@ -1109,11 +1162,13 @@ namespace
           }
           if ( c2->IsKind( STANDARD_TYPE(Geom2d_BSplineCurve )) ||
                c2->IsKind( STANDARD_TYPE(Geom2d_BezierCurve )))
-            return false;
+            edgeIsSafe = false;
         }
       }
+      if ( !edgeIsSafe && !noSafeTShapes.insert((const Standard_Transient*) e.TShape() ).second )
+        isSafe = false;
     }
-    return true;
+    return isSafe;
   }
   //================================================================================
   /*!
@@ -1223,6 +1278,7 @@ namespace
    */
   void Hexahedron::init( size_t i, size_t j, size_t k )
   {
+    _i = i; _j = j; _k = k;
     // set nodes of grid to nodes of the hexahedron and
     // count nodes at hexahedron corners located IN and ON geometry
     _nbCornerNodes = _nbBndNodes = 0;
@@ -1476,6 +1532,7 @@ namespace
         multiset< IntersectionPoint >::const_iterator ip = line._intPoints.begin();
         for ( ; ip != line._intPoints.end(); ++ip )
         {
+          if ( !ip->_node ) continue;
           lineInd.SetIndexOnLine( ip->_indexOnLine );
           for ( int iL = 0; iL < 4; ++iL ) // loop on 4 cells sharing a link
           {
@@ -1497,12 +1554,8 @@ namespace
               ++nbIntHex;
             }
             const int iLink = iL + iDir * 4;
-            hex->_hexLinks[iLink]._line = &line;
-            if ( ip->_node )
-            {
-              hex->_hexLinks[iLink]._intNodes.push_back( _Node( 0, &(*ip) ));
-              hex->_nbIntNodes++;
-            }
+            hex->_hexLinks[iLink]._intNodes.push_back( _Node( 0, &(*ip) ));
+            hex->_nbIntNodes++;
           }
         }
       }
@@ -1514,7 +1567,7 @@ namespace
     nbIntHex = 0;
     for ( size_t i = 0; i < intersectedHex.size(); ++i )
     {
-      Hexahedron * hex = intersectedHex[ i ];
+      Hexahedron * hex = intersectedHex[ i ];
       if ( hex )
       {
         intHexInd[ nbIntHex++ ] = i;
@@ -1542,6 +1595,14 @@ namespace
           --nbIntHex;
         }
       }
+      else if ( _nbCornerNodes > 3  && !hex )
+      {
+        // all intersection of hex with geometry are at grid nodes
+        hex = new Hexahedron( *this );
+        hex->init( i );
+        intHexInd.push_back(0);
+        intHexInd[ nbIntHex++ ] = i;
+      }
     }
 
     // add elements resulted from hexadron intersection
@@ -1613,23 +1674,29 @@ namespace
     const int ijk[3] = { _i, _j, _k };
     IntersectionPoint curIntPnt;
 
+    // consider a cell to be in a hole if all links in any direction
+    // comes OUT of geometry
     for ( int iDir = 0; iDir < 3; ++iDir )
     {
       const vector<double>& coords = _grid->_coords[ iDir ];
-      bool allLinksOut = true;
-      int linkID = iDir * 4;
-      for ( int i = 0; i < 4 && allLinksOut; ++i )
+      LineIndexer               li = _grid->GetLineIndexer( iDir );
+      li.SetIJK( _i,_j,_k );
+      size_t lineIndex[4] = { li.LineIndex  (),
+                              li.LineIndex10(),
+                              li.LineIndex01(),
+                              li.LineIndex11() };
+      bool allLinksOut = true, hasLinks = false;
+      for ( int iL = 0; iL < 4 && allLinksOut; ++iL ) // loop on 4 links parallel to iDir
       {
-        const _Link& link = _hexLinks[ linkID++ ];
-        if ( !link._line ) return false;
-        if ( link._splits.empty() ) continue;
+        const _Link& link = _hexLinks[ iL + 4*iDir ];
         // check transition of the first node of a link
         const IntersectionPoint* firstIntPnt = 0;
         if ( link._nodes[0]->Node() ) // 1st node is a hexa corner
         {
           curIntPnt._paramOnLine = coords[ ijk[ iDir ]] - coords[0];
+          const GridLine& line = _grid->_lines[ iDir ][ lineIndex[ iL ]];
           multiset< IntersectionPoint >::const_iterator ip =
-            link._line->_intPoints.upper_bound( curIntPnt );
+            line._intPoints.upper_bound( curIntPnt );
           --ip;
           firstIntPnt = &(*ip);
         }
@@ -1638,10 +1705,13 @@ namespace
           firstIntPnt = link._intNodes[0]._intPoint;
         }
 
-        if ( firstIntPnt && firstIntPnt->_transition == Trans_IN )
-          allLinksOut = false;
+        if ( firstIntPnt )
+        {
+          hasLinks = true;
+          allLinksOut = ( firstIntPnt->_transition == Trans_OUT );
+        }
       }
-      if ( allLinksOut )
+      if ( hasLinks && allLinksOut )
         return true;
     }
     return false;
@@ -1882,8 +1952,13 @@ bool StdMeshers_Cartesian_3D::Compute(SMESH_Mesh &         theMesh,
       if ( gridBox.IsOut( gp_Pnt( x0,y0,z0 )) ||
            gridBox.IsOut( gp_Pnt( x1,y1,z1 )))
         for ( size_t i = 0; i < facesItersectors.size(); ++i )
+        {
           if ( !facesItersectors[i].IsInGrid( gridBox ))
             return error("The grid doesn't enclose the geometry");
+#ifdef ELLIPSOLID_WORKAROUND
+          delete facesItersectors[i]._surfaceInt, facesItersectors[i]._surfaceInt = 0;
+#endif
+        }
     }
 
 #ifdef WITH_TBB
@@ -1892,8 +1967,7 @@ bool StdMeshers_Cartesian_3D::Compute(SMESH_Mesh &         theMesh,
       BRepBuilderAPI_Copy copier;
       for ( size_t i = 0; i < facesItersectors.size(); ++i )
       {
-        if ( !facesItersectors[i].IsThreadSafe() &&
-             !tshapes.insert((const Standard_Transient*) facesItersectors[i]._face.TShape() ).second )
+        if ( !facesItersectors[i].IsThreadSafe(tshapes) )
         {
           copier.Perform( facesItersectors[i]._face );
           facesItersectors[i]._face = TopoDS::Face( copier );