Salome HOME
[bos #20282] EDF 22320 - general compute fails
[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.cxx
index 5925964b06d03df0014bdc17f291ce85bb60122d..17c9d9047585186101e401967d437630aaf40b3e 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2020  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
@@ -36,6 +36,7 @@
 #include "SMESH_HypoFilter.hxx"
 #include "SMESH_Mesh.hxx"
 #include "SMESH_MeshAlgos.hxx"
+#include "SMESH_MeshEditor.hxx"
 #include "SMESH_ProxyMesh.hxx"
 #include "SMESH_subMesh.hxx"
 
@@ -73,7 +74,7 @@ using namespace std;
 
 namespace {
 
-  inline SMESH_TNodeXYZ XYZ(const SMDS_MeshNode* n) { return SMESH_TNodeXYZ(n); }
+  inline SMESH_NodeXYZ XYZ(const SMDS_MeshNode* n) { return SMESH_NodeXYZ(n); }
 
   enum { U_periodic = 1, V_periodic = 2 };
 }
@@ -150,8 +151,6 @@ bool SMESH_MesherHelper::IsQuadraticSubMesh(const TopoDS_Shape& aSh)
   // we can create quadratic elements only if all elements
   // created on sub-shapes of given shape are quadratic
   myCreateQuadratic = true;
-  mySeamShapeIds.clear();
-  myDegenShapeIds.clear();
   TopAbs_ShapeEnum subType( aSh.ShapeType()==TopAbs_FACE ? TopAbs_EDGE : TopAbs_FACE );
   if ( aSh.ShapeType()==TopAbs_COMPOUND )
   {
@@ -740,7 +739,8 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
               if ( !C2d.IsNull() ) {
                 double u = ( V == IthVertex( 0, edge )) ?  f : l;
                 uv = C2d->Value( u );
-                uvOK = true;
+                gp_Pnt p = GetSurface( F )->Value( uv );
+                uvOK = ( p.Distance( BRep_Tool::Pnt( V )) < getFaceMaxTol( F ));
                 break;
               }
             }
@@ -1035,6 +1035,16 @@ double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge&   E,
       int vertexID = n->getshapeId();
       const TopoDS_Vertex& V = TopoDS::Vertex(meshDS->IndexToShape(vertexID));
       param =  BRep_Tool::Parameter( V, E );
+
+      if ( inEdgeNode )
+      {
+        BRepAdaptor_Curve curve( E );
+        if ( curve.IsPeriodic() )
+        {
+          double uInEdge = GetNodeU( E, inEdgeNode );
+          param += ShapeAnalysis::AdjustByPeriod( param, uInEdge, curve.Period() );
+        }
+      }
     }
   }
   if ( check )
@@ -1662,30 +1672,8 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
   // get positions of the given nodes on shapes
   if ( pos.second == TopAbs_FACE )
   {
-    F = TopoDS::Face(meshDS->IndexToShape( faceID = pos.first ));
+    F = TopoDS::Face( meshDS->IndexToShape( faceID = pos.first ));
     uv[0] = GetNodeUV(F,n1,n2, force3d ? 0 : &uvOK[0]);
-    if (( !force3d ) &&
-        ( HasDegeneratedEdges() || GetSurface( F )->HasSingularities( 1e-7 )))
-    {
-      // IPAL52850 (degen VERTEX not at singularity)
-      // project middle point to a surface
-      SMESH_TNodeXYZ p1( n1 ), p2( n2 );
-      gp_Pnt pMid = 0.5 * ( p1 + p2 );
-      Handle(ShapeAnalysis_Surface) projector = GetSurface( F );
-      gp_Pnt2d uvMid;
-      if ( uvOK[0] )
-        uvMid = projector->NextValueOfUV( uv[0], pMid, BRep_Tool::Tolerance( F ));
-      else
-        uvMid = projector->ValueOfUV( pMid, getFaceMaxTol( F ));
-      if ( projector->Gap() * projector->Gap() < ( p1 - p2 ).SquareModulus() / 4 )
-      {
-        gp_Pnt pProj = projector->Value( uvMid );
-        n12  = meshDS->AddNode( pProj.X(), pProj.Y(), pProj.Z() );
-        meshDS->SetNodeOnFace( n12, faceID, uvMid.X(), uvMid.Y() );
-        myTLinkNodeMap.insert( make_pair ( link, n12 ));
-        return n12;
-      }
-    }
     uv[1] = GetNodeUV(F,n2,n1, force3d ? 0 : &uvOK[1]);
   }
   else if ( pos.second == TopAbs_EDGE )
@@ -1718,26 +1706,43 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
     // nodes, else - medium between corresponding 3d points
     if( ! F.IsNull() )
     {
-      //if ( uvOK[0] && uvOK[1] )
+      if ( IsDegenShape( n1->getshapeId() )) {
+        if ( myParIndex & U_periodic ) uv[0].SetCoord( 1, uv[1].Coord( 1 ));
+        else                           uv[0].SetCoord( 2, uv[1].Coord( 2 ));
+      }
+      else if ( IsDegenShape( n2->getshapeId() )) {
+        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] );
+      gp_Pnt P = S->Value( UV.X(), UV.Y() ).Transformed(loc);
+
+      SMESH_TNodeXYZ p1( n1 ), p2( n2 );
+      gp_Pnt pMid = 0.5 * ( p1 + p2 );
+      double distMid = pMid.SquareDistance( P );
+      double dist12  = ( p1 - p2 ).SquareModulus();
+      Handle(ShapeAnalysis_Surface) surfInfo = GetSurface( F );
+      if ( distMid > dist12 ||
+           HasDegeneratedEdges() ||
+           surfInfo->HasSingularities( 1e-7 ) )
       {
-        if ( IsDegenShape( n1->getshapeId() )) {
-          if ( myParIndex & U_periodic ) uv[0].SetCoord( 1, uv[1].Coord( 1 ));
-          else                           uv[0].SetCoord( 2, uv[1].Coord( 2 ));
-        }
-        else if ( IsDegenShape( n2->getshapeId() )) {
-          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] );
-        gp_Pnt P = S->Value( UV.X(), UV.Y() ).Transformed(loc);
-        n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
-        // if ( mySetElemOnShape ) node is not elem!
-        meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y());
-        myTLinkNodeMap.insert(make_pair(link,n12));
-        return n12;
+        // IPAL52850 (degen VERTEX not at singularity)
+        // project middle point to a surface
+        gp_Pnt2d uvMid;
+        if ( uvOK[0] )
+          uvMid = surfInfo->NextValueOfUV( uv[0], pMid, BRep_Tool::Tolerance( F ));
+        else
+          uvMid = surfInfo->ValueOfUV( pMid, getFaceMaxTol( F ));
+        if ( surfInfo->Gap() * surfInfo->Gap() < distMid )
+          P = surfInfo->Value( uvMid );
       }
+      n12 = meshDS->AddNode(P.X(), P.Y(), P.Z());
+      // if ( mySetElemOnShape ) node is not elem!
+      meshDS->SetNodeOnFace(n12, faceID, UV.X(), UV.Y());
+      myTLinkNodeMap.insert(make_pair(link,n12));
+      return n12;
     }
     else if ( !E.IsNull() )
     {
@@ -4640,6 +4645,7 @@ namespace { // Structures used by FixQuadraticElements()
             const SMDS_MeshElement* f = faceIt->next();
             if ( !faceSM->Contains( f ) ||
                  f->NbNodes() < 6       || // check quadratic triangles only
+                 f->NbNodes() > 7       ||
                  !checkedFaces.insert( f ).second )
               continue;
 
@@ -4851,7 +4857,12 @@ namespace { // Structures used by FixQuadraticElements()
                     gp_Pnt pMedium = SMESH_TNodeXYZ( linkIt->second );
                     double hMedium = faceNorm * gp_Vec( pOnFace0, pMedium ).XYZ();
                     double hVol    = faceNorm * gp_Vec( pOnFace0, pInSolid ).XYZ();
-                    isDistorted = ( Abs( hMedium ) > Abs( hVol * 0.75 ));
+                    if ( Abs( hMedium ) > Abs( hVol * 0.75 ))
+                    {
+                      SMESH_TNodeXYZ pI( nOnFace[i]), pJ( nOnFace[j]);
+                      double angle = gp_Vec( pI, pMedium ).Angle( gp_Vec( pI, pJ ));
+                      isDistorted  = ( angle > M_PI / 20 );
+                    }
                   }
                 }
               }