"Clipping on Structured mesh". Store and provide access to info on structured mesh

diff --git a/src/VISU_I/VISU_Result_i.cc b/src/VISU_I/VISU_Result_i.cc
index 5b79dc1d097fa88cd64cb903fd44de1d72683b1b..98a17f0d84f91e49c9422eb604901fd474022af8 100644 (file)
--- a/src/VISU_I/VISU_Result_i.cc
+++ b/src/VISU_I/VISU_Result_i.cc
@@ -52,12 +52,18 @@
 
 #include <memory>
 #include <fstream>     
+#include <set>
 
 //#include <vtkRenderer.h>
 //#include <vtkActorCollection.h>
 
 #include <vtkUnstructuredGridReader.h>
 #include <vtkUnstructuredGridWriter.h>
+#include <vtkCell.h>
+
+#include <gp_Pnt.hxx>
+#include <gp_Vec.hxx>
+#include <Bnd_Box.hxx>
 
 using namespace VISU;
 using namespace std;
@@ -209,7 +215,7 @@ CORBA::Boolean VISU::Result_i::BuildAll(){
   try{
     const VISU::TMeshMap& aMeshMap = myInput->GetMeshMap();
     VISU::TMeshMap::const_iterator aMeshMapIter = aMeshMap.begin();
-    for(; aMeshMapIter != aMeshMap.end(); aMeshMapIter++){
+    for(; aMeshMapIter != aMeshMap.end(); aMeshMapIter++) {
       const string& aMeshName = aMeshMapIter->first;
       const VISU::PMesh aMesh = aMeshMapIter->second;
       const VISU::TMeshOnEntityMap& aMeshOnEntityMap = aMesh->myMeshOnEntityMap;
@@ -680,3 +686,168 @@ VISU::Result_i::~Result_i() {
   }
   if(myInput) delete myInput;
 }
+
+//=======================================================================
+//function : GetAxisInfo
+//purpose  : 
+//=======================================================================
+
+const vector< float >* Result_i::GetAxisInfo(const string& theMeshName,
+                                             TAxis         theAxis,
+                                             gp_Dir&       thePlaneNormal)
+{
+  const vector< float >* components = NULL;
+
+  if ( theAxis < AXIS_X || theAxis > AXIS_Z ) {
+    MESSAGE(" Bad axis index " << theAxis );
+    return components;
+  }
+    
+  map< string, TGridInfo >::iterator name_info;
+  name_info = myMeshName2GridInfoMap.find( theMeshName );
+  TGridInfo * gInfo = 0;
+
+  if ( name_info != myMeshName2GridInfoMap.end() )
+  {
+    gInfo = & name_info->second;
+  }
+  else if ( myInput && IsPossible() && theAxis >= AXIS_X && theAxis <= AXIS_Z )
+  {
+    // check presence of theMeshName
+    const VISU::TMeshMap& meshMap = myInput->GetMeshMap();
+    if ( meshMap.find( theMeshName ) == meshMap.end() ) {
+      MESSAGE("No mesh named " << theMeshName );
+      return components;
+    }      
+    VISU_Convertor::TOutput* vtkMesh = myInput->GetMeshOnEntity (theMeshName,
+                                                                 CELL_ENTITY);
+    if ( !vtkMesh || vtkMesh->GetNumberOfCells() == 0 ) {
+      MESSAGE( "No cells in the mesh: " << theMeshName );
+      return components;
+    }
+
+    // define axis directions and min cell size in each direction
+    const int nbAxes = 3;
+    int iAx;
+    gp_Vec axDirs[ nbAxes ];
+    float minSize[3] = { FLT_MAX, FLT_MAX, FLT_MAX };
+    bool axesComputed = false;
+    for ( vtkIdType iCell = 0; iCell < vtkMesh->GetNumberOfCells(); ++iCell )
+    {
+      vtkCell* cell = vtkMesh->GetCell( iCell );
+      int nbPnt = cell->GetNumberOfPoints();
+      if ( nbPnt != 8 )
+        continue;
+      vtkPoints * points = cell->GetPoints();
+      float* coords[ 4 ];
+      coords[0] = points->GetPoint( 0 );
+      coords[1] = points->GetPoint( 1 );
+      coords[2] = points->GetPoint( 3 );
+      coords[3] = points->GetPoint( 4 );
+      gp_Pnt p0( coords[0][0], coords[0][1], coords[0][2] );
+      for ( iAx = 0; iAx < nbAxes; ++iAx )
+      {
+        float* coo = coords[ iAx + 1 ];
+        gp_Pnt p( coo[0], coo[1], coo[2] );
+        // min size
+        float size = p0.SquareDistance( p );
+        if ( size > FLT_MIN && size < minSize[ iAx ] )
+          minSize[ iAx ] = size;
+        // axis direction
+        if ( !axesComputed ) {
+          gp_Vec dir( p0, p );
+          if ( dir.SquareMagnitude() <= DBL_MIN )
+            break;
+          axDirs[ iAx ] = dir;
+        }
+      }
+      if ( iAx == nbAxes )
+        axesComputed = true;
+    }
+    if ( !axesComputed ) {
+      MESSAGE("No good hexahedrons in the mesh: " << theMeshName );
+      return components;
+    }
+
+    // compute axes dirs
+    gInfo = & myMeshName2GridInfoMap[ theMeshName ];
+    for ( iAx = 0; iAx < nbAxes; ++iAx )
+    {
+      int iPrev = ( iAx == 0 ) ? 2 : iAx - 1;
+      int iNext = ( iAx == 2 ) ? 0 : iAx + 1;
+      gInfo->myAxis[ iAx ] = axDirs[ iPrev ] ^ axDirs[ iNext ];
+    }
+
+    // get and sort intermediate component values - projections of nodes
+    // on axis direction; define bnd box
+    set< float > comps[ 3 ];
+    Bnd_Box box;
+    vtkPoints * points = vtkMesh->GetPoints();
+    vtkIdType iP, nbP = vtkMesh->GetNumberOfPoints();
+    for ( iP = 0; iP < nbP; ++iP )
+    {
+      float* coo = points->GetPoint( iP );
+      gp_Pnt p( coo[0], coo[1], coo[2] );
+      box.Add( p );
+      for ( iAx = 0; iAx < nbAxes; ++iAx ) {
+        const gp_Dir& dir = gInfo->myAxis[ iAx ];
+        float dot = dir.XYZ() * p.XYZ();
+        comps[ iAx ].insert( dot );
+      }
+    }
+
+    // find a range of projections of bnd box corners on each axis
+    float range[3], firstValue[3]; 
+    double x[2],y[2],z[2];
+    box.Get(x[0],y[0],z[0],x[1],y[1],z[1]);
+    for ( iAx = 0; iAx < nbAxes; ++iAx ) {
+      set< float > bndComps;
+      const gp_Dir& dir = gInfo->myAxis[ iAx ];
+      for ( int iX = 0; iX < 2; ++iX ) {
+        for ( int iY = 0; iY < 2; ++iY ) {
+          for ( int iZ = 0; iZ < 2; ++iZ ) {
+            gp_Pnt p( x[ iX ], y[ iY ], z[ iZ ] );
+            float dot = dir.XYZ() * p.XYZ();
+            bndComps.insert( dot );
+          }
+        }
+      }
+      firstValue[ iAx ] = *bndComps.begin();
+      range[ iAx ] = *bndComps.rbegin() - *bndComps.begin();
+    }
+      
+    // compute component values
+    for ( iAx = 0; iAx < nbAxes; ++iAx )
+    {
+      list< float > values;
+      int nbVals = 0;
+      set< float >& comp = comps[ iAx ];
+      set< float >::iterator val = comp.begin();
+      float bnd = -1., rng = range[ iAx ], first = firstValue[ iAx ];
+      float tol = 0.1 * sqrt( minSize[ iAx ]) / rng;
+      for ( ; val != comp.end(); ++val ) {
+        float value = ( *val - first ) / rng;
+        if ( value > bnd ) {
+          values.push_back( value );
+          bnd = value + tol;
+          nbVals++;
+        }
+      }
+      // store values in gInfo
+      vector< float >& myComp = gInfo->myComponets[ iAx ];
+      myComp.resize( nbVals );
+      list< float >::iterator v = values.begin();
+      for ( int i = 0; v != values.end(); ++v )
+        myComp[ i++ ] = *v;
+    }
+  }
+
+  // set return values
+  if ( gInfo )
+  {
+    thePlaneNormal = gInfo->myAxis[ theAxis ];
+    components = & gInfo->myComponets[ theAxis ];
+  }
+
+  return components;
+}

diff --git a/src/VISU_I/VISU_Result_i.hh b/src/VISU_I/VISU_Result_i.hh
index f113595164c81a6686aa3c7d6674e1d2563b43d4..48fbd4b1dc6fababdc73f56b436a332a70d04535 100644 (file)
--- a/src/VISU_I/VISU_Result_i.hh
+++ b/src/VISU_I/VISU_Result_i.hh
@@ -30,6 +30,10 @@
 #include "VISUConfig.hh"
 #include "SALOME_GenericObj_i.hh"
 
+#include <gp_Dir.hxx>
+
+#include <vector>
+
 class VISU_Convertor;
 
 namespace VISU{
@@ -94,6 +98,21 @@ namespace VISU{
     const SALOMEDS::SComponent_var& GetSComponent() const;
     std::string GetEntry(const std::string& theComment);
 
+    // Info on structured mesh contained in TInput
+  public:
+    typedef enum { AXIS_X = 0, AXIS_Y, AXIS_Z } TAxis;
+    const std::vector< float >* GetAxisInfo(const std::string& theMeshName,
+                                            TAxis              theAxis,
+                                            gp_Dir&            thePlaneNormal);
+    // Return i,j or k values and cutting plane normal for theAxis.
+    // In the case of any problems, return NULL pointer
+  private:
+    struct TGridInfo {
+      std::vector< float > myComponets[ 3 ];
+      gp_Dir               myAxis     [ 3 ];
+    };
+    map< string, TGridInfo > myMeshName2GridInfoMap;
+
   };
 
   Result_var FindResult(SALOMEDS::SObject_ptr theSObject);