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 decd8803a16d4513fce62fa5c6b669123eae5aea..784853234b4157551317d0220806b50bffa94f06 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2014  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
@@ -6,7 +6,7 @@
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
 // License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// version 2.1 of the License, or (at your option) any later version.
 //
 // This library is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
@@ -245,8 +245,8 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
     TopLoc_Location loc;
     Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc );
 
-    if ( surface->IsUPeriodic() || surface->IsVPeriodic() ||
-         surface->IsUClosed()   || surface->IsVClosed() )
+    // if ( surface->IsUPeriodic() || surface->IsVPeriodic() ||
+    //      surface->IsUClosed()   || surface->IsVClosed() )
     {
       //while ( surface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface )))
       //surface = Handle(Geom_RectangularTrimmedSurface)::DownCast( surface )->BasisSurface();
@@ -281,7 +281,7 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh)
         }
 
         // look for a degenerated edge
-        if ( BRep_Tool::Degenerated( edge )) {
+        if ( SMESH_Algo::isDegenerated( edge )) {
           myDegenShapeIds.insert( meshDS->ShapeToIndex( edge ));
           for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() )
             myDegenShapeIds.insert( meshDS->ShapeToIndex( v.Current() ));
@@ -478,7 +478,10 @@ bool SMESH_MesherHelper::toCheckPosOnShape(int shapeID ) const
 
 void SMESH_MesherHelper::setPosOnShapeValidity(int shapeID, bool ok ) const
 {
-  ((SMESH_MesherHelper*)this)->myNodePosShapesValidity.insert( make_pair( shapeID, ok));
+  std::map< int,bool >::iterator sh_ok = 
+    ((SMESH_MesherHelper*)this)->myNodePosShapesValidity.insert( make_pair( shapeID, ok)).first;
+  if ( !ok )
+    sh_ok->second = ok;
 }
 
 //=======================================================================
@@ -600,10 +603,8 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face&   F,
         for ( TopExp_Explorer vert(F,TopAbs_VERTEX); !uvOK && vert.More(); vert.Next() )
           uvOK = ( V == vert.Current() );
         if ( !uvOK ) {
-#ifdef _DEBUG_
           MESSAGE ( "SMESH_MesherHelper::GetNodeUV(); Vertex " << vertexID
-               << " not in face " << GetMeshDS()->ShapeToIndex( F ) );
-#endif
+                    << " not in face " << GetMeshDS()->ShapeToIndex( F ) );
           // get UV of a vertex closest to the node
           double dist = 1e100;
           gp_Pnt pn = XYZ( n );
@@ -663,9 +664,10 @@ bool SMESH_MesherHelper::CheckNodeUV(const TopoDS_Face&   F,
                                      const bool           force,
                                      double               distXYZ[4]) const
 {
-  int shapeID = n->getshapeId();
+  int  shapeID = n->getshapeId();
   bool infinit = ( Precision::IsInfinite( uv.X() ) || Precision::IsInfinite( uv.Y() ));
-  if ( force || toCheckPosOnShape( shapeID ) || infinit )
+  bool zero    = ( uv.X() == 0. && uv.Y() == 0. );
+  if ( force || toCheckPosOnShape( shapeID ) || infinit || zero )
   {
     // check that uv is correct
     TopLoc_Location loc;
@@ -809,6 +811,34 @@ gp_XY SMESH_MesherHelper::GetMiddleUV(const Handle(Geom_Surface)& surface,
   return applyIn2D( surf, p1, p2, & AverageUV );
 }
 
+//=======================================================================
+//function : GetCenterUV
+//purpose  : Return UV for the central node of a biquadratic triangle
+//=======================================================================
+
+gp_XY SMESH_MesherHelper::GetCenterUV(const gp_XY& uv1,
+                                      const gp_XY& uv2, 
+                                      const gp_XY& uv3, 
+                                      const gp_XY& uv12,
+                                      const gp_XY& uv23,
+                                      const gp_XY& uv31,
+                                      bool *       isBadTria/*=0*/)
+{
+  bool badTria;
+  gp_XY uvAvg = ( uv12 + uv23 + uv31 ) / 3.;
+
+  if      (( badTria = (( uvAvg - uv1 ) * ( uvAvg - uv23 ) > 0 )))
+    uvAvg = ( uv1 + uv23 ) / 2.;
+  else if (( badTria = (( uvAvg - uv2 ) * ( uvAvg - uv31 ) > 0 )))
+    uvAvg = ( uv2 + uv31 ) / 2.;
+  else if (( badTria = (( uvAvg - uv3 ) * ( uvAvg - uv12 ) > 0 )))
+    uvAvg = ( uv3 + uv12 ) / 2.;
+
+  if ( isBadTria )
+    *isBadTria = badTria;
+  return uvAvg;
+}
+
 //=======================================================================
 //function : GetNodeU
 //purpose  : Return node U on edge
@@ -817,7 +847,7 @@ gp_XY SMESH_MesherHelper::GetMiddleUV(const Handle(Geom_Surface)& surface,
 double SMESH_MesherHelper::GetNodeU(const TopoDS_Edge&   E,
                                     const SMDS_MeshNode* n,
                                     const SMDS_MeshNode* inEdgeNode,
-                                    bool*                check)
+                                    bool*                check) const
 {
   double param = Precision::Infinite();
 
@@ -870,8 +900,10 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
                                     const bool           force,
                                     double               distXYZ[4]) const
 {
-  int shapeID = n->getshapeId();
-  if ( force || toCheckPosOnShape( shapeID ))
+  int  shapeID = n->getshapeId();
+  bool infinit = Precision::IsInfinite( u );
+  bool zero    = ( u == 0. );
+  if ( force || toCheckPosOnShape( shapeID ) || infinit || zero )
   {
     TopLoc_Location loc; double f,l;
     Handle(Geom_Curve) curve = BRep_Tool::Curve( E,loc,f,l );
@@ -887,12 +919,17 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
     {
       gp_Pnt nodePnt = SMESH_TNodeXYZ( n );
       if ( !loc.IsIdentity() ) nodePnt.Transform( loc.Transformation().Inverted() );
-      gp_Pnt curvPnt = curve->Value( u );
-      double dist = nodePnt.Distance( curvPnt );
-      if ( distXYZ ) {
-        curvPnt.Transform( loc );
-        distXYZ[0] = dist;
-        distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z();
+      gp_Pnt curvPnt;
+      double dist = u;
+      if ( !infinit )
+      {
+        curvPnt = curve->Value( u );
+        dist    = nodePnt.Distance( curvPnt );
+        if ( distXYZ ) {
+          curvPnt.Transform( loc );
+          distXYZ[0] = dist;
+          distXYZ[1] = curvPnt.X(); distXYZ[2] = curvPnt.Y(); distXYZ[3]=curvPnt.Z();
+        }
       }
       if ( dist > tol )
       {
@@ -916,6 +953,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
         }
         Quantity_Parameter U = projector->LowerDistanceParameter();
         u = double( U );
+        MESSAGE(" f " << f << " l " << l << " u " << u);
         curvPnt = curve->Value( u );
         dist = nodePnt.Distance( curvPnt );
         if ( distXYZ ) {
@@ -940,6 +978,7 @@ bool SMESH_MesherHelper::CheckNodeU(const TopoDS_Edge&   E,
       }
       if (( u < f-tol || u > l+tol ) && force )
       {
+        MESSAGE("u < f-tol || u > l+tol  ; u " << u << " f " << f << " l " << l);
         // node is on vertex but is set on periodic but trimmed edge (issue 0020890)
         try
         {
@@ -1048,7 +1087,7 @@ SMESH_MesherHelper::GetMediumPos(const SMDS_MeshNode* n1,
 //purpose  : Return existing or create a new central node for a quardilateral
 //           quadratic face given its 8 nodes.
 //@param   : force3d - true means node creation in between the given nodes,
-//             else node position is found on a geometrical face if any.
+//           else node position is found on a geometrical face if any.
 //=======================================================================
 
 const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1,
@@ -1144,11 +1183,12 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1,
 
   gp_XY  uvAvg;
   gp_Pnt P;
+  bool toCheck = true;
   if ( !F.IsNull() && !force3d )
   {
     uvAvg = calcTFI (0.5, 0.5,
-                     GetNodeUV(F,n1,n3),  GetNodeUV(F,n2,n4),
-                     GetNodeUV(F,n3,n1),  GetNodeUV(F,n4,n2), 
+                     GetNodeUV(F,n1,n3,&toCheck), GetNodeUV(F,n2,n4,&toCheck),
+                     GetNodeUV(F,n3,n1,&toCheck), GetNodeUV(F,n4,n2,&toCheck), 
                      GetNodeUV(F,n12,n3), GetNodeUV(F,n23,n4),
                      GetNodeUV(F,n34,n2), GetNodeUV(F,n41,n2));
     TopLoc_Location loc;
@@ -1160,18 +1200,19 @@ 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
     {
-      uvAvg = (GetNodeUV(F,n1,n3) +
-               GetNodeUV(F,n2,n4) +
-               GetNodeUV(F,n3,n1) +
-               GetNodeUV(F,n4,n2)) / 4;
+      uvAvg = (GetNodeUV(F,n1,n3,&toCheck) +
+               GetNodeUV(F,n2,n4,&toCheck) +
+               GetNodeUV(F,n3,n1,&toCheck) +
+               GetNodeUV(F,n4,n2,&toCheck)) / 4;
       //CheckNodeUV( F, centralNode, uvAvg, 2*BRep_Tool::Tolerance( F ), /*force=*/true);
       meshDS->SetNodeOnFace( centralNode, faceID, uvAvg.X(), uvAvg.Y() );
     }
@@ -1193,7 +1234,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1,
 //purpose  : Return existing or create a new central node for a
 //           quadratic triangle given its 6 nodes.
 //@param   : force3d - true means node creation in between the given nodes,
-//             else node position is found on a geometrical face if any.
+//           else node position is found on a geometrical face if any.
 //=======================================================================
 
 const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1,
@@ -1278,21 +1319,30 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1,
   }
 
   TopoDS_Face F;
+  gp_XY       uvAvg;
+  bool        badTria=false;
+
   if ( shapeType == TopAbs_FACE )
   {
     F = TopoDS::Face( meshDS->IndexToShape( faceID ));
+    bool check;
+    gp_XY uv1  = GetNodeUV( F, n1, n23, &check );
+    gp_XY uv2  = GetNodeUV( F, n2, n31, &check );
+    gp_XY uv3  = GetNodeUV( F, n3, n12, &check );
+    gp_XY uv12 = GetNodeUV( F, n12, n3, &check );
+    gp_XY uv23 = GetNodeUV( F, n23, n1, &check );
+    gp_XY uv31 = GetNodeUV( F, n31, n2, &check );
+    uvAvg = GetCenterUV( uv1,uv2,uv3, uv12,uv23,uv31, &badTria );
+    if ( badTria )
+      force3d = false;
   }
 
-  // Create a node
+  // Create a central node
 
-  gp_XY  uvAvg;
   gp_Pnt P;
   if ( !F.IsNull() && !force3d )
   {
-    uvAvg = ( GetNodeUV(F,n12,n3) +
-              GetNodeUV(F,n23,n1) +
-              GetNodeUV(F,n31,n2) ) / 3.;
-    TopLoc_Location loc;
+    TopLoc_Location        loc;
     Handle( Geom_Surface ) S = BRep_Tool::Surface( F, loc );
     P = S->Value( uvAvg.X(), uvAvg.Y() ).Transformed( loc );
     centralNode = meshDS->AddNode( P.X(), P.Y(), P.Z() );
@@ -1308,9 +1358,6 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1,
 
     if ( !F.IsNull() ) // force3d
     {
-      uvAvg = ( GetNodeUV(F,n12,n3) +
-                GetNodeUV(F,n23,n1) +
-                GetNodeUV(F,n31,n2) ) / 3.;
       meshDS->SetNodeOnFace( centralNode, faceID, uvAvg.X(), uvAvg.Y() );
     }
     else if ( solidID > 0 )
@@ -1382,8 +1429,16 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1,
       return getMediumNodeOnComposedWire(n1,n2,force3d);
     }
     E = TopoDS::Edge(meshDS->IndexToShape( edgeID = pos.first ));
-    u[0] = GetNodeU(E,n1,n2, force3d ? 0 : &uvOK[0]);
-    u[1] = GetNodeU(E,n2,n1, force3d ? 0 : &uvOK[1]);
+    try {
+      u[0] = GetNodeU(E,n1,n2, force3d ? 0 : &uvOK[0]);
+      u[1] = GetNodeU(E,n2,n1, force3d ? 0 : &uvOK[1]);
+    }
+    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);
+    }
   }
 
   if ( !force3d & uvOK[0] && uvOK[1] )
@@ -1481,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 )
-  {
-    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() );
   }
+  // {
+  //   // 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 ( mySetElemOnShape ) node is not elem!
+  if ( !bestEdge.IsNull() )
   {
-    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 ));
 
@@ -2455,9 +2551,16 @@ bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace)
     {
       SMESH_TNodeXYZ nPnt[3];
       SMDS_ElemIteratorPtr nodesIt = elem->nodesIterator();
+      int iNodeOnFace = 0, iPosDim = SMDS_TOP_VERTEX;
       for ( int iN = 0; nodesIt->more() && iN < 3; ++iN) // loop on nodes
+      {
         nPnt[ iN ] = nodesIt->next();
-
+        if ( nPnt[ iN ]._node->GetPosition()->GetTypeOfPosition() > iPosDim )
+        {
+          iNodeOnFace = iN;
+          iPosDim = nPnt[ iN ]._node->GetPosition()->GetTypeOfPosition();
+        }
+      }
       // compute normal
       gp_Vec v01( nPnt[0], nPnt[1] ), v02( nPnt[0], nPnt[2] );
       if ( v01.SquareMagnitude() > RealSmall() &&
@@ -2465,7 +2568,7 @@ bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace)
       {
         Ne = v01 ^ v02;
         if (( normalOK = ( Ne.SquareMagnitude() > RealSmall() )))
-          uv = GetNodeUV( theFace, nPnt[0]._node, nPnt[2]._node, &normalOK );
+          uv = GetNodeUV( theFace, nPnt[iNodeOnFace]._node, 0, &normalOK );
       }
     }
   }
@@ -2612,6 +2715,76 @@ double SMESH_MesherHelper::MaxTolerance( const TopoDS_Shape& shape )
   return tol;
 }
 
+//================================================================================
+/*!
+ * \brief Return an angle between two EDGEs sharing a common VERTEX with reference
+ *        of the FACE normal
+ *  \return double - the angle (between -Pi and Pi), negative if the angle is concave,
+ *                   1e100 in case of failure
+ *  \waring Care about order of the EDGEs and their orientation to be as they are
+ *          within the FACE! Don't pass degenerated EDGEs neither!
+ */
+//================================================================================
+
+double SMESH_MesherHelper::GetAngle( const TopoDS_Edge & theE1,
+                                     const TopoDS_Edge & theE2,
+                                     const TopoDS_Face & theFace)
+{
+  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 );
+    if ( c1.IsNull() || c2.IsNull() )
+      return angle;
+    gp_XY uv = c2d1->Value( p1 ).XY();
+    gp_Vec du, dv; gp_Pnt p;
+    surf->D1( uv.X(), uv.Y(), p, du, dv );
+    gp_Vec vec1, vec2, vecRef = du ^ dv;
+    int  nbLoops = 0;
+    double p1tmp = p1;
+    while ( vecRef.SquareMagnitude() < std::numeric_limits<double>::min() )
+    {
+      double dp = ( l - f ) / 1000.;
+      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;
+      if ( ++nbLoops > 10 )
+      {
+#ifdef _DEBUG_
+        cout << "SMESH_MesherHelper::GetAngle(): Captured in a sigularity" << endl;
+#endif
+        return angle;
+      }
+    }
+    if ( theFace.Orientation() == TopAbs_REVERSED )
+      vecRef.Reverse();
+    c1->D1( p1, p, vec1 );
+    c2->D1( p2, p, vec2 );
+    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 );
+  }
+  catch (...)
+  {
+  }
+  return angle;
+}
+
 //================================================================================
 /*!
  * \brief Check if the first and last vertices of an edge are the same
@@ -2802,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)
 
@@ -3095,7 +3268,7 @@ namespace { // Structures used by FixQuadraticElements()
     chLink->SetFace( this );
     MSGBEG( *this );
 
-    // propagate from quadrangle to neighbour faces
+    // propagate from quadrangle to neighbour faces
     if ( link->MediumPos() >= pos ) {
       int nbLinkFaces = link->_faces.size();
       if ( nbLinkFaces == 4 || (/*nbLinkFaces < 4 && */link->OnBoundary())) {
@@ -3210,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();
@@ -3284,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();
@@ -3296,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
@@ -3307,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
@@ -3330,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;
     }
@@ -3420,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 );
     }
   }
   //================================================================================
@@ -3527,7 +3700,7 @@ namespace { // Structures used by FixQuadraticElements()
           if ( pInterLink == interLinks.end() ) continue; // not internal link
           interLink->Move( bndLink->_nodeMove );
           // treated internal links become new boundary ones
-          interLinks. erase( pInterLink );
+          interLinks.erase( pInterLink );
           newBndLinks->insert( interLink );
         }
       }
@@ -3784,7 +3957,7 @@ namespace { // Structures used by FixQuadraticElements()
       {
         // check if the EDGE needs checking
         const TopoDS_Edge& edge = TopoDS::Edge( edgeIt.Current() );
-        if ( BRep_Tool::Degenerated( edge ) )
+        if ( SMESH_Algo::isDegenerated( edge ) )
           continue;
         if ( theHelper.IsRealSeam( edge ) &&
              edge.Orientation() == TopAbs_REVERSED )
@@ -3847,7 +4020,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() < 6       || // check quadratic triangles only
                  !checkedFaces.insert( f ).second )
               continue;
 
@@ -4440,6 +4613,9 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
               gp_XY newUV   = applyIn2D( s, oldUV, gp_XY( move.X(),move.Y()), gp_XY_Added);
               gp_Pnt newPnt = s->Value( newUV.X(), newUV.Y());
               move = gp_Vec( XYZ((*link1)->_mediumNode), newPnt.Transformed(loc) );
+              if ( SMDS_FacePosition* nPos =
+                   dynamic_cast< SMDS_FacePosition* >((*link1)->_mediumNode->GetPosition()))
+                nPos->SetParameters( newUV.X(), newUV.Y() );
 #ifdef _DEBUG_
               if ( (XYZ((*link1)->node1()) - XYZ((*link1)->node2())).SquareModulus() <
                    move.SquareMagnitude())
@@ -4535,23 +4711,23 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
 
   // treat bi-quad triangles
   {
-    const SMDS_MeshNode* nodes[3]; // medium nodes
-    gp_XY uv[ 3 ]; // UV of medium nodes
+    vector< const SMDS_MeshNode* > nodes;
+    gp_XY uv[ 6 ];
     TIDSortedElemSet::iterator triIt = biQuadTris.begin();
     for ( ; triIt != biQuadTris.end(); ++triIt )
     {
       const SMDS_MeshElement* tria = *triIt;
-      // nodes
-      for ( int i = 3; i < 6; ++i )
-        nodes[ i-3 ] = tria->GetNode( i );
       // FACE
       const TopoDS_Shape& S = GetMeshDS()->IndexToShape( tria->getshapeId() );
       if ( S.IsNull() || S.ShapeType() != TopAbs_FACE ) continue;
       const TopoDS_Face& F = TopoDS::Face( S );
       Handle( Geom_Surface ) surf = BRep_Tool::Surface( F, loc );
       const double tol = BRep_Tool::Tolerance( F );
+
+      // nodes
+      nodes.assign( tria->begin_nodes(), tria->end_nodes() );
       // UV
-      for ( int i = 0; i < 3; ++i )
+      for ( int i = 0; i < 6; ++i )
       {
         uv[ i ] = GetNodeUV( F, nodes[i], nodes[(i+1)%3], &checkUV );
         // as this method is used after mesh generation, UV of nodes is not
@@ -4560,7 +4736,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
           CheckNodeUV( F, nodes[i], uv[ i ], 2*tol, /*force=*/true );
       }
       // move the central node
-      gp_XY uvCent = ( uv[0] + uv[1] + uv[2] ) / 3.;
+      gp_XY uvCent = GetCenterUV( uv[0], uv[1], uv[2], uv[3], uv[4], uv[5] );
       gp_Pnt p = surf->Value( uvCent.X(), uvCent.Y() ).Transformed( loc );
       GetMeshDS()->MoveNode( tria->GetNode(6), p.X(), p.Y(), p.Z() );
     }
@@ -4628,62 +4804,4 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError,
                              nCenterCoords.X(), nCenterCoords.Y(), nCenterCoords.Z());
     }
   }
-
-  // Issue 0020982
-  // Move the apex of pyramid together with the most curved link.
-  // TIDSortedNodeSet::iterator apexIt = apexOfPyramid.begin();
-  // for ( ; apexIt != apexOfPyramid.end(); ++apexIt )
-  // {
-  //   SMESH_TNodeXYZ apex = *apexIt;
-
-  //   gp_Vec maxMove( 0,0,0 );
-  //   double maxMoveSize2 = 0;
-
-  //   // shift of node index to get medium nodes between the base nodes
-  //   const int base2MediumShift = 5;
-
-  //   // find maximal movement of medium node
-  //   SMDS_ElemIteratorPtr volIt = apex._node->GetInverseElementIterator( SMDSAbs_Volume );
-  //   vector< const SMDS_MeshElement* > pyramids;
-  //   while ( volIt->more() )
-  //   {
-  //     const SMDS_MeshElement* pyram = volIt->next();
-  //     if ( pyram->GetEntityType() != SMDSEntity_Quad_Pyramid ) continue;
-  //     pyramids.push_back( pyram );
-
-  //     for ( int iBase = 0; iBase < apexIndex; ++iBase )
-  //     {
-  //       SMESH_TNodeXYZ medium = pyram->GetNode( iBase + base2MediumShift );
-  //       if ( medium._node->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE )
-  //       {
-  //         SMESH_TNodeXYZ n1 = pyram->GetNode( iBase );
-  //         SMESH_TNodeXYZ n2 = pyram->GetNode( ( iBase+1 ) % 4 );
-  //         gp_Pnt middle = 0.5 * ( n1 + n2 );
-  //         gp_Vec move( middle, medium );
-  //         double moveSize2 = move.SquareMagnitude();
-  //         if ( moveSize2 > maxMoveSize2 )
-  //           maxMove = move, maxMoveSize2 = moveSize2;
-  //       }
-  //     }
-  //   }
-
-  //   // move the apex
-  //   if ( maxMoveSize2 > 1e-20 )
-  //   {
-  //     apex += maxMove.XYZ();
-  //     GetMeshDS()->MoveNode( apex._node, apex.X(), apex.Y(), apex.Z());
-
-  //     // move medium nodes neighboring the apex to the middle
-  //     const int base2MediumShift_2 = 9;
-  //     for ( unsigned i = 0; i < pyramids.size(); ++i )
-  //       for ( int iBase = 0; iBase < apexIndex; ++iBase )
-  //       {
-  //         SMESH_TNodeXYZ         base = pyramids[i]->GetNode( iBase );
-  //         const SMDS_MeshNode* medium = pyramids[i]->GetNode( iBase + base2MediumShift_2 );
-  //         gp_XYZ middle = 0.5 * ( apex + base );
-  //         GetMeshDS()->MoveNode( medium, middle.X(), middle.Y(), middle.Z());
-  //       }
-  //   }
-  // }
 }
-