#16843 EDF 19340 - wrong quadratic mesh

src/SMESH/SMESH_Algo.cxx

@@ -1268,3 +1268,33 @@ bool SMESH_2D_Algo::FixInternalNodes(const SMESH_ProxyMesh& mesh,
   }
   return true;
 }
+
+//=======================================================================
+//function : IsApplicableToShape
+//purpose  : Return true if the algorithm can mesh a given shape
+//=======================================================================
+
+bool SMESH_1D_Algo::IsApplicableToShape(const TopoDS_Shape & shape, bool toCheckAll) const
+{
+  return ( !shape.IsNull() && TopExp_Explorer( shape, TopAbs_EDGE ).More() );
+}
+
+//=======================================================================
+//function : IsApplicableToShape
+//purpose  : Return true if the algorithm can mesh a given shape
+//=======================================================================
+
+bool SMESH_2D_Algo::IsApplicableToShape(const TopoDS_Shape & shape, bool toCheckAll) const
+{
+  return ( !shape.IsNull() && TopExp_Explorer( shape, TopAbs_FACE ).More() );
+}
+
+//=======================================================================
+//function : IsApplicableToShape
+//purpose  : Return true if the algorithm can mesh a given shape
+//=======================================================================
+
+bool SMESH_3D_Algo::IsApplicableToShape(const TopoDS_Shape & shape, bool toCheckAll) const
+{
+  return ( !shape.IsNull() && TopExp_Explorer( shape, TopAbs_SOLID ).More() );
+}

src/SMESH/SMESH_Algo.hxx

@@ -482,12 +482,28 @@ class SMESH_EXPORT SMESH_1D_Algo: public SMESH_Algo
 {
 public:
   SMESH_1D_Algo(int hypId, SMESH_Gen* gen);
+  /*!
+   * \brief Return true if the algorithm can mesh a given shape
+   *  \param [in] aShape - shape to check
+   *  \param [in] toCheckAll - if true, this check returns OK if all shapes are OK,
+   *              else, returns OK if at least one shape is OK
+   *  \retval bool - \c true by default
+   */
+  virtual bool IsApplicableToShape(const TopoDS_Shape & shape, bool toCheckAll) const;
 };
 
 class SMESH_EXPORT SMESH_2D_Algo: public SMESH_Algo
 {
 public:
   SMESH_2D_Algo(int hypId, SMESH_Gen* gen);
+  /*!
+   * \brief Return true if the algorithm can mesh a given shape
+   *  \param [in] aShape - shape to check
+   *  \param [in] toCheckAll - if true, this check returns OK if all shapes are OK,
+   *              else, returns OK if at least one shape is OK
+   *  \retval bool - \c true by default
+   */
+  virtual bool IsApplicableToShape(const TopoDS_Shape & shape, bool toCheckAll) const;
   /*!
    * \brief Method in which an algorithm generating a structured mesh
    *        fixes positions of in-face nodes after there movement
@@ -501,6 +517,14 @@ class SMESH_EXPORT SMESH_3D_Algo: public SMESH_Algo
 {
 public:
   SMESH_3D_Algo(int hypId, SMESH_Gen* gen);
+  /*!
+   * \brief Return true if the algorithm can mesh a given shape
+   *  \param [in] aShape - shape to check
+   *  \param [in] toCheckAll - if true, this check returns OK if all shapes are OK,
+   *              else, returns OK if at least one shape is OK
+   *  \retval bool - \c true by default
+   */
+  virtual bool IsApplicableToShape(const TopoDS_Shape & shape, bool toCheckAll) const;
 };
 
 #endif

src/SMESH/SMESH_MesherHelper.cxx

@@ -1675,28 +1675,6 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
   {
     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 )
@@ -1728,8 +1706,6 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
     // we try to create medium node using UV parameters of
     // 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 ));
@@ -1743,13 +1719,32 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
       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 ) )
+      {
+        // 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() )
     {
       double f,l;

src/SMESH_SWIG/smeshBuilder.py

@@ -4490,8 +4490,8 @@ class Mesh(metaclass = MeshMeta):
         Reorient faces contained in *the2DObject*.
 
         Parameters:
-                the2DObject: is a :class:`mesh, sub-mesh, group or filter <SMESH.SMESH_IDSource>` or list of IDs of 2D elements
+                the2DObject: a :class:`mesh, sub-mesh, group or filter <SMESH.SMESH_IDSource>` or list of IDs of 2D elements
-                theDirection: is a desired direction of normal of *theFace*.
+                theDirection: a desired direction of normal of *theFace*.
                         It can be either a GEOM vector or a list of coordinates [x,y,z].
                 theFaceOrPoint: defines a face of *the2DObject* whose normal will be
                         compared with theDirection. It can be either ID of face or a point

src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx

@@ -258,6 +258,7 @@ namespace
           meshDS->RemoveFreeNode( nodesToRemove[i], sm, /*fromGroups=*/false);
       }
     }
+    return;
   }
 
   //================================================================================
@@ -319,7 +320,7 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement*     Pr
                                                   set<const SMDS_MeshNode*> & nodesToMove)
 {
   // cout << endl << "Merge " << PrmI->GetID() << " " << PrmJ->GetID() << " "
-  //      << PrmI->GetNode(4) << PrmJ->GetNode(4) << endl;
+  //      << PrmI->GetNode(4)->GetID() << " " << PrmJ->GetNode(4)->GetID() << endl;
   const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove
   //int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume );
   SMESH_TNodeXYZ Pj( Nrem );
@@ -338,6 +339,9 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement*     Pr
   typedef SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > TStdElemIterator;
   TStdElemIterator itEnd;
 
+  typedef std::map< const SMDS_MeshNode*, const SMDS_MeshNode* > TNNMap;
+  TNNMap mediumReplaceMap;
+
   // find and remove coincided faces of merged pyramids
   vector< const SMDS_MeshElement* > inverseElems
     // copy inverse elements to avoid iteration on changing container
@@ -355,6 +359,11 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement*     Pr
     }
     if ( FJEqual )
     {
+      if ( FJEqual->NbNodes() == 6 ) // find medium nodes to replace
+      {
+        mediumReplaceMap.insert( std::make_pair( FJEqual->GetNode(3), FI->GetNode(5) ));
+        mediumReplaceMap.insert( std::make_pair( FJEqual->GetNode(5), FI->GetNode(3) ));
+      }
       removeTmpElement( FI );
       removeTmpElement( FJEqual );
       myRemovedTrias.insert( FI );
@@ -370,12 +379,27 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement*     Pr
     const SMDS_MeshElement* elem = inverseElems[i];
     nodes.assign( elem->begin_nodes(), elem->end_nodes() );
     nodes[ elem->GetType() == SMDSAbs_Volume ? PYRAM_APEX : TRIA_APEX ] = CommonNode;
+    if ( !mediumReplaceMap.empty() )
+      for ( size_t iN = elem->NbCornerNodes(); iN < nodes.size(); ++iN )
+      {
+        TNNMap::iterator n2n = mediumReplaceMap.find( nodes[iN] );
+        if ( n2n != mediumReplaceMap.end() )
+          nodes[iN] = n2n->second;
+      }
     GetMeshDS()->ChangeElementNodes( elem, &nodes[0], nodes.size());
   }
   ASSERT( Nrem->NbInverseElements() == 0 );
   GetMeshDS()->RemoveFreeNode( Nrem,
                                GetMeshDS()->MeshElements( Nrem->getshapeId()),
                                /*fromGroups=*/false);
+  if ( !mediumReplaceMap.empty() )
+    for ( TNNMap::iterator n2n = mediumReplaceMap.begin(); n2n != mediumReplaceMap.end(); ++n2n )
+    {
+      const SMDS_MeshNode* remNode = n2n->first;
+      if ( !remNode->IsNull() && remNode->NbInverseElements() == 0 )
+        GetMeshDS()->RemoveFreeNode( remNode, 0, /*fromGroups=*/false);
+    }
+  return;
 }
 
 //================================================================================
@@ -414,6 +438,7 @@ void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement*    PrmI
     for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm)
       MergeAdjacent( *prm, nodesToMove, true );
   }
+  return;
 }
 
 //================================================================================
@@ -522,7 +547,27 @@ static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
 
   Pint = orig + dir * t;
 
-  return ( t > 0.  &&  t < segLen );
+  bool hasInt = ( t > 0.  &&  t < segLen );
+
+  if ( hasInt && det < EPSILON ) // t is inaccurate, additionally check
+  {
+    gp_XYZ triNorm = edge1 ^ edge2;
+    gp_XYZ int0vec = Pint.XYZ() - vert0;
+    gp_XYZ in      = triNorm ^ edge1; // dir inside triangle from edge1
+    double dot     = int0vec * in;
+    if ( dot < 0 && dot / triNorm.Modulus() < -EPSILON )
+      return false;
+    in  = edge2 ^ triNorm;
+    dot = int0vec * in;
+    if ( dot < 0 && dot / triNorm.Modulus() < -EPSILON )
+      return false;
+    gp_XYZ int1vec = Pint.XYZ() - vert1;
+    in  = triNorm ^ ( vert2 - vert1 );
+    dot = int1vec * in;
+    if ( dot < 0 && dot / triNorm.Modulus() < -EPSILON )
+      return false;
+  }
+  return hasInt;
 }
 
 //=======================================================================
@@ -817,7 +862,6 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&         aMesh,
   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
   SMESH_MesherHelper helper(aMesh);
   helper.IsQuadraticSubMesh(aShape);
-  helper.SetElementsOnShape( true );
 
   if ( myElemSearcher ) delete myElemSearcher;
   vector< SMDS_ElemIteratorPtr > itVec;
@@ -877,10 +921,11 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&         aMesh,
             // degenerate face
             // add triangles to result map
             SMDS_MeshFace* NewFace;
+            helper.SetElementsOnShape( false );
             if(!isRev)
-              NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] );
+              NewFace = helper.AddFace( FNodes[0], FNodes[1], FNodes[2] );
             else
-              NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] );
+              NewFace = helper.AddFace( FNodes[0], FNodes[2], FNodes[1] );
             storeTmpElement( NewFace );
             trias.push_back ( NewFace );
             quads.push_back( face );
@@ -916,20 +961,25 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&         aMesh,
               if ( !LimitHeight( PCbest, PC, PN, FNodes, aMesh, face, /*UseApexRay=*/true, aShape ))
                 return false;
             }
-            // create node for PCbest
+            // create node at PCbest
+            helper.SetElementsOnShape( true );
             SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
 
+            // create a pyramid
+            SMDS_MeshVolume* aPyram;
+            if ( isRev )
+              aPyram = helper.AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode );
+            else
+              aPyram = helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
+            myPyramids.push_back(aPyram);
+
             // add triangles to result map
+            helper.SetElementsOnShape( false );
             for ( i = 0; i < 4; i++ )
             {
-              trias.push_back ( meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] ));
+              trias.push_back ( helper.AddFace( NewNode, FNodes[i], FNodes[i+1] ));
               storeTmpElement( trias.back() );
             }
-            // create a pyramid
-            if ( isRev ) swap( FNodes[1], FNodes[3]);
-            SMDS_MeshVolume* aPyram =
-              helper.AddVolume( FNodes[0], FNodes[1], FNodes[2], FNodes[3], NewNode );
-            myPyramids.push_back(aPyram);
 
             quads.push_back( face );
             hasNewTrias = true;
@@ -1018,7 +1068,6 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
   vector<const SMDS_MeshElement*> myPyramids;
   SMESH_MesherHelper helper(aMesh);
   helper.IsQuadraticSubMesh(aMesh.GetShapeToMesh());
-  helper.SetElementsOnShape( true );
 
   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
   SMESH_ProxyMesh::SubMesh* prxSubMesh = getProxySubMesh();
@@ -1116,10 +1165,11 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
           IsRev = true;
         }
       }
+      helper.SetElementsOnShape( false );
       if(!IsRev)
-        NewFace = meshDS->AddFace( FNodes[0], FNodes[1], FNodes[2] );
+        NewFace = helper.AddFace( FNodes[0], FNodes[1], FNodes[2] );
       else
-        NewFace = meshDS->AddFace( FNodes[0], FNodes[2], FNodes[1] );
+        NewFace = helper.AddFace( FNodes[0], FNodes[2], FNodes[1] );
       storeTmpElement( NewFace );
       prxSubMesh->AddElement( NewFace );
       continue;
@@ -1238,18 +1288,9 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
         return false;
 
       // create node for Papex
+      helper.SetElementsOnShape( true );
       SMDS_MeshNode* NewNode = helper.AddNode( Papex.X(), Papex.Y(), Papex.Z() );
 
-      // add triangles to result map
-      for ( i = 0; i < 4; i++) {
-        SMDS_MeshFace* NewFace;
-        if(isRev)
-          NewFace = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] );
-        else
-          NewFace = meshDS->AddFace( NewNode, FNodes[i+1], FNodes[i] );
-        storeTmpElement( NewFace );
-        prxSubMesh->AddElement( NewFace );
-      }
       // create a pyramid
       SMDS_MeshVolume* aPyram;
       if(isRev)
@@ -1257,6 +1298,18 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
       else
         aPyram = helper.AddVolume( FNodes[0], FNodes[3], FNodes[2], FNodes[1], NewNode );
       myPyramids.push_back(aPyram);
+
+      // add triangles to result map
+      helper.SetElementsOnShape( false );
+      for ( i = 0; i < 4; i++) {
+        SMDS_MeshFace* NewFace;
+        if(isRev)
+          NewFace = helper.AddFace( NewNode, FNodes[i], FNodes[i+1] );
+        else
+          NewFace = helper.AddFace( NewNode, FNodes[i+1], FNodes[i] );
+        storeTmpElement( NewFace );
+        prxSubMesh->AddElement( NewFace );
+      }
     }
   } // end loop on all faces