Salome HOME
Merge branch 'V7_dev'
[plugins/hexablockplugin.git] / src / HEXABLOCKPlugin / HEXABLOCKPlugin_FromSkin_3D.cxx
index ab54a39f2e3c381b0380e2ff4c237d65e94eb3a9..65569e37aa5b564462311534c3d7afdd082c3b97 100755 (executable)
@@ -1,20 +1,20 @@
-//  Copyright (C) 2007-2010  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2009-2016  CEA/DEN, EDF R&D
 //
-//  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.
+// 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, 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
-//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-//  Lesser General Public License for more details.
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
 //
-//  You should have received a copy of the GNU Lesser General Public
-//  License along with this library; if not, write to the Free Software
-//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
 //
-//  See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
 
 // File      : SMESH_HexaFromSkin_3D.cxx
 #include "SMDS_VolumeTool.hxx"
 #include "SMESH_Block.hxx"
 #include "SMESH_MesherHelper.hxx"
+#include "SMESH_MeshEditor.hxx"
+#include "SMESH_MeshAlgos.hxx"
 
 #include <gp_Ax2.hxx>
 
 //#include "utilities.h"
+#include <limits>
 
-// Define error message
+// Define error message and _MYDEBUG_ if needed
 #ifdef _DEBUG_
 #define BAD_MESH_ERR \
   error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh.\n" \
                       __FILE__ ":" )<<__LINE__)
+//#define _MYDEBUG_
 #else
 #define BAD_MESH_ERR \
   error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh"))
 #endif
 
 // Debug output
-#define _DUMP_(msg) //cout << msg << endl
+#ifdef _MYDEBUG_
+#define _DUMP_(msg) cout << msg << endl
+#else
+#define _DUMP_(msg)
+#endif
 
 
 #ifdef _DEBUG_
-static int MYDEBUG = 1;
+static int MYDEBUG = HEXA_NS::on_debug ();
 #else
-static int MYDEBUG = 1;
+static int MYDEBUG = 0;
 #endif
 
 
@@ -59,10 +67,12 @@ namespace
     {
       B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, NB_BLOCK_SIDES
     };
+#ifdef _MYDEBUG_
   const char* SBoxSides[] = //!< names of block sides
     {
       "BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED"
     };
+#endif    
   enum EQuadEdge //!< edges of quadrangle side
     {
       Q_BOTTOM = 0, Q_RIGHT, Q_TOP, Q_LEFT, NB_QUAD_SIDES
@@ -90,37 +100,6 @@ namespace
     return true;
   }
 
-  //================================================================================
-  /*!
-   * \brief Description of node used to detect corner nodes
-   */
-  struct _NodeDescriptor
-  {
-    int nbInverseFaces, nbNodesInInverseFaces;
-
-    _NodeDescriptor(const SMDS_MeshNode* n): nbInverseFaces(0), nbNodesInInverseFaces(0)
-    {
-      if ( n )
-      {
-        set<const SMDS_MeshNode*> nodes;
-        SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face );
-        while ( fIt->more() )
-        {
-          const SMDS_MeshElement* face = fIt->next();
-          nodes.insert( face->begin_nodes(), face->end_nodes() );
-          nbInverseFaces++;
-        }
-        nbNodesInInverseFaces = nodes.size();
-      }
-    }
-    bool operator==(const _NodeDescriptor& other) const
-    {
-      return
-        nbInverseFaces == other.nbInverseFaces &&
-        nbNodesInInverseFaces == other.nbNodesInInverseFaces;
-    }
-  };
-
   //================================================================================
   /*!
    * \brief return true if a node is at block corner
@@ -131,7 +110,19 @@ namespace
 
   bool isCornerNode( const SMDS_MeshNode* n )
   {
-    return n && n->NbInverseElements( SMDSAbs_Face ) % 2;
+    int nbF = n ? n->NbInverseElements( SMDSAbs_Face ) : 1;
+    if ( nbF % 2 )
+      return true;
+
+    set<const SMDS_MeshNode*> nodesInInverseFaces;
+    SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face );
+    while ( fIt->more() )
+    {
+      const SMDS_MeshElement* face = fIt->next();
+      nodesInInverseFaces.insert( face->begin_nodes(), face->end_nodes() );
+    }
+
+    return nodesInInverseFaces.size() != ( 6 + (nbF/2-1)*3 );
   }
 
   //================================================================================
@@ -213,9 +204,9 @@ namespace
     typedef void (*TFun)(int& x, int& y);
     TFun _xRevFun, _yRevFun, _swapFun;
     
-    static void lazy(int&, int&) {}
+    static void lazy   (int&, int&) {}
     static void reverse(int& x, int& size) { x = size - x - 1; }
-    static void swap(int& x, int& y) { std::swap(x,y); }
+    static void swap   (int& x, int& y) { std::swap(x,y); }
   };
   //================================================================================
   /*!
@@ -249,10 +240,10 @@ namespace
       return getNode( isXMax ? _index._xSize-1 : 0 , isYMax ? _index._ySize-1 : 0 );
     }
     const SMDS_MeshElement* getCornerFace(const SMDS_MeshNode* cornerNode) const;
-    //!< True if all blocks this side belongs to have beem found
+    //!< True if all blocks this side belongs to have been found
     bool isBound() const { return _nbBlocksExpected <= _nbBlocksFound; }
     //!< Return coordinates of node at XY
-    gp_XYZ getXYZ(int x, int y) const { return SMESH_MeshEditor::TNodeXYZ( getNode( x, y )); }
+    gp_XYZ getXYZ(int x, int y) const { return SMESH_TNodeXYZ( getNode( x, y )); }
     //!< Return gravity center of the four corners and the middle node
     gp_XYZ getGC() const
     {
@@ -281,7 +272,7 @@ namespace
     //!< return coordinates by XY
     gp_XYZ xyz(int x, int y) const
     {
-      return SMESH_MeshEditor::TNodeXYZ( _grid_access_(_side, _index( x, y )) );
+      return SMESH_TNodeXYZ( _grid_access_(_side, _index( x, y )) );
     }
     //!< safely return a node by XY
     const SMDS_MeshNode* node(int x, int y) const
@@ -363,10 +354,11 @@ namespace
                              vector<const SMDS_MeshNode*>& verRow1,
                              vector<const SMDS_MeshNode*>& verRow2,
                              bool alongN1N2 );
-    _OrientedBlockSide findBlockSide( EBoxSides startBlockSide,
-                                      EQuadEdge sharedSideEdge1,
-                                      EQuadEdge sharedSideEdge2,
-                                      bool      withGeometricAnalysis);
+    _OrientedBlockSide findBlockSide( EBoxSides           startBlockSide,
+                                      EQuadEdge           sharedSideEdge1,
+                                      EQuadEdge           sharedSideEdge2,
+                                      bool                withGeometricAnalysis,
+                                      set< _BlockSide* >& sidesAround);
     //!< update own data and data of the side bound to block
     void setSideBoundToBlock( _BlockSide& side )
     {
@@ -529,20 +521,25 @@ namespace
       EQuadEdge edgeOfAdj  [4] = { Q_BOTTOM, Q_LEFT, Q_BOTTOM, Q_LEFT };
       // first find all sides detectable w/o advanced analysis,
       // then repeat the search, which then may pass without advanced analysis
+      set< _BlockSide* > sidesAround;
       for ( int advAnalys = 0; advAnalys < 2; ++advAnalys )
       {
+        // try to find 4 sides adjacent to a FRONT side
         for ( int i = 0; (ok || !advAnalys) && i < NB_QUAD_SIDES; ++i )
-          if ( !block._side[i] ) // try to find 4 sides adjacent to front side
-            ok = block.setSide( i, findBlockSide( B_FRONT, edgeOfFront[i],edgeOfAdj[i],advAnalys));
+          if ( !block._side[i] )
+            ok = block.setSide( i, findBlockSide( B_FRONT, edgeOfFront[i], edgeOfAdj[i],
+                                                  advAnalys, sidesAround));
+        // try to find a BACK side by a TOP one
         if ( ok || !advAnalys)
-          if ( !block._side[B_BACK] && block._side[B_TOP] ) // try to find back side by top one
-            ok = block.setSide( B_BACK, findBlockSide( B_TOP, Q_TOP, Q_TOP, advAnalys ));
+          if ( !block._side[B_BACK] && block._side[B_TOP] )
+            ok = block.setSide( B_BACK, findBlockSide( B_TOP, Q_TOP, Q_TOP,
+                                                       advAnalys, sidesAround ));
         if ( !advAnalys ) ok = true;
       }
       ok = block.isValid();
       if ( ok )
       {
-        // check if just found block is same as one of previously found
+        // check if just found block is same as one of previously found blocks
         bool isSame = false;
         for ( int i = 1; i < _blocks.size() && !isSame; ++i )
           isSame = ( block._corners == _blocks[i-1]._corners );
@@ -646,8 +643,8 @@ namespace
       const SMDS_MeshNode* n1down = side.getNode( 0, y-1 );
       const SMDS_MeshNode* n2down = side.getNode( 1, y-1 );
       avoidSet.clear(); avoidSet.insert( firstQuad );
-      firstQuad = SMESH_MeshEditor::FindFaceInSet( n1down, n2down, emptySet, avoidSet,
-                                                   &i1down, &i2down);
+      firstQuad = SMESH_MeshAlgos::FindFaceInSet( n1down, n2down, emptySet, avoidSet,
+                                                  &i1down, &i2down);
       if ( !isQuadrangle( firstQuad ))
         return BAD_MESH_ERR;
 
@@ -659,8 +656,8 @@ namespace
       x = 1; 
       while ( ++x < nbX )
       {
-        const SMDS_MeshElement* quad = SMESH_MeshEditor::FindFaceInSet( n2up, n2down, emptySet,
-                                                                        avoidSet, &i2up, &i2down);
+        const SMDS_MeshElement* quad = SMESH_MeshAlgos::FindFaceInSet( n2up, n2down, emptySet,
+                                                                       avoidSet, &i2up, &i2down);
         if ( !isQuadrangle( quad ))
           return BAD_MESH_ERR;
 
@@ -682,16 +679,77 @@ namespace
     return ok;
   }
 
+  //================================================================================
+  /*!
+   * \brief Return true if it's possible to make a loop over corner2Sides starting
+   * from the startSide
+   */
+  //================================================================================
+
+  bool isClosedChainOfSides( _BlockSide*                                        startSide,
+                             map< const SMDS_MeshNode*, list< _BlockSide* > > & corner2Sides )
+  {
+    // get start and end nodes
+    const SMDS_MeshNode *n1 = 0, *n2 = 0, *n;
+    for ( int y = 0; y < 2; ++y )
+      for ( int x = 0; x < 2; ++x )
+      {
+        n = startSide->getCornerNode(x,y);
+        if ( !corner2Sides.count( n )) continue;
+        if ( n1 )
+          n2 = n;
+        else
+          n1 = n;
+      }
+    if ( !n2 ) return false;
+
+    map< const SMDS_MeshNode*, list< _BlockSide* > >::iterator
+      c2sides = corner2Sides.find( n1 );
+    if ( c2sides == corner2Sides.end() ) return false;
+
+    int nbChainLinks = 1;
+    n = n1;
+    _BlockSide* prevSide = startSide;
+    while ( n != n2 )
+    {
+      // get the next side sharing n
+      list< _BlockSide* > & sides = c2sides->second;
+      _BlockSide* nextSide = ( sides.back() == prevSide ? sides.front() : sides.back() );
+      if ( nextSide == prevSide ) return false;
+
+      // find the next corner of the nextSide being in corner2Sides
+      n1 = n;
+      n = 0;
+      for ( int y = 0; y < 2 && !n; ++y )
+        for ( int x = 0; x < 2; ++x )
+        {
+          n = nextSide->getCornerNode(x,y);
+          c2sides = corner2Sides.find( n );
+          if ( n == n1 || c2sides == corner2Sides.end() )
+            n = 0;
+          else
+            break;
+        }
+      if ( !n ) return false;
+
+      prevSide = nextSide;
+      nbChainLinks++;
+    }
+
+    return ( n == n2 && nbChainLinks == NB_QUAD_SIDES );
+  }
+
   //================================================================================
   /*!
    * \brief Try to find a block side adjacent to the given side by given edge
    */
   //================================================================================
 
-  _OrientedBlockSide _Skin::findBlockSide( EBoxSides startBlockSide,
-                                           EQuadEdge sharedSideEdge1,
-                                           EQuadEdge sharedSideEdge2,
-                                           bool      withGeometricAnalysis)
+  _OrientedBlockSide _Skin::findBlockSide( EBoxSides           startBlockSide,
+                                           EQuadEdge           sharedSideEdge1,
+                                           EQuadEdge           sharedSideEdge2,
+                                           bool                withGeometricAnalysis,
+                                           set< _BlockSide* >& sidesAround)
   {
     _Block& block = _blocks.back();
     _OrientedBlockSide& side1 = block._side[ startBlockSide ];
@@ -742,45 +800,75 @@ namespace
 
       if ( !foundSide )
       {
-        if ( !withGeometricAnalysis ) return 0;
-
-        // Select one of found sides most close to startBlockSide
-
-        gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()),  p2 (n2->X(),n2->Y(),n2->Z());
-        gp_Vec p1p2( p1, p2 );
-
-        const SMDS_MeshElement* face1 = side1->getCornerFace( n1 );
-        gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(n1)));
-        gp_Vec side1Dir( p1, p1Op );
-        gp_Ax2 pln( p1, p1p2, side1Dir ); // plane with normal p1p2 and X dir side1Dir
-        _DUMP_("  Select adjacent for "<< side1._side << " - side dir ("
-               << side1Dir.X() << ", " << side1Dir.Y() << ", " << side1Dir.Z() << ")" );
-
-        map < double , _BlockSide* > angleOfSide;
-        for (sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
+        if ( !withGeometricAnalysis )
         {
-          _BlockSide* sideI = *sideIt;
-          const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 );
-          gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(n1)));
-          gp_Vec sideIDir( p1, p1Op );
-          // compute angle of (sideIDir projection to pln) and (X dir of pln)
-          gp_Vec2d sideIDirProj( sideIDir * pln.XDirection(), sideIDir * pln.YDirection());
-          double angle = sideIDirProj.Angle( gp::DX2d() );
-          if ( angle < 0 ) angle += 2 * PI; // angle [0-2*PI]
-          angleOfSide.insert( make_pair( angle, sideI ));
-          _DUMP_("  "<< sideI << " - side dir ("
-                 << sideIDir.X() << ", " << sideIDir.Y() << ", " << sideIDir.Z() << ")"
-                 << " angle " << angle);
+          sidesAround.insert( sidesOnEdge.begin(), sidesOnEdge.end() );
+          return 0;
         }
         if ( nbLoadedSides == 1 )
         {
-          double angF = angleOfSide.begin()->first, angL = angleOfSide.rbegin()->first;
-          if ( angF > PI ) angF = 2*PI - angF;
-          if ( angL > PI ) angL = 2*PI - angL;
-          foundSide = angF < angL ? angleOfSide.begin()->second : angleOfSide.rbegin()->second;
+          // Issue 0021529. There are at least 2 sides by each edge and
+          // position of block gravity center is undefined.
+          // Find a side starting from which we can walk around the startBlockSide
+
+          // fill in corner2Sides
+          map< const SMDS_MeshNode*, list< _BlockSide* > > corner2Sides;
+          for ( sideIt = sidesAround.begin(); sideIt != sidesAround.end(); ++sideIt )
+          {
+            _BlockSide* sideI = *sideIt;
+            corner2Sides[ sideI->getCornerNode(0,0) ].push_back( sideI );
+            corner2Sides[ sideI->getCornerNode(1,0) ].push_back( sideI );
+            corner2Sides[ sideI->getCornerNode(0,1) ].push_back( sideI );
+            corner2Sides[ sideI->getCornerNode(1,1) ].push_back( sideI );
+          }
+          // remove corners of startBlockSide from corner2Sides
+          set<const SMDS_MeshNode*>::iterator nIt = block._corners.begin();
+          for ( ; nIt != block._corners.end(); ++nIt )
+            corner2Sides.erase( *nIt );
+
+          // select a side
+          for ( sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
+          {
+            if ( isClosedChainOfSides( *sideIt, corner2Sides ))
+            {
+              foundSide = *sideIt;
+              break;
+            }
+          }
+          if ( !foundSide )
+            return 0;
         }
         else
         {
+          // Select one of found sides most close to startBlockSide
+
+          gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()),  p2 (n2->X(),n2->Y(),n2->Z());
+          gp_Vec p1p2( p1, p2 );
+
+          const SMDS_MeshElement* face1 = side1->getCornerFace( n1 );
+          gp_XYZ p1Op = SMESH_TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(n1)));
+          gp_Vec side1Dir( p1, p1Op );
+          gp_Ax2 pln( p1, p1p2, side1Dir ); // plane with normal p1p2 and X dir side1Dir
+          _DUMP_("  Select adjacent for "<< side1._side << " - side dir ("
+                 << side1Dir.X() << ", " << side1Dir.Y() << ", " << side1Dir.Z() << ")" );
+
+          map < double , _BlockSide* > angleOfSide;
+          for (sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
+          {
+            _BlockSide* sideI = *sideIt;
+            const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 );
+            gp_XYZ p1Op = SMESH_TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(n1)));
+            gp_Vec sideIDir( p1, p1Op );
+            // compute angle of (sideIDir projection to pln) and (X dir of pln)
+            gp_Vec2d sideIDirProj( sideIDir * pln.XDirection(), sideIDir * pln.YDirection());
+            double angle = sideIDirProj.Angle( gp::DX2d() );
+            if ( angle < 0 ) angle += 2. * M_PI; // angle [0-2*PI]
+            angleOfSide.insert( make_pair( angle, sideI ));
+            _DUMP_("  "<< sideI << " - side dir ("
+                   << sideIDir.X() << ", " << sideIDir.Y() << ", " << sideIDir.Z() << ")"
+                   << " angle " << angle);
+          }
+
           gp_XYZ gc(0,0,0); // gravity center of already loaded block sides
           for (int i = 0; i < NB_BLOCK_SIDES; ++i )
             if ( block._side[ i ] )
@@ -854,17 +942,15 @@ namespace
       row2.push_back( n1 = oppositeNode( quad, i1 ));
     }
 
-    _NodeDescriptor nDesc( row1[1] );
-    if ( nDesc == _NodeDescriptor( row1[0] ))
-      return true; // row of 2 nodes
+    if ( isCornerNode( row1[1] ))
+      return true;
 
     // Find the rest nodes
     TIDSortedElemSet emptySet, avoidSet;
-    //while ( !isCornerNode( n2 ))
-    while ( nDesc == _NodeDescriptor( n2 ))
+    while ( !isCornerNode( n2 ) )
     {
       avoidSet.clear(); avoidSet.insert( quad );
-      quad = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2 );
+      quad = SMESH_MeshAlgos::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2 );
       if ( !isQuadrangle( quad ))
         return BAD_MESH_ERR;
 
@@ -960,9 +1046,9 @@ HEXA_NS::Hexa* _block2Hexa( const _Block& block,
 
 
     HEXA_NS::Hexa* hexa = NULL;
-    int nHexa = doc->countHexa();
+    int nHexa = doc->countUsedHexa();
     for (int j=0; j <nHexa; ++j ){
-      hexa = doc->getHexa(j);
+      hexa = doc->getUsedHexa(j);
       list<const SMDS_MeshNode* > nodeFromHexa;
       int nVx = hexa->countVertex();
       for ( int i=0; i <nVx; ++i ){
@@ -1005,14 +1091,14 @@ SMESH_HexaFromSkin_3D::SMESH_HexaFromSkin_3D(int hypId, int studyId, SMESH_Gen*
   :SMESH_3D_Algo(hypId, studyId, gen),
   _doc( doc )
 {
-  MESSAGE("SMESH_HexaFromSkin_3D::SMESH_HexaFromSkin_3D");
+  if(MYDEBUG) MESSAGE("SMESH_HexaFromSkin_3D::SMESH_HexaFromSkin_3D");
   _name = "HexaFromSkin_3D";
 }
 
 
 SMESH_HexaFromSkin_3D::~SMESH_HexaFromSkin_3D()
 {
-  MESSAGE("SMESH_HexaFromSkin_3D::~SMESH_HexaFromSkin_3D");
+  if(MYDEBUG) MESSAGE("SMESH_HexaFromSkin_3D::~SMESH_HexaFromSkin_3D");
 }
 
 //================================================================================
@@ -1082,7 +1168,7 @@ bool SMESH_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHel
     // ----------------------------
     // Add internal nodes of a box
     // ----------------------------
-    // projection points of internal nodes on box subshapes by which
+    // projection points of internal nodes on box sub-shapes by which
     // coordinates of internal nodes are computed
     vector<gp_XYZ> pointOnShape( SMESH_Block::ID_Shell );
 
@@ -1154,18 +1240,49 @@ bool SMESH_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHel
     // Add hexahedrons
     // ----------------
 
-    // find out orientation
-    const SMDS_MeshNode* n000 = block.getSide(B_BOTTOM).cornerNode( 0, 0 );
-    const SMDS_MeshNode* n100 = block.getSide(B_BOTTOM).cornerNode( 1, 0 );
-    const SMDS_MeshNode* n010 = block.getSide(B_BOTTOM).cornerNode( 0, 1 );
-    const SMDS_MeshNode* n110 = block.getSide(B_BOTTOM).cornerNode( 1, 1 );
-    const SMDS_MeshNode* n001 = block.getSide(B_TOP).cornerNode( 0, 0 );
-    const SMDS_MeshNode* n101 = block.getSide(B_TOP).cornerNode( 1, 0 );
-    const SMDS_MeshNode* n011 = block.getSide(B_TOP).cornerNode( 0, 1 );
-    const SMDS_MeshNode* n111 = block.getSide(B_TOP).cornerNode( 1, 1 );
-    SMDS_VolumeOfNodes probeVolume (n000,n010,n110,n100,
-                                    n001,n011,n111,n101);
-    bool isForw = SMDS_VolumeTool( &probeVolume ).IsForward();
+    // find out orientation by a least distorted hexahedron (issue 0020855);
+    // the last is defined by evaluating sum of face normals of 8 corner hexahedrons
+    double badness = numeric_limits<double>::max();
+    bool isForw = true;
+    for ( int xMax = 0; xMax < 2; ++xMax )
+      for ( int yMax = 0; yMax < 2; ++yMax )
+        for ( int zMax = 0; zMax < 2; ++zMax )
+        {
+          x = xMax ? xSize-1 : 1;
+          y = yMax ? ySize-1 : 1;
+          z = zMax ? zSize-1 : 1;
+          vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x-1, y-1 )];
+          vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x  , y-1 )];
+          vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x-1, y   )];
+          vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x  , y )];
+          
+          const SMDS_MeshNode* n000 = col00[z-1];
+          const SMDS_MeshNode* n100 = col10[z-1];
+          const SMDS_MeshNode* n010 = col01[z-1];
+          const SMDS_MeshNode* n110 = col11[z-1];
+          const SMDS_MeshNode* n001 = col00[z];
+          const SMDS_MeshNode* n101 = col10[z];
+          const SMDS_MeshNode* n011 = col01[z];
+          const SMDS_MeshNode* n111 = col11[z];
+          SMDS_VolumeOfNodes probeVolume (n000,n010,n110,n100,
+                                          n001,n011,n111,n101);
+          SMDS_VolumeTool volTool( &probeVolume );
+          double Nx=0.,Ny=0.,Nz=0.;
+          for ( int iFace = 0; iFace < volTool.NbFaces(); ++iFace )
+          {
+            double nx,ny,nz;
+            volTool.GetFaceNormal( iFace, nx,ny,nz );
+            Nx += nx;
+            Ny += ny;
+            Nz += nz;
+          }
+          double quality = Nx*Nx + Ny*Ny + Nz*Nz;
+          if ( quality < badness )
+          {
+            badness = quality;
+            isForw = volTool.IsForward();
+          }
+        }
 
     // add elements
     for ( x = 0; x < xSize-1; ++x ) {
@@ -1199,7 +1316,7 @@ bool SMESH_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHel
 bool SMESH_HexaFromSkin_3D::Compute( SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper,
       std::map<HEXA_NS::Hexa*, SMESH_HexaBlocks::SMESHVolumes>& volumesOnHexa,
       std::map<HEXA_NS::Vertex*, SMDS_MeshNode*> vertexNode )
-{
+          {
   if(MYDEBUG) MESSAGE("SMESH_HexaFromSkin_3D::Compute BEGIN");
   _Skin skin;
   int nbBlocks = skin.findBlocks(aMesh);
@@ -1234,7 +1351,7 @@ bool SMESH_HexaFromSkin_3D::Compute( SMESH_Mesh & aMesh, SMESH_MesherHelper* aHe
       for ( z = 0; z < zSize; ++z ) {
         column0[ z ] = block.getSide(B_FRONT).node( x, z );
         column1[ z ] = block.getSide(B_BACK) .node( x, z );
-      }
+          }
     }
     // fill node columns by left and right box sides
     for ( y = 1; y < ySize-1; ++y ) {
@@ -1275,7 +1392,7 @@ bool SMESH_HexaFromSkin_3D::Compute( SMESH_Mesh & aMesh, SMESH_MesherHelper* aHe
     pointOnShape[ SMESH_Block::ID_V111 ] = block.getSide(B_TOP).xyz( X, Y );
 
     for ( x = 1; x < xSize-1; ++x )
-    {
+          {
       gp_XYZ params; // normalized parameters of internal node within a unit box
       params.SetCoord( 1, x / double(X) );
       for ( y = 1; y < ySize-1; ++y )
@@ -1325,8 +1442,8 @@ bool SMESH_HexaFromSkin_3D::Compute( SMESH_Mesh & aMesh, SMESH_MesherHelper* aHe
           //cout << "Params: ( "<< params.X()<<", "<<params.Y()<<", "<<params.Z()<<" )"<<endl;
           //cout << "coords: ( "<< coords.X()<<", "<<coords.Y()<<", "<<coords.Z()<<" )"<<endl;
 #endif
+          }
         }
-      }
     }
     // ----------------
     // Add hexahedrons
@@ -1368,9 +1485,9 @@ bool SMESH_HexaFromSkin_3D::Compute( SMESH_Mesh & aMesh, SMESH_MesherHelper* aHe
             aHelper->AddVolume(col00[z],   col10[z],   col11[z],   col01[z],
                                col00[z+1], col10[z+1], col11[z+1], col01[z+1]);
             volumesOnBlock.push_back( newVolume );
-          }
       }
     }
+    }
 //     std::cout << "block i = " << i << std::endl;
     HEXA_NS::Hexa* currentHexa = _block2Hexa( block, _doc, vertexNode );
     if ( currentHexa != NULL ){
@@ -1391,14 +1508,6 @@ bool SMESH_HexaFromSkin_3D::Compute( SMESH_Mesh & aMesh, SMESH_MesherHelper* aHe
   return true;
 }
 
-
-
-
-
-
-
-
-
 //================================================================================
 /*!
  * \brief Evaluate nb of hexa
@@ -1470,3 +1579,4 @@ bool SMESH_HexaFromSkin_3D::Compute(SMESH_Mesh&, const TopoDS_Shape&)
   return error("Algorithm can't work with geometrical shapes");
 }
 
+