Salome HOME
[SALOME platform 0013410]: SubMesh not taken into account with Netgen 1D-2D et 1D...
[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.cxx
index 02627a0e6a9957d486a78ba8405f8a730dbd53f9..a02192ca30d9afaf3ad38fa836d68aaece48ef23 100644 (file)
 #include <BRepAdaptor_Surface.hxx>
 #include <BRepTools.hxx>
 #include <BRep_Tool.hxx>
+#include <BRepTools_WireExplorer.hxx>
 #include <Geom2d_Curve.hxx>
 #include <Geom_Curve.hxx>
 #include <Geom_Surface.hxx>
+#include <ShapeAnalysis.hxx>
+#include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
+#include <TopTools_ListIteratorOfListOfShape.hxx>
 #include <TopTools_MapOfShape.hxx>
+#include <TopoDS.hxx>
 #include <gp_Pnt2d.hxx>
-#include <ShapeAnalysis.hxx>
+
+#include <Standard_Failure.hxx>
+#include <Standard_ErrorHandler.hxx>
 
 #include <utilities.h>
 
@@ -69,6 +76,7 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
   // also we have to fill myNLinkNodeMap
   myCreateQuadratic = true;
   mySeamShapeIds.clear();
+  myDegenShapeIds.clear();
   TopAbs_ShapeEnum subType( aSh.ShapeType()==TopAbs_FACE ? TopAbs_EDGE : TopAbs_FACE );
   SMDSAbs_ElementType elemType( subType==TopAbs_FACE ? SMDSAbs_Face : SMDSAbs_Edge );
 
@@ -146,6 +154,7 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
 
   myShape = aSh;
   mySeamShapeIds.clear();
+  myDegenShapeIds.clear();
 
   if ( myShape.IsNull() ) {
     myShapeID  = -1;
@@ -161,8 +170,9 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
     BRepAdaptor_Surface surface( face );
     if ( surface.IsUPeriodic() || surface.IsVPeriodic() )
     {
-      // look for a seam edge
-      for ( TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next()) {
+      for (TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next())
+      {
+        // look for a seam edge
         const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
         if ( BRep_Tool::IsClosed( edge, face )) {
           // initialize myPar1, myPar2 and myParIndex
@@ -181,10 +191,20 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
               myPar2 = surface.LastVParameter();
             }
           }
-          // store shapes indices
-          mySeamShapeIds.insert( meshDS->ShapeToIndex( exp.Current() ));
-          for ( TopExp_Explorer v( exp.Current(), TopAbs_VERTEX ); v.More(); v.Next() )
-            mySeamShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
+          // store seam shape indices, negative if shape encounters twice
+          int edgeID = meshDS->ShapeToIndex( edge );
+          mySeamShapeIds.insert( IsSeamShape( edgeID ) ? -edgeID : edgeID );
+          for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() ) {
+            int vertexID = meshDS->ShapeToIndex( v.Current() );
+            mySeamShapeIds.insert( IsSeamShape( vertexID ) ? -vertexID : vertexID );
+          }
+        }
+
+        // look for a degenerated edge
+        if ( BRep_Tool::Degenerated( edge )) {
+          myDegenShapeIds.insert( meshDS->ShapeToIndex( edge ));
+          for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() )
+            myDegenShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
         }
       }
     }
@@ -298,20 +318,64 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
     int edgeID = Pos->GetShapeId();
     TopoDS_Edge E = TopoDS::Edge(meshDS->IndexToShape(edgeID));
     double f, l;
-    TopLoc_Location loc;
     Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(E, F, f, l);
     uv = C2d->Value( epos->GetUParameter() );
     // for a node on a seam edge select one of UVs on 2 pcurves
-    if ( n2 && mySeamShapeIds.find( edgeID ) != mySeamShapeIds.end() )
+    if ( n2 && IsSeamShape( edgeID ) )
       uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
   }
   else if(Pos->GetTypeOfPosition()==SMDS_TOP_VERTEX)
   {
-    int vertexID = n->GetPosition()->GetShapeId();
-    const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID));
-    uv = BRep_Tool::Parameters( V, F );
-    if ( n2 && mySeamShapeIds.find( vertexID ) != mySeamShapeIds.end() )
-      uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
+    if ( int vertexID = n->GetPosition()->GetShapeId() ) {
+      bool ok = true;
+      const TopoDS_Vertex& V = TopoDS::Vertex(GetMeshDS()->IndexToShape(vertexID));
+      try {
+        uv = BRep_Tool::Parameters( V, F );
+      }
+      catch (Standard_Failure& exc) {
+        ok = false;
+      }
+      if ( !ok ) {
+        for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !ok && vert.More(); vert.Next() )
+          ok = ( V == vert.Current() );
+        if ( !ok ) {
+#ifdef _DEBUG_
+          MESSAGE ( "SMESH_MesherHelper::GetNodeUV(); Vertex " << vertexID
+               << " not in face " << GetMeshDS()->ShapeToIndex( F ) );
+#endif
+          // get UV of a vertex closest to the node
+          double dist = 1e100;
+          gp_Pnt pn ( n->X(),n->Y(),n->Z() );
+          for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !ok && vert.More(); vert.Next() ) {
+            TopoDS_Vertex curV = TopoDS::Vertex( vert.Current() );
+            gp_Pnt p = BRep_Tool::Pnt( curV );
+            double curDist = p.SquareDistance( pn );
+            if ( curDist < dist ) {
+              dist = curDist;
+              uv = BRep_Tool::Parameters( curV, F );
+              if ( dist < DBL_MIN ) break;
+            }
+          }
+        }
+        else {
+          TopTools_ListIteratorOfListOfShape it( myMesh->GetAncestors( V ));
+          for ( ; it.More(); it.Next() ) {
+            if ( it.Value().ShapeType() == TopAbs_EDGE ) {
+              const TopoDS_Edge & edge = TopoDS::Edge( it.Value() );
+              double f,l;
+              Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface(edge, F, f, l);
+              if ( !C2d.IsNull() ) {
+                double u = ( V == TopExp::FirstVertex( edge ) ) ?  f : l;
+                uv = C2d->Value( u );
+                break;
+              }
+            }
+          }
+        }
+      }
+      if ( n2 && IsSeamShape( vertexID ) )
+        uv = GetUVOnSeam( uv, GetNodeUV( F, n2, 0 ));
+    }
   }
   return uv.XY();
 }
@@ -819,19 +883,19 @@ SMDS_MeshVolume* SMESH_MesherHelper::AddVolume(const SMDS_MeshNode* n1,
 }
 
 //=======================================================================
-  /*!
  * \brief Load nodes bound to face into a map of node columns
   * \param theParam2ColumnMap - map of node columns to fill
   * \param theFace - the face on which nodes are searched for
   * \param theBaseEdge - the edge nodes of which are columns' bases
   * \param theMesh - the mesh containing nodes
   * \retval bool - false if something is wrong
  
  * The key of the map is a normalized parameter of each
  * base node on theBaseEdge.
  * This method works in supposition that nodes on the face
  * forms a rectangular grid and elements can be quardrangles or triangles
  */
+/*!
+ * \brief Load nodes bound to face into a map of node columns
+ * \param theParam2ColumnMap - map of node columns to fill
+ * \param theFace - the face on which nodes are searched for
+ * \param theBaseEdge - the edge nodes of which are columns' bases
+ * \param theMesh - the mesh containing nodes
+ * \retval bool - false if something is wrong
+ * 
+ * The key of the map is a normalized parameter of each
+ * base node on theBaseEdge.
+ * This method works in supposition that nodes on the face
+ * forms a rectangular grid and elements can be quardrangles or triangles
+ */
 //=======================================================================
 
 bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
@@ -933,14 +997,14 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
 
   // load nodes from theBaseEdge
 
-  set<const SMDS_MeshNode*> loadedNodes;
+  std::set<const SMDS_MeshNode*> loadedNodes;
   const SMDS_MeshNode* nullNode = 0;
 
-  vector<const SMDS_MeshNode*> & nVecf = theParam2ColumnMap[ 0.];
+  std::vector<const SMDS_MeshNode*> & nVecf = theParam2ColumnMap[ 0.];
   nVecf.resize( vsize, nullNode );
   loadedNodes.insert( nVecf[ 0 ] = smVfb->GetNodes()->next() );
 
-  vector<const SMDS_MeshNode*> & nVecl = theParam2ColumnMap[ 1.];
+  std::vector<const SMDS_MeshNode*> & nVecl = theParam2ColumnMap[ 1.];
   nVecl.resize( vsize, nullNode );
   loadedNodes.insert( nVecl[ 0 ] = smVlb->GetNodes()->next() );
 
@@ -959,7 +1023,7 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
       return false;
     }
     double u = ( pos->GetUParameter() - f ) / range;
-    vector<const SMDS_MeshNode*> & nVec = theParam2ColumnMap[ u ];
+    std::vector<const SMDS_MeshNode*> & nVec = theParam2ColumnMap[ u ];
     nVec.resize( vsize, nullNode );
     loadedNodes.insert( nVec[ 0 ] = node );
   }
@@ -969,7 +1033,7 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
 
   // load nodes from e1
 
-  map< double, const SMDS_MeshNode*> sortedNodes; // sort by param on edge
+  std::map< double, const SMDS_MeshNode*> sortedNodes; // sort by param on edge
   nIt = sm1->GetNodes();
   while ( nIt->more() ) {
     node = nIt->next();
@@ -980,10 +1044,10 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
     if ( !pos ) {
       return false;
     }
-    sortedNodes.insert( make_pair( pos->GetUParameter(), node ));
+    sortedNodes.insert( std::make_pair( pos->GetUParameter(), node ));
   }
   loadedNodes.insert( nVecf[ vsize - 1 ] = smVft->GetNodes()->next() );
-  map< double, const SMDS_MeshNode*>::iterator u_n = sortedNodes.begin();
+  std::map< double, const SMDS_MeshNode*>::iterator u_n = sortedNodes.begin();
   int row  = rev1 ? vsize - 1 : 0;
   int dRow = rev1 ? -1 : +1;
   for ( ; u_n != sortedNodes.end(); u_n++ ) {
@@ -1079,3 +1143,49 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2ColumnMap,
 
   return true;
 }
+
+//=======================================================================
+/**
+ * Check mesh without geometry for: if all elements on this shape are quadratic,
+ * quadratic elements will be created.
+ * Used then generated 3D mesh without geometry.
+ */
+//=======================================================================
+
+SMESH_MesherHelper:: MType SMESH_MesherHelper::IsQuadraticMesh()
+{
+  int NbAllEdgsAndFaces=0;
+  int NbQuadFacesAndEdgs=0;
+  int NbFacesAndEdges=0;
+  //All faces and edges
+  NbAllEdgsAndFaces = myMesh->NbEdges() + myMesh->NbFaces();
+  
+  //Quadratic faces and edges
+  NbQuadFacesAndEdgs = myMesh->NbEdges(ORDER_QUADRATIC) + myMesh->NbFaces(ORDER_QUADRATIC);
+
+  //Linear faces and edges
+  NbFacesAndEdges = myMesh->NbEdges(ORDER_LINEAR) + myMesh->NbFaces(ORDER_LINEAR);
+  
+  if (NbAllEdgsAndFaces == NbQuadFacesAndEdgs) {
+    //Quadratic mesh
+    return SMESH_MesherHelper::QUADRATIC;
+  }
+  else if (NbAllEdgsAndFaces == NbFacesAndEdges) {
+    //Linear mesh
+    return SMESH_MesherHelper::LINEAR;
+  }
+  else
+    //Mesh with both type of elements
+    return SMESH_MesherHelper::COMP;
+}
+
+//=======================================================================
+/*!
+ * \brief Return an alternative parameter for a node on seam
+ */
+//=======================================================================
+
+double SMESH_MesherHelper::GetOtherParam(const double param) const
+{
+  return fabs(param-myPar1) < fabs(param-myPar2) ? myPar2 : myPar1;
+}