Salome HOME
Merge from BR_imps_2013 14/01/2014
[modules/smesh.git] / src / StdMeshers / StdMeshers_Prism_3D.cxx
index 048b3e7ad7c6fba5c5f0b4e12f0ff4b5f3957f99..33bd5e1dcfc2397a2370dff82b9401318a59898b 100644 (file)
@@ -107,7 +107,7 @@ namespace {
            algo->myProxyMesh->GetMesh() != helper->GetMesh() )
         algo->myProxyMesh.reset( new SMESH_ProxyMesh( *helper->GetMesh() ));
 
-      algo->myQuadStruct.reset();
+      algo->myQuadList.clear();
 
       if ( helper )
         algo->_quadraticMesh = helper->GetIsQuadratic();
@@ -166,15 +166,15 @@ namespace {
   //================================================================================
 
   bool setBottomEdge( const TopoDS_Edge&   botE,
-                      faceQuadStruct::Ptr& quad,
+                      FaceQuadStruct::Ptr& quad,
                       const TopoDS_Shape&  face)
   {
-    quad->side[ QUAD_TOP_SIDE  ]->Reverse();
-    quad->side[ QUAD_LEFT_SIDE ]->Reverse();
+    quad->side[ QUAD_TOP_SIDE  ].grid->Reverse();
+    quad->side[ QUAD_LEFT_SIDE ].grid->Reverse();
     int edgeIndex = 0;
     for ( size_t i = 0; i < quad->side.size(); ++i )
     {
-      StdMeshers_FaceSide* quadSide = quad->side[i];
+      StdMeshers_FaceSidePtr quadSide = quad->side[i];
       for ( int iE = 0; iE < quadSide->NbEdges(); ++iE )
         if ( botE.IsSame( quadSide->Edge( iE )))
         {
@@ -373,6 +373,23 @@ namespace {
     return nbRemoved;
   }
 
+  //================================================================================
+  /*!
+   * \brief Return and angle between two EDGEs
+   *  \return double - the angle normalized so that
+   * >~ 0  -> 2.0
+   *  PI/2 -> 1.0
+   *  PI   -> 0.0
+   * -PI/2 -> -1.0
+   * <~ 0  -> -2.0
+   */
+  //================================================================================
+
+  double normAngle(const TopoDS_Edge & E1, const TopoDS_Edge & E2, const TopoDS_Face & F)
+  {
+    return SMESH_MesherHelper::GetAngle( E1, E2, F ) / ( 0.5 * M_PI );
+  }
+
   //================================================================================
   /*!
    * Consider continuous straight EDGES as one side - mark them to unite
@@ -380,40 +397,80 @@ namespace {
   //================================================================================
 
   int countNbSides( const Prism_3D::TPrismTopo & thePrism,
-                    vector<int> &                nbUnitePerEdge )
+                    vector<int> &                nbUnitePerEdge,
+                    vector< double > &           edgeLength)
   {
     int nbEdges = thePrism.myNbEdgesInWires.front();  // nb outer edges
     int nbSides = nbEdges;
 
+    
     list< TopoDS_Edge >::const_iterator edgeIt = thePrism.myBottomEdges.begin();
     std::advance( edgeIt, nbEdges-1 );
     TopoDS_Edge   prevE = *edgeIt;
-    bool isPrevStraight = SMESH_Algo::isStraight( prevE );
+    // bool isPrevStraight = SMESH_Algo::IsStraight( prevE );
     int           iPrev = nbEdges - 1;
 
     int iUnite = -1; // the first of united EDGEs
 
+    // analyse angles between EDGEs
+    int nbCorners = 0;
+    vector< bool > isCorner( nbEdges );
     edgeIt = thePrism.myBottomEdges.begin();
     for ( int iE = 0; iE < nbEdges; ++iE, ++edgeIt )
     {
       const TopoDS_Edge&  curE = *edgeIt;
-      const bool isCurStraight = SMESH_Algo::isStraight( curE );
-      if ( isPrevStraight && isCurStraight && SMESH_Algo::IsContinuous( prevE, curE ))
-      {
-        if ( iUnite < 0 )
-          iUnite = iPrev;
-        nbUnitePerEdge[ iUnite ]++;
-        nbUnitePerEdge[ iE ] = -1;
-        --nbSides;
-      }
-      else
-      {
-        iUnite = -1;
-      }
-      prevE          = curE;
-      isPrevStraight = isCurStraight;
-      iPrev = iE;
+      edgeLength[ iE ] = SMESH_Algo::EdgeLength( curE );
+
+      // double normAngle = normAngle( prevE, curE, thePrism.myBottom );
+      // isCorner[ iE ] = false;
+      // if ( normAngle < 2.0 )
+      // {
+      //   if ( normAngle < 0.001 ) // straight or obtuse angle
+      //   {
+      //     // unite EDGEs in order not to put a corner of the unit quadrangle at this VERTEX
+      //     if ( iUnite < 0 )
+      //       iUnite = iPrev;
+      //     nbUnitePerEdge[ iUnite ]++;
+      //     nbUnitePerEdge[ iE ] = -1;
+      //     --nbSides;
+      //   }
+      //   else
+      //   {
+      //     isCorner[ iE ] = true;
+      //     nbCorners++;
+      //     iUnite = -1;
+      //   }
+      // }
+      // prevE = curE;
+    }
+
+    if ( nbCorners > 4 )
+    {
+      // define which of corners to put on a side of the unit quadrangle
     }
+    // edgeIt = thePrism.myBottomEdges.begin();
+    // for ( int iE = 0; iE < nbEdges; ++iE, ++edgeIt )
+    // {
+    //   const TopoDS_Edge&  curE = *edgeIt;
+    //   edgeLength[ iE ] = SMESH_Algo::EdgeLength( curE );
+
+    //   const bool isCurStraight = SMESH_Algo::IsStraight( curE );
+    //   if ( isPrevStraight && isCurStraight && SMESH_Algo::IsContinuous( prevE, curE ))
+    //   {
+    //     if ( iUnite < 0 )
+    //       iUnite = iPrev;
+    //     nbUnitePerEdge[ iUnite ]++;
+    //     nbUnitePerEdge[ iE ] = -1;
+    //     --nbSides;
+    //   }
+    //   else
+    //   {
+    //     iUnite = -1;
+    //   }
+    //   prevE          = curE;
+    //   isPrevStraight = isCurStraight;
+    //   iPrev = iE;
+    // }
     
     return nbSides;
   }
@@ -624,7 +681,7 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh
               continue; // already computed prism
             }
             // find a source FACE of the SOLID: it's a FACE sharing a bottom EDGE with wFace
-            const TopoDS_Edge& wEdge = (*wQuad)->side[ QUAD_TOP_SIDE ]->Edge(0);
+            const TopoDS_Edge& wEdge = (*wQuad)->side[ QUAD_TOP_SIDE ].grid->Edge(0);
             PShapeIteratorPtr faceIt = myHelper->GetAncestors( wEdge, *myHelper->GetMesh(),
                                                                TopAbs_FACE);
             while ( const TopoDS_Shape* f = faceIt->next() )
@@ -822,7 +879,7 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
     int nbKnownFaces;
     do {
       nbKnownFaces = faceMap.Extent();
-      StdMeshers_FaceSide *rightSide, *topSide; // sides of the quad
+      StdMeshers_FaceSidePtr rightSide, topSide; // sides of the quad
       for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
       {
         rightSide = thePrism.myWallQuads[i].back()->side[ QUAD_RIGHT_SIDE ];
@@ -854,8 +911,8 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
     {
       for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
       {
-        StdMeshers_FaceSide* topSide = thePrism.myWallQuads[i].back()->side[ QUAD_TOP_SIDE ];
-        const TopoDS_Edge &     topE = topSide->Edge( 0 );
+        StdMeshers_FaceSidePtr topSide = thePrism.myWallQuads[i].back()->side[ QUAD_TOP_SIDE ];
+        const TopoDS_Edge &       topE = topSide->Edge( 0 );
         if ( topSide->NbEdges() > 1 )
           return toSM( error(COMPERR_BAD_SHAPE, TCom("Side face #") <<
                              shapeID( thePrism.myWallQuads[i].back()->face )
@@ -901,8 +958,8 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
   // Check that the top FACE shares all the top EDGEs
   for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
   {
-    StdMeshers_FaceSide* topSide = thePrism.myWallQuads[i].back()->side[ QUAD_TOP_SIDE ];
-    const TopoDS_Edge &     topE = topSide->Edge( 0 );
+    StdMeshers_FaceSidePtr topSide = thePrism.myWallQuads[i].back()->side[ QUAD_TOP_SIDE ];
+    const TopoDS_Edge &       topE = topSide->Edge( 0 );
     if ( !myHelper->IsSubShape( topE, thePrism.myTop ))
       return toSM( error( TCom("Wrong source face (#") << shapeID( thePrism.myBottom )));
   }
@@ -925,7 +982,7 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
   if ( !computeWalls( thePrism ))
     return false;
 
-  // Analyse mesh and geometry to find block sub-shapes and submeshes
+  // Analyse mesh and geometry to find all block sub-shapes and submeshes
   if ( !myBlock.Init( myHelper, thePrism ))
     return toSM( error( myBlock.GetError()));
 
@@ -1148,7 +1205,7 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism)
     int wgt = 0; // "weight"
     for ( ; quad != thePrism.myWallQuads[iW].end(); ++quad )
     {
-      StdMeshers_FaceSide* lftSide = (*quad)->side[ QUAD_LEFT_SIDE ];
+      StdMeshers_FaceSidePtr lftSide = (*quad)->side[ QUAD_LEFT_SIDE ];
       for ( int i = 0; i < lftSide->NbEdges(); ++i )
       {
         ++wgt;
@@ -1167,7 +1224,7 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism)
       quad = thePrism.myWallQuads[iW].begin();
       for ( ; quad != thePrism.myWallQuads[iW].end(); ++quad )
         for ( int i = 0; i < NB_QUAD_SIDES; ++i )
-          (*quad)->side[ i ]->SetIgnoreMediumNodes( true );
+          (*quad)->side[ i ].grid->SetIgnoreMediumNodes( true );
     }
   }
 
@@ -1180,8 +1237,8 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism)
     Prism_3D::TQuadList::const_iterator quad = quads.begin();
     for ( ; quad != quads.end(); ++quad )
     {
-      StdMeshers_FaceSide* rgtSide = (*quad)->side[ QUAD_RIGHT_SIDE ]; // tgt
-      StdMeshers_FaceSide* lftSide = (*quad)->side[ QUAD_LEFT_SIDE ];  // src
+      StdMeshers_FaceSidePtr rgtSide = (*quad)->side[ QUAD_RIGHT_SIDE ]; // tgt
+      StdMeshers_FaceSidePtr lftSide = (*quad)->side[ QUAD_LEFT_SIDE ];  // src
       bool swapLeftRight = ( lftSide->NbSegments( /*update=*/true ) == 0 &&
                              rgtSide->NbSegments( /*update=*/true )  > 0 );
       if ( swapLeftRight )
@@ -1316,13 +1373,13 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism)
       // to compute stuctured quad mesh on wall FACEs
       // ---------------------------------------------------
       {
-        const TopoDS_Edge& botE = (*quad)->side[ QUAD_BOTTOM_SIDE ]->Edge(0);
-        const TopoDS_Edge& topE = (*quad)->side[ QUAD_TOP_SIDE    ]->Edge(0);
+        const TopoDS_Edge& botE = (*quad)->side[ QUAD_BOTTOM_SIDE ].grid->Edge(0);
+        const TopoDS_Edge& topE = (*quad)->side[ QUAD_TOP_SIDE    ].grid->Edge(0);
         SMESH_subMesh*    botSM = mesh->GetSubMesh( botE );
         SMESH_subMesh*    topSM = mesh->GetSubMesh( topE );
         SMESH_subMesh*    srcSM = botSM;
         SMESH_subMesh*    tgtSM = topSM;
-        if ( !srcSM->IsMeshComputed() && topSM->IsMeshComputed() )
+        if ( !srcSM->IsMeshComputed() && tgtSM->IsMeshComputed() )
           std::swap( srcSM, tgtSM );
 
         if ( !srcSM->IsMeshComputed() )
@@ -1333,10 +1390,31 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism)
         }
         srcSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
 
+        if ( tgtSM->IsMeshComputed() &&
+             tgtSM->GetSubMeshDS()->NbNodes() != srcSM->GetSubMeshDS()->NbNodes() )
+        {
+          // the top EDGE is computed differently than the bottom one,
+          // try to clear a wrong mesh
+          bool isAdjFaceMeshed = false;
+          PShapeIteratorPtr fIt = myHelper->GetAncestors( tgtSM->GetSubShape(),
+                                                          *mesh, TopAbs_FACE );
+          while ( const TopoDS_Shape* f = fIt->next() )
+            if (( isAdjFaceMeshed = mesh->GetSubMesh( *f )->IsMeshComputed() ))
+              break;
+          if ( isAdjFaceMeshed )
+            return toSM( error( TCom("Different nb of segment on logically horizontal edges #")
+                                << shapeID( botE ) << " and #"
+                                << shapeID( topE ) << ": "
+                                << tgtSM->GetSubMeshDS()->NbElements() << " != "
+                                << srcSM->GetSubMeshDS()->NbElements() ));
+          tgtSM->ComputeStateEngine( SMESH_subMesh::CLEAN );
+        }
         if ( !tgtSM->IsMeshComputed() )
         {
           // compute nodes on VERTEXes
-          tgtSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE );
+          SMESH_subMeshIteratorPtr smIt = tgtSM->getDependsOnIterator(/*includeSelf=*/false);
+          while ( smIt->more() )
+            smIt->next()->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE );
           // project segments
           DBGOUT( "COMPUTE H edge (proj) " << tgtSM->GetId());
           projector1D->myHyp.SetSourceEdge( TopoDS::Edge( srcSM->GetSubShape() ));
@@ -2264,7 +2342,7 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper*         helper,
   vector<int> nbUnitePerEdge( nbEdges, 0 ); // -1 means "joined to a previous"
 
   // consider continuous straight EDGEs as one side
-  const int nbSides = countNbSides( thePrism, nbUnitePerEdge );
+  const int nbSides = countNbSides( thePrism, nbUnitePerEdge, edgeLength );
 
   list< TopoDS_Edge >::const_iterator edgeIt = thePrism.myBottomEdges.begin();
   for ( iE = 0; iE < nbEdges; ++iE, ++edgeIt )
@@ -2274,7 +2352,7 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper*         helper,
     Prism_3D::TQuadList::const_iterator quad = thePrism.myWallQuads[ iE ].begin();
     for ( ; quad != thePrism.myWallQuads[ iE ].end(); ++quad )
     {
-      const TopoDS_Edge& quadBot = (*quad)->side[ QUAD_BOTTOM_SIDE ]->Edge( 0 );
+      const TopoDS_Edge& quadBot = (*quad)->side[ QUAD_BOTTOM_SIDE ].grid->Edge( 0 );
       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 ));
@@ -2283,8 +2361,6 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper*         helper,
     SHOWYXZ("p2 F "   <<iE, gpXYZ(faceColumns.rbegin()->second.front() ));
     SHOWYXZ("V First "<<iE, BRep_Tool::Pnt( TopExp::FirstVertex(*edgeIt,true )));
 
-    edgeLength[ iE ] = SMESH_Algo::EdgeLength( *edgeIt );
-
     if ( nbSides < NB_WALL_FACES ) // fill map used to split faces
       len2edgeMap.insert( make_pair( edgeLength[ iE ], iE )); // sort edges by length
   }
@@ -2297,7 +2373,7 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper*         helper,
     Prism_3D::TQuadList::const_iterator quad = thePrism.myWallQuads[ iE ].begin();
     for ( ; quad != thePrism.myWallQuads[ iE ].end(); ++quad )
     {
-      const TopoDS_Edge& quadBot = (*quad)->side[ QUAD_BOTTOM_SIDE ]->Edge( 0 );
+      const TopoDS_Edge& quadBot = (*quad)->side[ QUAD_BOTTOM_SIDE ].grid->Edge( 0 );
       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 ));
@@ -2524,8 +2600,8 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper*         helper,
     tFace.Set( ID_TOP_FACE, new BRepAdaptor_Surface( thePrism.myTop ), topPcurves, isForward );
     SMESH_Block::Insert( thePrism.myTop, ID_TOP_FACE, myShapeIDMap );
   }
-  // faceGridToPythonDump( SMESH_Block::ID_Fxy0 );
-  // faceGridToPythonDump( SMESH_Block::ID_Fxy1 );
+  //faceGridToPythonDump( SMESH_Block::ID_Fxy0, 50 );
+  //faceGridToPythonDump( SMESH_Block::ID_Fxy1 );
 
   // Fill map ShapeIndex to TParam2ColumnMap
   // ----------------------------------------
@@ -2729,7 +2805,8 @@ bool StdMeshers_PrismAsBlock::IsForwardEdge(SMESHDS_Mesh*           meshDS,
 //purpose  : Prints a script creating a normal grid on the prism side
 //=======================================================================
 
-void StdMeshers_PrismAsBlock::faceGridToPythonDump(const SMESH_Block::TShapeID face)
+void StdMeshers_PrismAsBlock::faceGridToPythonDump(const SMESH_Block::TShapeID face,
+                                                   const int                   nb)
 {
 #ifdef _DEBUG_
   gp_XYZ pOnF[6] = { gp_XYZ(0,0,0), gp_XYZ(0,0,1),
@@ -2739,7 +2816,7 @@ void StdMeshers_PrismAsBlock::faceGridToPythonDump(const SMESH_Block::TShapeID f
   cout << "mesh = smesh.Mesh( 'Face " << face << "')" << endl;
   SMESH_Block::TFace& f = myFace[ face - ID_FirstF ];
   gp_XYZ params = pOnF[ face - ID_FirstF ];
-  const int nb = 10; // nb face rows
+  //const int nb = 10; // nb face rows
   for ( int j = 0; j <= nb; ++j )
   {
     params.SetCoord( f.GetVInd(), double( j )/ nb );