Salome HOME
IPAL0052448: Too thin viscous layers are constructed
[modules/smesh.git] / src / StdMeshers / StdMeshers_Quadrangle_2D.cxx
index e24cc375c85fdd1394fc0b583116eb6322b301bb..f720468f7bd3bfce81d71f3dee855186610ac29c 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
 #include "StdMeshers_QuadrangleParams.hxx"
 #include "StdMeshers_ViscousLayers2D.hxx"
 
+#include <BRepBndLib.hxx>
 #include <BRepClass_FaceClassifier.hxx>
 #include <BRep_Tool.hxx>
+#include <Bnd_Box.hxx>
 #include <GeomAPI_ProjectPointOnSurf.hxx>
 #include <Geom_Surface.hxx>
 #include <NCollection_DefineArray2.hxx>
@@ -88,6 +90,7 @@ StdMeshers_Quadrangle_2D::StdMeshers_Quadrangle_2D (int hypId, int studyId,
     myTrianglePreference(false),
     myTriaVertexID(-1),
     myNeedSmooth(false),
+    myCheckOri(false),
     myParams( NULL ),
     myQuadType(QUAD_STANDARD),
     myHelper( NULL )
@@ -213,15 +216,20 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh&         aMesh,
   const TopoDS_Face& F = TopoDS::Face(aShape);
   aMesh.GetSubMesh( F );
 
+  // do not initialize my fields before this as StdMeshers_ViscousLayers2D
+  // can call Compute() recursively
+  SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F );
+  if ( !proxyMesh )
+    return false;
+
+  myProxyMesh = proxyMesh;
+
   SMESH_MesherHelper helper (aMesh);
   myHelper = &helper;
 
-  myProxyMesh = StdMeshers_ViscousLayers2D::Compute( aMesh, F );
-  if ( !myProxyMesh )
-    return false;
-
   _quadraticMesh = myHelper->IsQuadraticSubMesh(aShape);
   myNeedSmooth = false;
+  myCheckOri   = false;
 
   FaceQuadStruct::Ptr quad = CheckNbEdges( aMesh, F, /*considerMesh=*/true );
   if (!quad)
@@ -288,6 +296,9 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh&         aMesh,
   if ( res == COMPUTE_OK && myNeedSmooth )
     smooth( quad );
 
+  if ( res == COMPUTE_OK )
+    res = check();
+
   return ( res == COMPUTE_OK );
 }
 
@@ -317,17 +328,25 @@ bool StdMeshers_Quadrangle_2D::computeTriangles(SMESH_Mesh&         aMesh,
   if ( !setNormalizedGrid( quad ))
     return false;
 
-  if ( quad->nbNodeOut( QUAD_BOTTOM_SIDE ))
-  {
-    splitQuad( quad, 0, 1 );
-  }
   if ( quad->nbNodeOut( QUAD_TOP_SIDE    ))
   {
     splitQuad( quad, 0, quad->jSize-2 );
   }
+  if ( quad->nbNodeOut( QUAD_BOTTOM_SIDE )) // this should not happen
+  {
+    splitQuad( quad, 0, 1 );
+  }
   FaceQuadStruct::Ptr newQuad = myQuadList.back();
   if ( quad != newQuad ) // split done
   {
+    { // update left side limit till where to make triangles
+      FaceQuadStruct::Ptr botQuad = // a bottom part
+        ( quad->side[ QUAD_LEFT_SIDE ].from == 0 ) ? quad : newQuad;
+      if ( botQuad->nbNodeOut( QUAD_LEFT_SIDE ) > 0 )
+        botQuad->side[ QUAD_LEFT_SIDE ].to += botQuad->nbNodeOut( QUAD_LEFT_SIDE );
+      else if ( botQuad->nbNodeOut( QUAD_RIGHT_SIDE ) > 0 )
+        botQuad->side[ QUAD_RIGHT_SIDE ].to += botQuad->nbNodeOut( QUAD_RIGHT_SIDE );
+    }
     // make quad be a greatest one
     if ( quad->side[ QUAD_LEFT_SIDE ].NbPoints() == 2 ||
          quad->side[ QUAD_RIGHT_SIDE ].NbPoints() == 2  )
@@ -340,12 +359,33 @@ bool StdMeshers_Quadrangle_2D::computeTriangles(SMESH_Mesh&         aMesh,
   {
     splitQuad( quad, quad->iSize-2, 0 );
   }
-  if ( quad->nbNodeOut( QUAD_LEFT_SIDE    ))
+  if ( quad->nbNodeOut( QUAD_LEFT_SIDE  ))
   {
     splitQuad( quad, 1, 0 );
+
+    if ( quad->nbNodeOut( QUAD_TOP_SIDE ))
+    {
+      newQuad = myQuadList.back();
+      if ( newQuad == quad ) // too narrow to split
+      {
+        // update left side limit till where to make triangles
+        quad->side[ QUAD_LEFT_SIDE ].to--;
+      }
+      else
+      {
+        FaceQuadStruct::Ptr leftQuad =
+          ( quad->side[ QUAD_BOTTOM_SIDE ].from == 0 ) ? quad : newQuad;
+        leftQuad->nbNodeOut( QUAD_TOP_SIDE ) = 0;
+      }
+    }
   }
 
-  return computeQuadDominant( aMesh, aFace );
+  if ( ! computeQuadDominant( aMesh, aFace ))
+    return false;
+
+  // try to fix zero-area triangles near straight-angle corners
+
+  return true;
 }
 
 //================================================================================
@@ -455,7 +495,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh&         aMesh,
   int nbright = (int) uv_e1.size();
   int nbleft  = (int) uv_e3.size();
 
-  if (quad->nbNodeOut(0) && nbvertic == 2)
+  if (quad->nbNodeOut(0) && nbvertic == 2) // this should not occure
   {
     // Down edge is out
     // 
@@ -572,17 +612,59 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh&         aMesh,
       iup = nbhoriz - 1;
 
       int stop = 0;
-      // if left edge is out, we will stop at a second node
-      //if (quad->nbNodeOut(3)) stop++;
-      if ( quad->nbNodeOut( QUAD_RIGHT_SIDE ))
-        quad->UVPt( nbhoriz-1, 0 ).node = uv_e1[ nbright-2 ].node;
-      if ( quad->nbNodeOut( QUAD_LEFT_SIDE ))
-        quad->UVPt( 0, 0 ).node = uv_e3[ nbleft-2 ].node;
+      if ( quad->side[3].grid->Edge(0).IsNull() ) // left side is simulated one
+      {
+        // quad divided at I but not at J, as nbvertic==nbright==2
+        stop++; // we stop at a second node
+      }
+      else
+      {
+        if ( quad->nbNodeOut( QUAD_RIGHT_SIDE ))
+          quad->UVPt( nbhoriz-1, 0 ).node = uv_e1[ nbright-2 ].node;
+        if ( quad->nbNodeOut( QUAD_LEFT_SIDE ))
+          quad->UVPt( 0, 0 ).node = uv_e3[ nbleft-2 ].node;
 
+        if ( nbright > 2 ) // there was a split at J
+          quad->nbNodeOut( QUAD_LEFT_SIDE ) = 0;
+      }
+      const SMDS_MeshNode *a, *b, *c, *d;
+      i = nbup - 1;
+      // avoid creating zero-area triangles near a straight-angle corner
+      {
+        a = uv_e2[i].node;
+        b = uv_e2[i-1].node;
+        c = uv_e1[nbright-2].node;
+        SMESH_TNodeXYZ pa( a ), pb( b ), pc( c );
+        double area = 0.5 * (( pb - pa ) ^ ( pc - pa )).Modulus();
+        if ( Abs( area ) < 1e-20 )
+        {
+          --g;
+          d = quad->UVPt( g, nbvertic-2 ).node;
+          if ( myTrianglePreference )
+          {
+            if ( SMDS_MeshFace* face = myHelper->AddFace(a, d, c))
+              meshDS->SetMeshElementOnShape(face, geomFaceID);
+          }
+          else
+          {
+            if ( SMDS_MeshFace* face = myHelper->AddFace(a, b, d, c))
+            {
+              meshDS->SetMeshElementOnShape(face, geomFaceID);
+              SMESH_ComputeErrorPtr& err = aMesh.GetSubMesh( aFace )->GetComputeError();
+              if ( !err || err->IsOK() || err->myName < COMPERR_WARNING )
+              {
+                err.reset( new SMESH_ComputeError( COMPERR_WARNING,
+                                                   "Bad quality quad created"));
+                err->myBadElements.push_back( face );
+              }
+            }
+            --i;
+          }
+        }
+      }
       // for each node of the up edge find nearest node
       // in the first row of the regular grid and link them
-      for (i = nbup - 1; i > stop; i--) {
-        const SMDS_MeshNode *a, *b, *c, *d;
+      for ( ; i > stop; i--) {
         a = uv_e2[i].node;
         b = uv_e2[i - 1].node;
         gp_Pnt pb (b->X(), b->Y(), b->Z());
@@ -649,7 +731,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh&         aMesh,
   }
 
   // right or left boundary quadrangles
-  if (quad->nbNodeOut( QUAD_RIGHT_SIDE ) && nbhoriz == 2)
+  if (quad->nbNodeOut( QUAD_RIGHT_SIDE ) && nbhoriz == 2) // this should not occure
   {
     int g = 0; // last processed node in the grid
     int stop = nbright - 1;
@@ -725,11 +807,44 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh&         aMesh,
 //      MESSAGE("left edge is out");
       int g = nbvertic - 1; // last processed node in the grid
       int stop = 0;
-      i = nbleft - 1;
-      if (quad->side[3].from != stop ) stop++;
-      if (quad->side[3].to   != i    ) i--;
-      for (; i > stop; i--) {
-        const SMDS_MeshNode *a, *b, *c, *d;
+      i = quad->side[ QUAD_LEFT_SIDE ].to-1; // nbleft - 1;
+
+      const SMDS_MeshNode *a, *b, *c, *d;
+      // avoid creating zero-area triangles near a straight-angle corner
+      {
+        a = uv_e3[i].node;
+        b = uv_e3[i-1].node;
+        c = quad->UVPt( 1, g ).node;
+        SMESH_TNodeXYZ pa( a ), pb( b ), pc( c );
+        double area = 0.5 * (( pb - pa ) ^ ( pc - pa )).Modulus();
+        if ( Abs( area ) < 1e-20 )
+        {
+          --g;
+          d = quad->UVPt( 1, g ).node;
+          if ( myTrianglePreference )
+          {
+            if ( SMDS_MeshFace* face = myHelper->AddFace(a, d, c))
+              meshDS->SetMeshElementOnShape(face, geomFaceID);
+          }
+          else
+          {
+            if ( SMDS_MeshFace* face = myHelper->AddFace(a, b, d, c))
+            {
+              meshDS->SetMeshElementOnShape(face, geomFaceID);
+              SMESH_ComputeErrorPtr& err = aMesh.GetSubMesh( aFace )->GetComputeError();
+              if ( !err || err->IsOK() || err->myName < COMPERR_WARNING )
+              {
+                err.reset( new SMESH_ComputeError( COMPERR_WARNING,
+                                                   "Bad quality quad created"));
+                err->myBadElements.push_back( face );
+              }
+            }
+            --i;
+          }
+        }
+      }
+      for (; i > stop; i--) // loop on nodes on the left side
+      {
         a = uv_e3[i].node;
         b = uv_e3[i - 1].node;
         gp_Pnt pb (b->X(), b->Y(), b->Z());
@@ -739,12 +854,13 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh&         aMesh,
         if (i == stop + 1) { // down bondary reached
           c = quad->uv_grid[nbhoriz*jlow + 1].node;
           near = jlow;
-        } else {
+        }
+        else {
           double mind = RealLast();
           for (int k = g; k >= jlow; k--) {
             const SMDS_MeshNode *nk;
             if (k > jup)
-              nk = uv_e2[1].node;
+              nk = quad->uv_grid[nbhoriz*jup + 1].node; //uv_e2[1].node;
             else
               nk = quad->uv_grid[nbhoriz*k + 1].node;
             gp_Pnt pnk (nk->X(), nk->Y(), nk->Z());
@@ -765,11 +881,10 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh&         aMesh,
         }
         else { // make quadrangle
           if (near + 1 > jup)
-            d = uv_e2[1].node;
+            d = quad->uv_grid[nbhoriz*jup + 1].node; //uv_e2[1].node;
           else
             d = quad->uv_grid[nbhoriz*(near + 1) + 1].node;
-          //SMDS_MeshFace* face = meshDS->AddFace(a, b, c, d);
-          if (!myTrianglePreference){
+          if (!myTrianglePreference) {
             SMDS_MeshFace* face = myHelper->AddFace(a, b, c, d);
             if (face) meshDS->SetMeshElementOnShape(face, geomFaceID);
           }
@@ -781,7 +896,7 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh&         aMesh,
             for (int k = near + 1; k < g; k++) {
               c = quad->uv_grid[nbhoriz*k + 1].node;
               if (k + 1 > jup)
-                d = uv_e2[1].node;
+                d = quad->uv_grid[nbhoriz*jup + 1].node; //uv_e2[1].node;
               else
                 d = quad->uv_grid[nbhoriz*(k + 1) + 1].node;
               SMDS_MeshFace* face = myHelper->AddFace(a, c, d);
@@ -892,6 +1007,38 @@ bool StdMeshers_Quadrangle_2D::Evaluate(SMESH_Mesh&         aMesh,
   return true;
 }
 
+//================================================================================
+/*!
+ * \brief Return true if the algorithm can mesh this 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
+ */
+//================================================================================
+
+bool StdMeshers_Quadrangle_2D::IsApplicable( const TopoDS_Shape & aShape, bool toCheckAll )
+{
+  int nbFoundFaces = 0;
+  for (TopExp_Explorer exp( aShape, TopAbs_FACE ); exp.More(); exp.Next(), ++nbFoundFaces )
+  {
+    const TopoDS_Shape& aFace = exp.Current();
+    int nbWire = SMESH_MesherHelper::Count( aFace, TopAbs_WIRE, false );
+    if ( nbWire != 1 ) {
+      if ( toCheckAll ) return false;
+      continue;
+    }
+
+    int nbNoDegenEdges = 0;
+    TopExp_Explorer eExp( aFace, TopAbs_EDGE );
+    for ( ; eExp.More() && nbNoDegenEdges < 3; eExp.Next() ) {
+      if ( !SMESH_Algo::isDegenerated( TopoDS::Edge( eExp.Current() )))
+        ++nbNoDegenEdges;
+    }
+    if ( toCheckAll  && nbNoDegenEdges <  3 ) return false;
+    if ( !toCheckAll && nbNoDegenEdges >= 3 ) return true;
+  }
+  return ( toCheckAll && nbFoundFaces != 0 );
+}
 
 //================================================================================
 /*!
@@ -931,7 +1078,7 @@ FaceQuadStruct::Ptr StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh &
   if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD );
   const bool ignoreMediumNodes = _quadraticMesh;
 
-  // verify 1 wire only, with 4 edges
+  // verify 1 wire only
   list< TopoDS_Edge > edges;
   list< int > nbEdgesInWire;
   int nbWire = SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire);
@@ -987,45 +1134,40 @@ FaceQuadStruct::Ptr StdMeshers_Quadrangle_2D::CheckNbEdges(SMESH_Mesh &
     {
       list< TopoDS_Edge > sideEdges;
       TopoDS_Vertex nextSideV = corners[( iSide + 1 - nbUsedDegen ) % corners.size() ];
-      while ( edgeIt != edges.end() &&
-              !nextSideV.IsSame( myHelper->IthVertex( 0, *edgeIt )))
+      bool nextSideVReached = false;
+      do
       {
-        if ( SMESH_Algo::isDegenerated( *edgeIt ) )
+        const TopoDS_Edge& edge = *edgeIt;
+        nextSideVReached = nextSideV.IsSame( myHelper->IthVertex( 1, edge ));
+        if ( SMESH_Algo::isDegenerated( edge ))
         {
-          if ( myNeedSmooth )
-          {
-            ++edgeIt; // no side on the degenerated EDGE
-          }
-          else
+          if ( !myNeedSmooth ) // need to make a side on a degen edge
           {
             if ( sideEdges.empty() )
             {
+              sideEdges.push_back( edge );
               ++nbUsedDegen;
-              sideEdges.push_back( *edgeIt++ ); // a degenerated side
-              break;
+              nextSideVReached = true;
             }
             else
             {
-              break; // do not append a degenerated EDGE to a regular side
+              break;
             }
           }
         }
         else
         {
-          sideEdges.push_back( *edgeIt++ );
+          sideEdges.push_back( edge );
         }
+        ++edgeIt;
       }
+      while ( edgeIt != edges.end() && !nextSideVReached );
+
       if ( !sideEdges.empty() )
       {
-        quad->side.push_back( StdMeshers_FaceSide::New( F, sideEdges, &aMesh, iSide < QUAD_TOP_SIDE,
-                                                        ignoreMediumNodes, myProxyMesh ));
-        ++iSide;
-      }
-      else if ( !SMESH_Algo::isDegenerated( *edgeIt ) && // closed EDGE
-                myHelper->IthVertex( 0, *edgeIt ).IsSame( myHelper->IthVertex( 1, *edgeIt )))
-      {
-        quad->side.push_back( StdMeshers_FaceSide::New( F, *edgeIt++, &aMesh, iSide < QUAD_TOP_SIDE,
-                                                        ignoreMediumNodes, myProxyMesh));
+        quad->side.push_back
+          ( StdMeshers_FaceSide::New( F, sideEdges, &aMesh, iSide < QUAD_TOP_SIDE,
+                                      ignoreMediumNodes, myProxyMesh ));
         ++iSide;
       }
       if ( quad->side.size() == 4 )
@@ -3639,6 +3781,18 @@ namespace // data for smoothing
     double d = v1 ^ v2;
     return d > 1e-100;
   }
+  //================================================================================
+  /*!
+   * \brief Returns area of a triangle
+   */
+  //================================================================================
+
+  double getArea( const gp_UV uv1, const gp_UV uv2, const gp_UV uv3 )
+  {
+    gp_XY v1 = uv1 - uv2, v2 = uv3 - uv2;
+    double a = v2 ^ v1;
+    return a;
+  }
 }
 
 //================================================================================
@@ -3689,7 +3843,7 @@ void StdMeshers_Quadrangle_2D::updateDegenUV(FaceQuadStruct::Ptr quad)
       uv1.v = uv2.v = 0.5 * ( uv1.v + uv2.v );
     }
 
-  else if ( quad->side.size() == 4 && myQuadType == QUAD_STANDARD)
+  else if ( quad->side.size() == 4 /*&& myQuadType == QUAD_STANDARD*/)
 
     // Set number of nodes on a degenerated side to be same as on an opposite side
     // ----------------------------------------------------------------------------
@@ -3726,6 +3880,8 @@ void StdMeshers_Quadrangle_2D::smooth (FaceQuadStruct::Ptr quad)
 
   // Get nodes to smooth
 
+  // TODO: do not smooth fixed nodes
+
   typedef map< const SMDS_MeshNode*, TSmoothNode, TIDCompare > TNo2SmooNoMap;
   TNo2SmooNoMap smooNoMap;
 
@@ -3879,6 +4035,145 @@ void StdMeshers_Quadrangle_2D::smooth (FaceQuadStruct::Ptr quad)
   }
 }
 
+//================================================================================
+/*!
+ * \brief Checks validity of generated faces
+ */
+//================================================================================
+
+bool StdMeshers_Quadrangle_2D::check()
+{
+  const bool isOK = true;
+  if ( !myCheckOri || myQuadList.empty() || !myQuadList.front() || !myHelper )
+    return isOK;
+
+  TopoDS_Face      geomFace = TopoDS::Face( myHelper->GetSubShape() );
+  SMESHDS_Mesh*    meshDS   = myHelper->GetMeshDS();
+  SMESHDS_SubMesh* fSubMesh = meshDS->MeshElements( geomFace );
+  bool toCheckUV;
+  if ( geomFace.Orientation() >= TopAbs_INTERNAL ) geomFace.Orientation( TopAbs_FORWARD );
+
+  // Get a reference orientation sign
+
+  double okSign;
+  {
+    TError err;
+    TSideVector wireVec =
+      StdMeshers_FaceSide::GetFaceWires( geomFace, *myHelper->GetMesh(), true, err );
+    StdMeshers_FaceSidePtr wire = wireVec[0];
+
+    // find a right angle VERTEX
+    int iVertex;
+    double maxAngle = -1e100;
+    for ( int i = 0; i < wire->NbEdges(); ++i )
+    {
+      int iPrev = myHelper->WrapIndex( i-1, wire->NbEdges() );
+      const TopoDS_Edge& e1 = wire->Edge( iPrev );
+      const TopoDS_Edge& e2 = wire->Edge( i );
+      double angle = myHelper->GetAngle( e1, e2, geomFace, wire->FirstVertex( i ));
+      if (( maxAngle < angle ) &&
+          ( 5.* M_PI/180 < angle && angle < 175.* M_PI/180  ))
+      {
+        maxAngle = angle;
+        iVertex = i;
+      }
+    }
+    if ( maxAngle < -2*M_PI ) return isOK;
+
+    // get a sign of 2D area of a corner face
+
+    int iPrev = myHelper->WrapIndex( iVertex-1, wire->NbEdges() );
+    const TopoDS_Edge& e1 = wire->Edge( iPrev );
+    const TopoDS_Edge& e2 = wire->Edge( iVertex );
+
+    gp_Vec2d v1, v2; gp_Pnt2d p;
+    double u[2];
+    {
+      bool rev = ( e1.Orientation() == TopAbs_REVERSED );
+      Handle(Geom2d_Curve) c = BRep_Tool::CurveOnSurface( e1, geomFace, u[0], u[1] );
+      c->D1( u[ !rev ], p, v1 );
+      if ( !rev )
+        v1.Reverse();
+    }
+    {
+      bool rev = ( e2.Orientation() == TopAbs_REVERSED );
+      Handle(Geom2d_Curve) c = BRep_Tool::CurveOnSurface( e2, geomFace, u[0], u[1] );
+      c->D1( u[ rev ], p, v2 );
+      if ( rev )
+        v2.Reverse();
+    }
+
+    okSign = v2 ^ v1;
+
+    if ( maxAngle < 0 )
+      okSign *= -1;
+  }
+
+  // Look for incorrectly oriented faces
+
+  std::list<const SMDS_MeshElement*> badFaces;
+
+  const SMDS_MeshNode* nn [ 8 ]; // 8 is just for safety
+  gp_UV                uv [ 8 ];
+  SMDS_ElemIteratorPtr fIt = fSubMesh->GetElements();
+  while ( fIt->more() ) // loop on faces bound to a FACE
+  {
+    const SMDS_MeshElement* f = fIt->next();
+
+    const int nbN = f->NbCornerNodes();
+    for ( int i = 0; i < nbN; ++i )
+      nn[ i ] = f->GetNode( i );
+
+    const SMDS_MeshNode* nInFace = 0;
+    if ( myHelper->HasSeam() )
+      for ( int i = 0; i < nbN && !nInFace; ++i )
+        if ( !myHelper->IsSeamShape( nn[i]->getshapeId() ))
+          nInFace = nn[i];
+
+    for ( int i = 0; i < nbN; ++i )
+      uv[ i ] = myHelper->GetNodeUV( geomFace, nn[i], nInFace, &toCheckUV );
+
+    switch ( nbN ) {
+    case 4:
+    {
+      double sign1 = getArea( uv[0], uv[1], uv[2] );
+      double sign2 = getArea( uv[0], uv[2], uv[3] );
+      if ( sign1 * sign2 < 0 )
+      {
+        sign2 = getArea( uv[1], uv[2], uv[3] );
+        sign1 = getArea( uv[1], uv[3], uv[0] );
+        if ( sign1 * sign2 < 0 )
+          continue; // this should not happen
+      }
+      if ( sign1 * okSign < 0 )
+        badFaces.push_back ( f );
+      break;
+    }
+    case 3:
+    {
+      double sign = getArea( uv[0], uv[1], uv[2] );
+      if ( sign * okSign < 0 )
+        badFaces.push_back ( f );
+      break;
+    }
+    default:;
+    }
+  }
+
+  if ( !badFaces.empty() )
+  {
+    SMESH_subMesh* fSM = myHelper->GetMesh()->GetSubMesh( geomFace );
+    SMESH_ComputeErrorPtr& err = fSM->GetComputeError();
+    err.reset ( new SMESH_ComputeError( COMPERR_ALGO_FAILED,
+                                        "Inverted elements generated"));
+    err->myBadElements.swap( badFaces );
+
+    return !isOK;
+  }
+
+  return isOK;
+}
+
 /*//================================================================================
 /*!
  * \brief Finds vertices at the most sharp face corners
@@ -3929,7 +4224,7 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face&          theFace,
     TopoDS_Vertex v = helper.IthVertex( 0, *edge );
     if ( !theConsiderMesh || SMESH_Algo::VertexNode( v, helper.GetMeshDS() ))
     {
-      double angle = SMESH_MesherHelper::GetAngle( prevE, *edge, theFace );
+      double angle = SMESH_MesherHelper::GetAngle( prevE, *edge, theFace, v );
       vertexByAngle.insert( make_pair( angle, v ));
       angleByVertex.Bind( v, angle );
     }
@@ -3941,7 +4236,8 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face&          theFace,
   TopoDS_Shape triaVertex = helper.GetMeshDS()->IndexToShape( myTriaVertexID );
   if ( !triaVertex.IsNull() &&
        triaVertex.ShapeType() == TopAbs_VERTEX &&
-       helper.IsSubShape( triaVertex, theFace ))
+       helper.IsSubShape( triaVertex, theFace ) &&
+       ( vertexByAngle.size() != 4 || vertexByAngle.begin()->first < 5 * M_PI/180. ))
     nbCorners = 3;
   else
     triaVertex.Nullify();
@@ -3980,7 +4276,7 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face&          theFace,
   if ( nbCorners == 3 )
     vMap.Add( triaVertex );
   multimap<double, TopoDS_Vertex>::reverse_iterator a2v = vertexByAngle.rbegin();
-  for ( ; a2v != vertexByAngle.rend() && vMap.Extent() < nbCorners; ++a2v )
+  for ( int iC = 0; a2v != vertexByAngle.rend() && iC < nbCorners; ++a2v, ++iC )
     vMap.Add( (*a2v).second );
 
   // check if there are possible variations in choosing corners
@@ -3992,6 +4288,9 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face&          theFace,
     isThereVariants  = ( lostAngle * 1.1 >= lastAngle );
   }
 
+  myCheckOri = ( vertexByAngle.size() > nbCorners ||
+                 vertexByAngle.begin()->first < 5.* M_PI/180 );
+
   // make theWire begin from a corner vertex or triaVertex
   if ( nbCorners == 3 )
     while ( !triaVertex.IsSame( ( helper.IthVertex( 0, theWire.front() ))) ||
@@ -4194,6 +4493,9 @@ bool StdMeshers_Quadrangle_2D::getEnforcedUV()
   surf->Bounds( u1,u2,v1,v2 );
   GeomAPI_ProjectPointOnSurf project;
   project.Init(surf, u1,u2, v1,v2, tol );
+  Bnd_Box bbox;
+  BRepBndLib::Add( face, bbox );
+  double farTol = 0.01 * sqrt( bbox.SquareExtent() );
 
   for ( size_t iP = 0; iP < points.size(); ++iP )
   {
@@ -4206,7 +4508,7 @@ bool StdMeshers_Quadrangle_2D::getEnforcedUV()
            << points[ iP ].X() << ", "<< points[ iP ].Y() << ", "<< points[ iP ].Z() << " )");
       continue;
     }
-    if ( project.LowerDistance() > tol*1000 )
+    if ( project.LowerDistance() > farTol )
     {
       if ( isStrictCheck && iP < nbPoints )
         return error
@@ -4499,7 +4801,7 @@ int StdMeshers_Quadrangle_2D::splitQuad(FaceQuadStruct::Ptr quad, int I, int J)
   myQuadList.push_back( FaceQuadStruct::Ptr( newQuad ));
 
   vector<UVPtStruct> points;
-  if ( I > 0 )
+  if ( I > 0 && I <= quad->iSize-2 )
   {
     points.reserve( quad->jSize );
     for ( int jP = 0; jP < quad->jSize; ++jP )
@@ -4538,7 +4840,7 @@ int StdMeshers_Quadrangle_2D::splitQuad(FaceQuadStruct::Ptr quad, int I, int J)
 
     return QUAD_LEFT_SIDE;
   }
-  else if ( J > 0 ) //// split horizontally, a new quad is below an old one
+  else if ( J > 0  && J <= quad->jSize-2 ) //// split horizontally, a new quad is below an old one
   {
     points.reserve( quad->iSize );
     for ( int iP = 0; iP < quad->iSize; ++iP )
@@ -4577,6 +4879,9 @@ int StdMeshers_Quadrangle_2D::splitQuad(FaceQuadStruct::Ptr quad, int I, int J)
 
     return QUAD_TOP_SIDE;
   }
+
+  myQuadList.pop_back();
+  return -1;
 }
 
 //================================================================================
@@ -5033,6 +5338,7 @@ FaceQuadStruct::Side& FaceQuadStruct::Side::operator=(const Side& otherSide)
         oSide->contacts[iOC].other_side = this;
       }
   }
+  return *this;
 }
 
 //================================================================================
@@ -5079,6 +5385,7 @@ bool FaceQuadStruct::Side::Reverse(bool keepGrid)
       grid->Reverse();
     }
   }
+  return (bool)grid;
 }
 
 //================================================================================