Salome HOME
52568: Quadrangle (Mapping) differently meshes equal L-shaped faces.
[modules/smesh.git] / src / StdMeshers / StdMeshers_Quadrangle_2D.cxx
index 016a5098db2c27e8f35e096d8077e519502b4483..4d1321497a6d08eb6ad1fe25b8c4b43227dc3647 100644 (file)
@@ -67,9 +67,7 @@
 #ifndef StdMeshers_Array2OfNode_HeaderFile
 #define StdMeshers_Array2OfNode_HeaderFile
 typedef const SMDS_MeshNode* SMDS_MeshNodePtr;
-DEFINE_BASECOLLECTION (StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr)
-DEFINE_ARRAY2(StdMeshers_Array2OfNode,
-              StdMeshers_BaseCollectionNodePtr, SMDS_MeshNodePtr)
+typedef NCollection_Array2<SMDS_MeshNodePtr> StdMeshers_Array2OfNode;
 #endif
 
 using namespace std;
@@ -201,7 +199,9 @@ bool StdMeshers_Quadrangle_2D::CheckHypothesis
     }
   }
 
-  return isOk;
+  error( StdMeshers_ViscousLayers2D::CheckHypothesis( aMesh, aShape, aStatus ));
+
+  return aStatus == HYP_OK;
 }
 
 //=============================================================================
@@ -277,7 +277,7 @@ bool StdMeshers_Quadrangle_2D::Compute (SMESH_Mesh&         aMesh,
                "two opposite sides should have same number of segments, "
                "but actual number of segments is different on all sides. "
                "'Standard' transion has been used.");
-      else
+      else if ( ! ( n1 == n3 && n2 == n4 ))
         error( COMPERR_WARNING,
                "To use 'Reduced' transition, "
                "two opposite sides should have an even difference in number of segments. "
@@ -339,7 +339,7 @@ bool StdMeshers_Quadrangle_2D::computeTriangles(SMESH_Mesh&         aMesh,
   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 )
@@ -359,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;
 }
 
 //================================================================================
@@ -591,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());
@@ -745,8 +808,43 @@ bool StdMeshers_Quadrangle_2D::computeQuadDominant(SMESH_Mesh&         aMesh,
       int g = nbvertic - 1; // last processed node in the grid
       int stop = 0;
       i = quad->side[ QUAD_LEFT_SIDE ].to-1; // nbleft - 1;
-      for (; i > stop; i--) {
-        const SMDS_MeshNode *a, *b, *c, *d;
+
+      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());
@@ -756,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());
@@ -782,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);
           }
@@ -798,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);
@@ -1036,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 )
@@ -3750,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
     // ----------------------------------------------------------------------------
@@ -3977,8 +4070,9 @@ bool StdMeshers_Quadrangle_2D::check()
       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 );
-      if ( maxAngle < angle && angle < 0.9 * M_PI )
+      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;
@@ -4130,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 );
     }
@@ -4182,20 +4276,21 @@ 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
-  bool isThereVariants = false;
+  bool haveVariants = false;
   if ( vertexByAngle.size() > nbCorners )
   {
     double lostAngle = a2v->first;
     double lastAngle = ( --a2v, a2v->first );
-    isThereVariants  = ( lostAngle * 1.1 >= lastAngle );
+    haveVariants  = ( lostAngle * 1.1 >= lastAngle );
   }
 
+  const double angleTol = 5.* M_PI/180;
   myCheckOri = ( vertexByAngle.size() > nbCorners ||
-                 vertexByAngle.begin()->first < 5.* M_PI/180 );
+                 vertexByAngle.begin()->first < angleTol );
 
   // make theWire begin from a corner vertex or triaVertex
   if ( nbCorners == 3 )
@@ -4212,9 +4307,10 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face&          theFace,
   vector< double >      angles;
   vector< TopoDS_Edge > edgeVec;
   vector< int >         cornerInd, nbSeg;
-  angles.reserve( vertexByAngle.size() );
+  int nbSegTot = 0;
+  angles .reserve( vertexByAngle.size() );
   edgeVec.reserve( vertexByAngle.size() );
-  nbSeg.reserve( vertexByAngle.size() );
+  nbSeg  .reserve( vertexByAngle.size() );
   cornerInd.reserve( nbCorners );
   for ( edge = theWire.begin(); edge != theWire.end(); ++edge )
   {
@@ -4227,105 +4323,219 @@ int StdMeshers_Quadrangle_2D::getCorners(const TopoDS_Face&          theFace,
       theVertices.push_back( v );
       cornerInd.push_back( angles.size() );
     }
-    angles.push_back( angleByVertex.IsBound( v ) ? angleByVertex( v ) : -M_PI );
+    angles .push_back( angleByVertex.IsBound( v ) ? angleByVertex( v ) : -M_PI );
     edgeVec.push_back( *edge );
-    if ( theConsiderMesh && isThereVariants )
+    if ( theConsiderMesh && haveVariants )
     {
       if ( SMESHDS_SubMesh* sm = helper.GetMeshDS()->MeshElements( *edge ))
         nbSeg.push_back( sm->NbNodes() + 1 );
       else
         nbSeg.push_back( 0 );
+      nbSegTot += nbSeg.back();
     }
   }
 
-  // refine the result vector - make sides elual by length if
+  // refine the result vector - make sides equal by length if
   // there are several equal angles
-  if ( isThereVariants )
+  if ( haveVariants )
   {
     if ( nbCorners == 3 )
       angles[0] = 2 * M_PI; // not to move the base triangle VERTEX
 
-    set< int > refinedCorners;
+    // here we refer to VERTEX'es and EDGEs by indices in angles and edgeVec vectors
+    typedef int TGeoIndex;
+
+    // for each vertex find a vertex till which there are nbSegHalf segments
+    const int nbSegHalf = ( nbSegTot % 2 || nbCorners == 3 ) ? 0 : nbSegTot / 2;
+    vector< TGeoIndex > halfDivider( angles.size(), -1 );
+    int nbHalfDividers = 0;
+    if ( nbSegHalf )
+    {
+      // get min angle of corners
+      double minAngle = 10.;
+      for ( size_t iC = 0; iC < cornerInd.size(); ++iC )
+        minAngle = Min( minAngle, angles[ cornerInd[ iC ]]);
+
+      // find halfDivider's
+      for ( TGeoIndex iV1 = 0; iV1 < TGeoIndex( angles.size() ); ++iV1 )
+      {
+        int nbSegs = 0;
+        TGeoIndex iV2 = iV1;
+        do {
+          nbSegs += nbSeg[ iV2 ];
+          iV2 = helper.WrapIndex( iV2 + 1, nbSeg.size() );
+        } while ( nbSegs < nbSegHalf );
+
+        if ( nbSegs == nbSegHalf &&
+             angles[ iV1 ] + angleTol >= minAngle &&
+             angles[ iV2 ] + angleTol >= minAngle )
+        {
+          halfDivider[ iV1 ] = iV2;
+          ++nbHalfDividers;
+        }
+      }
+    }
+
+    set< TGeoIndex > refinedCorners, treatedCorners;
     for ( size_t iC = 0; iC < cornerInd.size(); ++iC )
     {
-      int iV = cornerInd[iC];
-      if ( !refinedCorners.insert( iV ).second )
+      TGeoIndex iV = cornerInd[iC];
+      if ( !treatedCorners.insert( iV ).second )
         continue;
-      list< int > equalVertices;
-      equalVertices.push_back( iV );
+      list< TGeoIndex > equVerts; // inds of vertices that can become corners
+      equVerts.push_back( iV );
       int nbC[2] = { 0, 0 };
       // find equal angles backward and forward from the iV-th corner vertex
       for ( int isFwd = 0; isFwd < 2; ++isFwd )
       {
-        int     dV = isFwd ? +1 : -1;
-        int iCNext = helper.WrapIndex( iC + dV, cornerInd.size() );
-        int iVNext = helper.WrapIndex( iV + dV, angles.size() );
+        int           dV = isFwd ? +1 : -1;
+        int       iCNext = helper.WrapIndex( iC + dV, cornerInd.size() );
+        TGeoIndex iVNext = helper.WrapIndex( iV + dV, angles.size() );
         while ( iVNext != iV )
         {
-          bool equal = Abs( angles[iV] - angles[iVNext] ) < 0.1 * angles[iV];
+          bool equal = Abs( angles[iV] - angles[iVNext] ) < angleTol;
           if ( equal )
-            equalVertices.insert( isFwd ? equalVertices.end() : equalVertices.begin(), iVNext );
+            equVerts.insert( isFwd ? equVerts.end() : equVerts.begin(), iVNext );
           if ( iVNext == cornerInd[ iCNext ])
           {
             if ( !equal )
+            {
+              if ( angles[iV] < angles[iVNext] )
+                refinedCorners.insert( iVNext );
               break;
+            }
             nbC[ isFwd ]++;
-            refinedCorners.insert( cornerInd[ iCNext ] );
+            treatedCorners.insert( cornerInd[ iCNext ] );
             iCNext = helper.WrapIndex( iCNext + dV, cornerInd.size() );
           }
           iVNext = helper.WrapIndex( iVNext + dV, angles.size() );
         }
+        if ( iVNext == iV )
+          break; // all angles equal
       }
+
+      const bool allCornersSame = ( nbC[0] == 3 );
+      if ( allCornersSame && nbHalfDividers > 0 )
+      {
+        // select two halfDivider's as corners
+        TGeoIndex hd1, hd2 = -1;
+        int iC2;
+        for ( iC2 = 0; iC2 < cornerInd.size() && hd2 < 0; ++iC2 )
+        {
+          hd1 = cornerInd[ iC2 ];
+          hd2 = halfDivider[ hd1 ];
+          if ( std::find( equVerts.begin(), equVerts.end(), hd2 ) == equVerts.end() )
+            hd2 = -1; // hd2-th vertex can't become a corner
+          else
+            break;
+        }
+        if ( hd2 >= 0 )
+        {
+          angles[ hd1 ] = 2 * M_PI; // make hd1-th vertex no more "equal"
+          angles[ hd2 ] = 2 * M_PI;
+          refinedCorners.insert( hd1 );
+          refinedCorners.insert( hd2 );
+          treatedCorners = refinedCorners;
+          // update cornerInd
+          equVerts.push_front( equVerts.back() );
+          equVerts.push_back( equVerts.front() );
+          list< TGeoIndex >::iterator hdPos =
+            std::find( equVerts.begin(), equVerts.end(), hd2 );
+          if ( hdPos == equVerts.end() ) break;
+          cornerInd[ helper.WrapIndex( iC2 + 0, cornerInd.size()) ] = hd1;
+          cornerInd[ helper.WrapIndex( iC2 + 1, cornerInd.size()) ] = *( --hdPos );
+          cornerInd[ helper.WrapIndex( iC2 + 2, cornerInd.size()) ] = hd2;
+          cornerInd[ helper.WrapIndex( iC2 + 3, cornerInd.size()) ] = *( ++hdPos, ++hdPos );
+
+          theVertices[ 0 ] = helper.IthVertex( 0, edgeVec[ cornerInd[0] ]);
+          theVertices[ 1 ] = helper.IthVertex( 0, edgeVec[ cornerInd[1] ]);
+          theVertices[ 2 ] = helper.IthVertex( 0, edgeVec[ cornerInd[2] ]);
+          theVertices[ 3 ] = helper.IthVertex( 0, edgeVec[ cornerInd[3] ]);
+          iC = -1;
+          continue;
+        }
+      }
+
       // move corners to make sides equal by length
-      int nbEqualV  = equalVertices.size();
+      int nbEqualV  = equVerts.size();
       int nbExcessV = nbEqualV - ( 1 + nbC[0] + nbC[1] );
-      if ( nbExcessV > 0 )
+      if ( nbExcessV > 0 ) // there is nbExcessV vertices that can become corners
       {
-        // calculate normalized length of each side enclosed between neighbor equalVertices
-        vector< double > curLengths;
+        // calculate normalized length of each "side" enclosed between neighbor equVerts
+        vector< double > accuLength;
         double totalLen = 0;
-        vector< int > evVec( equalVertices.begin(), equalVertices.end() );
-        int   iEV = 0;
-        int    iE = cornerInd[ helper.WrapIndex( iC - nbC[0] - 1, cornerInd.size() )];
-        int iEEnd = cornerInd[ helper.WrapIndex( iC + nbC[1] + 1, cornerInd.size() )];
-        while ( curLengths.size() < nbEqualV + 1 )
+        vector< TGeoIndex > evVec( equVerts.begin(), equVerts.end() );
+        int          iEV = 0;
+        TGeoIndex    iE = cornerInd[ helper.WrapIndex( iC - nbC[0] - 1, cornerInd.size() )];
+        TGeoIndex iEEnd = cornerInd[ helper.WrapIndex( iC + nbC[1] + 1, cornerInd.size() )];
+        while ( accuLength.size() < nbEqualV + int( !allCornersSame ) )
         {
-          curLengths.push_back( totalLen );
+          // accumulate length of edges before iEV-th equal vertex
+          accuLength.push_back( totalLen );
           do {
-            curLengths.back() += SMESH_Algo::EdgeLength( edgeVec[ iE ]);
+            accuLength.back() += SMESH_Algo::EdgeLength( edgeVec[ iE ]);
             iE = helper.WrapIndex( iE + 1, edgeVec.size());
-            if ( iEV < evVec.size() && iE == evVec[ iEV++ ] )
-              break;
+            if ( iEV < evVec.size() && iE == evVec[ iEV ] ) {
+              iEV++;
+              break; // equal vertex reached
+            }
           }
           while( iE != iEEnd );
-          totalLen = curLengths.back();
+          totalLen = accuLength.back();
         }
-        curLengths.resize( equalVertices.size() );
-        for ( size_t iS = 0; iS < curLengths.size(); ++iS )
-          curLengths[ iS ] /= totalLen;
+        accuLength.resize( equVerts.size() );
+        for ( size_t iS = 0; iS < accuLength.size(); ++iS )
+          accuLength[ iS ] /= totalLen;
 
-        // find equalVertices most close to the ideal sub-division of all sides
+        // find equVerts most close to the ideal sub-division of all sides
         int iBestEV = 0;
         int iCorner = helper.WrapIndex( iC - nbC[0], cornerInd.size() );
-        int nbSides = 2 + nbC[0] + nbC[1];
+        int nbSides = Min( nbCorners, 2 + nbC[0] + nbC[1] );
         for ( int iS = 1; iS < nbSides; ++iS, ++iBestEV )
         {
           double idealLen = iS / double( nbSides );
-          double d, bestDist = 1.;
-          for ( iEV = iBestEV; iEV < curLengths.size(); ++iEV )
-            if (( d = Abs( idealLen - curLengths[ iEV ])) < bestDist )
+          double d, bestDist = 2.;
+          for ( iEV = iBestEV; iEV < accuLength.size(); ++iEV )
+          {
+            d = Abs( idealLen - accuLength[ iEV ]);
+
+            // take into account presence of a coresponding halfDivider
+            const double cornerWgt = 0.5  / nbSides;
+            const double vertexWgt = 0.25 / nbSides;
+            TGeoIndex hd = halfDivider[ evVec[ iEV ]];
+            if ( hd < 0 )
+              d += vertexWgt;
+            else if( refinedCorners.count( hd ))
+              d -= cornerWgt;
+            else
+              d -= vertexWgt;
+
+            // choose vertex with the best d
+            if ( d < bestDist )
             {
               bestDist = d;
               iBestEV  = iEV;
             }
+          }
           if ( iBestEV > iS-1 + nbExcessV )
             iBestEV = iS-1 + nbExcessV;
           theVertices[ iCorner ] = helper.IthVertex( 0, edgeVec[ evVec[ iBestEV ]]);
+          refinedCorners.insert( evVec[ iBestEV ]);
           iCorner = helper.WrapIndex( iCorner + 1, cornerInd.size() );
         }
+
+      } // if ( nbExcessV > 0 )
+      else
+      {
+        refinedCorners.insert( cornerInd[ iC ]);
       }
-    }
-  }
+    } // loop on cornerInd
+
+    // make theWire begin from the cornerInd[0]-th EDGE
+    while ( !theWire.front().IsSame( edgeVec[ cornerInd[0] ]))
+      theWire.splice( theWire.begin(), theWire, --theWire.end() );
+
+  } // if ( haveVariants )
 
   return nbCorners;
 }
@@ -4707,7 +4917,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 )
@@ -4746,7 +4956,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 )