23142: EDF 11419 SMESH: Details about extrusion methods

Debug of [HYDRO module - Feature #523] river, channel, embankment meshing

+ Wrong meshing progress if an algo calls SMESH_Gen::Compute()

src/SMESH/SMESH_Gen.cxx

@@ -71,7 +71,6 @@ SMESH_Gen::SMESH_Gen()
   SMDS_Mesh::_meshList.clear();
   MESSAGE(SMDS_Mesh::_meshList.size());
   _compute_canceled = false;
-  _sm_current = NULL;
   //vtkDebugLeaks::SetExitError(0);
 }
 
@@ -181,9 +180,9 @@ bool SMESH_Gen::Compute(SMESH_Mesh &          aMesh,
       {
         if (_compute_canceled)
           return false;
-        _sm_current = smToCompute;
+        setCurrentSubMesh( smToCompute );
         smToCompute->ComputeStateEngine( computeEvent );
-        _sm_current = NULL;
+        setCurrentSubMesh( NULL );
       }
 
       // we check all the sub-meshes here and detect if any of them failed to compute
@@ -268,9 +267,9 @@ bool SMESH_Gen::Compute(SMESH_Mesh &          aMesh,
         {
           if (_compute_canceled)
             return false;
-          _sm_current = smToCompute;
+          setCurrentSubMesh( smToCompute );
           smToCompute->ComputeStateEngine( computeEvent );
-          _sm_current = NULL;
+          setCurrentSubMesh( NULL );
           if ( aShapesId )
             aShapesId->insert( smToCompute->GetId() );
         }
@@ -355,9 +354,9 @@ bool SMESH_Gen::Compute(SMESH_Mesh &          aMesh,
 
           if (_compute_canceled)
             return false;
-          _sm_current = sm;
+          setCurrentSubMesh( sm );
           sm->ComputeStateEngine( computeEvent );
-          _sm_current = NULL;
+          setCurrentSubMesh( NULL );
           if ( aShapesId )
             aShapesId->insert( sm->GetId() );
         }
@@ -400,8 +399,9 @@ void SMESH_Gen::PrepareCompute(SMESH_Mesh &          aMesh,
                                const TopoDS_Shape &  aShape)
 {
   _compute_canceled = false;
-  _sm_current = NULL;
+  resetCurrentSubMesh();
 }
+
 //=============================================================================
 /*!
  * Cancel Compute a mesh
@@ -411,10 +411,43 @@ void SMESH_Gen::CancelCompute(SMESH_Mesh &          aMesh,
                               const TopoDS_Shape &  aShape)
 {
   _compute_canceled = true;
-  if(_sm_current)
+  if ( const SMESH_subMesh* sm = GetCurrentSubMesh() )
-    {
+  {
-      _sm_current->ComputeStateEngine( SMESH_subMesh::COMPUTE_CANCELED );
+    const_cast< SMESH_subMesh* >( sm )->ComputeStateEngine( SMESH_subMesh::COMPUTE_CANCELED );
-    }
+  }
+  resetCurrentSubMesh();
+}
+
+//================================================================================
+/*!
+ * \brief Returns a sub-mesh being currently computed
+ */
+//================================================================================
+
+const SMESH_subMesh* SMESH_Gen::GetCurrentSubMesh() const
+{
+  return _sm_current.empty() ? 0 : _sm_current.back();
+}
+
+//================================================================================
+/*!
+ * \brief Sets a sub-mesh being currently computed.
+ *
+ * An algorithm can call Compute() for a sub-shape, hence we keep a stack of sub-meshes
+ */
+//================================================================================
+
+void SMESH_Gen::setCurrentSubMesh(SMESH_subMesh* sm)
+{
+  if ( sm )
+    _sm_current.push_back( sm );
+  else
+    _sm_current.pop_back();
+}
+
+void SMESH_Gen::resetCurrentSubMesh()
+{
+  _sm_current.clear();
 }
 
 //=============================================================================

src/SMESH/SMESH_Gen.hxx

@@ -87,7 +87,7 @@ public:
   void CancelCompute(::SMESH_Mesh &        aMesh,
                      const TopoDS_Shape &  aShape);
 
-  const SMESH_subMesh* GetCurrentSubMesh() const { return _sm_current; }
+  const SMESH_subMesh* GetCurrentSubMesh() const;
 
   /*!
    * \brief evaluates size of prospective mesh on a shape 
@@ -173,8 +173,11 @@ private:
   // default number of segments
   int _nbSegments;
 
-  volatile bool  _compute_canceled;
+  void setCurrentSubMesh(SMESH_subMesh* sm);
-  SMESH_subMesh* _sm_current;
+  void resetCurrentSubMesh();
+
+  volatile bool               _compute_canceled;
+  std::list< SMESH_subMesh* > _sm_current;
 };
 
 #endif

src/SMESH/SMESH_MeshEditor.cxx

@@ -5002,14 +5002,17 @@ void SMESH_MeshEditor::makeWalls (TNodeOfNodeListMap &     mapNewNodes,
     SMDSAbs_ElementType highType = SMDSAbs_Edge; // count most complex elements only
     while ( eIt->more() && nbInitElems < 2 ) {
       const SMDS_MeshElement* e = eIt->next();
-      SMDSAbs_ElementType type = e->GetType();
+      SMDSAbs_ElementType  type = e->GetType();
-      if ( type == SMDSAbs_Volume || type < highType ) continue;
+      if ( type == SMDSAbs_Volume ||
+           type < highType ||
+           !elemSet.count(e))
+        continue;
       if ( type > highType ) {
         nbInitElems = 0;
-        highType = type;
+        highType    = type;
       }
       el = e;
-      nbInitElems += elemSet.count(el);
+      ++nbInitElems;
     }
     if ( nbInitElems == 1 ) {
       bool NotCreateEdge = el && el->IsMediumNode(node);

src/SMESHUtils/SMESH_MAT2d.cxx

@@ -47,6 +47,7 @@
 #include <TopoDS_Wire.hxx>
 
 #ifdef _DEBUG_
+//#define _MYDEBUG_
 #include "SMESH_File.hxx"
 #include "SMESH_Comment.hxx"
 #endif
@@ -76,6 +77,8 @@ namespace
 
     size_t index( const vector< InPoint >& inPoints ) const { return this - &inPoints[0]; }
     bool operator==( const InPoint& other ) const { return _a == other._a && _b == other._b; }
+    bool operator==( const TVDVertex* v ) const { return ( Abs( _a - v->x() ) < 1. &&
+                                                           Abs( _b - v->y() ) < 1. ); }
   };
   // -------------------------------------------------------------------------------------
 
@@ -105,10 +108,12 @@ namespace
   // check  if a TVDEdge begins at my end or ends at my start
   inline bool InSegment::isConnected( const TVDEdge* edge )
   {
-    return ((Abs( edge->vertex0()->x() - _p1->_a ) < 1.&&
+    return (( edge->vertex0() && edge->vertex1() )
-             Abs( edge->vertex0()->y() - _p1->_b ) < 1.  ) ||
+            &&
-            (Abs( edge->vertex1()->x() - _p0->_a ) < 1.&&
+            ((Abs( edge->vertex0()->x() - _p1->_a ) < 1.&&
-             Abs( edge->vertex1()->y() - _p0->_b ) < 1.  ));
+              Abs( edge->vertex0()->y() - _p1->_b ) < 1.  ) ||
+             (Abs( edge->vertex1()->x() - _p0->_a ) < 1.&&
+              Abs( edge->vertex1()->y() - _p0->_b ) < 1.  )));
   }
 
   // check if a MA TVDEdge is outside of a domain
@@ -147,7 +152,7 @@ namespace
   // }
 
   // -------------------------------------------------------------------------------------
-#ifdef _DEBUG_
+#ifdef _MYDEBUG_
   // writes segments into a txt file readable by voronoi_visualizer
   void inSegmentsToFile( vector< InSegment>& inSegments)
   {
@@ -155,6 +160,7 @@ namespace
       return;
     const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAdebug.txt";
     SMESH_File file(fileName, false );
+    file.remove();
     file.openForWriting();
     SMESH_Comment text;
     text << "0\n"; // nb points
@@ -180,7 +186,7 @@ namespace
     if ( !edge->vertex1() )
       cout << ") -> ( INF, INF";
     else
-      cout << ") -> (" << edge->vertex1()->x() << ", " << edge->vertex1()->y();
+      cout << ") -> ( " << edge->vertex1()->x() << ", " << edge->vertex1()->y();
     cout << ")\t cell=" << edge->cell()
          << " iBnd=" << edge->color()
          << " twin=" << edge->twin()
@@ -253,6 +259,7 @@ namespace boost {
 namespace
 {
   const int theNoBrachID = 0; // std::numeric_limits<int>::max();
+  double theScale[2]; // scale used in bndSegsToMesh()
 
   // -------------------------------------------------------------------------------------
   /*!
@@ -349,7 +356,64 @@ namespace
 
   //================================================================================
   /*!
-   * \brief Computes length of of TVDEdge
+   * \brief debug: to visually check found MA edges
+   */
+  //================================================================================
+
+  void bndSegsToMesh( const vector< BndSeg >& bndSegs )
+  {
+#ifdef _MYDEBUG_
+    if ( !getenv("bndSegsToMesh")) return;
+    map< const TVDVertex *, int > v2Node;
+    map< const TVDVertex *, int >::iterator v2n;
+    set< const TVDEdge* > addedEdges;
+
+    const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAedges.py";
+    SMESH_File file(fileName, false );
+    file.remove();
+    file.openForWriting();
+    SMESH_Comment text;
+    text << "import salome, SMESH\n";
+    text << "salome.salome_init()\n";
+    text << "from salome.smesh import smeshBuilder\n";
+    text << "smesh = smeshBuilder.New(salome.myStudy)\n";
+    text << "m=smesh.Mesh()\n";
+    for ( size_t i = 0; i < bndSegs.size(); ++i )
+    {
+      if ( !bndSegs[i]._edge )
+        text << "# " << i << " NULL edge";
+      else if ( !bndSegs[i]._edge->vertex0() ||
+                !bndSegs[i]._edge->vertex1() )
+        text << "# " << i << " INFINITE edge";
+      else if ( addedEdges.insert( bndSegs[i]._edge ).second &&
+                addedEdges.insert( bndSegs[i]._edge->twin() ).second )
+      {
+        v2n = v2Node.insert( make_pair( bndSegs[i]._edge->vertex0(), v2Node.size() + 1 )).first;
+        int n0 = v2n->second;
+        if ( n0 == v2Node.size() )
+          text << "n" << n0 << " = m.AddNode( "
+               << bndSegs[i]._edge->vertex0()->x() / theScale[0] << ", "
+               << bndSegs[i]._edge->vertex0()->y() / theScale[1] << ", 0 )\n";
+
+        v2n = v2Node.insert( make_pair( bndSegs[i]._edge->vertex1(), v2Node.size() + 1 )).first;
+        int n1 = v2n->second;
+        if ( n1 == v2Node.size() )
+          text << "n" << n1 << " = m.AddNode( "
+               << bndSegs[i]._edge->vertex1()->x() / theScale[0] << ", "
+               << bndSegs[i]._edge->vertex1()->y() / theScale[1] << ", 0 )\n";
+
+        text << "e" << i << " = m.AddEdge([ n" << n0 << ", n" << n1 << " ])\n";
+      }
+    }
+    text << "\n";
+    file.write( text.c_str(), text.size() );
+    cout << "Write " << fileName << endl;
+#endif
+  }
+
+  //================================================================================
+  /*!
+   * \brief Computes length of a TVDEdge
    */
   //================================================================================
 
@@ -593,6 +657,10 @@ namespace
         }
       }
     }
+    // debug
+    theScale[0] = scale[0];
+    theScale[1] = scale[1];
+
     return true;
   }
 
@@ -729,9 +797,12 @@ namespace
           continue;
         inPntChecked[ pInd ] = true;
 
-        const TVDEdge* edge =  // a TVDEdge passing through an end of inSeg
+        const TVDEdge* maE = is2nd ? maEdges.front() : maEdges.back();
-          is2nd ? maEdges.front()->prev() : maEdges.back()->next();
+        if ( inPnt == ( is2nd ? maE->vertex0() : maE->vertex1() ))
-        while ( true )
+          continue;
+        const TVDEdge* edge =  // a secondary TVDEdge connecting inPnt and maE
+          is2nd ? maE->prev() : maE->next();
+        while ( inSeg.isConnected( edge ))
         {
           if ( edge->is_primary() ) break; // this should not happen
           const TVDEdge* edge2 = edge->twin(); // we are in a neighbor cell, add MA edges to inPnt
@@ -825,6 +896,8 @@ namespace
     for ( size_t i = 0; i < bndSegs.size(); ++i )
       bndSegs[i].setIndexToEdge( i );
 
+    bndSegsToMesh( bndSegs ); // debug: visually check found MA edges
+
 
     // Find TVDEdge's of Branches and associate them with bndSegs
 
@@ -839,7 +912,7 @@ namespace
     size_t i1st = 0;
     while ( i1st < bndSegs.size() && !bndSegs[i1st].hasOppositeEdge( noEdgeID ))
       ++i1st;
-    bndSegs[i1st].setBranch( branchID, bndSegs ); // set to the i-th and the opposite bndSeg
+    bndSegs[i1st].setBranch( branchID, bndSegs ); // set to the i-th and to the opposite bndSeg
     branchEdges[ branchID ].push_back( bndSegs[i1st]._edge );
 
     for ( size_t i = i1st+1; i < bndSegs.size(); ++i )
@@ -866,7 +939,7 @@ namespace
           endType.insert( make_pair( bndSegs[i]._edge->vertex1(),
                                      SMESH_MAT2d::BE_BRANCH_POINT ));
       }
-      bndSegs[i].setBranch( branchID, bndSegs ); // set to i-th and the opposite bndSeg
+      bndSegs[i].setBranch( branchID, bndSegs ); // set to i-th and to the opposite bndSeg
       if ( bndSegs[i].hasOppositeEdge( noEdgeID ))
         branchEdges[ bndSegs[i].branchID() ].push_back( bndSegs[i]._edge );
     }
@@ -1066,7 +1139,7 @@ namespace
 
       iSeg = iSegEnd;
 
-    } // loop on all bndSegs
+    } // loop on all bndSegs to construct Boundary
 
     // Initialize branches
 

src/SMESH_SWIG/smeshBuilder.py

@@ -3839,7 +3839,7 @@ class Mesh:
 
 
     ## Generates new elements by extrusion of the elements with given ids
-    #  @param IDsOfElements the list of elements ids for extrusion
+    #  @param IDsOfElements the list of ids of elements or nodes for extrusion
     #  @param StepVector vector or DirStruct or 3 vector components, defining
     #         the direction and value of extrusion for one step (the total extrusion
     #         length will be NbOfSteps * ||StepVector||)
@@ -3855,8 +3855,8 @@ class Mesh:
         return self.ExtrusionSweepObjects(n,e,f, StepVector, NbOfSteps, MakeGroups)
 
     ## Generates new elements by extrusion along the normal to a discretized surface or wire
-    #  @param Elements elements to extrude - a list including ids, groups, sub-meshes or a mesh
+    #  @param Elements elements to extrude - a list including ids, groups, sub-meshes or a mesh.
-    #         Only faces can be extruded so far. Sub-mesh should be a sub-mesh on geom faces.
+    #         Only faces can be extruded so far. A sub-mesh should be a sub-mesh on geom faces.
     #  @param StepSize length of one extrusion step (the total extrusion
     #         length will be \a NbOfSteps * \a StepSize ).
     #  @param NbOfSteps number of extrusion steps.
@@ -3892,15 +3892,15 @@ class Mesh:
         return self.editor.ExtrusionByNormal(Elements, StepSize, NbOfSteps,
                                              ByAverageNormal, UseInputElemsOnly, MakeGroups, Dim)
 
-    ## Generates new elements by extrusion of the elements which belong to the object
+    ## Generates new elements by extrusion of the elements or nodes which belong to the object
-    #  @param theObject the object which elements should be processed.
+    #  @param theObject the object whose elements or nodes should be processed.
-    #                   It can be a mesh, a sub mesh or a group.
+    #                   It can be a mesh, a sub-mesh or a group.
     #  @param StepVector vector or DirStruct or 3 vector components, defining
     #         the direction and value of extrusion for one step (the total extrusion
     #         length will be NbOfSteps * ||StepVector||)
     #  @param NbOfSteps the number of steps
     #  @param MakeGroups forces the generation of new groups from existing ones
-    #  @param  IsNodes is True if elements to extrude are nodes
+    #  @param IsNodes is True if elements to extrude are nodes
     #  @return list of created groups (SMESH_GroupBase) if MakeGroups=True, empty list otherwise
     #  @ingroup l2_modif_extrurev
     def ExtrusionSweepObject(self, theObject, StepVector, NbOfSteps, MakeGroups=False, IsNodes=False):
@@ -3909,9 +3909,9 @@ class Mesh:
         else      : e,f, = theObject,theObject
         return self.ExtrusionSweepObjects(n,e,f, StepVector, NbOfSteps, MakeGroups)
 
-    ## Generates new elements by extrusion of the elements which belong to the object
+    ## Generates new elements by extrusion of edges which belong to the object
-    #  @param theObject object which elements should be processed.
+    #  @param theObject object whose 1D elements should be processed.
-    #                   It can be a mesh, a sub mesh or a group.
+    #                   It can be a mesh, a sub-mesh or a group.
     #  @param StepVector vector or DirStruct or 3 vector components, defining
     #         the direction and value of extrusion for one step (the total extrusion
     #         length will be NbOfSteps * ||StepVector||)
@@ -3922,9 +3922,9 @@ class Mesh:
     def ExtrusionSweepObject1D(self, theObject, StepVector, NbOfSteps, MakeGroups=False):
         return self.ExtrusionSweepObjects([],theObject,[], StepVector, NbOfSteps, MakeGroups)
 
-    ## Generates new elements by extrusion of the elements which belong to the object
+    ## Generates new elements by extrusion of faces which belong to the object
-    #  @param theObject object which elements should be processed.
+    #  @param theObject object whose 2D elements should be processed.
-    #                   It can be a mesh, a sub mesh or a group.
+    #                   It can be a mesh, a sub-mesh or a group.
     #  @param StepVector vector or DirStruct or 3 vector components, defining
     #         the direction and value of extrusion for one step (the total extrusion
     #         length will be NbOfSteps * ||StepVector||)
@@ -4000,7 +4000,7 @@ class Mesh:
 
     ## Generates new elements by extrusion of the given elements
     #  The path of extrusion must be a meshed edge.
-    #  @param Base mesh or group, or submesh, or list of ids of elements for extrusion
+    #  @param Base mesh or group, or sub-mesh, or list of ids of elements for extrusion
     #  @param Path - 1D mesh or 1D sub-mesh, along which proceeds the extrusion
     #  @param NodeStart the start node from Path. Defines the direction of extrusion
     #  @param HasAngles allows the shape to be rotated around the path
@@ -4062,7 +4062,7 @@ class Mesh:
 
     ## Generates new elements by extrusion of the elements which belong to the object
     #  The path of extrusion must be a meshed edge.
-    #  @param theObject the object which elements should be processed.
+    #  @param theObject the object whose elements should be processed.
     #                   It can be a mesh, a sub-mesh or a group.
     #  @param PathMesh mesh containing a 1D sub-mesh on the edge, along which the extrusion proceeds
     #  @param PathShape shape(edge) defines the sub-mesh for the path
@@ -4089,10 +4089,10 @@ class Mesh:
         if MakeGroups: return gr,er
         return er
 
-    ## Generates new elements by extrusion of the elements which belong to the object
+    ## Generates new elements by extrusion of mesh segments which belong to the object
     #  The path of extrusion must be a meshed edge.
-    #  @param theObject the object which elements should be processed.
+    #  @param theObject the object whose 1D elements should be processed.
-    #                   It can be a mesh, a sub mesh or a group.
+    #                   It can be a mesh, a sub-mesh or a group.
     #  @param PathMesh mesh containing a 1D sub-mesh on the edge, along which the extrusion proceeds
     #  @param PathShape shape(edge) defines the sub-mesh for the path
     #  @param NodeStart the first or the last node on the edge. Defines the direction of extrusion
@@ -4118,10 +4118,10 @@ class Mesh:
         if MakeGroups: return gr,er
         return er
 
-    ## Generates new elements by extrusion of the elements which belong to the object
+    ## Generates new elements by extrusion of faces which belong to the object
     #  The path of extrusion must be a meshed edge.
-    #  @param theObject the object which elements should be processed.
+    #  @param theObject the object whose 2D elements should be processed.
-    #                   It can be a mesh, a sub mesh or a group.
+    #                   It can be a mesh, a sub-mesh or a group.
     #  @param PathMesh mesh containing a 1D sub-mesh on the edge, along which the extrusion proceeds
     #  @param PathShape shape(edge) defines the sub-mesh for the path
     #  @param NodeStart the first or the last node on the edge. Defines the direction of extrusion

src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx

@@ -1259,8 +1259,10 @@ namespace
           }
           list< const SMDS_MeshNode* >& mergeNodes = theSinuFace._nodesToMerge[ existingNode ];
 
-          TIterator u2NPprev = sameU2NP.front(); u2NPprev--;
+          TIterator u2NPprev = sameU2NP.front();
-          TIterator u2NPnext = sameU2NP.back() ; u2NPnext++;
+          TIterator u2NPnext = sameU2NP.back() ;
+          if ( u2NPprev->first > 0. ) --u2NPprev;
+          if ( u2NPnext->first < 1. ) ++u2NPprev;
 
           set< int >::iterator edgeID = edgeInds.begin();
           for ( ; edgeID != edgeInds.end(); ++edgeID )
@@ -1933,7 +1935,7 @@ bool StdMeshers_QuadFromMedialAxis_1D2D::Compute(SMESH_Mesh&         theMesh,
 
   if ( getSinuousEdges( helper, sinuFace ))
   {
-    _progress = 0.2;
+    _progress = 0.4;
 
     double minSegLen = getMinSegLen( helper, sinuFace._sinuEdges );
     SMESH_MAT2d::MedialAxis ma( F, sinuFace._sinuEdges, minSegLen, /*ignoreCorners=*/true );
@@ -1946,7 +1948,7 @@ bool StdMeshers_QuadFromMedialAxis_1D2D::Compute(SMESH_Mesh&         theMesh,
     if ( ! divideMA( helper, ma, sinuFace, _regular1D, minSegLen, maParams ))
       return error(COMPERR_BAD_SHAPE);
 
-    _progress = 0.4;
+    _progress = 0.8;
     if ( _hyp2D )
       _regular1D->SetRadialDistribution( _hyp2D );
 
@@ -1954,12 +1956,12 @@ bool StdMeshers_QuadFromMedialAxis_1D2D::Compute(SMESH_Mesh&         theMesh,
          !computeShortEdges( helper, sinuFace._shortSide[1], _regular1D, _hyp2D, 1 ))
       return error("Failed to mesh short edges");
 
-    _progress = 0.6;
+    _progress = 0.85;
 
     if ( !computeSinuEdges( helper, minSegLen, ma, maParams, sinuFace, _regular1D ))
       return error("Failed to mesh sinuous edges");
 
-    _progress = 0.8;
+    _progress = 0.9;
 
     bool ok = computeQuads( helper, sinuFace._quad );