From 986c507e02a020ea13a0975d3a45769fda8f6ada Mon Sep 17 00:00:00 2001
From: eap <eap@opencascade.com>
Date: Wed, 18 Aug 2010 08:00:09 +0000
Subject: [PATCH] 0020855: [CEA] Problem with ijk algo      find out
 orientation by a least distorted hexahedron

---
 src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx | 56 +++++++++++++++----
 1 file changed, 44 insertions(+), 12 deletions(-)

diff --git a/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx b/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx
index 6e8448e0a..b4172bf27 100644
--- a/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx
+++ b/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx
@@ -31,6 +31,7 @@
 #include <gp_Ax2.hxx>
 
 //#include "utilities.h"
+#include <limits>
 
 // Define error message
 #ifdef _DEBUG_
@@ -1080,18 +1081,49 @@ bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper*
     // 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 ) {