Salome HOME
0021422: EDF 1963 SMESH: Viscous layer algorithm fails in some cases
[modules/smesh.git] / src / StdMeshers / StdMeshers_Prism_3D.cxx
index ef99f758657ec4e757d55f45b97d4ad364254635..151926415fb5de5b57d312af62301607b0c192f0 100644 (file)
@@ -1,24 +1,25 @@
-//  Copyright (C) 2007-2008  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
 // File      : StdMeshers_Prism_3D.cxx
 // Module    : SMESH
@@ -213,6 +214,94 @@ namespace {
 
     return gp_Ax2( O, Z, X);
   }
+
+  //================================================================================
+  /*!
+   * \brief Removes submeshes meshed with regular grid from given list
+   *  \retval int - nb of removed submeshes
+   */
+  //================================================================================
+
+  int removeQuasiQuads(list< SMESH_subMesh* >& notQuadSubMesh)
+  {
+    int oldNbSM = notQuadSubMesh.size();
+    SMESHDS_Mesh* mesh = notQuadSubMesh.front()->GetFather()->GetMeshDS();
+    list< SMESH_subMesh* >::iterator smIt = notQuadSubMesh.begin();
+#define __NEXT_SM { ++smIt; continue; }
+    while ( smIt != notQuadSubMesh.end() )
+    {
+      SMESH_subMesh* faceSm = *smIt;
+      SMESHDS_SubMesh* faceSmDS = faceSm->GetSubMeshDS();
+      int nbQuads = faceSmDS->NbElements();
+      if ( nbQuads == 0 ) __NEXT_SM;
+
+      // get oredered edges
+      list< TopoDS_Edge > orderedEdges;
+      list< int >         nbEdgesInWires;
+      TopoDS_Vertex       V000;
+      int nbWires = SMESH_Block::GetOrderedEdges( TopoDS::Face( faceSm->GetSubShape() ),
+                                                  V000, orderedEdges, nbEdgesInWires );
+      if ( nbWires != 1 || nbEdgesInWires.front() <= 4 )
+        __NEXT_SM;
+
+      // get nb of segements on edges
+      list<int> nbSegOnEdge;
+      list< TopoDS_Edge >::iterator edge = orderedEdges.begin();
+      for ( ; edge != orderedEdges.end(); ++edge )
+      {
+        if ( SMESHDS_SubMesh* edgeSmDS = mesh->MeshElements( *edge ))
+          nbSegOnEdge.push_back( edgeSmDS->NbElements() );
+        else
+          nbSegOnEdge.push_back(0);
+      }
+
+      // unite nbSegOnEdge of continues edges
+      int nbEdges = nbEdgesInWires.front();
+      list<int>::iterator nbSegIt = nbSegOnEdge.begin();
+      for ( edge = orderedEdges.begin(); edge != orderedEdges.end(); )
+      {
+        const TopoDS_Edge& e1 = *edge++;
+        const TopoDS_Edge& e2 = ( edge == orderedEdges.end() ? orderedEdges.front() : *edge );
+        if ( SMESH_Algo::IsContinuous( e1, e2 ))
+        {
+          // common vertex of continues edges must be shared by two 2D mesh elems of geom face
+          TopoDS_Vertex vCommon = TopExp::LastVertex( e1, true );
+          const SMDS_MeshNode* vNode = SMESH_Algo::VertexNode( vCommon, mesh );
+          int nbF = 0;
+          if ( vNode )
+          {
+            SMDS_ElemIteratorPtr fIt = vNode->GetInverseElementIterator(SMDSAbs_Face);
+            while ( fIt->more() )
+              nbF += faceSmDS->Contains( fIt->next() );
+          }
+          list<int>::iterator nbSegIt1 = nbSegIt++;
+          if ( !vNode || nbF == 2 ) // !vNode - two edges can be meshed as one
+          {
+            // unite
+            if ( nbSegIt == nbSegOnEdge.end() ) nbSegIt = nbSegOnEdge.begin();
+            *nbSegIt += *nbSegIt1;
+            nbSegOnEdge.erase( nbSegIt1 );
+            --nbEdges;
+          }
+        }
+        else
+        {
+          ++nbSegIt;
+        }
+      }
+      vector<int> nbSegVec( nbSegOnEdge.begin(), nbSegOnEdge.end());
+      if ( nbSegVec.size() == 4 &&
+           nbSegVec[0] == nbSegVec[2] &&
+           nbSegVec[1] == nbSegVec[3] &&
+           nbSegVec[0] * nbSegVec[1] == nbQuads
+           )
+        smIt = notQuadSubMesh.erase( smIt );
+      else
+        __NEXT_SM;
+    }
+
+    return oldNbSM - notQuadSubMesh.size();
+  }
 }
 
 //=======================================================================
@@ -319,7 +408,7 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh
   // Projections on the top and bottom faces are taken from nodes existing
   // on these faces; find correspondence between bottom and top nodes
   myBotToColumnMap.clear();
-  if ( !assocOrProjBottom2Top() ) // it also fill myBotToColumnMap
+  if ( !assocOrProjBottom2Top() ) // it also fills myBotToColumnMap
     return false;
 
 
@@ -389,7 +478,7 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh
         if ( !myBlock.ComputeParameters( topCoords, topParams, ID_TOP_FACE, topParams ))
           return error(TCom("Can't compute normalized parameters ")
                        << "for node " << column.back()->GetID()
-                       << " on the face #"<< column.back()->GetPosition()->GetShapeId() );
+                       << " on the face #"<< column.back()->getshapeId() );
       }
 
       // vertical loop
@@ -466,6 +555,10 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh
     AddPrisms( columns, myHelper );
 
   } // loop on bottom mesh faces
+
+  // clear data
+  myBotToColumnMap.clear();
+  myBlock.Clear();
         
   return true;
 }
@@ -717,6 +810,8 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top()
        botSMDS->NbElements() != topSMDS->NbElements() ||
        botSMDS->NbNodes()    != topSMDS->NbNodes())
   {
+    MESSAGE("nb elem bot " << botSMDS->NbElements() << " top " << topSMDS->NbElements());
+    MESSAGE("nb node bot " << botSMDS->NbNodes() << " top " << topSMDS->NbNodes());
     if ( myBlock.HasNotQuadElemOnTop() )
       return error(TCom("Mesh on faces #") << botSM->GetId()
                    <<" and #"<< topSM->GetId() << " seems different" );
@@ -762,19 +857,8 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top()
     const SMDS_MeshNode* topNode = bN_tN->second;
     if ( botNode->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE )
       continue; // wall columns are contained in myBlock
-    // compute bottom node params
-    TNode bN( botNode );
-//     if ( zSize > 2 ) {
-//       gp_XYZ paramHint(-1,-1,-1);
-//       if ( prevTNode.IsNeighbor( bN ))
-//         paramHint = prevTNode.GetParams();
-//       if ( !myBlock.ComputeParameters( bN.GetCoords(), bN.ChangeParams(),
-//                                        ID_BOT_FACE, paramHint ))
-//         return error(TCom("Can't compute normalized parameters for node ")
-//                      << botNode->GetID() << " on the face #"<< botSM->GetId() );
-//       prevTNode = bN;
-//     }
     // create node column
+    TNode bN( botNode );
     TNode2ColumnMap::iterator bN_col = 
       myBotToColumnMap.insert( make_pair ( bN, TNodeColumn() )).first;
     TNodeColumn & column = bN_col->second;
@@ -965,9 +1049,19 @@ StdMeshers_PrismAsBlock::StdMeshers_PrismAsBlock()
 
 StdMeshers_PrismAsBlock::~StdMeshers_PrismAsBlock()
 {
+  Clear();
+}
+void StdMeshers_PrismAsBlock::Clear()
+{
+  myHelper = 0;
+  myShapeIDMap.Clear();
+  myError.reset();
+
   if ( mySide ) {
     delete mySide; mySide = 0;
   }
+  myParam2ColumnMaps.clear();
+  myShapeIndex2ColumnMap.clear();
 }
 
 //================================================================================
@@ -1077,7 +1171,7 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper,
   }
 
   // ----------------------------------------------------------------------
-  // Analyse faces mesh and topology: choose the bottom submesh.
+  // Analyse mesh and topology of faces: choose the bottom submesh.
   // If there are not quadrangle geom faces, they are top and bottom ones.
   // Not quadrangle geom faces must be only on top and bottom.
   // ----------------------------------------------------------------------
@@ -1090,14 +1184,24 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper,
   bool hasNotQuad = ( nbNotQuad || nbNotQuadMeshed );
 
   // detect bad cases
-  if ( nbNotQuad > 0 && nbNotQuad != 2 )
-    return error(COMPERR_BAD_SHAPE,
-                 TCom("More than 2 not quadrilateral faces: ")
-                 <<nbNotQuad);
   if ( nbNotQuadMeshed > 2 )
+  {
     return error(COMPERR_BAD_INPUT_MESH,
                  TCom("More than 2 faces with not quadrangle elements: ")
                  <<nbNotQuadMeshed);
+  }
+  int nbQuasiQuads = 0;
+  if ( nbNotQuad > 0 && nbNotQuad != 2 )
+  {
+    // Issue 0020843 - one of side faces is quasi-quadrilateral.
+    // Remove from notQuadGeomSubMesh faces meshed with regular grid
+    nbQuasiQuads = removeQuasiQuads( notQuadGeomSubMesh );
+    nbNotQuad -= nbQuasiQuads;
+    if ( nbNotQuad > 0 && nbNotQuad != 2 )
+      return error(COMPERR_BAD_SHAPE,
+                   TCom("More than 2 not quadrilateral faces: ")
+                   <<nbNotQuad);
+  }
 
   // get found submeshes
   if ( hasNotQuad )
@@ -1120,6 +1224,7 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper,
   }
 
   myNotQuadOnTop = ( nbNotQuadMeshed > 1 );
+  MESSAGE("myNotQuadOnTop " << myNotQuadOnTop << " nbNotQuadMeshed " << nbNotQuadMeshed);
  
   // ----------------------------------------------------------
 
@@ -1213,24 +1318,47 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper,
 
   // Get ordered bottom edges
   list< TopoDS_Edge > orderedEdges;
-  list< int >         nbVertexInWires;
+  list< int >         nbEInW;
   SMESH_Block::GetOrderedEdges( TopoDS::Face( botSM->GetSubShape().Reversed() ),
-                                V000, orderedEdges, nbVertexInWires );
-//   if ( nbVertexInWires.size() != 1 )
+                                V000, orderedEdges, nbEInW );
+//   if ( nbEInW.size() != 1 )
 //     RETURN_BAD_RESULT("Wrong prism geometry");
 
   // Get Wall faces corresponding to the ordered bottom edges
   list< TopoDS_Face > wallFaces;
-  if ( !GetWallFaces( Mesh(), shape3D, botSM->GetSubShape(), orderedEdges, wallFaces))
+  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 )
+    {
+      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
   // --------------------------------------------------------
 
   myParam2ColumnMaps.clear();
   myParam2ColumnMaps.resize( orderedEdges.size() ); // total nb edges
 
-  int iE, nbEdges = nbVertexInWires.front(); // nb outer edges
+  int iE, nbEdges = nbEInW.front(); // nb outer edges
   vector< double > edgeLength( nbEdges );
   map< double, int > len2edgeMap;
 
@@ -1276,11 +1404,11 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper,
     // columns for vertices
     // 1
     const SMDS_MeshNode* n0 = faceColumns.begin()->second.front();
-    id = n0->GetPosition()->GetShapeId();
+    id = n0->getshapeId();
     myShapeIndex2ColumnMap[ id ] = make_pair( & faceColumns, isForward );
     // 2
     const SMDS_MeshNode* n1 = faceColumns.rbegin()->second.front();
-    id = n1->GetPosition()->GetShapeId();
+    id = n1->getshapeId();
     myShapeIndex2ColumnMap[ id ] = make_pair( & faceColumns, isForward );
 //     SHOWYXZ("\np1 F "<<iE, gpXYZ(faceColumns.begin()->second.front() ));
 //     SHOWYXZ("p2 F "<<iE, gpXYZ(faceColumns.rbegin()->second.front() ));
@@ -1490,11 +1618,11 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper,
 
       // columns for vertices
       const SMDS_MeshNode* n0 = cols->begin()->second.front();
-      id = n0->GetPosition()->GetShapeId();
+      id = n0->getshapeId();
       myShapeIndex2ColumnMap[ id ] = make_pair( cols, isForward );
 
       const SMDS_MeshNode* n1 = cols->rbegin()->second.front();
-      id = n1->GetPosition()->GetShapeId();
+      id = n1->getshapeId();
       myShapeIndex2ColumnMap[ id ] = make_pair( cols, !isForward );
     }
   }
@@ -1521,7 +1649,7 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper,
 
 const TNodeColumn* StdMeshers_PrismAsBlock::GetNodeColumn(const SMDS_MeshNode* node) const
 {
-  int sID = node->GetPosition()->GetShapeId();
+  int sID = node->getshapeId();
 
   map<int, pair< TParam2ColumnMap*, bool > >::const_iterator col_frw =
     myShapeIndex2ColumnMap.find( sID );
@@ -1560,16 +1688,18 @@ bool StdMeshers_PrismAsBlock::GetLayersTransformation(vector<gp_Trsf> & trsf) co
     list< TopoDS_Edge >::iterator edgeIt = orderedEdges.begin();
     for ( int iE = 0; iE < nbEdgesInWires.front(); ++iE, ++edgeIt )
     {
-      const TParam2ColumnMap& u2colMap =
+      if ( BRep_Tool::Degenerated( *edgeIt )) continue;
+      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 );
       }
@@ -1583,7 +1713,7 @@ bool StdMeshers_PrismAsBlock::GetLayersTransformation(vector<gp_Trsf> & trsf) co
     Bnd_B3d bndBox;
     for ( int i = 0; i < columns.size(); ++i )
       bndBox.Add( gpXYZ( columns[i]->front() ));
-    tol2 = bndBox.SquareExtent() * 4 * 1e-4;
+    tol2 = bndBox.SquareExtent() * 1e-5;
   }
 
   // Compute transformations
@@ -1633,7 +1763,7 @@ bool StdMeshers_PrismAsBlock::IsForwardEdge(SMESHDS_Mesh*           meshDS,
                                             const int               sideFaceID)
 {
   bool isForward = false;
-  if ( TAssocTool::IsClosedEdge( bottomEdge ))
+  if ( SMESH_MesherHelper::IsClosedEdge( bottomEdge ))
   {
     isForward = ( bottomEdge.Orientation() == TopAbs_FORWARD );
   }
@@ -1661,31 +1791,49 @@ bool StdMeshers_PrismAsBlock::IsForwardEdge(SMESHDS_Mesh*           meshDS,
  */
 //================================================================================
 
-bool StdMeshers_PrismAsBlock::GetWallFaces( SMESH_Mesh*                     mesh,
-                                            const TopoDS_Shape &            mainShape,
-                                            const TopoDS_Shape &            bottomFace,
-                                            const std::list< TopoDS_Edge >& bottomEdges,
-                                            std::list< TopoDS_Face >&       wallFaces)
+bool StdMeshers_PrismAsBlock::GetWallFaces( SMESH_Mesh*               mesh,
+                                            const TopoDS_Shape &      mainShape,
+                                            const TopoDS_Shape &      bottomFace,
+                                            std::list< TopoDS_Edge >& bottomEdges,
+                                            std::list< int > &        nbEInW,
+                                            std::list< TopoDS_Face >& wallFaces)
 {
   wallFaces.clear();
 
   TopTools_IndexedMapOfShape faceMap;
   TopExp::MapShapes( mainShape, TopAbs_FACE, faceMap );
 
-  list< TopoDS_Edge >::const_iterator edge = bottomEdges.begin();
-  for ( ; edge != bottomEdges.end(); ++edge )
+  list< TopoDS_Edge >::iterator edge = bottomEdges.begin();
+  std::list< int >::iterator nbE = nbEInW.begin();
+  int iE = 0;
+  while ( edge != bottomEdges.end() )
   {
-    TopTools_ListIteratorOfListOfShape ancestIt = mesh->GetAncestors( *edge );
-    for ( ; ancestIt.More(); ancestIt.Next() )
+    ++iE;
+    if ( BRep_Tool::Degenerated( *edge ))
     {
-      const TopoDS_Shape& ancestor = ancestIt.Value();
-      if ( ancestor.ShapeType() == TopAbs_FACE && // face
-           !bottomFace.IsSame( ancestor ) &&      // not bottom
-           faceMap.FindIndex( ancestor ))         // belongs to the prism
+      edge = bottomEdges.erase( edge );
+      --iE;
+      --(*nbE);
+    }
+    else
+    {
+      PShapeIteratorPtr fIt = myHelper->GetAncestors( *edge, *mesh, TopAbs_FACE );
+      while ( fIt->more() )
       {
-        wallFaces.push_back( TopoDS::Face( ancestor ));
-        break;
+        const TopoDS_Shape* face = fIt->next();
+        if ( !bottomFace.IsSame( *face ) &&      // not bottom
+             faceMap.FindIndex( *face ))         // belongs to the prism
+        {
+          wallFaces.push_back( TopoDS::Face( *face ));
+          break;
+        }
       }
+      ++edge;
+    }
+    if ( iE == *nbE )
+    {
+      iE = 0;
+      ++nbE;
     }
   }
   return ( wallFaces.size() == bottomEdges.size() );
@@ -2054,16 +2202,9 @@ TopoDS_Edge StdMeshers_PrismAsBlock::TSideFace::GetEdge(const int iEdge) const
   TopoDS_Shape V2 = myHelper->GetSubShapeByNode( node, meshDS );
   if ( 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();
 }