Salome HOME
IPAL 0051561: Hexotic algorithm fail
[modules/smesh.git] / src / StdMeshers / StdMeshers_Prism_3D.cxx
index 8998a294986d56db759ef0af39a559e399e594e7..85eaef56619c46b5d460884730c24d96960fecd0 100644 (file)
@@ -1,26 +1,25 @@
-//  Copyright (C) 2007-2010  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2012  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
 // File      : StdMeshers_Prism_3D.cxx
 // Module    : SMESH
 // Created   : Fri Oct 20 11:37:07 2006
@@ -385,7 +384,7 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh
 
   myHelper->IsQuadraticSubMesh( theShape );
 
-  // Analyse mesh and geomerty to find block subshapes and submeshes
+  // Analyse mesh and geomerty to find block sub-shapes and submeshes
   if ( !myBlock.Init( myHelper, theShape ))
     return error( myBlock.GetError());
 
@@ -689,9 +688,6 @@ bool StdMeshers_Prism_3D::Evaluate(SMESH_Mesh& theMesh,
 void StdMeshers_Prism_3D::AddPrisms( vector<const TNodeColumn*> & columns,
                                      SMESH_MesherHelper*          helper)
 {
-  SMESHDS_Mesh * meshDS = helper->GetMeshDS();
-  int shapeID = helper->GetSubShapeID();
-
   int nbNodes = columns.size();
   int nbZ     = columns[0]->size();
   if ( nbZ < 2 ) return;
@@ -702,87 +698,103 @@ void StdMeshers_Prism_3D::AddPrisms( vector<const TNodeColumn*> & columns,
   int z = 1;
   switch ( nbNodes ) {
   case 3: {
-    const SMDS_MeshNode* botNodes[3] = { (*columns[0])[z-1],
-                                         (*columns[1])[z-1],
-                                         (*columns[2])[z-1] };
-    const SMDS_MeshNode* topNodes[3] = { (*columns[0])[z],
-                                         (*columns[1])[z],
-                                         (*columns[2])[z] };
-    SMDS_VolumeOfNodes tmpVol ( botNodes[0], botNodes[1], botNodes[2],
-                                topNodes[0], topNodes[1], topNodes[2]);
-    vTool.Set( &tmpVol );
+    SMDS_VolumeOfNodes tmpPenta ( (*columns[0])[z-1], // bottom
+                                  (*columns[1])[z-1],
+                                  (*columns[2])[z-1],
+                                  (*columns[0])[z],   // top
+                                  (*columns[1])[z],
+                                  (*columns[2])[z] );
+    vTool.Set( &tmpPenta );
     isForward  = vTool.IsForward();
     break;
   }
   case 4: {
-    const SMDS_MeshNode* botNodes[4] = { (*columns[0])[z-1], (*columns[1])[z-1],
-                                         (*columns[2])[z-1], (*columns[3])[z-1] };
-    const SMDS_MeshNode* topNodes[4] = { (*columns[0])[z], (*columns[1])[z],
-                                         (*columns[2])[z], (*columns[3])[z] };
-    SMDS_VolumeOfNodes tmpVol ( botNodes[0], botNodes[1], botNodes[2], botNodes[3],
-                                topNodes[0], topNodes[1], topNodes[2], topNodes[3]);
-    vTool.Set( &tmpVol );
+    SMDS_VolumeOfNodes tmpHex( (*columns[0])[z-1], (*columns[1])[z-1], // bottom
+                               (*columns[2])[z-1], (*columns[3])[z-1],
+                               (*columns[0])[z],   (*columns[1])[z],   // top
+                               (*columns[2])[z],   (*columns[3])[z] );
+    vTool.Set( &tmpHex );
     isForward  = vTool.IsForward();
     break;
   }
+  default:
+    const int di = (nbNodes+1) / 3;
+    SMDS_VolumeOfNodes tmpVol ( (*columns[0]   )[z-1],
+                                (*columns[di]  )[z-1],
+                                (*columns[2*di])[z-1],
+                                (*columns[0]   )[z],
+                                (*columns[di]  )[z],
+                                (*columns[2*di])[z] );
+    vTool.Set( &tmpVol );
+    isForward  = vTool.IsForward();
   }
 
   // vertical loop on columns
-  for ( z = 1; z < nbZ; ++z )
-  {
-    SMDS_MeshElement* vol = 0;
-    switch ( nbNodes ) {
 
-    case 3: {
-      const SMDS_MeshNode* botNodes[3] = { (*columns[0])[z-1],
-                                           (*columns[1])[z-1],
-                                           (*columns[2])[z-1] };
-      const SMDS_MeshNode* topNodes[3] = { (*columns[0])[z],
-                                           (*columns[1])[z],
-                                           (*columns[2])[z] };
-      if ( isForward )
-        vol = helper->AddVolume( botNodes[0], botNodes[1], botNodes[2],
-                                 topNodes[0], topNodes[1], topNodes[2]);
-      else
-        vol = helper->AddVolume( topNodes[0], topNodes[1], topNodes[2],
-                                 botNodes[0], botNodes[1], botNodes[2]);
-      break;
-      }
-    case 4: {
-      const SMDS_MeshNode* botNodes[4] = { (*columns[0])[z-1], (*columns[1])[z-1],
-                                           (*columns[2])[z-1], (*columns[3])[z-1] };
-      const SMDS_MeshNode* topNodes[4] = { (*columns[0])[z], (*columns[1])[z],
-                                           (*columns[2])[z], (*columns[3])[z] };
-      if ( isForward )
-        vol = helper->AddVolume( botNodes[0], botNodes[1], botNodes[2], botNodes[3],
-                                 topNodes[0], topNodes[1], topNodes[2], topNodes[3]);
-      else
-        vol = helper->AddVolume( topNodes[0], topNodes[1], topNodes[2], topNodes[3],
-                                 botNodes[0], botNodes[1], botNodes[2], botNodes[3]);
-      break;
-      }
-    default:
-      // polyhedron
-      vector<const SMDS_MeshNode*> nodes( 2*nbNodes + 4*nbNodes);
-      vector<int> quantities( 2 + nbNodes, 4 );
-      quantities[0] = quantities[1] = nbNodes;
-      columns.resize( nbNodes + 1 );
-      columns[ nbNodes ] = columns[ 0 ];
+  helper->SetElementsOnShape( true );
+
+  switch ( nbNodes ) {
+
+  case 3: { // ---------- pentahedra
+    const int i1 = isForward ? 1 : 2;
+    const int i2 = isForward ? 2 : 1;
+    for ( z = 1; z < nbZ; ++z )
+      helper->AddVolume( (*columns[0 ])[z-1], // bottom
+                         (*columns[i1])[z-1],
+                         (*columns[i2])[z-1],
+                         (*columns[0 ])[z],   // top
+                         (*columns[i1])[z],
+                         (*columns[i2])[z] );
+    break;
+  }
+  case 4: { // ---------- hexahedra
+    const int i1 = isForward ? 1 : 3;
+    const int i3 = isForward ? 3 : 1;
+    for ( z = 1; z < nbZ; ++z )
+      helper->AddVolume( (*columns[0])[z-1], (*columns[i1])[z-1], // bottom
+                         (*columns[2])[z-1], (*columns[i3])[z-1],
+                         (*columns[0])[z],   (*columns[i1])[z],     // top
+                         (*columns[2])[z],   (*columns[i3])[z] );
+    break;
+  }
+  case 6: { // ---------- octahedra
+    const int iBase1 = isForward ? -1 : 0;
+    const int iBase2 = isForward ?  0 :-1;
+    for ( z = 1; z < nbZ; ++z )
+      helper->AddVolume( (*columns[0])[z+iBase1], (*columns[1])[z+iBase1], // bottom or top
+                         (*columns[2])[z+iBase1], (*columns[3])[z+iBase1],
+                         (*columns[4])[z+iBase1], (*columns[5])[z+iBase1],
+                         (*columns[0])[z+iBase2], (*columns[1])[z+iBase2], // top or bottom
+                         (*columns[2])[z+iBase2], (*columns[3])[z+iBase2],
+                         (*columns[4])[z+iBase2], (*columns[5])[z+iBase2] );
+    break;
+  }
+  default: // ---------- polyhedra
+    vector<int> quantities( 2 + nbNodes, 4 );
+    quantities[0] = quantities[1] = nbNodes;
+    columns.resize( nbNodes + 1 );
+    columns[ nbNodes ] = columns[ 0 ];
+    const int i1 = isForward ? 1 : 3;
+    const int i3 = isForward ? 3 : 1;
+    const int iBase1 = isForward ? -1 : 0;
+    const int iBase2 = isForward ?  0 :-1;
+    vector<const SMDS_MeshNode*> nodes( 2*nbNodes + 4*nbNodes);
+    for ( z = 1; z < nbZ; ++z )
+    {
       for ( int i = 0; i < nbNodes; ++i ) {
-        nodes[ i         ] = (*columns[ i ])[z-1]; // bottom
-        nodes[ i+nbNodes ] = (*columns[ i ])[z  ]; // top
+        nodes[ i             ] = (*columns[ i ])[z+iBase1]; // bottom or top
+        nodes[ 2*nbNodes-i-1 ] = (*columns[ i ])[z+iBase2]; // top or bottom
         // side
-        int di = 2*nbNodes + 4*i - 1;
-        nodes[ di   ] = (*columns[i  ])[z-1];
-        nodes[ di+1 ] = (*columns[i+1])[z-1];
-        nodes[ di+2 ] = (*columns[i+1])[z  ];
-        nodes[ di+3 ] = (*columns[i  ])[z  ];
+        int di = 2*nbNodes + 4*i;
+        nodes[ di+0 ] = (*columns[i  ])[z  ];
+        nodes[ di+i1] = (*columns[i+1])[z  ];
+        nodes[ di+2 ] = (*columns[i+1])[z-1];
+        nodes[ di+i3] = (*columns[i  ])[z-1];
       }
-      vol = meshDS->AddPolyhedralVolume( nodes, quantities );
+      helper->AddPolyhedralVolume( nodes, quantities );
     }
-    if ( vol && shapeID > 0 )
-      meshDS->SetMeshElementOnShape( vol, shapeID );
-  }
+
+  } // switch ( nbNodes )
 }
 
 //================================================================================
@@ -984,7 +996,7 @@ bool StdMeshers_Prism_3D::projectBottomToTop()
 
 //================================================================================
 /*!
- * \brief Set projection coordinates of a node to a face and it's subshapes
+ * \brief Set projection coordinates of a node to a face and it's sub-shapes
  * \param faceID - the face given by in-block ID
  * \param params - node normalized parameters
  * \retval bool - is a success
@@ -1329,15 +1341,27 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper,
   if ( !GetWallFaces( Mesh(), shape3D, botSM->GetSubShape(), orderedEdges, nbEInW, wallFaces))
     return error(COMPERR_BAD_SHAPE, "Can't find side faces");
 
+  // check that the found top and bottom faces are opposite
+  {
+    for (TopExp_Explorer edge(botSM->GetSubShape(), TopAbs_EDGE); edge.More(); edge.Next())
+      if ( helper->IsSubShape( edge.Current(), topSM->GetSubShape() ))
+        return error(notQuadGeomSubMesh.empty() ? COMPERR_BAD_INPUT_MESH : COMPERR_BAD_SHAPE,
+                     "Non-quadrilateral faces are not opposite");
+  }
+
   // Protect from a distorted block (test 3D_mesh_HEXA3D/B7 on 32bit platform)
   // check that all wall faces have an edge common with the top face
   {
     list< TopoDS_Face >::iterator faceIt = wallFaces.begin();
     for ( ; faceIt != wallFaces.end(); ++faceIt )
-      for (TopExp_Explorer edge(*faceIt, TopAbs_EDGE); edge.More(); edge.Next()) {
-        if ( !helper->IsSubShape( edge.Current(), topSM->GetSubShape() ))
-          return error(COMPERR_BAD_SHAPE);
-      }
+    {
+      bool hasCommon = false;
+      for (TopExp_Explorer edge(*faceIt, TopAbs_EDGE); !hasCommon && edge.More(); edge.Next())
+        if ( helper->IsSubShape( edge.Current(), topSM->GetSubShape() ))
+          hasCommon = true;
+      if ( !hasCommon )
+        return error(COMPERR_BAD_SHAPE);
+    }
   }
 
   // Find columns of wall nodes and calculate edges' lengths
@@ -1448,9 +1472,12 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper,
         int nbSplit = i_nb->second;
         vector< double > params;
         splitParams( nbSplit, &myParam2ColumnMaps[ iE ], params );
-        bool isForward = ( edgeIt->Orientation() == TopAbs_FORWARD );
+        const bool isForward =
+          StdMeshers_PrismAsBlock::IsForwardEdge( myHelper->GetMeshDS(),
+                                                  myParam2ColumnMaps[iE],
+                                                  *edgeIt, SMESH_Block::ID_Fx0z );
         for ( int i = 0; i < nbSplit; ++i ) {
-          double f = ( isForward ? params[ i ] : params[ nbSplit - i-1 ]);
+          double f = ( isForward ? params[ i ]   : params[ nbSplit - i-1 ]);
           double l = ( isForward ? params[ i+1 ] : params[ nbSplit - i ]);
           TSideFace* comp = new TSideFace( myHelper, wallFaceIds[ iSide ],
                                            *faceIt, *edgeIt,
@@ -1677,16 +1704,17 @@ bool StdMeshers_PrismAsBlock::GetLayersTransformation(vector<gp_Trsf> & trsf) co
     for ( int iE = 0; iE < nbEdgesInWires.front(); ++iE, ++edgeIt )
     {
       if ( BRep_Tool::Degenerated( *edgeIt )) continue;
-      const TParam2ColumnMap& u2colMap =
+      const TParam2ColumnMap* u2colMap =
         GetParam2ColumnMap( myHelper->GetMeshDS()->ShapeToIndex( *edgeIt ), isReverse );
+      if ( !u2colMap ) return false;
       isReverse = ( edgeIt->Orientation() == TopAbs_REVERSED );
-      double f = u2colMap.begin()->first, l = u2colMap.rbegin()->first;
+      double f = u2colMap->begin()->first, l = u2colMap->rbegin()->first;
       if ( isReverse ) swap ( f, l );
       const int nbCol = 5;
       for ( int i = 0; i < nbCol; ++i )
       {
         double u = f + i/double(nbCol) * ( l - f );
-        const TNodeColumn* col = & getColumn( u2colMap, u )->second;
+        const TNodeColumn* col = & getColumn( u2colMap, u )->second;
         if ( columns.empty() || col != columns.back() )
           columns.push_back( col );
       }
@@ -1740,7 +1768,7 @@ bool StdMeshers_PrismAsBlock::GetLayersTransformation(vector<gp_Trsf> & trsf) co
   * \param columnsMap - node columns map of side face
   * \param bottomEdge - the bootom edge
   * \param sideFaceID - side face in-block ID
-  * \retval bool - true if orientation coinside with in-block froward orientation
+  * \retval bool - true if orientation coinside with in-block forward orientation
  */
 //================================================================================
 
@@ -2187,27 +2215,20 @@ TopoDS_Edge StdMeshers_PrismAsBlock::TSideFace::GetEdge(const int iEdge) const
   // find edge by 2 vertices
   TopoDS_Shape V1 = edge;
   TopoDS_Shape V2 = myHelper->GetSubShapeByNode( node, meshDS );
-  if ( V2.ShapeType() == TopAbs_VERTEX && !V2.IsSame( V1 ))
+  if ( !V2.IsNull() && V2.ShapeType() == TopAbs_VERTEX && !V2.IsSame( V1 ))
   {
-    TopTools_ListIteratorOfListOfShape ancestIt =
-      myHelper->GetMesh()->GetAncestors( V1 );
-    for ( ; ancestIt.More(); ancestIt.Next() )
-    {
-      const TopoDS_Shape & ancestor = ancestIt.Value();
-      if ( ancestor.ShapeType() == TopAbs_EDGE )
-        for ( TopExp_Explorer e( ancestor, TopAbs_VERTEX ); e.More(); e.Next() )
-          if ( V2.IsSame( e.Current() ))
-            return TopoDS::Edge( ancestor );
-    }
+    TopoDS_Shape ancestor = myHelper->GetCommonAncestor( V1, V2, *myHelper->GetMesh(), TopAbs_EDGE);
+    if ( !ancestor.IsNull() )
+      return TopoDS::Edge( ancestor );
   }
   return TopoDS_Edge();
 }
 
 //================================================================================
 /*!
- * \brief Fill block subshapes
+ * \brief Fill block sub-shapes
   * \param shapeMap - map to fill in
-  * \retval int - nb inserted subshapes
+  * \retval int - nb inserted sub-shapes
  */
 //================================================================================