Salome HOME
22568: [CEA 1151] ConvertToQuadratic does not work if a face was not assigned with...
[modules/smesh.git] / src / SMESH / SMESH_MesherHelper.cxx
index 6a97d47246cc5928909c9a7037536ebf849d245b..784853234b4157551317d0220806b50bffa94f06 100644 (file)
@@ -1200,10 +1200,11 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1,
   }
   else // ( force3d || F.IsNull() )
   {
-    P = ( SMESH_TNodeXYZ( n1 ) +
-          SMESH_TNodeXYZ( n2 ) +
-          SMESH_TNodeXYZ( n3 ) +
-          SMESH_TNodeXYZ( n4 ) ) / 4;
+    P = calcTFI (0.5, 0.5,
+                 SMESH_TNodeXYZ(n1),  SMESH_TNodeXYZ(n2),
+                 SMESH_TNodeXYZ(n3),  SMESH_TNodeXYZ(n4), 
+                 SMESH_TNodeXYZ(n12), SMESH_TNodeXYZ(n23),
+                 SMESH_TNodeXYZ(n34), SMESH_TNodeXYZ(n41));
     centralNode = meshDS->AddNode( P.X(), P.Y(), P.Z() );
 
     if ( !F.IsNull() ) // force3d
@@ -1435,6 +1436,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
     catch ( Standard_Failure& f )
     {
       // issue 22502 / a node is on VERTEX not belonging to E
+      // issue 22568 / both nodes are on non-connected VERTEXes
       return getMediumNodeOnComposedWire(n1,n2,force3d);
     }
   }
@@ -1534,67 +1536,108 @@ const SMDS_MeshNode* SMESH_MesherHelper::getMediumNodeOnComposedWire(const SMDS_
                                                                      const SMDS_MeshNode* n2,
                                                                      bool                 force3d)
 {
-  gp_Pnt middle = 0.5 * XYZ(n1) + 0.5 * XYZ(n2);
+  SMESH_TNodeXYZ p1( n1 ), p2( n2 );
+  gp_Pnt      middle = 0.5 * p1 + 0.5 * p2;
   SMDS_MeshNode* n12 = AddNode( middle.X(), middle.Y(), middle.Z() );
 
   // To find position on edge and 3D position for n12,
   // project <middle> to 2 edges and select projection most close to <middle>
 
-  double u = 0, distMiddleProj = Precision::Infinite(), distXYZ[4];
-  int iOkEdge = 0;
-  TopoDS_Edge edges[2];
+  TopoDS_Edge bestEdge;
+  double u = 0, distMiddleProj = Precision::Infinite(), distXYZ[4], f,l;
+
+  // get shapes under the nodes
+  TopoDS_Shape shape[2];
+  int nbShapes = 0;
   for ( int is2nd = 0; is2nd < 2; ++is2nd )
   {
-    // get an edge
     const SMDS_MeshNode* n = is2nd ? n2 : n1;
-    TopoDS_Shape shape = GetSubShapeByNode( n, GetMeshDS() );
-    if ( shape.IsNull() || shape.ShapeType() != TopAbs_EDGE )
-      continue;
+    TopoDS_Shape S = GetSubShapeByNode( n, GetMeshDS() );
+    if ( !S.IsNull() )
+      shape[ nbShapes++ ] = S;
+  }
+  // get EDGEs
+  vector< TopoDS_Shape > edges;
+  for ( int iS = 0; iS < nbShapes; ++iS )
+  {
+    switch ( shape[iS].ShapeType() ) {
+    case TopAbs_EDGE:
+    {
+      edges.push_back( shape[iS] );
+      break;
+    }
+    case TopAbs_VERTEX:
+    {
+      TopoDS_Shape edge;
+      if ( nbShapes == 2 && iS==0 && shape[1-iS].ShapeType() == TopAbs_VERTEX )
+        edge = GetCommonAncestor( shape[iS], shape[1-iS], *myMesh, TopAbs_EDGE );
 
-    // project to get U of projection and distance from middle to projection
-    TopoDS_Edge edge = edges[ is2nd ] = TopoDS::Edge( shape );
-    double node2MiddleDist = middle.Distance( XYZ(n) );
-    double foundU = GetNodeU( edge, n );
-    CheckNodeU( edge, n12, foundU, 2*BRep_Tool::Tolerance(edge), /*force=*/true, distXYZ );
-    if ( distXYZ[0] < node2MiddleDist )
+      if ( edge.IsNull() )
+      {
+        PShapeIteratorPtr eIt = GetAncestors( shape[iS], *myMesh, TopAbs_EDGE );
+        while( const TopoDS_Shape* e = eIt->next() )
+          edges.push_back( *e );
+      }
+      break;
+    }
+    case TopAbs_FACE:
     {
-      distMiddleProj = distXYZ[0];
-      u = foundU;
-      iOkEdge = is2nd;
+      if ( nbShapes == 1 || shape[1-iS].ShapeType() < TopAbs_EDGE )
+        for ( TopExp_Explorer e( shape[iS], TopAbs_EDGE ); e.More(); e.Next() )
+          edges.push_back( e.Current() );
+      break;
+    }
+    default:
+      continue;
     }
   }
-  if ( Precision::IsInfinite( distMiddleProj ))
+  // project to get U of projection and distance from middle to projection
+  for ( size_t iE = 0; iE < edges.size(); ++iE )
   {
-    // both projections failed; set n12 on the edge of n1 with U of a common vertex
-    TopoDS_Vertex vCommon;
-    if ( TopExp::CommonVertex( edges[0], edges[1], vCommon ))
-      u = BRep_Tool::Parameter( vCommon, edges[0] );
-    else
+    const TopoDS_Edge& edge = TopoDS::Edge( edges[ iE ]);
+    distXYZ[0] = distMiddleProj;
+    double testU = 0;
+    CheckNodeU( edge, n12, testU, 2 * BRep_Tool::Tolerance(edge), /*force=*/true, distXYZ );
+    if ( distXYZ[0] < distMiddleProj )
     {
-      double f,l, u0 = GetNodeU( edges[0], n1 );
-      BRep_Tool::Range( edges[0],f,l );
-      u = ( fabs(u0-f) < fabs(u0-l) ) ? f : l;
+      distMiddleProj = distXYZ[0];
+      u = testU;
+      bestEdge = edge;
     }
-    iOkEdge = 0;
-    distMiddleProj = 0;
   }
-
-  // move n12 to position of a successfull projection
-  double tol = BRep_Tool::Tolerance(edges[ iOkEdge ]);
-  if ( !force3d && distMiddleProj > 2*tol )
+  // {
+  //   // both projections failed; set n12 on the edge of n1 with U of a common vertex
+  //   TopoDS_Vertex vCommon;
+  //   if ( TopExp::CommonVertex( edges[0], edges[1], vCommon ))
+  //     u = BRep_Tool::Parameter( vCommon, edges[0] );
+  //   else
+  //   {
+  //     double f,l, u0 = GetNodeU( edges[0], n1 );
+  //     BRep_Tool::Range( edges[0],f,l );
+  //     u = ( fabs(u0-f) < fabs(u0-l) ) ? f : l;
+  //   }
+  //   iOkEdge = 0;
+  //   distMiddleProj = 0;
+  // }
+
+  if ( !bestEdge.IsNull() )
   {
-    TopLoc_Location loc; double f,l;
-    Handle(Geom_Curve) curve = BRep_Tool::Curve( edges[iOkEdge],loc,f,l );
-    gp_Pnt p = curve->Value( u );
-    GetMeshDS()->MoveNode( n12, p.X(), p.Y(), p.Z() );
-  }
-
-  //if ( mySetElemOnShape ) node is not elem!
-  {
-    int edgeID = GetMeshDS()->ShapeToIndex( edges[iOkEdge] );
-    if ( edgeID != n12->getshapeId() )
-      GetMeshDS()->UnSetNodeOnShape( n12 );
-    GetMeshDS()->SetNodeOnEdge(n12, edgeID, u);
+    // move n12 to position of a successfull projection
+    //double tol = BRep_Tool::Tolerance(edges[ iOkEdge ]);
+    if ( !force3d /*&& distMiddleProj > 2*tol*/ )
+    {
+      TopLoc_Location loc;
+      Handle(Geom_Curve) curve = BRep_Tool::Curve( bestEdge,loc,f,l );
+      gp_Pnt p = curve->Value( u ).Transformed( loc );
+      GetMeshDS()->MoveNode( n12, p.X(), p.Y(), p.Z() );
+    }
+    //if ( mySetElemOnShape ) node is not elem!
+    {
+      int edgeID = GetMeshDS()->ShapeToIndex( bestEdge );
+      if ( edgeID != n12->getshapeId() )
+        GetMeshDS()->UnSetNodeOnShape( n12 );
+      GetMeshDS()->SetNodeOnEdge(n12, edgeID, u);
+    }
   }
   myTLinkNodeMap.insert( make_pair( SMESH_TLink(n1,n2), n12 ));
 
@@ -2932,7 +2975,7 @@ namespace { // Structures used by FixQuadraticElements()
 //=======================================================================
 
 #define __DMP__(txt) \
-  //cout << txt
+  // cout << txt
 #define MSG(txt) __DMP__(txt<<endl)
 #define MSGBEG(txt) __DMP__(txt)
 
@@ -3340,27 +3383,11 @@ namespace { // Structures used by FixQuadraticElements()
 
   gp_Vec QFace::LinkNorm(const int i, SMESH_MesherHelper* /*uvHelper*/) const
   {
-    gp_Vec norm, vecOut;
-//     if ( uvHelper ) {
-//       TopoDS_Face face = TopoDS::Face( uvHelper->GetSubShape());
-//       const SMDS_MeshNode* inFaceNode = uvHelper->GetNodeUVneedInFaceNode() ? GetNodeInFace() : 0;
-//       gp_XY uv1 = uvHelper->GetNodeUV( face, _sides[i]->node1(), inFaceNode );
-//       gp_XY uv2 = uvHelper->GetNodeUV( face, _sides[i]->node2(), inFaceNode );
-//       norm.SetCoord( uv1.Y() - uv2.Y(), uv2.X() - uv1.X(), 0 );
-
-//       const QLink* otherLink = _sides[(i + 1) % _sides.size()];
-//       const SMDS_MeshNode* otherNode =
-//         otherLink->node1() == _sides[i]->node1() ? otherLink->node2() : otherLink->node1();
-//       gp_XY pIn = uvHelper->GetNodeUV( face, otherNode, inFaceNode );
-//       vecOut.SetCoord( uv1.X() - pIn.X(), uv1.Y() - pIn.Y(), 0 );
-//     }
-//     else {
-      norm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2()));
-      gp_XYZ pIn = ( XYZ( _sides[0]->node1() ) +
-                     XYZ( _sides[0]->node2() ) +
-                     XYZ( _sides[1]->node1() )) / 3.;
-      vecOut.SetXYZ( _sides[i]->MiddlePnt() - pIn );
-      //}
+    gp_Vec   norm = _normal ^ gp_Vec( XYZ(_sides[i]->node1()), XYZ(_sides[i]->node2()));
+    gp_XYZ    pIn = ( _sides[ (i+1)%3 ]->MiddlePnt() +
+                      _sides[ (i+2)%3 ]->MiddlePnt() ) / 2.;
+    gp_Vec vecOut = ( _sides[i]->MiddlePnt() - pIn );
+
     if ( norm * vecOut < 0 )
       norm.Reverse();
     double mag2 = norm.SquareMagnitude();
@@ -3414,10 +3441,26 @@ namespace { // Structures used by FixQuadraticElements()
     int iL1 = (iL + 1) % 3, iL2 = (iL + 2) % 3; // indices of the two other links of triangle
     TLinkInSet link1 = theLinks.find( _sides[iL1] );
     TLinkInSet link2 = theLinks.find( _sides[iL2] );
-    if ( link1 == theLinks.end() || link2 == theLinks.end() )
-      return thePrevLen;
-    const QFace* f1 = link1->NextFace( this ); // adjacent faces
-    const QFace* f2 = link2->NextFace( this );
+
+    const QFace *f1 = 0, *f2 = 0;  // adjacent faces
+    bool isBndLink1 = true, isBndLink2 = true;
+    if ( link1 != theLinks.end() && link2 != theLinks.end() )
+    {
+      f1 = link1->NextFace( this );
+      f2 = link2->NextFace( this );
+
+      isBndLink1 = ( theLink->MediumPos() > (*link1)->MediumPos() );
+      isBndLink2 = ( theLink->MediumPos() > (*link2)->MediumPos() );
+      if ( theStep == theFirstStep ) // (issue 22541) quad-dominant mesh
+      {
+        if ( !isBndLink1 && !f1 )
+          f1 = (*link1)->GetContinuesFace( this ); // get a quadrangle face
+        if ( !isBndLink2 && !f2 )
+          f2 = (*link2)->GetContinuesFace( this );
+      }
+    }
+    else if ( _sides.size() < 4 )
+      return thePrevLen;      
 
     // propagate to adjacent faces till limit step or boundary
     double len1 = thePrevLen + (theLink->MiddlePnt() - _sides[iL1]->MiddlePnt()).Modulus();
@@ -3426,7 +3469,7 @@ namespace { // Structures used by FixQuadraticElements()
     gp_Vec linkDir2(0,0,0);
     try {
       OCC_CATCH_SIGNALS;
-      if ( f1 && theLink->MediumPos() <= (*link1)->MediumPos() )
+      if ( f1 && !isBndLink1 )
         len1 = f1->MoveByBoundary
           ( *link1, theRefVec, theLinks, theFaceHelper, len1, theStep-1, &linkDir1, theSign);
       else
@@ -3437,7 +3480,7 @@ namespace { // Structures used by FixQuadraticElements()
     }
     try {
       OCC_CATCH_SIGNALS;
-      if ( f2 && theLink->MediumPos() <= (*link2)->MediumPos() )
+      if ( f2 && !isBndLink2 )
         len2 = f2->MoveByBoundary
           ( *link2, theRefVec, theLinks, theFaceHelper, len2, theStep-1, &linkDir2, theSign);
       else
@@ -3460,7 +3503,7 @@ namespace { // Structures used by FixQuadraticElements()
 
       MSG(string(theStep,'.')<<" Move "<< theLink->_mediumNode->GetID()<<
           " by " << refProj * ( 1 - r ) << " following " <<
-          (choose1 ? *link1->_qlink : *link2->_qlink));
+          (choose1 ? *link1->_qlink : *link2->_qlink)); // warning: link1 can be invalid
 
       if ( theLinkNorm ) *theLinkNorm = linkNorm;
     }
@@ -3550,7 +3593,7 @@ namespace { // Structures used by FixQuadraticElements()
     }
     else if ( _faces.size() > 1 ) // not found, set NULL by the first face
     {
-      _faces.insert( ++_faces.begin(), 0 );
+      _faces.insert( ++_faces.begin(), (QFace*) 0 );
     }
   }
   //================================================================================