Salome HOME
0021336: EDF 1717 SMESH: New algorithm "body fitting" cartesian unstructured
[modules/smesh.git] / src / StdMeshers / StdMeshers_Hexa_3D.cxx
index 917f48e70527fed8395cc3382faf7490570f5b45..ceb29fa71f0e67e7e93b76e5bc025194eafb13fd 100644 (file)
@@ -1,23 +1,23 @@
-//  Copyright (C) 2007-2010  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2011  CEA/DEN, EDF R&D, OPEN CASCADE
 //
-//  Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
-//  CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+// Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
 //
-//  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.
 //
-//  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
 //
 
 //  SMESH SMESH : implementaion of SMESH idl descriptions
 #include "SMESH_MesherHelper.hxx"
 #include "SMESH_subMesh.hxx"
 
-#include "SMDS_MeshElement.hxx"
 #include "SMDS_MeshNode.hxx"
-#include "SMDS_FacePosition.hxx"
-#include "SMDS_VolumeTool.hxx"
-#include "SMDS_VolumeOfNodes.hxx"
 
 #include <TopExp.hxx>
 #include <TopExp_Explorer.hxx>
-#include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
-#include <TopTools_ListIteratorOfListOfShape.hxx>
-#include <TopTools_ListOfShape.hxx>
 #include <TopTools_SequenceOfShape.hxx>
 #include <TopTools_MapOfShape.hxx>
 #include <TopoDS.hxx>
-#include <gp_Pnt2d.hxx>
 
 #include "utilities.h"
 #include "Utils_ExceptHandlers.hxx"
@@ -157,24 +149,23 @@ namespace
 {
   //=============================================================================
 
+  typedef boost::shared_ptr< FaceQuadStruct > FaceQuadStructPtr;
+
   // symbolic names of box sides
   enum EBoxSides{ B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_NB_SIDES };
 
-  // indices of FACE's of box sides in terms of SMESH_Block::TShapeID enum
-  enum ESideIndex{ I_BOTTOM    = SMESH_Block::ID_Fxy0,
-                   I_RIGHT     = SMESH_Block::ID_F1yz,
-                   I_TOP       = SMESH_Block::ID_Fxy1,
-                   I_LEFT      = SMESH_Block::ID_F0yz,
-                   I_FRONT     = SMESH_Block::ID_Fx0z,
-                   I_BACK      = SMESH_Block::ID_Fx1z,
-                   I_UNDEFINED = SMESH_Block::ID_NONE
-  };
+  // symbolic names of sides of quadrangle
+  enum EQuadSides{ Q_BOTTOM=0, Q_RIGHT, Q_TOP, Q_LEFT, Q_NB_SIDES };
+
   //=============================================================================
   /*!
    * \brief Container of nodes of structured mesh on a qudrangular geom FACE
    */
   struct _FaceGrid
   {
+    // face sides
+    FaceQuadStructPtr _quad;
+
     // map of (node parameter on EDGE) to (column (vector) of nodes)
     TParam2ColumnMap _u2nodesMap;
 
@@ -183,7 +174,6 @@ namespace
 
     // geometry of a cube side
     TopoDS_Face _sideF;
-    TopoDS_Edge _baseE;
 
     const SMDS_MeshNode* GetNode(int iCol, int iRow) const
     {
@@ -191,7 +181,7 @@ namespace
     }
     gp_XYZ GetXYZ(int iCol, int iRow) const
     {
-      return SMESH_MeshEditor::TNodeXYZ( GetNode( iCol, iRow ));
+      return SMESH_TNodeXYZ( GetNode( iCol, iRow ));
     }
   };
 
@@ -206,6 +196,79 @@ namespace
     int size() const { return _xSize * _ySize; }
     int operator()(const int x, const int y) const { return y * _xSize + x; }
   };
+
+  //================================================================================
+  /*!
+   * \brief Appends a range of node columns from a map to another map
+   */
+  template< class TMapIterator >
+  void append( TParam2ColumnMap& toMap, TMapIterator from, TMapIterator to )
+  {
+    const SMDS_MeshNode* lastNode = toMap.rbegin()->second[0];
+    const SMDS_MeshNode* firstNode = from->second[0];
+    if ( lastNode == firstNode )
+      from++;
+    double u = toMap.rbegin()->first;
+    for (; from != to; ++from )
+    {
+      u += 1;
+      TParam2ColumnMap::iterator u2nn = toMap.insert( toMap.end(), make_pair ( u, TNodeColumn()));
+      u2nn->second.swap( from->second );
+    }
+  }
+
+  //================================================================================
+  /*!
+   * \brief Finds FaceQuadStruct having a side equal to a given one and rearranges
+   *  the found FaceQuadStruct::side to have the given side at a Q_BOTTOM place
+   */
+  FaceQuadStructPtr getQuadWithBottom( StdMeshers_FaceSide* side,
+                                       FaceQuadStructPtr    quad[ 6 ])
+  {
+    FaceQuadStructPtr foundQuad;
+    for ( int i = 1; i < 6; ++i )
+    {
+      if ( !quad[i] ) continue;
+      for ( unsigned iS = 0; iS < quad[i]->side.size(); ++iS )
+      {
+        const StdMeshers_FaceSide* side2 = quad[i]->side[iS];
+        if (( side->FirstVertex().IsSame( side2->FirstVertex() ) ||
+              side->FirstVertex().IsSame( side2->LastVertex() ))
+            &&
+            ( side->LastVertex().IsSame( side2->FirstVertex() ) ||
+              side->LastVertex().IsSame( side2->LastVertex() ))
+            )
+        {
+          if ( iS != Q_BOTTOM )
+          {
+            vector< StdMeshers_FaceSide*> newSides;
+            for ( unsigned j = iS; j < quad[i]->side.size(); ++j )
+              newSides.push_back( quad[i]->side[j] );
+            for ( unsigned j = 0; j < iS; ++j )
+              newSides.push_back( quad[i]->side[j] );
+            quad[i]->side.swap( newSides );
+          }
+          foundQuad.swap(quad[i]);
+          return foundQuad;
+        }
+      }
+    }
+    return foundQuad;
+  }
+  //================================================================================
+  /*!
+   * \brief Returns true if the 1st base node of sideGrid1 belongs to sideGrid2
+   */
+  //================================================================================
+
+  bool beginsAtSide( const _FaceGrid& sideGrid1, const _FaceGrid& sideGrid2 )
+  {
+    const SMDS_MeshNode* n00 = (sideGrid1._u2nodesMap.begin()->second)[0];
+    const TNodeColumn& col0  = sideGrid2._u2nodesMap.begin()->second;
+    const TNodeColumn& col1  = sideGrid2._u2nodesMap.rbegin()->second;
+    return ( n00 == col0.front() || n00 == col0.back() ||
+             n00 == col1.front() || n00 == col1.back() );
+  }
 }
 
 //=============================================================================
@@ -226,20 +289,18 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
   MESSAGE("StdMeshers_Hexa_3D::Compute");
   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
 
-  // 1) Shape verification
+  // Shape verification
   // ----------------------
 
   // shape must be a solid (or a shell) with 6 faces
-  TopoDS_Shell shell;
   TopExp_Explorer exp(aShape,TopAbs_SHELL);
   if ( !exp.More() )
     return error(COMPERR_BAD_SHAPE, "No SHELL in the geometry");
-  shell = TopoDS::Shell( exp.Current());
   if ( exp.Next(), exp.More() )
     return error(COMPERR_BAD_SHAPE, "More than one SHELL in the geometry");
 
   TopTools_IndexedMapOfShape FF;
-  TopExp::MapShapes( shell, TopAbs_FACE, FF);
+  TopExp::MapShapes( aShape, TopAbs_FACE, FF);
   if ( FF.Extent() != 6)
   {
     static StdMeshers_CompositeHexa_3D compositeHexa(_gen->GetANewId(), 0, _gen);
@@ -248,13 +309,38 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
     return true;
   }
 
-  // Find sub-shapes of a cube
-  TopTools_IndexedMapOfOrientedShape cubeSubShapes;
-  TopoDS_Vertex V;
-  if ( !SMESH_Block::FindBlockShapes( shell, V,V, cubeSubShapes ))
-    return error(COMPERR_BAD_SHAPE, "Not a 6 sides cube");
+  // Find sides of a cube
+  // ---------------------
+  
+  FaceQuadStructPtr quad[ 6 ];
+  StdMeshers_Quadrangle_2D quadAlgo( _gen->GetANewId(), GetStudyId(), _gen);
+  for ( int i = 0; i < 6; ++i )
+  {
+    if ( !( quad[i] = FaceQuadStructPtr( quadAlgo.CheckNbEdges( aMesh, FF( i+1 )))))
+      return error( quadAlgo.GetComputeError() );
+    if ( quad[i]->side.size() != 4 )
+      return error( COMPERR_BAD_SHAPE, "Not a quadrangular box side" );
+  }
+
+  _FaceGrid aCubeSide[ 6 ];
+
+  swap( aCubeSide[B_BOTTOM]._quad, quad[0] );
+  swap( aCubeSide[B_BOTTOM]._quad->side[ Q_RIGHT],// direct the normal of bottom quad inside cube
+        aCubeSide[B_BOTTOM]._quad->side[ Q_LEFT ] );
+
+  aCubeSide[B_FRONT]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_BOTTOM], quad );
+  aCubeSide[B_RIGHT]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_RIGHT ], quad );
+  aCubeSide[B_BACK ]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_TOP   ], quad );
+  aCubeSide[B_LEFT ]._quad = getQuadWithBottom( aCubeSide[B_BOTTOM]._quad->side[Q_LEFT  ], quad );
+  if ( aCubeSide[B_FRONT ]._quad )
+    aCubeSide[B_TOP]._quad = getQuadWithBottom( aCubeSide[B_FRONT ]._quad->side[Q_TOP ], quad );
 
-  // 2) make viscous layers
+  for ( int i = 1; i < 6; ++i )
+    if ( !aCubeSide[i]._quad )
+      return error( COMPERR_BAD_SHAPE );
+
+  // Make viscous layers
+  // --------------------
 
   SMESH_ProxyMesh::Ptr proxymesh;
   if ( _viscousLayersHyp )
@@ -264,61 +350,15 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
       return false;
   }
 
-  // 3) Check presence of regular grid mesh on FACEs of the cube
-  // ------------------------------------------------------------
-
-  // indices of FACEs of cube sides within cubeSubShapes
-  const int sideIndex[6] = { I_BOTTOM, I_RIGHT, I_TOP, I_LEFT, I_FRONT, I_BACK };
-  // indices of base EDGEs of cube sides within cubeSubShapes, corresponding to sideIndex
-  const int baseEdgeIndex[6] = {
-    SMESH_Block::ID_Ex00, // bottom side
-    SMESH_Block::ID_E1y0, // right side
-    SMESH_Block::ID_Ex01, // top side
-    SMESH_Block::ID_E0y0, // left side
-    SMESH_Block::ID_Ex00, // front side
-    SMESH_Block::ID_Ex10  // back side
-  };
-
-  // Load mesh of cube sides
+  // Check if there are triangles on cube sides
+  // -------------------------------------------
 
-  _FaceGrid aCubeSide[ 6 ] ;
-
-  // tool creating quadratic elements if needed
-  SMESH_MesherHelper helper (aMesh);
-  _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
-
-  for ( int i = 0; i < 6; ++i )
+  if ( aMesh.NbTriangles() > 0 )
   {
-    aCubeSide[i]._sideF = TopoDS::Face( cubeSubShapes( sideIndex[i] ));
-    aCubeSide[i]._baseE = TopoDS::Edge( cubeSubShapes( baseEdgeIndex[i] ));
-
-    // assure correctness of node positions on _baseE
-    if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( aCubeSide[i]._baseE ))
-    {
-      bool ok;
-      helper.SetSubShape( aCubeSide[i]._baseE );
-      SMDS_ElemIteratorPtr eIt = smDS->GetElements();
-      while ( eIt->more() )
-      {
-        const SMDS_MeshElement* e = eIt->next();
-        helper.GetNodeU( aCubeSide[i]._baseE, e->GetNode(0), e->GetNode(1), &ok);
-        helper.GetNodeU( aCubeSide[i]._baseE, e->GetNode(1), e->GetNode(0), &ok);
-      }
-    }
-
-    // load grid
-    if ( !helper.LoadNodeColumns( aCubeSide[i]._u2nodesMap,
-                                  aCubeSide[i]._sideF,
-                                  aCubeSide[i]._baseE, meshDS, proxymesh.get() ))
-    {
-      SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape, proxymesh.get());
-      return error( err );
-    }
-
-    // check if there are triangles on aCubeSide[i]._sideF
-    if ( aMesh.NbTriangles() > 0 )
+    for ( int i = 0; i < 6; ++i )
     {
-      if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( aCubeSide[i]._sideF ))
+      const TopoDS_Face& sideF = aCubeSide[i]._quad->face;
+      if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( sideF ))
       {
         bool isAllQuad = true;
         SMDS_ElemIteratorPtr fIt = smDS->GetElements();
@@ -336,27 +376,75 @@ bool StdMeshers_Hexa_3D::Compute(SMESH_Mesh &         aMesh,
     }
   }
 
-  // Orient loaded grids of cube sides along axis of the unitary cube coord system
+  // Check presence of regular grid mesh on FACEs of the cube
+  // ------------------------------------------------------------
+
+  // tool creating quadratic elements if needed
+  SMESH_MesherHelper helper (aMesh);
+  _quadraticMesh = helper.IsQuadraticSubMesh(aShape);
+
   for ( int i = 0; i < 6; ++i )
   {
-    bool reverse = false;
-    if ( helper.GetSubShapeOri( shell.Oriented( TopAbs_FORWARD ),
-                                aCubeSide[i]._sideF ) == TopAbs_REVERSED )
-      reverse = !reverse;
+    const TopoDS_Face& F = aCubeSide[i]._quad->face;
+    StdMeshers_FaceSide* baseQuadSide = aCubeSide[i]._quad->side[ Q_BOTTOM ];
+    list<TopoDS_Edge> baseEdges( baseQuadSide->Edges().begin(), baseQuadSide->Edges().end() );
 
-    if ( helper.GetSubShapeOri( aCubeSide[i]._sideF.Oriented( TopAbs_FORWARD ),
-                                aCubeSide[i]._baseE ) == TopAbs_REVERSED )
-      reverse = !reverse;
+    // assure correctness of node positions on baseE:
+    // helper.GetNodeU() will fix positions if they are wrong
+    for ( int iE = 0; iE < baseQuadSide->NbEdges(); ++iE )
+    {
+      const TopoDS_Edge& baseE = baseQuadSide->Edge( iE );
+      if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( baseE ))
+      {
+        bool ok;
+        helper.SetSubShape( baseE );
+        SMDS_ElemIteratorPtr eIt = smDS->GetElements();
+        while ( eIt->more() )
+        {
+          const SMDS_MeshElement* e = eIt->next();
+          // expect problems on a composite side
+          try { helper.GetNodeU( baseE, e->GetNode(0), e->GetNode(1), &ok); }
+          catch (...) {}
+          try { helper.GetNodeU( baseE, e->GetNode(1), e->GetNode(0), &ok); }
+          catch (...) {}
+        }
+      }
+    }
 
-    if ( sideIndex[i] == I_BOTTOM ||
-         sideIndex[i] == I_LEFT   ||
-         sideIndex[i] == I_BACK )
-      reverse = !reverse;
+    // load grid
+    bool ok =
+      helper.LoadNodeColumns( aCubeSide[i]._u2nodesMap, F, baseEdges, meshDS, proxymesh.get());
+    if ( ok )
+    {
+      // check if the loaded grid corresponds to nb of quadrangles on the FACE
+      const SMESHDS_SubMesh* faceSubMesh =
+        proxymesh ? proxymesh->GetSubMesh( F ) : meshDS->MeshElements( F );
+      const int nbQuads = faceSubMesh->NbElements();
+      const int nbHor = aCubeSide[i]._u2nodesMap.size() - 1;
+      const int nbVer = aCubeSide[i]._u2nodesMap.begin()->second.size() - 1;
+      ok = ( nbQuads == nbHor * nbVer );
+    }
+    if ( !ok )
+    {
+      SMESH_ComputeErrorPtr err = ComputePentahedralMesh(aMesh, aShape, proxymesh.get());
+      return error( err );
+    }
+  }
 
+  // Orient loaded grids of cube sides along axis of the unitary cube coord system
+  bool isReverse[6];
+  isReverse[B_BOTTOM] = beginsAtSide( aCubeSide[B_BOTTOM], aCubeSide[B_RIGHT ] );
+  isReverse[B_TOP   ] = beginsAtSide( aCubeSide[B_TOP   ], aCubeSide[B_RIGHT ] );
+  isReverse[B_FRONT ] = beginsAtSide( aCubeSide[B_FRONT ], aCubeSide[B_RIGHT ] );
+  isReverse[B_BACK  ] = beginsAtSide( aCubeSide[B_BACK  ], aCubeSide[B_RIGHT ] );
+  isReverse[B_LEFT  ] = beginsAtSide( aCubeSide[B_LEFT  ], aCubeSide[B_BACK  ] );
+  isReverse[B_RIGHT ] = beginsAtSide( aCubeSide[B_RIGHT ], aCubeSide[B_BACK  ] );
+  for ( int i = 0; i < 6; ++i )
+  {
     aCubeSide[i]._columns.resize( aCubeSide[i]._u2nodesMap.size() );
 
     int iFwd = 0, iRev = aCubeSide[i]._columns.size()-1;
-    int* pi = reverse ? &iRev : &iFwd;
+    int* pi = isReverse[i] ? &iRev : &iFwd;
     TParam2ColumnMap::iterator u2nn = aCubeSide[i]._u2nodesMap.begin();
     for ( ; iFwd < aCubeSide[i]._columns.size(); --iRev, ++iFwd, ++u2nn )
       aCubeSide[i]._columns[ *pi ].swap( u2nn->second );