Salome HOME
23269: Meshing a composite block with IJK
[modules/smesh.git] / src / StdMeshers / StdMeshers_Prism_3D.cxx
index 1485332b3a83ff81669513bfe1d46dc3f9620669..d60f1397e9d3ee62ab209b2c9d57cd49b1700628 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2015  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2016  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
@@ -62,6 +62,7 @@
 #include <gp_Ax3.hxx>
 
 #include <limits>
+#include <numeric>
 
 using namespace std;
 
@@ -977,7 +978,7 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
     {
       iE = 0;
       ++nbE;
-      int nbQuadPrev = nbQuadsPerWire.empty() ? 0 : nbQuadsPerWire.back();
+      int nbQuadPrev = std::accumulate( nbQuadsPerWire.begin(), nbQuadsPerWire.end(), 0 );
       nbQuadsPerWire.push_back( thePrism.myWallQuads.size() - nbQuadPrev );
     }
   }
@@ -1115,7 +1116,7 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
       ( ! botSM->GetAlgo() ||
         ! _gen->Compute( *botSM->GetFather(), botSM->GetSubShape(), /*shapeOnly=*/true )))
     return error( COMPERR_BAD_INPUT_MESH,
-                  TCom( "No mesher defined to compute the face #")
+                  TCom( "No mesher defined to compute the base face #")
                   << shapeID( thePrism.myBottom ));
 
   // Make all side FACEs of thePrism meshed with quads
@@ -1291,6 +1292,7 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
   if ( !smDS ) return toSM( error(COMPERR_BAD_INPUT_MESH, "Null submesh"));
 
   // loop on bottom mesh faces
+  vector< const TNodeColumn* > columns;
   SMDS_ElemIteratorPtr faceIt = smDS->GetElements();
   while ( faceIt->more() )
   {
@@ -1300,7 +1302,7 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
 
     // find node columns for each node
     int nbNodes = face->NbCornerNodes();
-    vector< const TNodeColumn* > columns( nbNodes );
+    columns.resize( nbNodes );
     for ( int i = 0; i < nbNodes; ++i )
     {
       const SMDS_MeshNode* n = face->GetNode( i );
@@ -1317,7 +1319,8 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
       }
     }
     // create prisms
-    AddPrisms( columns, myHelper );
+    if ( !AddPrisms( columns, myHelper ))
+      return toSM( error("Different 'vertical' discretization"));
 
   } // loop on bottom mesh faces
 
@@ -1327,7 +1330,7 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
 
   // update state of sub-meshes (mostly in order to erase improper errors)
   SMESH_subMesh* sm = myHelper->GetMesh()->GetSubMesh( thePrism.myShape3D );
-  SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/false);
+  SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
   while ( smIt->more() )
   {
     sm = smIt->next();
@@ -1815,17 +1818,20 @@ bool StdMeshers_Prism_3D::Evaluate(SMESH_Mesh&         theMesh,
  */
 //================================================================================
 
-void StdMeshers_Prism_3D::AddPrisms( vector<const TNodeColumn*> & columns,
+bool StdMeshers_Prism_3D::AddPrisms( vector<const TNodeColumn*> & columns,
                                      SMESH_MesherHelper*          helper)
 {
-  int nbNodes = columns.size();
-  int nbZ     = columns[0]->size();
-  if ( nbZ < 2 ) return;
+  size_t nbNodes = columns.size();
+  size_t nbZ     = columns[0]->size();
+  if ( nbZ < 2 ) return false;
+  for ( size_t i = 1; i < nbNodes; ++i )
+    if ( columns[i]->size() != nbZ )
+      return false;
 
   // find out orientation
   bool isForward = true;
   SMDS_VolumeTool vTool;
-  int z = 1;
+  size_t z = 1;
   switch ( nbNodes ) {
   case 3: {
     SMDS_VolumeOfNodes tmpPenta ( (*columns[0])[z-1], // bottom
@@ -1911,7 +1917,7 @@ void StdMeshers_Prism_3D::AddPrisms( vector<const TNodeColumn*> & columns,
     vector<const SMDS_MeshNode*> nodes( 2*nbNodes + 4*nbNodes);
     for ( z = 1; z < nbZ; ++z )
     {
-      for ( int i = 0; i < nbNodes; ++i ) {
+      for ( size_t i = 0; i < nbNodes; ++i ) {
         nodes[ i             ] = (*columns[ i ])[z+iBase1]; // bottom or top
         nodes[ 2*nbNodes-i-1 ] = (*columns[ i ])[z+iBase2]; // top or bottom
         // side
@@ -1925,6 +1931,8 @@ void StdMeshers_Prism_3D::AddPrisms( vector<const TNodeColumn*> & columns,
     }
 
   } // switch ( nbNodes )
+
+  return true;
 }
 
 //================================================================================
@@ -2501,20 +2509,24 @@ namespace // utils used by StdMeshers_Prism_3D::IsApplicable()
       _iL( SMESH_MesherHelper::WrapIndex( iE-1, nbE ) + shift ),
       _iR( SMESH_MesherHelper::WrapIndex( iE+1, nbE ) + shift )
     {
+      //_edge.Orientation( TopAbs_FORWARD ); // for operator==() to work
     }
     EdgeWithNeighbors() {}
+    bool IsInternal() const { return !_edge.IsNull() && _edge.Orientation() == TopAbs_INTERNAL; }
   };
-  struct PrismSide
+  // PrismSide contains all FACEs linking a bottom EDGE with a top one. 
+  struct PrismSide 
   {
-    TopoDS_Face                 _face;
-    TopTools_IndexedMapOfShape *_faces; // pointer because its copy constructor is private
-    TopoDS_Edge                 _topEdge;
-    vector< EdgeWithNeighbors >*_edges;
-    int                         _iBotEdge;
-    vector< bool >              _isCheckedEdge;
+    TopoDS_Face                 _face;    // a currently treated upper FACE
+    TopTools_IndexedMapOfShape *_faces;   // all FACEs (pointer because of a private copy constructor)
+    TopoDS_Edge                 _topEdge; // a current top EDGE
+    vector< EdgeWithNeighbors >*_edges;   // all EDGEs of _face
+    int                         _iBotEdge;       // index of _topEdge within _edges
+    vector< bool >              _isCheckedEdge;  // mark EDGEs whose two owner FACEs found
     int                         _nbCheckedEdges; // nb of EDGEs whose location is defined
-    PrismSide                  *_leftSide;
+    PrismSide                  *_leftSide;       // neighbor sides
     PrismSide                  *_rightSide;
+    bool                        _isInternal; // whether this side raises from an INTERNAL EDGE
     void SetExcluded() { _leftSide = _rightSide = NULL; }
     bool IsExcluded() const { return !_leftSide; }
     const TopoDS_Edge& Edge( int i ) const
@@ -2527,10 +2539,15 @@ namespace // utils used by StdMeshers_Prism_3D::IsApplicable()
         if ( E.IsSame( Edge( i ))) return i;
       return -1;
     }
-    bool IsSideFace( const TopoDS_Shape& face ) const
+    bool IsSideFace( const TopoDS_Shape& face, const bool checkNeighbors ) const
     {
       if ( _faces->Contains( face )) // avoid returning true for a prism top FACE
         return ( !_face.IsNull() || !( face.IsSame( _faces->FindKey( _faces->Extent() ))));
+
+      if ( checkNeighbors )
+        return (( _leftSide  && _leftSide->IsSideFace ( face, false )) ||
+                ( _rightSide && _rightSide->IsSideFace( face, false )));
+
       return false;
     }
   };
@@ -2542,15 +2559,19 @@ namespace // utils used by StdMeshers_Prism_3D::IsApplicable()
                  vector< EdgeWithNeighbors > & edges,
                  const bool                    noHolesAllowed)
   {
+    TopoDS_Face f = face;
+    if ( f.Orientation() != TopAbs_FORWARD &&
+         f.Orientation() != TopAbs_REVERSED )
+      f.Orientation( TopAbs_FORWARD );
     list< TopoDS_Edge > ee;
     list< int >         nbEdgesInWires;
-    int nbW = SMESH_Block::GetOrderedEdges( face, ee, nbEdgesInWires );
+    int nbW = SMESH_Block::GetOrderedEdges( f, ee, nbEdgesInWires );
     if ( nbW > 1 && noHolesAllowed )
       return false;
 
     int iE, nbTot = 0;
-    list< TopoDS_Edge >::iterator e = ee.begin();
-    list< int >::iterator       nbE = nbEdgesInWires.begin();
+    list< TopoDS_Edge >::iterator   e = ee.begin();
+    list< int         >::iterator nbE = nbEdgesInWires.begin();
     for ( ; nbE != nbEdgesInWires.end(); ++nbE )
       for ( iE = 0; iE < *nbE; ++e, ++iE )
         if ( SMESH_Algo::isDegenerated( *e ))
@@ -2559,10 +2580,6 @@ namespace // utils used by StdMeshers_Prism_3D::IsApplicable()
           --(*nbE);
           --iE;
         }
-        else
-        {
-          e->Orientation( TopAbs_FORWARD ); // for operator==() to work
-        }
 
     edges.clear();
     e = ee.begin();
@@ -2572,6 +2589,35 @@ namespace // utils used by StdMeshers_Prism_3D::IsApplicable()
         edges.push_back( EdgeWithNeighbors( *e, iE, *nbE, nbTot ));
       nbTot += *nbE;
     }
+
+    // IPAL53099. Set correct neighbors to INTERNAL EDGEs, which can be connected to
+    // EDGEs of the outer WIRE but this fact can't be detected by their order.
+    if ( nbW > 1 )
+    {
+      int iFirst = 0, iLast;
+      for ( nbE = nbEdgesInWires.begin(); nbE != nbEdgesInWires.end(); ++nbE )
+      {
+        iLast = iFirst + *nbE - 1;
+        TopoDS_Vertex vv[2] = { SMESH_MesherHelper::IthVertex( 0, edges[ iFirst ]._edge ),
+                                SMESH_MesherHelper::IthVertex( 1, edges[ iLast  ]._edge ) };
+        bool isConnectOk = ( vv[0].IsSame( vv[1] ));
+        if ( !isConnectOk )
+        {
+          // look for an EDGE of the outer WIRE connected to vv
+          TopoDS_Vertex v0, v1;
+          for ( iE = 0; iE < nbEdgesInWires.front(); ++iE )
+          {
+            v0 = SMESH_MesherHelper::IthVertex( 0, edges[ iE ]._edge );
+            v1 = SMESH_MesherHelper::IthVertex( 1, edges[ iE ]._edge );
+            if ( vv[0].IsSame( v0 ) || vv[0].IsSame( v1 ))
+              edges[ iFirst ]._iL = iE;
+            if ( vv[1].IsSame( v0 ) || vv[1].IsSame( v1 ))
+              edges[ iLast ]._iR = iE;
+          }
+        }
+        iFirst += *nbE;
+      }
+    }
     return edges.size();
   }
   //--------------------------------------------------------------------------------
@@ -2671,7 +2717,7 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
 
     typedef vector< EdgeWithNeighbors > TEdgeWithNeighborsVec;
     vector< TEdgeWithNeighborsVec > faceEdgesVec( allFaces.Extent() + 1 );
-    const size_t nbEdgesMax = facesOfEdge.Extent() * 2; // there can be seam EDGES
+    const size_t nbEdgesMax = facesOfEdge.Extent() * 2; // there can be seam EDGEs
     TopTools_IndexedMapOfShape* facesOfSide = new TopTools_IndexedMapOfShape[ nbEdgesMax ];
     SMESHUtils::ArrayDeleter<TopTools_IndexedMapOfShape> delFacesOfSide( facesOfSide );
 
@@ -2694,11 +2740,12 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
       sides.resize( botEdges.size() );
       for ( size_t iS = 0; iS < botEdges.size(); ++iS )
       {
-        sides[ iS ]._topEdge   = botEdges[ iS ]._edge;
-        sides[ iS ]._face      = botF;
-        sides[ iS ]._leftSide  = & sides[ botEdges[ iS ]._iR ];
-        sides[ iS ]._rightSide = & sides[ botEdges[ iS ]._iL ];
-        sides[ iS ]._faces = & facesOfSide[ iS ];
+        sides[ iS ]._topEdge    = botEdges[ iS ]._edge;
+        sides[ iS ]._face       = botF;
+        sides[ iS ]._leftSide   = & sides[ botEdges[ iS ]._iR ];
+        sides[ iS ]._rightSide  = & sides[ botEdges[ iS ]._iL ];
+        sides[ iS ]._isInternal = botEdges[ iS ].IsInternal();
+        sides[ iS ]._faces      = & facesOfSide[ iS ];
         sides[ iS ]._faces->Clear();
       }
 
@@ -2727,8 +2774,9 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
                 if ( side._isCheckedEdge[ iE ] ) continue;
                 const TopoDS_Edge&      vertE = side.Edge( iE );
                 const TopoDS_Shape& neighborF = getAnotherFace( side._face, vertE, facesOfEdge );
-                bool             isEdgeShared = adjSide->IsSideFace( neighborF );
-                if ( isEdgeShared )               // vertE is shared with adjSide
+                bool isEdgeShared = (( adjSide->IsSideFace( neighborF, side._isInternal )) ||
+                                     ( adjSide == &side && neighborF.IsSame( side._face )) );
+                if ( isEdgeShared ) // vertE is shared with adjSide
                 {
                   isAdvanced = true;
                   side._isCheckedEdge[ iE ] = true;
@@ -2769,20 +2817,19 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
             {
               stop = true;
             }
-            else if ( side._leftSide != & side ) // not closed side face
+            else if ( side._leftSide != & side && // not closed side face
+                      side._leftSide->_faces->Contains( f ))
             {
-              if ( side._leftSide->_faces->Contains( f ))
-              {
-                stop = true; // probably f is the prism top face
-                side._leftSide->_face.Nullify();
-                side._leftSide->_topEdge.Nullify();
-              }
-              if ( side._rightSide->_faces->Contains( f ))
-              {
-                stop = true; // probably f is the prism top face
-                side._rightSide->_face.Nullify();
-                side._rightSide->_topEdge.Nullify();
-              }
+              stop = true; // probably f is the prism top face
+              side._leftSide->_face.Nullify();
+              side._leftSide->_topEdge.Nullify();
+            }
+            else if ( side._rightSide != & side &&
+                      side._rightSide->_faces->Contains( f ))
+            {
+              stop = true; // probably f is the prism top face
+              side._rightSide->_face.Nullify();
+              side._rightSide->_topEdge.Nullify();
             }
             if ( stop )
             {
@@ -2807,6 +2854,10 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA
               side._isCheckedEdge[ side._iBotEdge ] = true;
               side._nbCheckedEdges = 1; // bottom EDGE is known
             }
+            else // probably a triangular top face found
+            {
+              side._face.Nullify();
+            }
             side._topEdge.Nullify();
             isOK = ( !side._edges->empty() || side._faces->Extent() > 1 );
 
@@ -3247,6 +3298,9 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper*         helper,
       if ( !myHelper->LoadNodeColumns( faceColumns, (*quad)->face, quadBot, meshDS ))
         return error(COMPERR_BAD_INPUT_MESH, TCom("Can't find regular quadrangle mesh ")
                      << "on a side face #" << MeshDS()->ShapeToIndex( (*quad)->face ));
+
+      if ( !faceColumns.empty() && (int)faceColumns.begin()->second.size() != VerticalSize() )
+        return error(COMPERR_BAD_INPUT_MESH, "Different 'vertical' discretization");
     }
     // edge columns
     int id = MeshDS()->ShapeToIndex( *edgeIt );