Salome HOME
IPAL52479: Mixed linear/quadratic mesh is created
[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.cxx
index 784853234b4157551317d0220806b50bffa94f06..31f2d88b7b98a25c5d8339de5f334be10b009709 100644 (file)
@@ -242,34 +242,39 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
   for ( TopExp_Explorer eF( aSh, TopAbs_FACE ); eF.More(); eF.Next() )
   {
     const TopoDS_Face& face = TopoDS::Face( eF.Current() );
-    TopLoc_Location loc;
-    Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
 
     // if ( surface->IsUPeriodic() || surface->IsVPeriodic() ||
     //      surface->IsUClosed()   || surface->IsVClosed() )
     {
       //while ( surface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface )))
       //surface = Handle(Geom_RectangularTrimmedSurface)::DownCast( surface )->BasisSurface();
-      GeomAdaptor_Surface surf( surface );
 
       for (TopExp_Explorer exp( face, TopAbs_EDGE ); exp.More(); exp.Next())
       {
         // look for a seam edge
-        const TopoDS_Edge& edge = TopoDS::Edge( exp.Current() );
+        TopoDS_Edge edge = TopoDS::Edge( exp.Current() );
         if ( BRep_Tool::IsClosed( edge, face )) {
           // initialize myPar1, myPar2 and myParIndex
           gp_Pnt2d uv1, uv2;
           BRep_Tool::UVPoints( edge, face, uv1, uv2 );
           if ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Abs( uv1.Coord(2) - uv2.Coord(2) ))
           {
+            double u1 = uv1.Coord(1);
+            edge.Reverse();
+            BRep_Tool::UVPoints( edge, face, uv1, uv2 );
+            double u2 = uv1.Coord(1);
             myParIndex |= U_periodic;
-            myPar1[0] = surf.FirstUParameter();
-            myPar2[0] = surf.LastUParameter();
+            myPar1[0] = Min( u1, u2 );
+            myPar2[0] = Max( u1, u2 );
           }
           else {
+            double v1 = uv1.Coord(2);
+            edge.Reverse();
+            BRep_Tool::UVPoints( edge, face, uv1, uv2 );
+            double v2 = uv1.Coord(2);
             myParIndex |= V_periodic;
-            myPar1[1] = surf.FirstVParameter();
-            myPar2[1] = surf.LastVParameter();
+            myPar1[1] = Min( v1, v2 );
+            myPar2[1] = Max( v1, v2 );
           }
           // store seam shape indices, negative if shape encounters twice
           int edgeID = meshDS->ShapeToIndex( edge );
@@ -287,13 +292,15 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
             myDegenShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
         }
       }
-      if ( !myDegenShapeIds.empty() && !myParIndex ) {
-        if ( surface->IsUPeriodic() || surface->IsUClosed() ) {
+      if ( !myDegenShapeIds.empty() && !myParIndex )
+      {
+        BRepAdaptor_Surface surf( face, false );
+        if ( surf.IsUPeriodic() || surf.IsUClosed() ) {
           myParIndex |= U_periodic;
           myPar1[0] = surf.FirstUParameter();
           myPar2[0] = surf.LastUParameter();
         }
-        else if ( surface->IsVPeriodic() || surface->IsVClosed() ) {
+        else if ( surf.IsVPeriodic() || surf.IsVClosed() ) {
           myParIndex |= V_periodic;
           myPar1[1] = surf.FirstVParameter();
           myPar2[1] = surf.LastVParameter();
@@ -373,10 +380,13 @@ void SMESH_MesherHelper::AddTLinkNode(const SMDS_MeshNode* n1,
  */
 //================================================================================
 
-void SMESH_MesherHelper::AddTLinks(const SMDS_MeshEdge* edge)
+bool SMESH_MesherHelper::AddTLinks(const SMDS_MeshEdge* edge)
 {
-  if ( edge->IsQuadratic() )
+  if ( edge && edge->IsQuadratic() )
     AddTLinkNode(edge->GetNode(0), edge->GetNode(1), edge->GetNode(2));
+  else
+    return false;
+  return true;
 }
 
 //================================================================================
@@ -385,8 +395,9 @@ void SMESH_MesherHelper::AddTLinks(const SMDS_MeshEdge* edge)
  */
 //================================================================================
 
-void SMESH_MesherHelper::AddTLinks(const SMDS_MeshFace* f)
+bool SMESH_MesherHelper::AddTLinks(const SMDS_MeshFace* f)
 {
+  bool isQuad = true;
   if ( !f->IsPoly() )
     switch ( f->NbNodes() ) {
     case 7:
@@ -410,7 +421,9 @@ void SMESH_MesherHelper::AddTLinks(const SMDS_MeshFace* f)
       AddTLinkNode(f->GetNode(2),f->GetNode(3),f->GetNode(6));
       AddTLinkNode(f->GetNode(3),f->GetNode(0),f->GetNode(7)); break;
     default:;
+      isQuad = false;
     }
+  return isQuad;
 }
 
 //================================================================================
@@ -419,7 +432,7 @@ void SMESH_MesherHelper::AddTLinks(const SMDS_MeshFace* f)
  */
 //================================================================================
 
-void SMESH_MesherHelper::AddTLinks(const SMDS_MeshVolume* volume)
+bool SMESH_MesherHelper::AddTLinks(const SMDS_MeshVolume* volume)
 {
   if ( volume->IsQuadratic() )
   {
@@ -453,7 +466,9 @@ void SMESH_MesherHelper::AddTLinks(const SMDS_MeshVolume* volume)
                          nFCenter ));
       }
     }
+    return true;
   }
+  return false;
 }
 
 //================================================================================
@@ -513,8 +528,8 @@ gp_Pnt2d SMESH_MesherHelper::GetUVOnSeam( const gp_Pnt2d& uv1, const gp_Pnt2d& u
       double p1 = uv1.Coord( i );
       double dp1 = Abs( p1-myPar1[i-1]), dp2 = Abs( p1-myPar2[i-1]);
       if ( myParIndex == i ||
-           dp1 < ( myPar2[i-1] - myPar2[i-1] ) / 100. ||
-           dp2 < ( myPar2[i-1] - myPar2[i-1] ) / 100. )
+           dp1 < ( myPar2[i-1] - myPar1[i-1] ) / 100. ||
+           dp2 < ( myPar2[i-1] - myPar1[i-1] ) / 100. )
       {
         double p2 = uv2.Coord( i );
         double p1Alt = ( dp1 < dp2 ) ? myPar2[i-1] : myPar1[i-1];
@@ -544,7 +559,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
   {
     // node has position on face
     const SMDS_FacePosition* fpos =
-      static_cast<const SMDS_FacePosition*>(n->GetPosition());
+      static_cast<const SMDS_FacePosition*>( Pos );
     uv.SetCoord(fpos->GetUParameter(),fpos->GetVParameter());
     if ( check )
       uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*MaxTolerance( F ));
@@ -555,7 +570,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
     // corresponding edge from face, get pcurve for this
     // edge and retrieve value from this pcurve
     const SMDS_EdgePosition* epos =
-      static_cast<const SMDS_EdgePosition*>(n->GetPosition());
+      static_cast<const SMDS_EdgePosition*>( Pos );
     int edgeID = n->getshapeId();
     TopoDS_Edge E = TopoDS::Edge(GetMeshDS()->IndexToShape(edgeID));
     double f, l, u = epos->GetUParameter();
@@ -1457,7 +1472,6 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
           if ( myParIndex & U_periodic ) uv[1].SetCoord( 1, uv[0].Coord( 1 ));
           else                           uv[1].SetCoord( 2, uv[0].Coord( 2 ));
         }
-
         TopLoc_Location loc;
         Handle(Geom_Surface) S = BRep_Tool::Surface(F,loc);
         gp_XY UV = GetMiddleUV( S, uv[0], uv[1] );
@@ -2322,11 +2336,12 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2
           while ( ++u_n != sortedBaseNN.end() && !isNodeInSubMesh( u_n->second, faceSubMesh ));
           sortedBaseNN.erase( sortedBaseNN.begin(), u_n );
         }
-        if ( u_n = --sortedBaseNN.end(), !isNodeInSubMesh( u_n->second, faceSubMesh ))
-        {
-          while ( u_n != sortedBaseNN.begin() && !isNodeInSubMesh( (--u_n)->second, faceSubMesh ));
-          sortedBaseNN.erase( ++u_n, sortedBaseNN.end() );
-        }
+        if ( !sortedBaseNN.empty() )
+          if ( u_n = --sortedBaseNN.end(), !isNodeInSubMesh( u_n->second, faceSubMesh ))
+          {
+            while ( u_n != sortedBaseNN.begin() && !isNodeInSubMesh( (--u_n)->second, faceSubMesh ));
+            sortedBaseNN.erase( ++u_n, sortedBaseNN.end() );
+          }
         if ( sortedBaseNN.empty() ) continue;
       }
 
@@ -2343,7 +2358,7 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap &            theParam2
         u2nn->second.push_back( u_n->second );
       }
     }
-    if ( theParam2ColumnMap.empty() )
+    if ( theParam2ColumnMap.size() < 2 )
       return false;
   }
 
@@ -2726,23 +2741,22 @@ double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape )
  */
 //================================================================================
 
-double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1,
-                                     const TopoDS_Edge & theE2,
-                                     const TopoDS_Face & theFace)
+double SMESH_MesherHelper::GetAngle( const TopoDS_Edge &   theE1,
+                                     const TopoDS_Edge &   theE2,
+                                     const TopoDS_Face &   theFace,
+                                     const TopoDS_Vertex & theCommonV,
+                                     gp_Vec*               theFaceNormal)
 {
   double angle = 1e100;
   try
   {
-    TopoDS_Vertex vCommon;
-    if ( !TopExp::CommonVertex( theE1, theE2, vCommon ))
-      return angle;
     double f,l;
     Handle(Geom_Curve)     c1 = BRep_Tool::Curve( theE1, f,l );
     Handle(Geom_Curve)     c2 = BRep_Tool::Curve( theE2, f,l );
     Handle(Geom2d_Curve) c2d1 = BRep_Tool::CurveOnSurface( theE1, theFace, f,l );
     Handle(Geom_Surface) surf = BRep_Tool::Surface( theFace );
-    double                 p1 = BRep_Tool::Parameter( vCommon, theE1 );
-    double                 p2 = BRep_Tool::Parameter( vCommon, theE2 );
+    double                 p1 = BRep_Tool::Parameter( theCommonV, theE1 );
+    double                 p2 = BRep_Tool::Parameter( theCommonV, theE2 );
     if ( c1.IsNull() || c2.IsNull() )
       return angle;
     gp_XY uv = c2d1->Value( p1 ).XY();
@@ -2751,10 +2765,10 @@ double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1,
     gp_Vec vec1, vec2, vecRef = du ^ dv;
     int  nbLoops = 0;
     double p1tmp = p1;
-    while ( vecRef.SquareMagnitude() < std::numeric_limits<double>::min() )
+    while ( vecRef.SquareMagnitude() < 1e-25 )
     {
       double dp = ( l - f ) / 1000.;
-      p1tmp += dp * (( Abs( p1 - f ) > Abs( p1 - l )) ? +1. : -1.);
+      p1tmp += dp * (( Abs( p1 - f ) > Abs( p1 - l )) ? -1. : +1.);
       uv = c2d1->Value( p1tmp ).XY();
       surf->D1( uv.X(), uv.Y(), p, du, dv );
       vecRef = du ^ dv;
@@ -2768,16 +2782,33 @@ double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1,
     }
     if ( theFace.Orientation() == TopAbs_REVERSED )
       vecRef.Reverse();
+    if ( theFaceNormal ) *theFaceNormal = vecRef;
+
     c1->D1( p1, p, vec1 );
     c2->D1( p2, p, vec2 );
-    TopoDS_Face F = theFace;
-    if ( F.Orientation() == TopAbs_INTERNAL )
-      F.Orientation( TopAbs_FORWARD );
+    // TopoDS_Face F = theFace;
+    // if ( F.Orientation() == TopAbs_INTERNAL )
+    //   F.Orientation( TopAbs_FORWARD );
     if ( theE1.Orientation() /*GetSubShapeOri( F, theE1 )*/ == TopAbs_REVERSED )
       vec1.Reverse();
     if ( theE2.Orientation() /*GetSubShapeOri( F, theE2 )*/ == TopAbs_REVERSED )
       vec2.Reverse();
     angle = vec1.AngleWithRef( vec2, vecRef );
+
+    if ( Abs ( angle ) >= 0.99 * M_PI )
+    {
+      BRep_Tool::Range( theE1, f, l );
+      p1 += 1e-7 * ( p1-f < l-p1 ? +1. : -1. );
+      c1->D1( p1, p, vec1 );
+      if ( theE1.Orientation() == TopAbs_REVERSED )
+        vec1.Reverse();
+      BRep_Tool::Range( theE2, f, l );
+      p2 += 1e-7 * ( p2-f < l-p2 ? +1. : -1. );
+      c2->D1( p2, p, vec2 );
+      if ( theE2.Orientation() == TopAbs_REVERSED )
+        vec2.Reverse();
+      angle = vec1.AngleWithRef( vec2, vecRef );
+    }
   }
   catch (...)
   {
@@ -4136,9 +4167,11 @@ namespace { // Structures used by FixQuadraticElements()
       //BRepClass3d_SolidClassifier solidClassifier( shape );
 
       TIDSortedElemSet checkedVols, movedNodes;
-      for ( faceIt.ReInit(); faceIt.More(); faceIt.Next() ) // loop on FACEs of a SOLID
+      //for ( faceIt.ReInit(); faceIt.More(); faceIt.Next() ) // loop on FACEs of a SOLID
+      for ( size_t iF = 0; iF < concaveFaces.size(); ++iF ) // loop on concave FACEs
       {
-        const TopoDS_Shape& face = faceIt.Current();
+        //const TopoDS_Shape& face = faceIt.Current();
+        const TopoDS_Shape& face = concaveFaces[ iF ];
         SMESHDS_SubMesh*  faceSM = meshDS->MeshElements( face );
         if ( !faceSM ) continue;
 
@@ -4153,11 +4186,16 @@ namespace { // Structures used by FixQuadraticElements()
           if ( !vertexSM ) continue;
           nodeIt = vertexSM->GetNodes();
         }
+        // get ids of sub-shapes of the FACE
+        set< int > subIDs;
+        SMESH_subMeshIteratorPtr smIt =
+          theHelper.GetMesh()->GetSubMesh( face )->getDependsOnIterator(/*includeSelf=*/true);
+        while ( smIt->more() )
+          subIDs.insert( smIt->next()->GetId() );
 
         // find suspicious volumes adjacent to the FACE
         vector< const SMDS_MeshNode* >    nOnFace( 4 );
         const SMDS_MeshNode*              nInSolid;
-        //vector< const SMDS_MeshElement* > intersectedFaces;
         while ( nodeIt->more() )
         {
           const SMDS_MeshNode* n     = nodeIt->next();
@@ -4179,8 +4217,10 @@ namespace { // Structures used by FixQuadraticElements()
               n = *volNode;
               if ( n->GetPosition()->GetDim() == 3 )
                 nInSolid = n;
-              else
+              else if ( subIDs.count( n->getshapeId() ))
                 nOnFace.push_back( n );
+              else
+                nInSolid = n;
             }
             if ( !nInSolid || nOnFace.size() != nbN - 1 )
               continue;