fix file header (comments)

[modules/smesh.git] / src / StdMeshers / StdMeshers_HexaFromSkin_3D.cxx

diff --git a/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx b/src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx
index 6e8448e0a2092d0bf0226dc03dc77a8393f55df2..2516de0848f44709e8680ebc5c087e7be3e6d09f 100644 (file)
--- 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_
@@ -53,10 +54,10 @@ namespace
     {
       B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, NB_BLOCK_SIDES
     };
-  const char* SBoxSides[] = //!< names of block sides
-    {
-      "BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED"
-    };
+//   const char* SBoxSides[] = //!< names of block sides -- needed for DEBUG only
+//     {
+//       "BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED"
+//     };
   enum EQuadEdge //!< edges of quadrangle side
     {
       Q_BOTTOM = 0, Q_RIGHT, Q_TOP, Q_LEFT, NB_QUAD_SIDES
@@ -246,7 +247,7 @@ namespace
     //!< True if all blocks this side belongs to have beem 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
     {
@@ -275,7 +276,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
@@ -744,7 +745,7 @@ namespace
         gp_Vec p1p2( p1, p2 );
 
         const SMDS_MeshElement* face1 = side1->getCornerFace( n1 );
-        gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(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 ("
@@ -755,7 +756,7 @@ namespace
         {
           _BlockSide* sideI = *sideIt;
           const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 );
-          gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(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());
@@ -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 ) {