src/StdMeshers/StdMeshers_Prism_3D.cxx

   1 // Copyright (C) 2007-2020  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 // File      : StdMeshers_Prism_3D.cxx
  24 // Module    : SMESH
  25 // Created   : Fri Oct 20 11:37:07 2006
  26 // Author    : Edward AGAPOV (eap)
  27 //
  28 #include "StdMeshers_Prism_3D.hxx"
  29
  30 #include "SMDS_EdgePosition.hxx"
  31 #include "SMDS_VolumeOfNodes.hxx"
  32 #include "SMDS_VolumeTool.hxx"
  33 #include "SMESH_Comment.hxx"
  34 #include "SMESH_Gen.hxx"
  35 #include "SMESH_HypoFilter.hxx"
  36 #include "SMESH_MeshEditor.hxx"
  37 #include "SMESH_MesherHelper.hxx"
  38 #include "StdMeshers_FaceSide.hxx"
  39 #include "StdMeshers_ProjectionSource1D.hxx"
  40 #include "StdMeshers_ProjectionSource2D.hxx"
  41 #include "StdMeshers_ProjectionUtils.hxx"
  42 #include "StdMeshers_Projection_1D.hxx"
  43 #include "StdMeshers_Projection_1D2D.hxx"
  44 #include "StdMeshers_Quadrangle_2D.hxx"
  45
  46 #include "utilities.h"
  47
  48 #include <BRepAdaptor_CompCurve.hxx>
  49 #include <BRep_Tool.hxx>
  50 #include <Bnd_B3d.hxx>
  51 #include <Geom2dAdaptor_Curve.hxx>
  52 #include <Geom2d_Line.hxx>
  53 #include <GeomLib_IsPlanarSurface.hxx>
  54 #include <Geom_Curve.hxx>
  55 #include <Standard_ErrorHandler.hxx>
  56 #include <TopExp.hxx>
  57 #include <TopExp_Explorer.hxx>
  58 #include <TopTools_ListIteratorOfListOfShape.hxx>
  59 #include <TopTools_ListOfShape.hxx>
  60 #include <TopTools_MapOfShape.hxx>
  61 #include <TopTools_SequenceOfShape.hxx>
  62 #include <TopoDS.hxx>
  63 #include <gp_Ax2.hxx>
  64 #include <gp_Ax3.hxx>
  65
  66 #include <limits>
  67 #include <numeric>
  68
  69 using namespace std;
  70
  71 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
  72 #define gpXYZ(n) SMESH_TNodeXYZ(n)
  73
  74 #ifdef _DEBUG_
  75 #define DBGOUT(msg) //cout << msg << endl;
  76 #define SHOWYXZ(msg, xyz)                                               \
  77   //{ gp_Pnt p (xyz); cout << msg << " ("<< p.X() << "; " <<p.Y() << "; " <<p.Z() << ") " <<endl; }
  78 #else
  79 #define DBGOUT(msg)
  80 #define SHOWYXZ(msg, xyz)
  81 #endif
  82
  83 namespace NSProjUtils = StdMeshers_ProjectionUtils;
  84
  85 typedef SMESH_Comment TCom;
  86
  87 enum { ID_BOT_FACE = SMESH_Block::ID_Fxy0,
  88        ID_TOP_FACE = SMESH_Block::ID_Fxy1,
  89        BOTTOM_EDGE = 0, TOP_EDGE, V0_EDGE, V1_EDGE, // edge IDs in face
  90        NB_WALL_FACES = 4 }; //
  91
  92 namespace
  93 {
  94   //=======================================================================
  95   /*!
  96    * \brief Auxiliary mesh
  97    */
  98   struct TmpMesh: public SMESH_Mesh
  99   {
 100     TmpMesh() {
 101       _isShapeToMesh = (_id = 0);
 102       _myMeshDS  = new SMESHDS_Mesh( _id, true );
 103     }
 104   };
 105   //=======================================================================
 106   /*!
 107    * \brief Quadrangle algorithm
 108    */
 109   class TQuadrangleAlgo : public StdMeshers_Quadrangle_2D
 110   {
 111     typedef NCollection_DataMap< TopoDS_Face, FaceQuadStruct::Ptr > TFace2QuadMap;
 112     TFace2QuadMap myFace2QuadMap;
 113
 114     TQuadrangleAlgo(SMESH_Gen* gen)
 115       : StdMeshers_Quadrangle_2D( gen->GetANewId(), gen)
 116     {
 117     }
 118   public:
 119
 120     //================================================================================
 121     // Clear data of TQuadrangleAlgo at destruction
 122     struct Cleaner
 123     {
 124       TQuadrangleAlgo* myAlgo;
 125
 126       Cleaner(TQuadrangleAlgo* algo): myAlgo( algo ){}
 127       ~Cleaner() { myAlgo->reset(); }
 128     };
 129
 130     //================================================================================
 131     // Return TQuadrangleAlgo singleton
 132     static TQuadrangleAlgo* instance( SMESH_Algo*         fatherAlgo,
 133                                       SMESH_MesherHelper* helper=0)
 134     {
 135       static TQuadrangleAlgo* algo = new TQuadrangleAlgo( fatherAlgo->GetGen() );
 136       if ( helper &&
 137            algo->myProxyMesh &&
 138            algo->myProxyMesh->GetMesh() != helper->GetMesh() )
 139         algo->myProxyMesh.reset( new SMESH_ProxyMesh( *helper->GetMesh() ));
 140
 141       algo->myQuadList.clear();
 142       algo->myHelper = 0;
 143
 144       if ( helper )
 145         algo->_quadraticMesh = helper->GetIsQuadratic();
 146
 147       return algo;
 148     }
 149
 150     //================================================================================
 151     // Clear collected data
 152     void reset()
 153     {
 154       myFace2QuadMap.Clear();
 155       StdMeshers_Quadrangle_2D::myQuadList.clear();
 156       StdMeshers_Quadrangle_2D::myHelper = nullptr;
 157     }
 158
 159     //================================================================================
 160     /*!
 161      * \brief Return FaceQuadStruct if a given FACE can be meshed by StdMeshers_Quadrangle_2D
 162      */
 163     FaceQuadStruct::Ptr CheckNbEdges(SMESH_Mesh&         theMesh,
 164                                      const TopoDS_Shape& theShape )
 165     {
 166       const TopoDS_Face& face = TopoDS::Face( theShape );
 167       if ( myFace2QuadMap.IsBound( face ))
 168         return myFace2QuadMap.Find( face );
 169
 170       FaceQuadStruct::Ptr &  resultQuad = * myFace2QuadMap.Bound( face, FaceQuadStruct::Ptr() );
 171
 172       FaceQuadStruct::Ptr quad =
 173         StdMeshers_Quadrangle_2D::CheckNbEdges( theMesh, face, /*considerMesh=*/false, myHelper );
 174       if ( quad )
 175       {
 176         // check if the quadrangle mesh would be valid
 177
 178         // check existing 1D mesh
 179         // int nbSegments[4], i = 0;
 180         // for ( FaceQuadStruct::Side & side : quad->side )
 181         //   nbSegments[ i++ ] = side.grid->NbSegments();
 182         // if ( nbSegments[0] > 0 && nbSegments[2] > 0 && nbSegments[0] != nbSegments[2] ||
 183         //      nbSegments[1] > 0 && nbSegments[3] > 0 && nbSegments[1] != nbSegments[3] )
 184         //   return resultQuad;
 185
 186         int nbEdges = 0;
 187         for ( FaceQuadStruct::Side & side : quad->side )
 188           nbEdges += side.grid->NbEdges();
 189         if ( nbEdges == 4 )
 190           return resultQuad = quad;
 191
 192         TmpMesh mesh;
 193         mesh.ShapeToMesh( face );
 194         SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
 195         SMESH_MesherHelper helper( mesh );
 196         helper.SetSubShape( face );
 197         helper.SetElementsOnShape( true );
 198
 199         // create nodes on all VERTEX'es
 200         for ( TopExp_Explorer vert( face, TopAbs_VERTEX ); vert.More(); vert.Next() )
 201           mesh.GetSubMesh( vert.Current() )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
 202
 203         FaceQuadStruct::Ptr tmpQuad( new FaceQuadStruct() );
 204         tmpQuad->side.resize( 4 );
 205
 206         // divide quad sides into halves at least
 207         const SMDS_MeshNode* node;
 208         for ( int iDir = 0; iDir < 2; ++iDir )
 209         {
 210           StdMeshers_FaceSidePtr sides[2] = { quad->side[iDir], quad->side[iDir+2] };
 211           std::map< double, const SMDS_MeshNode* > nodes[2];
 212           for ( int iS : { 0, 1 } )
 213           {
 214             node = SMESH_Algo::VertexNode( sides[iS]->FirstVertex(), meshDS );
 215             nodes[iS].insert( std::make_pair( 0, node ));
 216             double curLen = 0;
 217             for ( int iE = 1; iE < sides[iS]->NbEdges(); ++iE )
 218             {
 219               curLen += sides[iS]->EdgeLength( iE - 1 );
 220               double u = curLen / sides[iS]->Length();
 221               node = SMESH_Algo::VertexNode( sides[iS]->FirstVertex( iE ), meshDS );
 222               nodes[iS  ].insert( std::make_pair( u, node ));
 223               nodes[1-iS].insert( std::make_pair( u, nullptr ));
 224             }
 225             nodes[iS].insert( std::make_pair( 0.5, nullptr ));
 226             node = SMESH_Algo::VertexNode( sides[iS]->LastVertex(), meshDS );
 227             nodes[iS].insert( std::make_pair( 1, node ));
 228           }
 229
 230           for ( int iS : { 0, 1 } )
 231           {
 232             UVPtStructVec sideNodes;
 233             sideNodes.reserve( nodes[ iS ].size() );
 234             for ( auto & u_node : nodes[ iS ])
 235             {
 236               if ( !u_node.second )
 237               {
 238                 gp_Pnt p = sides[iS]->Value3d( u_node.first );
 239                 u_node.second = meshDS->AddNode( p.X(), p.Y(), p.Z() );
 240                 TopoDS_Edge edge;
 241                 double param = sides[iS]->Parameter( u_node.first, edge );
 242                 meshDS->SetNodeOnEdge( u_node.second, edge, param );
 243               }
 244               sideNodes.push_back( u_node.second );
 245               sideNodes.back().SetUV( helper.GetNodeUV( face, u_node.second ));
 246             }
 247             tmpQuad->side[ iS ? iDir+2 : iDir ] = StdMeshers_FaceSide::New( sideNodes, face );
 248           }
 249         }
 250         StdMeshers_Quadrangle_2D::myCheckOri = true;
 251         StdMeshers_Quadrangle_2D::myQuadList.clear();
 252         StdMeshers_Quadrangle_2D::myQuadList.push_back( tmpQuad );
 253         StdMeshers_Quadrangle_2D::myHelper = &helper;
 254         if ( StdMeshers_Quadrangle_2D::computeQuadDominant( mesh, face, tmpQuad ) &&
 255              StdMeshers_Quadrangle_2D::check())
 256         {
 257           resultQuad = quad;
 258         }
 259         StdMeshers_Quadrangle_2D::myQuadList.clear();
 260         StdMeshers_Quadrangle_2D::myHelper = nullptr;
 261       }
 262       return resultQuad;
 263     }
 264   };
 265
 266   //=======================================================================
 267   /*!
 268    * \brief Algorithm projecting 1D mesh
 269    */
 270   struct TProjction1dAlgo : public StdMeshers_Projection_1D
 271   {
 272     StdMeshers_ProjectionSource1D myHyp;
 273
 274     TProjction1dAlgo(SMESH_Gen* gen)
 275       : StdMeshers_Projection_1D( gen->GetANewId(), gen),
 276         myHyp( gen->GetANewId(), gen)
 277     {
 278       StdMeshers_Projection_1D::_sourceHypo = & myHyp;
 279     }
 280     static TProjction1dAlgo* instance( SMESH_Algo* fatherAlgo )
 281     {
 282       static TProjction1dAlgo* algo = new TProjction1dAlgo( fatherAlgo->GetGen() );
 283       return algo;
 284     }
 285   };
 286   //=======================================================================
 287   /*!
 288    * \brief Algorithm projecting 2D mesh
 289    */
 290   struct TProjction2dAlgo : public StdMeshers_Projection_1D2D
 291   {
 292     StdMeshers_ProjectionSource2D myHyp;
 293
 294     TProjction2dAlgo(SMESH_Gen* gen)
 295       : StdMeshers_Projection_1D2D( gen->GetANewId(), gen),
 296         myHyp( gen->GetANewId(), gen)
 297     {
 298       StdMeshers_Projection_2D::_sourceHypo = & myHyp;
 299     }
 300     static TProjction2dAlgo* instance( SMESH_Algo* fatherAlgo )
 301     {
 302       static TProjction2dAlgo* algo = new TProjction2dAlgo( fatherAlgo->GetGen() );
 303       return algo;
 304     }
 305     const NSProjUtils::TNodeNodeMap& GetNodesMap()
 306     {
 307       return _src2tgtNodes;
 308     }
 309     void SetEventListener( SMESH_subMesh* tgtSubMesh )
 310     {
 311       NSProjUtils::SetEventListener( tgtSubMesh,
 312                                      _sourceHypo->GetSourceFace(),
 313                                      _sourceHypo->GetSourceMesh() );
 314     }
 315   };
 316   //=======================================================================
 317   /*!
 318    * \brief Returns already computed EDGEs
 319    */
 320   void getPrecomputedEdges( SMESH_MesherHelper&    theHelper,
 321                             const TopoDS_Shape&    theShape,
 322                             vector< TopoDS_Edge >& theEdges)
 323   {
 324     theEdges.clear();
 325
 326     SMESHDS_Mesh* meshDS = theHelper.GetMeshDS();
 327     SMESHDS_SubMesh* sm;
 328
 329     TopTools_IndexedMapOfShape edges;
 330     TopExp::MapShapes( theShape, TopAbs_EDGE, edges );
 331     for ( int iE = 1; iE <= edges.Extent(); ++iE )
 332     {
 333       const TopoDS_Shape edge = edges( iE );
 334       if (( ! ( sm = meshDS->MeshElements( edge ))) ||
 335           ( sm->NbElements() == 0 ))
 336         continue;
 337
 338       // there must not be FACEs meshed with triangles and sharing a computed EDGE
 339       // as the precomputed EDGEs are used for propagation other to 'vertical' EDGEs
 340       bool faceFound = false;
 341       PShapeIteratorPtr faceIt =
 342         theHelper.GetAncestors( edge, *theHelper.GetMesh(), TopAbs_FACE );
 343       while ( const TopoDS_Shape* face = faceIt->next() )
 344
 345         if (( sm = meshDS->MeshElements( *face )) &&
 346             ( sm->NbElements() > 0 ) &&
 347             ( !theHelper.IsSameElemGeometry( sm, SMDSGeom_QUADRANGLE ) ))
 348         {
 349           faceFound = true;
 350           break;
 351         }
 352       if ( !faceFound )
 353         theEdges.push_back( TopoDS::Edge( edge ));
 354     }
 355   }
 356
 357   //================================================================================
 358   /*!
 359    * \brief Make \a botE be the BOTTOM_SIDE of \a quad.
 360    *        Return false if the BOTTOM_SIDE is composite
 361    */
 362   //================================================================================
 363
 364   bool setBottomEdge( const TopoDS_Edge&   botE,
 365                       FaceQuadStruct::Ptr& quad,
 366                       const TopoDS_Shape&  face)
 367   {
 368     quad->side[ QUAD_TOP_SIDE  ].grid->Reverse();
 369     quad->side[ QUAD_LEFT_SIDE ].grid->Reverse();
 370     int edgeIndex = 0;
 371     bool isComposite = false;
 372     for ( size_t i = 0; i < quad->side.size(); ++i )
 373     {
 374       StdMeshers_FaceSidePtr quadSide = quad->side[i];
 375       for ( int iE = 0; iE < quadSide->NbEdges(); ++iE )
 376         if ( botE.IsSame( quadSide->Edge( iE )))
 377         {
 378           if ( quadSide->NbEdges() > 1 )
 379             isComposite = true; //return false;
 380           edgeIndex = i;
 381           i = quad->side.size(); // to quit from the outer loop
 382           break;
 383         }
 384     }
 385     if ( edgeIndex != QUAD_BOTTOM_SIDE )
 386       quad->shift( quad->side.size() - edgeIndex, /*keepUnitOri=*/false );
 387
 388     quad->face = TopoDS::Face( face );
 389
 390     return !isComposite;
 391   }
 392
 393   //================================================================================
 394   /*!
 395    * \brief Return iterator pointing to node column for the given parameter
 396    * \param columnsMap - node column map
 397    * \param parameter - parameter
 398    * \retval TParam2ColumnMap::iterator - result
 399    *
 400    * it returns closest left column
 401    */
 402   //================================================================================
 403
 404   TParam2ColumnIt getColumn( const TParam2ColumnMap* columnsMap,
 405                              const double            parameter )
 406   {
 407     TParam2ColumnIt u_col = columnsMap->upper_bound( parameter );
 408     if ( u_col != columnsMap->begin() )
 409       --u_col;
 410     return u_col; // return left column
 411   }
 412
 413   //================================================================================
 414   /*!
 415    * \brief Return nodes around given parameter and a ratio
 416    * \param column - node column
 417    * \param param - parameter
 418    * \param node1 - lower node
 419    * \param node2 - upper node
 420    * \retval double - ratio
 421    */
 422   //================================================================================
 423
 424   double getRAndNodes( const TNodeColumn*     column,
 425                        const double           param,
 426                        const SMDS_MeshNode* & node1,
 427                        const SMDS_MeshNode* & node2)
 428   {
 429     if ( param >= 1.0 || column->size() == 1) {
 430       node1 = node2 = column->back();
 431       return 0;
 432     }
 433
 434     int i = int( param * ( column->size() - 1 ));
 435     double u0 = double( i )/ double( column->size() - 1 );
 436     double r = ( param - u0 ) * ( column->size() - 1 );
 437
 438     node1 = (*column)[ i ];
 439     node2 = (*column)[ i + 1];
 440     return r;
 441   }
 442
 443   //================================================================================
 444   /*!
 445    * \brief Compute boundary parameters of face parts
 446     * \param nbParts - nb of parts to split columns into
 447     * \param columnsMap - node columns of the face to split
 448     * \param params - computed parameters
 449    */
 450   //================================================================================
 451
 452   void splitParams( const int               nbParts,
 453                     const TParam2ColumnMap* columnsMap,
 454                     vector< double > &      params)
 455   {
 456     params.clear();
 457     params.reserve( nbParts + 1 );
 458     TParam2ColumnIt last_par_col = --columnsMap->end();
 459     double par = columnsMap->begin()->first; // 0.
 460     double parLast = last_par_col->first;
 461     params.push_back( par );
 462     for ( int i = 0; i < nbParts - 1; ++ i )
 463     {
 464       double partSize = ( parLast - par ) / double ( nbParts - i );
 465       TParam2ColumnIt par_col = getColumn( columnsMap, par + partSize );
 466       if ( par_col->first == par ) {
 467         ++par_col;
 468         if ( par_col == last_par_col ) {
 469           while ( i < nbParts - 1 )
 470             params.push_back( par + partSize * i++ );
 471           break;
 472         }
 473       }
 474       par = par_col->first;
 475       params.push_back( par );
 476     }
 477     params.push_back( parLast ); // 1.
 478   }
 479
 480   //================================================================================
 481   /*!
 482    * \brief Return coordinate system for z-th layer of nodes
 483    */
 484   //================================================================================
 485
 486   gp_Ax2 getLayerCoordSys(const int                           z,
 487                           const vector< const TNodeColumn* >& columns,
 488                           int&                                xColumn)
 489   {
 490     // gravity center of a layer
 491     gp_XYZ O(0,0,0);
 492     int vertexCol = -1;
 493     for ( size_t i = 0; i < columns.size(); ++i )
 494     {
 495       O += gpXYZ( (*columns[ i ])[ z ]);
 496       if ( vertexCol < 0 &&
 497            columns[ i ]->front()->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX )
 498         vertexCol = i;
 499     }
 500     O /= columns.size();
 501
 502     // Z axis
 503     gp_Vec Z(0,0,0);
 504     int iPrev = columns.size()-1;
 505     for ( size_t i = 0; i < columns.size(); ++i )
 506     {
 507       gp_Vec v1( O, gpXYZ( (*columns[ iPrev ])[ z ]));
 508       gp_Vec v2( O, gpXYZ( (*columns[ i ]    )[ z ]));
 509       Z += v1 ^ v2;
 510       iPrev = i;
 511     }
 512
 513     if ( vertexCol >= 0 )
 514     {
 515       O = gpXYZ( (*columns[ vertexCol ])[ z ]);
 516     }
 517     if ( xColumn < 0 || xColumn >= (int) columns.size() )
 518     {
 519       // select a column for X dir
 520       double maxDist = 0;
 521       for ( size_t i = 0; i < columns.size(); ++i )
 522       {
 523         double dist = ( O - gpXYZ((*columns[ i ])[ z ])).SquareModulus();
 524         if ( dist > maxDist )
 525         {
 526           xColumn = i;
 527           maxDist = dist;
 528         }
 529       }
 530     }
 531
 532     // X axis
 533     gp_Vec X( O, gpXYZ( (*columns[ xColumn ])[ z ]));
 534
 535     return gp_Ax2( O, Z, X);
 536   }
 537
 538   //================================================================================
 539   /*!
 540    * \brief Removes submeshes that are or can be meshed with regular grid from given list
 541    *  \retval int - nb of removed submeshes
 542    */
 543   //================================================================================
 544
 545   int removeQuasiQuads(list< SMESH_subMesh* >&   notQuadSubMesh,
 546                        SMESH_MesherHelper*       helper,
 547                        TQuadrangleAlgo*          quadAlgo)
 548   {
 549     int nbRemoved = 0;
 550     //SMESHDS_Mesh* mesh = notQuadSubMesh.front()->GetFather()->GetMeshDS();
 551     list< SMESH_subMesh* >::iterator smIt = notQuadSubMesh.begin();
 552     while ( smIt != notQuadSubMesh.end() )
 553     {
 554       SMESH_subMesh* faceSm = *smIt;
 555       SMESHDS_SubMesh* faceSmDS = faceSm->GetSubMeshDS();
 556       int nbQuads = faceSmDS ? faceSmDS->NbElements() : 0;
 557       bool toRemove;
 558       if ( nbQuads > 0 )
 559         toRemove = helper->IsStructured( faceSm );
 560       else
 561         toRemove = ( quadAlgo->CheckNbEdges( *helper->GetMesh(),
 562                                              faceSm->GetSubShape() ) != NULL );
 563       nbRemoved += toRemove;
 564       if ( toRemove )
 565         smIt = notQuadSubMesh.erase( smIt );
 566       else
 567         ++smIt;
 568     }
 569
 570     return nbRemoved;
 571   }
 572
 573   //================================================================================
 574   /*!
 575    * \brief Return and angle between two EDGEs
 576    *  \return double - the angle normalized so that
 577    * >~ 0  -> 2.0
 578    *  PI/2 -> 1.0
 579    *  PI   -> 0.0
 580    * -PI/2 -> -1.0
 581    * <~ 0  -> -2.0
 582    */
 583   //================================================================================
 584
 585   // double normAngle(const TopoDS_Edge & E1, const TopoDS_Edge & E2, const TopoDS_Face & F)
 586   // {
 587   //   return SMESH_MesherHelper::GetAngle( E1, E2, F ) / ( 0.5 * M_PI );
 588   // }
 589
 590   //================================================================================
 591   /*!
 592    * Consider continuous straight EDGES as one side - mark them to unite
 593    */
 594   //================================================================================
 595
 596   int countNbSides( const Prism_3D::TPrismTopo & thePrism,
 597                     vector<int> &                nbUnitePerEdge,
 598                     vector< double > &           edgeLength)
 599   {
 600     int nbEdges = thePrism.myNbEdgesInWires.front();  // nb outer edges
 601     int nbSides = nbEdges;
 602
 603
 604     list< TopoDS_Edge >::const_iterator edgeIt = thePrism.myBottomEdges.begin();
 605     std::advance( edgeIt, nbEdges-1 );
 606     TopoDS_Edge   prevE = *edgeIt;
 607     // bool isPrevStraight = SMESH_Algo::IsStraight( prevE );
 608     // int           iPrev = nbEdges - 1;
 609
 610     // int iUnite = -1; // the first of united EDGEs
 611
 612     // analyse angles between EDGEs
 613     int nbCorners = 0;
 614     vector< bool > isCorner( nbEdges );
 615     edgeIt = thePrism.myBottomEdges.begin();
 616     for ( int iE = 0; iE < nbEdges; ++iE, ++edgeIt )
 617     {
 618       const TopoDS_Edge&  curE = *edgeIt;
 619       edgeLength[ iE ] = SMESH_Algo::EdgeLength( curE );
 620
 621       // double normAngle = normAngle( prevE, curE, thePrism.myBottom );
 622       // isCorner[ iE ] = false;
 623       // if ( normAngle < 2.0 )
 624       // {
 625       //   if ( normAngle < 0.001 ) // straight or obtuse angle
 626       //   {
 627       //     // unite EDGEs in order not to put a corner of the unit quadrangle at this VERTEX
 628       //     if ( iUnite < 0 )
 629       //       iUnite = iPrev;
 630       //     nbUnitePerEdge[ iUnite ]++;
 631       //     nbUnitePerEdge[ iE ] = -1;
 632       //     --nbSides;
 633       //   }
 634       //   else
 635       //   {
 636       //     isCorner[ iE ] = true;
 637       //     nbCorners++;
 638       //     iUnite = -1;
 639       //   }
 640       // }
 641       // prevE = curE;
 642     }
 643
 644     if ( nbCorners > 4 )
 645     {
 646       // define which of corners to put on a side of the unit quadrangle
 647     }
 648     // edgeIt = thePrism.myBottomEdges.begin();
 649     // for ( int iE = 0; iE < nbEdges; ++iE, ++edgeIt )
 650     // {
 651     //   const TopoDS_Edge&  curE = *edgeIt;
 652     //   edgeLength[ iE ] = SMESH_Algo::EdgeLength( curE );
 653
 654     //   const bool isCurStraight = SMESH_Algo::IsStraight( curE );
 655     //   if ( isPrevStraight && isCurStraight && SMESH_Algo::IsContinuous( prevE, curE ))
 656     //   {
 657     //     if ( iUnite < 0 )
 658     //       iUnite = iPrev;
 659     //     nbUnitePerEdge[ iUnite ]++;
 660     //     nbUnitePerEdge[ iE ] = -1;
 661     //     --nbSides;
 662     //   }
 663     //   else
 664     //   {
 665     //     iUnite = -1;
 666     //   }
 667     //   prevE          = curE;
 668     //   isPrevStraight = isCurStraight;
 669     //   iPrev = iE;
 670     // }
 671
 672     return nbSides;
 673   }
 674
 675   //================================================================================
 676   /*!
 677    * \brief Count EDGEs ignoring degenerated ones
 678    */
 679   //================================================================================
 680
 681   int CountEdges( const TopoDS_Face& face )
 682   {
 683     int nbE = 0;
 684     for ( TopExp_Explorer edgeExp( face, TopAbs_EDGE ); edgeExp.More(); edgeExp.Next() )
 685       if ( !SMESH_Algo::isDegenerated( TopoDS::Edge( edgeExp.Current() )))
 686         ++nbE;
 687
 688     return nbE;
 689   }
 690
 691   //================================================================================
 692   /*!
 693    * \brief Set/get wire index to FaceQuadStruct
 694    */
 695   //================================================================================
 696
 697   void setWireIndex( TFaceQuadStructPtr& quad, int iWire )
 698   {
 699     quad->iSize = iWire;
 700   }
 701   int getWireIndex( const TFaceQuadStructPtr& quad )
 702   {
 703     return quad->iSize;
 704   }
 705
 706   //================================================================================
 707   /*!
 708    * \brief Print Python commands adding given points to a mesh
 709    */
 710   //================================================================================
 711
 712   void pointsToPython(const std::vector<gp_XYZ>& p)
 713   {
 714 #ifdef _DEBUG_
 715     for ( size_t i = SMESH_Block::ID_V000; i < p.size(); ++i )
 716     {
 717       cout << "mesh.AddNode( " << p[i].X() << ", "<< p[i].Y() << ", "<< p[i].Z() << ") # " << i <<" " ;
 718       SMESH_Block::DumpShapeID( i, cout ) << endl;
 719     }
 720 #endif
 721   }
 722
 723 } // namespace
 724
 725 //=======================================================================
 726 //function : StdMeshers_Prism_3D
 727 //purpose  :
 728 //=======================================================================
 729
 730 StdMeshers_Prism_3D::StdMeshers_Prism_3D(int hypId, SMESH_Gen* gen)
 731   :SMESH_3D_Algo(hypId, gen)
 732 {
 733   _name                    = "Prism_3D";
 734   _shapeType               = (1 << TopAbs_SOLID); // 1 bit per shape type
 735   _onlyUnaryInput          = false; // mesh all SOLIDs at once
 736   _requireDiscreteBoundary = false; // mesh FACEs and EDGEs by myself
 737   _supportSubmeshes        = true;  // "source" FACE must be meshed by other algo
 738   _neededLowerHyps[ 1 ]    = true;  // suppress warning on hiding a global 1D algo
 739   _neededLowerHyps[ 2 ]    = true;  // suppress warning on hiding a global 2D algo
 740
 741   //myProjectTriangles       = false;
 742   mySetErrorToSM           = true;  // to pass an error to a sub-mesh of a current solid or not
 743   myPrevBottomSM           = 0;     // last treated bottom sub-mesh with a suitable algorithm
 744 }
 745
 746 //================================================================================
 747 /*!
 748  * \brief Destructor
 749  */
 750 //================================================================================
 751
 752 StdMeshers_Prism_3D::~StdMeshers_Prism_3D()
 753 {
 754   pointsToPython( std::vector<gp_XYZ>() ); // avoid warning: pointsToPython defined but not used
 755 }
 756
 757 //=======================================================================
 758 //function : CheckHypothesis
 759 //purpose  :
 760 //=======================================================================
 761
 762 bool StdMeshers_Prism_3D::CheckHypothesis(SMESH_Mesh&                          aMesh,
 763                                           const TopoDS_Shape&                  aShape,
 764                                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 765 {
 766   // no hypothesis
 767   aStatus = SMESH_Hypothesis::HYP_OK;
 768   return true;
 769 }
 770
 771 //=======================================================================
 772 //function : Compute
 773 //purpose  : Compute mesh on a COMPOUND of SOLIDs
 774 //=======================================================================
 775
 776 bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theShape)
 777 {
 778   SMESH_MesherHelper helper( theMesh );
 779   myHelper = &helper;
 780   myPrevBottomSM = 0;
 781   TQuadrangleAlgo::Cleaner( TQuadrangleAlgo::instance( this ));
 782
 783   int nbSolids = helper.Count( theShape, TopAbs_SOLID, /*skipSame=*/false );
 784   if ( nbSolids < 1 )
 785     return true;
 786
 787   TopTools_IndexedDataMapOfShapeListOfShape faceToSolids;
 788   TopExp::MapShapesAndAncestors( theShape, TopAbs_FACE, TopAbs_SOLID, faceToSolids );
 789
 790   // look for meshed FACEs ("source" FACEs) that must be prism bottoms
 791   list< TopoDS_Face > meshedFaces, notQuadMeshedFaces, notQuadFaces;
 792   const bool meshHasQuads = ( theMesh.NbQuadrangles() > 0 );
 793   for ( int iF = 1; iF <= faceToSolids.Extent(); ++iF )
 794   {
 795     const TopoDS_Face& face = TopoDS::Face( faceToSolids.FindKey( iF ));
 796     SMESH_subMesh*   faceSM = theMesh.GetSubMesh( face );
 797     if ( !faceSM->IsEmpty() )
 798     {
 799       if ( !meshHasQuads ||
 800            !helper.IsSameElemGeometry( faceSM->GetSubMeshDS(), SMDSGeom_QUADRANGLE ) ||
 801            !helper.IsStructured( faceSM )
 802            )
 803         notQuadMeshedFaces.push_front( face );
 804       else if ( myHelper->Count( face, TopAbs_EDGE, /*ignoreSame=*/false ) != 4 )
 805         meshedFaces.push_front( face );
 806       else
 807         meshedFaces.push_back( face );
 808     }
 809     // not add not quadrilateral FACE as we can't compute it
 810     // else if ( !quadAlgo->CheckNbEdges( theMesh, face ))
 811     // // not add not quadrilateral FACE as it can be a prism side
 812     // // else if ( myHelper->Count( face, TopAbs_EDGE, /*ignoreSame=*/false ) != 4 )
 813     // {
 814     //   notQuadFaces.push_back( face );
 815     // }
 816   }
 817   // notQuadFaces are of medium priority, put them before ordinary meshed faces
 818   meshedFaces.splice( meshedFaces.begin(), notQuadFaces );
 819   // notQuadMeshedFaces are of highest priority, put them before notQuadFaces
 820   meshedFaces.splice( meshedFaces.begin(), notQuadMeshedFaces );
 821
 822   Prism_3D::TPrismTopo prism;
 823   myPropagChains = 0;
 824   bool selectBottom = meshedFaces.empty();
 825
 826   if ( nbSolids == 1 )
 827   {
 828     TopoDS_Shape solid = TopExp_Explorer( theShape, TopAbs_SOLID ).Current();
 829     if ( !meshedFaces.empty() )
 830       prism.myBottom = meshedFaces.front();
 831     return ( initPrism( prism, solid, selectBottom ) &&
 832              compute( prism ));
 833   }
 834
 835   // find propagation chains from already computed EDGEs
 836   vector< TopoDS_Edge > computedEdges;
 837   getPrecomputedEdges( helper, theShape, computedEdges );
 838   myPropagChains = new TopTools_IndexedMapOfShape[ computedEdges.size() + 1 ];
 839   SMESHUtils::ArrayDeleter< TopTools_IndexedMapOfShape > pcDel( myPropagChains );
 840   for ( size_t i = 0, nb = 0; i < computedEdges.size(); ++i )
 841   {
 842     StdMeshers_ProjectionUtils::GetPropagationEdge( &theMesh, TopoDS_Edge(),
 843                                                     computedEdges[i], myPropagChains + nb );
 844     if ( myPropagChains[ nb ].Extent() < 2 ) // an empty map is a termination sign
 845       myPropagChains[ nb ].Clear();
 846     else
 847       nb++;
 848   }
 849
 850   TopTools_MapOfShape meshedSolids;
 851   list< Prism_3D::TPrismTopo > meshedPrism;
 852   list< TopoDS_Face > suspectSourceFaces;
 853   TopTools_ListIteratorOfListOfShape solidIt;
 854
 855   while ( meshedSolids.Extent() < nbSolids )
 856   {
 857     if ( _computeCanceled )
 858       return toSM( error( SMESH_ComputeError::New(COMPERR_CANCELED)));
 859
 860     // compute prisms having avident computed source FACE
 861     while ( !meshedFaces.empty() )
 862     {
 863       TopoDS_Face face = meshedFaces.front();
 864       meshedFaces.pop_front();
 865       TopTools_ListOfShape& solidList = faceToSolids.ChangeFromKey( face );
 866       while ( !solidList.IsEmpty() )
 867       {
 868         TopoDS_Shape solid = solidList.First();
 869         solidList.RemoveFirst();
 870         if ( meshedSolids.Add( solid ))
 871         {
 872           prism.Clear();
 873           prism.myBottom = face;
 874           if ( !initPrism( prism, solid, selectBottom ) ||
 875                !compute( prism ))
 876             return false;
 877
 878           SMESHDS_SubMesh* smDS = theMesh.GetMeshDS()->MeshElements( prism.myTop );
 879           if ( !myHelper->IsSameElemGeometry( smDS, SMDSGeom_QUADRANGLE ) ||
 880                !myHelper->IsStructured( theMesh.GetSubMesh( prism.myTop )))
 881           {
 882             meshedFaces.push_front( prism.myTop );
 883           }
 884           else
 885           {
 886             suspectSourceFaces.push_back( prism.myTop );
 887           }
 888           meshedPrism.push_back( prism );
 889         }
 890       }
 891     }
 892     if ( meshedSolids.Extent() == nbSolids )
 893       break;
 894
 895     // below in the loop we try to find source FACEs somehow
 896
 897     // project mesh from source FACEs of computed prisms to
 898     // prisms sharing wall FACEs
 899     list< Prism_3D::TPrismTopo >::iterator prismIt = meshedPrism.begin();
 900     for ( ; prismIt != meshedPrism.end(); ++prismIt )
 901     {
 902       for ( size_t iW = 0; iW < prismIt->myWallQuads.size(); ++iW )
 903       {
 904         Prism_3D::TQuadList::iterator wQuad = prismIt->myWallQuads[iW].begin();
 905         for ( ; wQuad != prismIt->myWallQuads[iW].end(); ++ wQuad )
 906         {
 907           const TopoDS_Face& wFace = (*wQuad)->face;
 908           TopTools_ListOfShape& solidList = faceToSolids.ChangeFromKey( wFace );
 909           solidIt.Initialize( solidList );
 910           while ( solidIt.More() )
 911           {
 912             const TopoDS_Shape& solid = solidIt.Value();
 913             if ( meshedSolids.Contains( solid )) {
 914               solidList.Remove( solidIt );
 915               continue; // already computed prism
 916             }
 917             if ( myHelper->IsBlock( solid )) {
 918               solidIt.Next();
 919               continue; // too trivial
 920             }
 921             // find a source FACE of the SOLID: it's a FACE sharing a bottom EDGE with wFace
 922             const TopoDS_Edge& wEdge = (*wQuad)->side[ QUAD_TOP_SIDE ].grid->Edge(0);
 923             PShapeIteratorPtr faceIt = myHelper->GetAncestors( wEdge, *myHelper->GetMesh(),
 924                                                                TopAbs_FACE);
 925             while ( const TopoDS_Shape* f = faceIt->next() )
 926             {
 927               const TopoDS_Face& candidateF = TopoDS::Face( *f );
 928               if ( candidateF.IsSame( wFace )) continue;
 929               // select a source FACE: prismIt->myBottom or prismIt->myTop
 930               TopoDS_Face sourceF = prismIt->myBottom;
 931               for ( TopExp_Explorer v( prismIt->myTop, TopAbs_VERTEX ); v.More(); v.Next() )
 932                 if ( myHelper->IsSubShape( v.Current(), candidateF )) {
 933                   sourceF = prismIt->myTop;
 934                   break;
 935                 }
 936               prism.Clear();
 937               prism.myBottom = candidateF;
 938               mySetErrorToSM = false;
 939               if ( !myHelper->IsSubShape( candidateF, prismIt->myShape3D ) &&
 940                    myHelper ->IsSubShape( candidateF, solid ) &&
 941                    !myHelper->GetMesh()->GetSubMesh( candidateF )->IsMeshComputed() &&
 942                    initPrism( prism, solid, /*selectBottom=*/false ) &&
 943                    !myHelper->GetMesh()->GetSubMesh( prism.myTop )->IsMeshComputed() &&
 944                    !myHelper->GetMesh()->GetSubMesh( prism.myBottom )->IsMeshComputed() &&
 945                    project2dMesh( sourceF, prism.myBottom ))
 946               {
 947                 mySetErrorToSM = true;
 948                 if ( !compute( prism ))
 949                   return false;
 950                 SMESHDS_SubMesh* smDS = theMesh.GetMeshDS()->MeshElements( prism.myTop );
 951                 if ( !myHelper->IsSameElemGeometry( smDS, SMDSGeom_QUADRANGLE ))
 952                 {
 953                   meshedFaces.push_front( prism.myTop );
 954                   meshedFaces.push_front( prism.myBottom );
 955                   selectBottom = false;
 956                 }
 957                 meshedPrism.push_back( prism );
 958                 meshedSolids.Add( solid );
 959               }
 960               InitComputeError();
 961             }
 962             mySetErrorToSM = true;
 963             InitComputeError();
 964             if ( meshedSolids.Contains( solid ))
 965               solidList.Remove( solidIt );
 966             else
 967               solidIt.Next();
 968           }
 969         }
 970       }
 971       if ( !meshedFaces.empty() )
 972         break; // to compute prisms with avident sources
 973     }
 974
 975     if ( meshedFaces.empty() )
 976     {
 977       meshedFaces.splice( meshedFaces.end(), suspectSourceFaces );
 978       selectBottom = true;
 979     }
 980
 981     // find FACEs with local 1D hyps, which has to be computed by now,
 982     // or at least any computed FACEs
 983     if ( meshedFaces.empty() )
 984     {
 985       int prevNbFaces = 0;
 986       for ( int iF = 1; iF <= faceToSolids.Extent(); ++iF )
 987       {
 988         const TopoDS_Face&               face = TopoDS::Face( faceToSolids.FindKey( iF ));
 989         const TopTools_ListOfShape& solidList = faceToSolids.FindFromKey( face );
 990         if ( solidList.IsEmpty() ) continue;
 991         SMESH_subMesh*                 faceSM = theMesh.GetSubMesh( face );
 992         if ( !faceSM->IsEmpty() )
 993         {
 994           int nbFaces = faceSM->GetSubMeshDS()->NbElements();
 995           if ( prevNbFaces < nbFaces )
 996           {
 997             if ( !meshedFaces.empty() ) meshedFaces.pop_back();
 998             meshedFaces.push_back( face ); // lower priority
 999             selectBottom = true;
1000             prevNbFaces = nbFaces;
1001           }
1002         }
1003         else
1004         {
1005           bool allSubMeComputed = true;
1006           SMESH_subMeshIteratorPtr smIt = faceSM->getDependsOnIterator(false,true);
1007           while ( smIt->more() && allSubMeComputed )
1008             allSubMeComputed = smIt->next()->IsMeshComputed();
1009           if ( allSubMeComputed )
1010           {
1011             faceSM->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1012             if ( !faceSM->IsEmpty() ) {
1013               meshedFaces.push_front( face ); // higher priority
1014               selectBottom = true;
1015               break;
1016             }
1017             else {
1018               faceSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1019             }
1020           }
1021         }
1022       }
1023     }
1024
1025
1026     // TODO. there are other ways to find out the source FACE:
1027     // propagation, topological similarity, etc...
1028
1029     // simply try to mesh all not meshed SOLIDs
1030     if ( meshedFaces.empty() )
1031     {
1032       for ( TopExp_Explorer solid( theShape, TopAbs_SOLID ); solid.More(); solid.Next() )
1033       {
1034         mySetErrorToSM = false;
1035         prism.Clear();
1036         if ( !meshedSolids.Contains( solid.Current() ) &&
1037              initPrism( prism, solid.Current() ))
1038         {
1039           mySetErrorToSM = true;
1040           if ( !compute( prism ))
1041             return false;
1042           meshedFaces.push_front( prism.myTop );
1043           meshedFaces.push_front( prism.myBottom );
1044           meshedPrism.push_back( prism );
1045           meshedSolids.Add( solid.Current() );
1046           selectBottom = true;
1047         }
1048         mySetErrorToSM = true;
1049       }
1050     }
1051
1052     if ( meshedFaces.empty() ) // set same error to 10 not-computed solids
1053     {
1054       SMESH_ComputeErrorPtr err = SMESH_ComputeError::New
1055         ( COMPERR_BAD_INPUT_MESH, "No meshed source face found", this );
1056
1057       const int maxNbErrors = 10; // limit nb errors not to overload the Compute dialog
1058       TopExp_Explorer solid( theShape, TopAbs_SOLID );
1059       for ( int i = 0; ( i < maxNbErrors && solid.More() ); ++i, solid.Next() )
1060         if ( !meshedSolids.Contains( solid.Current() ))
1061         {
1062           SMESH_subMesh* sm = theMesh.GetSubMesh( solid.Current() );
1063           sm->GetComputeError() = err;
1064         }
1065       return error( err );
1066     }
1067   }
1068   return error( COMPERR_OK );
1069 }
1070
1071 //================================================================================
1072 /*!
1073  * \brief Find wall faces by bottom edges
1074  */
1075 //================================================================================
1076
1077 bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
1078                                         const int              totalNbFaces)
1079 {
1080   thePrism.myWallQuads.clear();
1081
1082   SMESH_Mesh* mesh = myHelper->GetMesh();
1083
1084   TQuadrangleAlgo* quadAlgo = TQuadrangleAlgo::instance( this, myHelper );
1085
1086   TopTools_MapOfShape faceMap;
1087   TopTools_IndexedDataMapOfShapeListOfShape edgeToFaces;
1088   TopExp::MapShapesAndAncestors( thePrism.myShape3D,
1089                                  TopAbs_EDGE, TopAbs_FACE, edgeToFaces );
1090
1091   // ------------------------------
1092   // Get the 1st row of wall FACEs
1093   // ------------------------------
1094
1095   list< TopoDS_Edge >::iterator edge = thePrism.myBottomEdges.begin();
1096   std::list< int >::iterator     nbE = thePrism.myNbEdgesInWires.begin();
1097   std::list< int > nbQuadsPerWire;
1098   int iE = 0, iWire = 0;
1099   while ( edge != thePrism.myBottomEdges.end() )
1100   {
1101     ++iE;
1102     if ( SMESH_Algo::isDegenerated( *edge ))
1103     {
1104       edge = thePrism.myBottomEdges.erase( edge );
1105       --iE;
1106       --(*nbE);
1107     }
1108     else
1109     {
1110       bool hasWallFace = false;
1111       TopTools_ListIteratorOfListOfShape faceIt( edgeToFaces.FindFromKey( *edge ));
1112       for ( ; faceIt.More(); faceIt.Next() )
1113       {
1114         const TopoDS_Face& face = TopoDS::Face( faceIt.Value() );
1115         if ( !thePrism.myBottom.IsSame( face ))
1116         {
1117           hasWallFace = true;
1118           Prism_3D::TQuadList quadList( 1, quadAlgo->CheckNbEdges( *mesh, face ));
1119           if ( !quadList.back() )
1120             return toSM( error(TCom("Side face #") << shapeID( face )
1121                                << " not meshable with quadrangles"));
1122           bool isCompositeBase = ! setBottomEdge( *edge, quadList.back(), face );
1123           if ( isCompositeBase )
1124           {
1125             // it's OK if all EDGEs of the bottom side belongs to the bottom FACE
1126             StdMeshers_FaceSidePtr botSide = quadList.back()->side[ QUAD_BOTTOM_SIDE ];
1127             for ( int iE = 0; iE < botSide->NbEdges(); ++iE )
1128               if ( !myHelper->IsSubShape( botSide->Edge(iE), thePrism.myBottom ))
1129                 return toSM( error(TCom("Composite 'horizontal' edges are not supported")));
1130           }
1131           if ( faceMap.Add( face ))
1132           {
1133             setWireIndex( quadList.back(), iWire ); // for use in makeQuadsForOutInProjection()
1134             thePrism.myWallQuads.push_back( quadList );
1135           }
1136           break;
1137         }
1138       }
1139       if ( hasWallFace )
1140       {
1141         ++edge;
1142       }
1143       else // seam edge (IPAL53561)
1144       {
1145         edge = thePrism.myBottomEdges.erase( edge );
1146         --iE;
1147         --(*nbE);
1148       }
1149     }
1150     if ( iE == *nbE )
1151     {
1152       iE = 0;
1153       ++iWire;
1154       ++nbE;
1155       int nbQuadPrev = std::accumulate( nbQuadsPerWire.begin(), nbQuadsPerWire.end(), 0 );
1156       nbQuadsPerWire.push_back( thePrism.myWallQuads.size() - nbQuadPrev );
1157     }
1158   }
1159
1160   // -------------------------
1161   // Find the rest wall FACEs
1162   // -------------------------
1163
1164   // Compose a vector of indixes of right neighbour FACE for each wall FACE
1165   // that is not so evident in case of several WIREs in the bottom FACE
1166   thePrism.myRightQuadIndex.clear();
1167   for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
1168   {
1169     thePrism.myRightQuadIndex.push_back( i+1 ); // OK for all but the last EDGE of a WIRE
1170   }
1171   list< int >::iterator nbQinW = nbQuadsPerWire.begin();
1172   for ( int iLeft = 0; nbQinW != nbQuadsPerWire.end(); ++nbQinW )
1173   {
1174     thePrism.myRightQuadIndex[ iLeft + *nbQinW - 1 ] = iLeft; // for the last EDGE of a WIRE
1175     iLeft += *nbQinW;
1176   }
1177
1178   while ( totalNbFaces - faceMap.Extent() > 2 )
1179   {
1180     // find wall FACEs adjacent to each of wallQuads by the right side EDGE
1181     int nbKnownFaces;
1182     do {
1183       nbKnownFaces = faceMap.Extent();
1184       StdMeshers_FaceSidePtr rightSide, topSide; // sides of the quad
1185       for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
1186       {
1187         rightSide = thePrism.myWallQuads[i].back()->side[ QUAD_RIGHT_SIDE ];
1188         for ( int iE = 0; iE < rightSide->NbEdges(); ++iE ) // rightSide can be composite
1189         {
1190           const TopoDS_Edge & rightE = rightSide->Edge( iE );
1191           TopTools_ListIteratorOfListOfShape face( edgeToFaces.FindFromKey( rightE ));
1192           for ( ; face.More(); face.Next() )
1193             if ( faceMap.Add( face.Value() ))
1194             {
1195               // a new wall FACE encountered, store it in thePrism.myWallQuads
1196               const int iRight = thePrism.myRightQuadIndex[i];
1197               topSide = thePrism.myWallQuads[ iRight ].back()->side[ QUAD_TOP_SIDE ];
1198               const TopoDS_Edge&   newBotE = topSide->Edge(0);
1199               const TopoDS_Shape& newWallF = face.Value();
1200               thePrism.myWallQuads[ iRight ].push_back( quadAlgo->CheckNbEdges( *mesh, newWallF ));
1201               if ( !thePrism.myWallQuads[ iRight ].back() )
1202                 return toSM( error(TCom("Side face #") << shapeID( newWallF ) <<
1203                                    " not meshable with quadrangles"));
1204               if ( ! setBottomEdge( newBotE, thePrism.myWallQuads[ iRight ].back(), newWallF ))
1205                 return toSM( error(TCom("Composite 'horizontal' edges are not supported")));
1206             }
1207         }
1208       }
1209     } while ( nbKnownFaces != faceMap.Extent() );
1210
1211     // find wall FACEs adjacent to each of thePrism.myWallQuads by the top side EDGE
1212     if ( totalNbFaces - faceMap.Extent() > 2 )
1213     {
1214       const int nbFoundWalls = faceMap.Extent();
1215       for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
1216       {
1217         StdMeshers_FaceSidePtr topSide = thePrism.myWallQuads[i].back()->side[ QUAD_TOP_SIDE ];
1218         const TopoDS_Edge &       topE = topSide->Edge( 0 );
1219         if ( topSide->NbEdges() > 1 )
1220           return toSM( error(COMPERR_BAD_SHAPE, TCom("Side face #") <<
1221                              shapeID( thePrism.myWallQuads[i].back()->face )
1222                              << " has a composite top edge"));
1223         TopTools_ListIteratorOfListOfShape faceIt( edgeToFaces.FindFromKey( topE ));
1224         for ( ; faceIt.More(); faceIt.Next() )
1225           if ( faceMap.Add( faceIt.Value() ))
1226           {
1227             // a new wall FACE encountered, store it in wallQuads
1228             thePrism.myWallQuads[ i ].push_back( quadAlgo->CheckNbEdges( *mesh, faceIt.Value() ));
1229             if ( !thePrism.myWallQuads[ i ].back() )
1230               return toSM( error(TCom("Side face #") << shapeID( faceIt.Value() ) <<
1231                                  " not meshable with quadrangles"));
1232             if ( ! setBottomEdge( topE, thePrism.myWallQuads[ i ].back(), faceIt.Value() ))
1233               return toSM( error(TCom("Composite 'horizontal' edges are not supported")));
1234             if ( totalNbFaces - faceMap.Extent() == 2 )
1235             {
1236               i = thePrism.myWallQuads.size(); // to quit from the outer loop
1237               break;
1238             }
1239           }
1240       }
1241       if ( nbFoundWalls == faceMap.Extent() )
1242         return toSM( error("Failed to find wall faces"));
1243
1244     }
1245   } // while ( totalNbFaces - faceMap.Extent() > 2 )
1246
1247   // ------------------
1248   // Find the top FACE
1249   // ------------------
1250
1251   if ( thePrism.myTop.IsNull() )
1252   {
1253     // now only top and bottom FACEs are not in the faceMap
1254     faceMap.Add( thePrism.myBottom );
1255     for ( TopExp_Explorer f( thePrism.myShape3D, TopAbs_FACE ); f.More(); f.Next() )
1256       if ( !faceMap.Contains( f.Current() )) {
1257         thePrism.myTop = TopoDS::Face( f.Current() );
1258         break;
1259       }
1260     if ( thePrism.myTop.IsNull() )
1261       return toSM( error("Top face not found"));
1262   }
1263
1264   // Check that the top FACE shares all the top EDGEs
1265   for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
1266   {
1267     StdMeshers_FaceSidePtr topSide = thePrism.myWallQuads[i].back()->side[ QUAD_TOP_SIDE ];
1268     const TopoDS_Edge &       topE = topSide->Edge( 0 );
1269     if ( !myHelper->IsSubShape( topE, thePrism.myTop ))
1270       return toSM( error( TCom("Wrong source face: #") << shapeID( thePrism.myBottom )));
1271   }
1272
1273   return true;
1274 }
1275
1276 //=======================================================================
1277 //function : compute
1278 //purpose  : Compute mesh on a SOLID
1279 //=======================================================================
1280
1281 bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
1282 {
1283   myHelper->IsQuadraticSubMesh( thePrism.myShape3D );
1284   if ( _computeCanceled )
1285     return toSM( error( SMESH_ComputeError::New(COMPERR_CANCELED)));
1286
1287   // Assure the bottom is meshed
1288   if ( !computeBase( thePrism ))
1289     return false;
1290
1291   // Make all side FACEs of thePrism meshed with quads
1292   if ( !computeWalls( thePrism ))
1293     return false;
1294
1295   // Analyse mesh and geometry to find all block sub-shapes and submeshes
1296   // (after fixing IPAL52499 myBlock is used as a holder of boundary nodes
1297   // and for 2D projection in hard cases where StdMeshers_Projection_2D fails;
1298   // location of internal nodes is usually computed by StdMeshers_Sweeper)
1299   if ( !myBlock.Init( myHelper, thePrism ))
1300     return toSM( error( myBlock.GetError()));
1301
1302   SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
1303
1304   int volumeID = meshDS->ShapeToIndex( thePrism.myShape3D );
1305
1306   // Try to get gp_Trsf to get all nodes from bottom ones
1307   vector<gp_Trsf> trsf;
1308   gp_Trsf bottomToTopTrsf;
1309   // if ( !myBlock.GetLayersTransformation( trsf, thePrism ))
1310   //   trsf.clear();
1311   // else if ( !trsf.empty() )
1312   //   bottomToTopTrsf = trsf.back();
1313
1314   // To compute coordinates of a node inside a block using "block approach",
1315   // it is necessary to know
1316   // 1. normalized parameters of the node by which
1317   // 2. coordinates of node projections on all block sub-shapes are computed
1318
1319   // So we fill projections on vertices at once as they are same for all nodes
1320   myShapeXYZ.resize( myBlock.NbSubShapes() );
1321   for ( int iV = SMESH_Block::ID_FirstV; iV < SMESH_Block::ID_FirstE; ++iV ) {
1322     myBlock.VertexPoint( iV, myShapeXYZ[ iV ]);
1323     SHOWYXZ("V point " <<iV << " ", myShapeXYZ[ iV ]);
1324   }
1325
1326   // Projections on the top and bottom faces are taken from nodes existing
1327   // on these faces; find correspondence between bottom and top nodes
1328   myUseBlock = false; // is set to true if projection is done using "block approach"
1329   myBotToColumnMap.clear();
1330   if ( !assocOrProjBottom2Top( bottomToTopTrsf, thePrism ) ) // it also fills myBotToColumnMap
1331     return false;
1332
1333   // If all "vertical" EDGEs are straight, then all nodes of an internal node column
1334   // are located on a line connecting the top node and the bottom node.
1335   bool isStrightColunm = allVerticalEdgesStraight( thePrism );
1336   if ( isStrightColunm )
1337     myUseBlock = false;
1338
1339   // Create nodes inside the block
1340
1341   if ( !myUseBlock )
1342   {
1343     // use transformation (issue 0020680, IPAL0052499) or a "straight line" approach
1344     StdMeshers_Sweeper sweeper;
1345     sweeper.myHelper  = myHelper;
1346     sweeper.myBotFace = thePrism.myBottom;
1347     sweeper.myTopFace = thePrism.myTop;
1348
1349     // load boundary nodes into sweeper
1350     bool dummy;
1351     std::set< const SMDS_MeshNode* > usedEndNodes;
1352     list< TopoDS_Edge >::const_iterator edge = thePrism.myBottomEdges.begin();
1353     for ( ; edge != thePrism.myBottomEdges.end(); ++edge )
1354     {
1355       int edgeID = meshDS->ShapeToIndex( *edge );
1356       TParam2ColumnMap* u2col = const_cast<TParam2ColumnMap*>
1357         ( myBlock.GetParam2ColumnMap( edgeID, dummy ));
1358
1359       TParam2ColumnMap::iterator u2colIt = u2col->begin(), u2colEnd = u2col->end();
1360       const SMDS_MeshNode* n0 = u2colIt->second[0];
1361       const SMDS_MeshNode* n1 = u2col->rbegin()->second[0];
1362       if ( !usedEndNodes.insert ( n0 ).second ) ++u2colIt;
1363       if ( !usedEndNodes.insert ( n1 ).second ) --u2colEnd;
1364
1365       for ( ; u2colIt != u2colEnd; ++u2colIt )
1366         sweeper.myBndColumns.push_back( & u2colIt->second );
1367     }
1368     // load node columns inside the bottom FACE
1369     sweeper.myIntColumns.reserve( myBotToColumnMap.size() );
1370     TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin();
1371     for ( ; bot_column != myBotToColumnMap.end(); ++bot_column )
1372       sweeper.myIntColumns.push_back( & bot_column->second );
1373
1374     myHelper->SetElementsOnShape( true );
1375
1376     if ( !isStrightColunm )
1377     {
1378       double tol = getSweepTolerance( thePrism );
1379       bool allowHighBndError = !isSimpleBottom( thePrism );
1380       myUseBlock = !sweeper.ComputeNodesByTrsf( tol, allowHighBndError );
1381     }
1382     else if ( sweeper.CheckSameZ() )
1383     {
1384       myUseBlock = !sweeper.ComputeNodesOnStraightSameZ();
1385     }
1386     else
1387     {
1388       myUseBlock = !sweeper.ComputeNodesOnStraight();
1389     }
1390     myHelper->SetElementsOnShape( false );
1391   }
1392
1393   if ( myUseBlock ) // use block approach
1394   {
1395     // loop on nodes inside the bottom face
1396     Prism_3D::TNode prevBNode;
1397     TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin();
1398     for ( ; bot_column != myBotToColumnMap.end(); ++bot_column )
1399     {
1400       const Prism_3D::TNode& tBotNode = bot_column->first; // bottom TNode
1401       if ( tBotNode.GetPositionType() != SMDS_TOP_FACE &&
1402            myBlock.HasNodeColumn( tBotNode.myNode ))
1403         continue; // node is not inside the FACE
1404
1405       // column nodes; middle part of the column are zero pointers
1406       TNodeColumn& column = bot_column->second;
1407
1408       // check if a column is already computed using non-block approach
1409       size_t i;
1410       for ( i = 0; i < column.size(); ++i )
1411         if ( !column[ i ])
1412           break;
1413       if ( i == column.size() )
1414         continue; // all nodes created
1415
1416       gp_XYZ botParams, topParams;
1417       if ( !tBotNode.HasParams() )
1418       {
1419         // compute bottom node parameters
1420         gp_XYZ paramHint(-1,-1,-1);
1421         if ( prevBNode.IsNeighbor( tBotNode ))
1422           paramHint = prevBNode.GetParams();
1423         if ( !myBlock.ComputeParameters( tBotNode.GetCoords(), tBotNode.ChangeParams(),
1424                                          ID_BOT_FACE, paramHint ))
1425           return toSM( error(TCom("Can't compute normalized parameters for node ")
1426                              << tBotNode.myNode->GetID() << " on the face #"
1427                              << myBlock.SubMesh( ID_BOT_FACE )->GetId() ));
1428         prevBNode = tBotNode;
1429
1430         botParams = topParams = tBotNode.GetParams();
1431         topParams.SetZ( 1 );
1432
1433         // compute top node parameters
1434         if ( column.size() > 2 ) {
1435           gp_Pnt topCoords = gpXYZ( column.back() );
1436           if ( !myBlock.ComputeParameters( topCoords, topParams, ID_TOP_FACE, topParams ))
1437             return toSM( error(TCom("Can't compute normalized parameters ")
1438                                << "for node " << column.back()->GetID()
1439                                << " on the face #"<< column.back()->getshapeId() ));
1440         }
1441       }
1442       else // top nodes are created by projection using parameters
1443       {
1444         botParams = topParams = tBotNode.GetParams();
1445         topParams.SetZ( 1 );
1446       }
1447
1448       myShapeXYZ[ ID_BOT_FACE ] = tBotNode.GetCoords();
1449       myShapeXYZ[ ID_TOP_FACE ] = gpXYZ( column.back() );
1450
1451       // vertical loop
1452       TNodeColumn::iterator columnNodes = column.begin();
1453       for ( int z = 0; columnNodes != column.end(); ++columnNodes, ++z)
1454       {
1455         const SMDS_MeshNode* & node = *columnNodes;
1456         if ( node ) continue; // skip bottom or top node
1457
1458         // params of a node to create
1459         double rz = (double) z / (double) ( column.size() - 1 );
1460         gp_XYZ params = botParams * ( 1 - rz ) + topParams * rz;
1461
1462         // set coords on all faces and nodes
1463         const int nbSideFaces = 4;
1464         int sideFaceIDs[nbSideFaces] = { SMESH_Block::ID_Fx0z,
1465                                          SMESH_Block::ID_Fx1z,
1466                                          SMESH_Block::ID_F0yz,
1467                                          SMESH_Block::ID_F1yz };
1468         for ( int iF = 0; iF < nbSideFaces; ++iF )
1469           if ( !setFaceAndEdgesXYZ( sideFaceIDs[ iF ], params, z ))
1470             return false;
1471
1472         // compute coords for a new node
1473         gp_XYZ coords;
1474         if ( !SMESH_Block::ShellPoint( params, myShapeXYZ, coords ))
1475           return toSM( error("Can't compute coordinates by normalized parameters"));
1476
1477         // if ( !meshDS->MeshElements( volumeID ) ||
1478         //      meshDS->MeshElements( volumeID )->NbNodes() == 0 )
1479         //   pointsToPython(myShapeXYZ);
1480         SHOWYXZ("TOPFacePoint ",myShapeXYZ[ ID_TOP_FACE]);
1481         SHOWYXZ("BOT Node "<< tBotNode.myNode->GetID(),gpXYZ(tBotNode.myNode));
1482         SHOWYXZ("ShellPoint ",coords);
1483
1484         // create a node
1485         node = meshDS->AddNode( coords.X(), coords.Y(), coords.Z() );
1486         meshDS->SetNodeInVolume( node, volumeID );
1487
1488         if ( _computeCanceled )
1489           return false;
1490       }
1491     } // loop on bottom nodes
1492   }
1493
1494   // Create volumes
1495
1496   SMESHDS_SubMesh* smDS = myBlock.SubMeshDS( ID_BOT_FACE );
1497   if ( !smDS ) return toSM( error(COMPERR_BAD_INPUT_MESH, "Null submesh"));
1498
1499   // loop on bottom mesh faces
1500   vector< const TNodeColumn* > columns;
1501   SMDS_ElemIteratorPtr faceIt = smDS->GetElements();
1502   while ( faceIt->more() )
1503   {
1504     const SMDS_MeshElement* face = faceIt->next();
1505     if ( !face || face->GetType() != SMDSAbs_Face )
1506       continue;
1507
1508     // find node columns for each node
1509     int nbNodes = face->NbCornerNodes();
1510     columns.resize( nbNodes );
1511     for ( int i = 0; i < nbNodes; ++i )
1512     {
1513       const SMDS_MeshNode* n = face->GetNode( i );
1514       columns[ i ] = NULL;
1515
1516       if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE )
1517         columns[ i ] = myBlock.GetNodeColumn( n );
1518
1519       if ( !columns[ i ] )
1520       {
1521         TNode2ColumnMap::iterator bot_column = myBotToColumnMap.find( n );
1522         if ( bot_column == myBotToColumnMap.end() )
1523           return toSM( error(TCom("No side nodes found above node ") << n->GetID() ));
1524         columns[ i ] = & bot_column->second;
1525       }
1526     }
1527     // create prisms
1528     if ( !AddPrisms( columns, myHelper ))
1529       return toSM( error("Different 'vertical' discretization"));
1530
1531   } // loop on bottom mesh faces
1532
1533   // clear data
1534   myBotToColumnMap.clear();
1535   myBlock.Clear();
1536
1537   // update state of sub-meshes (mostly in order to erase improper errors)
1538   SMESH_subMesh* sm = myHelper->GetMesh()->GetSubMesh( thePrism.myShape3D );
1539   SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true);
1540   while ( smIt->more() )
1541   {
1542     sm = smIt->next();
1543     sm->GetComputeError().reset();
1544     sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1545   }
1546
1547   return true;
1548 }
1549
1550 //=======================================================================
1551 //function : computeBase
1552 //purpose  : Compute the base face of a prism
1553 //=======================================================================
1554
1555 bool StdMeshers_Prism_3D::computeBase(const Prism_3D::TPrismTopo& thePrism)
1556 {
1557   SMESH_Mesh*     mesh = myHelper->GetMesh();
1558   SMESH_subMesh* botSM = mesh->GetSubMesh( thePrism.myBottom );
1559   if (( botSM->IsEmpty() ) &&
1560       ( ! botSM->GetAlgo() ||
1561         ! _gen->Compute( *botSM->GetFather(), botSM->GetSubShape(), /*shapeOnly=*/true )))
1562   {
1563     // find any applicable algorithm assigned to any FACE of the main shape
1564     std::vector< TopoDS_Shape > faces;
1565     if ( myPrevBottomSM &&
1566          myPrevBottomSM->GetAlgo()->IsApplicableToShape( thePrism.myBottom, /*all=*/false ))
1567       faces.push_back( myPrevBottomSM->GetSubShape() );
1568
1569     TopExp_Explorer faceIt( mesh->GetShapeToMesh(), TopAbs_FACE );
1570     for ( ; faceIt.More(); faceIt.Next() )
1571       faces.push_back( faceIt.Current() );
1572
1573     faces.push_back( TopoDS_Shape() ); // to try quadrangle algorithm
1574
1575     SMESH_Algo* algo = 0;
1576     for ( size_t i = 0; i < faces.size() &&  botSM->IsEmpty(); ++i )
1577     {
1578       if ( faces[i].IsNull() ) algo = TQuadrangleAlgo::instance( this, myHelper );
1579       else                     algo = mesh->GetSubMesh( faces[i] )->GetAlgo();
1580       if ( algo && algo->IsApplicableToShape( thePrism.myBottom, /*all=*/false ))
1581       {
1582         // try to compute the bottom FACE
1583         if ( algo->NeedDiscreteBoundary() )
1584         {
1585           // compute sub-shapes
1586           SMESH_subMeshIteratorPtr smIt = botSM->getDependsOnIterator(false,false);
1587           bool subOK = true;
1588           while ( smIt->more() && subOK )
1589           {
1590             SMESH_subMesh* sub = smIt->next();
1591             sub->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1592             subOK = sub->IsMeshComputed();
1593           }
1594           if ( !subOK )
1595             continue;
1596         }
1597         try {
1598           OCC_CATCH_SIGNALS;
1599           algo->InitComputeError();
1600           algo->Compute( *mesh, botSM->GetSubShape() );
1601         }
1602         catch (...) {
1603         }
1604       }
1605     }
1606   }
1607
1608   if ( botSM->IsEmpty() )
1609     return error( COMPERR_BAD_INPUT_MESH,
1610                   TCom( "No mesher defined to compute the base face #")
1611                   << shapeID( thePrism.myBottom ));
1612
1613   if ( botSM->GetAlgo() )
1614     myPrevBottomSM = botSM;
1615
1616   return true;
1617 }
1618
1619 //=======================================================================
1620 //function : computeWalls
1621 //purpose  : Compute 2D mesh on walls FACEs of a prism
1622 //=======================================================================
1623
1624 bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism)
1625 {
1626   SMESH_Mesh*     mesh = myHelper->GetMesh();
1627   SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
1628   DBGOUT( endl << "COMPUTE Prism " << meshDS->ShapeToIndex( thePrism.myShape3D ));
1629
1630   TProjction1dAlgo* projector1D = TProjction1dAlgo::instance( this );
1631   TQuadrangleAlgo*     quadAlgo = TQuadrangleAlgo::instance( this, myHelper );
1632
1633   // SMESH_HypoFilter hyp1dFilter( SMESH_HypoFilter::IsAlgo(),/*not=*/true);
1634   // hyp1dFilter.And( SMESH_HypoFilter::HasDim( 1 ));
1635   // hyp1dFilter.And( SMESH_HypoFilter::IsMoreLocalThan( thePrism.myShape3D, *mesh ));
1636
1637   // Discretize equally 'vertical' EDGEs
1638   // -----------------------------------
1639   // find source FACE sides for projection: either already computed ones or
1640   // the 'most composite' ones
1641   const size_t nbWalls = thePrism.myWallQuads.size();
1642   vector< int > wgt( nbWalls, 0 ); // "weight" of a wall
1643   for ( size_t iW = 0; iW != nbWalls; ++iW )
1644   {
1645     Prism_3D::TQuadList::const_iterator quad = thePrism.myWallQuads[iW].begin();
1646     for ( ; quad != thePrism.myWallQuads[iW].end(); ++quad )
1647     {
1648       StdMeshers_FaceSidePtr lftSide = (*quad)->side[ QUAD_LEFT_SIDE ];
1649       lftSide->Reverse(); // to go up
1650       for ( int i = 0; i < lftSide->NbEdges(); ++i )
1651       {
1652         ++wgt[ iW ];
1653         const TopoDS_Edge& E = lftSide->Edge(i);
1654         if ( mesh->GetSubMesh( E )->IsMeshComputed() )
1655         {
1656           wgt[ iW ] += 100;
1657           wgt[ myHelper->WrapIndex( iW+1, nbWalls)] += 10;
1658           wgt[ myHelper->WrapIndex( iW-1, nbWalls)] += 10;
1659         }
1660         // else if ( mesh->GetHypothesis( E, hyp1dFilter, true )) // local hypothesis!
1661         //   wgt += 100;
1662       }
1663     }
1664     // in quadratic mesh, pass ignoreMediumNodes to quad sides
1665     if ( myHelper->GetIsQuadratic() )
1666     {
1667       quad = thePrism.myWallQuads[iW].begin();
1668       for ( ; quad != thePrism.myWallQuads[iW].end(); ++quad )
1669         for ( int i = 0; i < NB_QUAD_SIDES; ++i )
1670           (*quad)->side[ i ].grid->SetIgnoreMediumNodes( true );
1671     }
1672   }
1673   multimap< int, int > wgt2quad;
1674   for ( size_t iW = 0; iW != nbWalls; ++iW )
1675     wgt2quad.insert( make_pair( wgt[ iW ], iW ));
1676
1677   // artificial quads to do outer <-> inner wall projection
1678   std::map< int, FaceQuadStruct > iW2oiQuads;
1679   std::map< int, FaceQuadStruct >::iterator w2oiq;
1680   makeQuadsForOutInProjection( thePrism, wgt2quad, iW2oiQuads );
1681
1682   // Project 'vertical' EDGEs, from left to right
1683   multimap< int, int >::reverse_iterator w2q = wgt2quad.rbegin();
1684   for ( ; w2q != wgt2quad.rend(); ++w2q )
1685   {
1686     const int iW = w2q->second;
1687     const Prism_3D::TQuadList&         quads = thePrism.myWallQuads[ iW ];
1688     Prism_3D::TQuadList::const_iterator quad = quads.begin();
1689     for ( ; quad != quads.end(); ++quad )
1690     {
1691       StdMeshers_FaceSidePtr rgtSide = (*quad)->side[ QUAD_RIGHT_SIDE ]; // tgt
1692       StdMeshers_FaceSidePtr lftSide = (*quad)->side[ QUAD_LEFT_SIDE ];  // src
1693       bool swapLeftRight = ( lftSide->NbSegments( /*update=*/true ) == 0 &&
1694                              rgtSide->NbSegments( /*update=*/true )  > 0 );
1695       if ( swapLeftRight )
1696         std::swap( lftSide, rgtSide );
1697
1698       bool isArtificialQuad = (( w2oiq = iW2oiQuads.find( iW )) != iW2oiQuads.end() );
1699       if ( isArtificialQuad )
1700       {
1701         // reset sides to perform the outer <-> inner projection
1702         FaceQuadStruct& oiQuad = w2oiq->second;
1703         rgtSide = oiQuad.side[ QUAD_RIGHT_SIDE ];
1704         lftSide = oiQuad.side[ QUAD_LEFT_SIDE ];
1705         iW2oiQuads.erase( w2oiq );
1706       }
1707
1708       // assure that all the source (left) EDGEs are meshed
1709       int nbSrcSegments = 0;
1710       for ( int i = 0; i < lftSide->NbEdges(); ++i )
1711       {
1712         if ( isArtificialQuad )
1713         {
1714           nbSrcSegments = lftSide->NbPoints()-1;
1715           continue;
1716         }
1717         const TopoDS_Edge& srcE = lftSide->Edge(i);
1718         SMESH_subMesh*    srcSM = mesh->GetSubMesh( srcE );
1719         if ( !srcSM->IsMeshComputed() ) {
1720           DBGOUT( "COMPUTE V edge " << srcSM->GetId() );
1721           TopoDS_Edge prpgSrcE = findPropagationSource( srcE );
1722           if ( !prpgSrcE.IsNull() ) {
1723             srcSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE );
1724             projector1D->myHyp.SetSourceEdge( prpgSrcE );
1725             projector1D->Compute( *mesh, srcE );
1726             srcSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1727           }
1728           else {
1729             srcSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE );
1730             srcSM->ComputeStateEngine       ( SMESH_subMesh::COMPUTE );
1731           }
1732           if ( !srcSM->IsMeshComputed() )
1733             return toSM( error( "Can't compute 1D mesh" ));
1734         }
1735         nbSrcSegments += srcSM->GetSubMeshDS()->NbElements();
1736       }
1737       // check target EDGEs
1738       int nbTgtMeshed = 0, nbTgtSegments = 0;
1739       vector< bool > isTgtEdgeComputed( rgtSide->NbEdges() );
1740       for ( int i = 0; i < rgtSide->NbEdges(); ++i )
1741       {
1742         const TopoDS_Edge& tgtE = rgtSide->Edge(i);
1743         SMESH_subMesh*    tgtSM = mesh->GetSubMesh( tgtE );
1744         if ( !( isTgtEdgeComputed[ i ] = tgtSM->IsMeshComputed() )) {
1745           tgtSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE );
1746           tgtSM->ComputeStateEngine       ( SMESH_subMesh::COMPUTE );
1747         }
1748         if ( tgtSM->IsMeshComputed() ) {
1749           ++nbTgtMeshed;
1750           nbTgtSegments += tgtSM->GetSubMeshDS()->NbElements();
1751         }
1752       }
1753       if ( rgtSide->NbEdges() == nbTgtMeshed ) // all tgt EDGEs meshed
1754       {
1755         if ( nbTgtSegments != nbSrcSegments )
1756         {
1757           bool badMeshRemoved = false;
1758           // remove just computed segments
1759           for ( int i = 0; i < rgtSide->NbEdges(); ++i )
1760             if ( !isTgtEdgeComputed[ i ])
1761             {
1762               const TopoDS_Edge& tgtE = rgtSide->Edge(i);
1763               SMESH_subMesh*    tgtSM = mesh->GetSubMesh( tgtE );
1764               tgtSM->ComputeStateEngine( SMESH_subMesh::CLEAN );
1765               badMeshRemoved = true;
1766               nbTgtMeshed--;
1767             }
1768           if ( !badMeshRemoved )
1769           {
1770             for ( int i = 0; i < lftSide->NbEdges(); ++i )
1771               addBadInputElements( meshDS->MeshElements( lftSide->Edge( i )));
1772             for ( int i = 0; i < rgtSide->NbEdges(); ++i )
1773               addBadInputElements( meshDS->MeshElements( rgtSide->Edge( i )));
1774             return toSM( error( TCom("Different nb of segment on logically vertical edges #")
1775                                 << shapeID( lftSide->Edge(0) ) << " and #"
1776                                 << shapeID( rgtSide->Edge(0) ) << ": "
1777                                 << nbSrcSegments << " != " << nbTgtSegments ));
1778           }
1779         }
1780         else // if ( nbTgtSegments == nbSrcSegments )
1781         {
1782           continue;
1783         }
1784       }
1785       // Compute 'vertical projection'
1786       if ( nbTgtMeshed == 0 )
1787       {
1788         // compute nodes on target VERTEXes
1789         const UVPtStructVec&  srcNodeStr = lftSide->GetUVPtStruct();
1790         if ( srcNodeStr.size() == 0 )
1791           return toSM( error( TCom("Invalid node positions on edge #") <<
1792                               lftSide->EdgeID(0) ));
1793         vector< SMDS_MeshNode* > newNodes( srcNodeStr.size() );
1794         for ( int is2ndV = 0; is2ndV < 2; ++is2ndV )
1795         {
1796           const TopoDS_Edge& E = rgtSide->Edge( is2ndV ? rgtSide->NbEdges()-1 : 0 );
1797           TopoDS_Vertex      v = myHelper->IthVertex( is2ndV, E );
1798           mesh->GetSubMesh( v )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1799           const SMDS_MeshNode* n = SMESH_Algo::VertexNode( v, meshDS );
1800           newNodes[ is2ndV ? newNodes.size()-1 : 0 ] = (SMDS_MeshNode*) n;
1801         }
1802
1803         // compute nodes on target EDGEs
1804         DBGOUT( "COMPUTE V edge (proj) " << shapeID( lftSide->Edge(0)));
1805         //rgtSide->Reverse(); // direct it same as the lftSide
1806         myHelper->SetElementsOnShape( false ); // myHelper holds the prism shape
1807         TopoDS_Edge tgtEdge;
1808         for ( size_t iN = 1; iN < srcNodeStr.size()-1; ++iN ) // add nodes
1809         {
1810           gp_Pnt       p = rgtSide->Value3d  ( srcNodeStr[ iN ].normParam );
1811           double       u = rgtSide->Parameter( srcNodeStr[ iN ].normParam, tgtEdge );
1812           newNodes[ iN ] = meshDS->AddNode( p.X(), p.Y(), p.Z() );
1813           meshDS->SetNodeOnEdge( newNodes[ iN ], tgtEdge, u );
1814         }
1815         for ( size_t iN = 1; iN < srcNodeStr.size(); ++iN ) // add segments
1816         {
1817           // find an EDGE to set a new segment
1818           std::pair<int, TopAbs_ShapeEnum> id2type =
1819             myHelper->GetMediumPos( newNodes[ iN-1 ], newNodes[ iN ] );
1820           if ( id2type.second != TopAbs_EDGE )
1821           {
1822             // new nodes are on different EDGEs; put one of them on VERTEX
1823             const int      edgeIndex = rgtSide->EdgeIndex( srcNodeStr[ iN-1 ].normParam );
1824             const double vertexParam = rgtSide->LastParameter( edgeIndex );
1825             TopoDS_Vertex     vertex = rgtSide->LastVertex( edgeIndex );
1826             const SMDS_MeshNode*  vn = SMESH_Algo::VertexNode( vertex, meshDS );
1827             const gp_Pnt           p = BRep_Tool::Pnt( vertex );
1828             const int         isPrev = ( Abs( srcNodeStr[ iN-1 ].normParam - vertexParam ) <
1829                                          Abs( srcNodeStr[ iN   ].normParam - vertexParam ));
1830             meshDS->UnSetNodeOnShape( newNodes[ iN-isPrev ] );
1831             meshDS->SetNodeOnVertex ( newNodes[ iN-isPrev ], vertex );
1832             meshDS->MoveNode        ( newNodes[ iN-isPrev ], p.X(), p.Y(), p.Z() );
1833             id2type.first = newNodes[ iN-(1-isPrev) ]->getshapeId();
1834             if ( vn )
1835             {
1836               SMESH_MeshEditor::TListOfListOfNodes lln( 1, list< const SMDS_MeshNode* >() );
1837               lln.back().push_back ( vn );
1838               lln.back().push_front( newNodes[ iN-isPrev ] ); // to keep
1839               SMESH_MeshEditor( mesh ).MergeNodes( lln );
1840             }
1841           }
1842           SMDS_MeshElement* newEdge = myHelper->AddEdge( newNodes[ iN-1 ], newNodes[ iN ] );
1843           meshDS->SetMeshElementOnShape( newEdge, id2type.first );
1844         }
1845         myHelper->SetElementsOnShape( true );
1846         for ( int i = 0; i < rgtSide->NbEdges(); ++i ) // update state of sub-meshes
1847         {
1848           const TopoDS_Edge& E = rgtSide->Edge( i );
1849           SMESH_subMesh* tgtSM = mesh->GetSubMesh( E );
1850           tgtSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1851         }
1852
1853         // to continue projection from the just computed side as a source
1854         if ( !swapLeftRight && rgtSide->NbEdges() > 1 && w2q->second == iW )
1855         {
1856           std::pair<int,int> wgt2quadKeyVal( w2q->first + 1, thePrism.myRightQuadIndex[ iW ]);
1857           wgt2quad.insert( wgt2quadKeyVal ); // it will be skipped by ++w2q
1858           wgt2quad.insert( wgt2quadKeyVal );
1859           w2q = wgt2quad.rbegin();
1860         }
1861       }
1862       else
1863       {
1864         // HOPE assigned hypotheses are OK, so that equal nb of segments will be generated
1865         //return toSM( error("Partial projection not implemented"));
1866       }
1867     } // loop on quads of a composite wall side
1868   } // loop on the ordered wall sides
1869
1870
1871
1872   for ( size_t iW = 0; iW != thePrism.myWallQuads.size(); ++iW )
1873   {
1874     Prism_3D::TQuadList::const_iterator quad = thePrism.myWallQuads[iW].begin();
1875     for ( ; quad != thePrism.myWallQuads[iW].end(); ++quad )
1876     {
1877       const TopoDS_Face& face = (*quad)->face;
1878       SMESH_subMesh*      fSM = mesh->GetSubMesh( face );
1879       if ( ! fSM->IsMeshComputed() )
1880       {
1881         // Top EDGEs must be projections from the bottom ones
1882         // to compute structured quad mesh on wall FACEs
1883         // ---------------------------------------------------
1884         const TopoDS_Edge& botE = (*quad)->side[ QUAD_BOTTOM_SIDE ].grid->Edge(0);
1885         const TopoDS_Edge& topE = (*quad)->side[ QUAD_TOP_SIDE    ].grid->Edge(0);
1886         SMESH_subMesh*    botSM = mesh->GetSubMesh( botE );
1887         SMESH_subMesh*    topSM = mesh->GetSubMesh( topE );
1888         SMESH_subMesh*    srcSM = botSM;
1889         SMESH_subMesh*    tgtSM = topSM;
1890         srcSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1891         tgtSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1892         if ( !srcSM->IsMeshComputed() && tgtSM->IsMeshComputed() )
1893           std::swap( srcSM, tgtSM );
1894
1895         if ( !srcSM->IsMeshComputed() )
1896         {
1897           DBGOUT( "COMPUTE H edge " << srcSM->GetId());
1898           srcSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE ); // nodes on VERTEXes
1899           srcSM->ComputeStateEngine( SMESH_subMesh::COMPUTE );        // segments on the EDGE
1900         }
1901
1902         if ( tgtSM->IsMeshComputed() &&
1903              tgtSM->GetSubMeshDS()->NbNodes() != srcSM->GetSubMeshDS()->NbNodes() )
1904         {
1905           // the top EDGE is computed differently than the bottom one,
1906           // try to clear a wrong mesh
1907           bool isAdjFaceMeshed = false;
1908           PShapeIteratorPtr fIt = myHelper->GetAncestors( tgtSM->GetSubShape(),
1909                                                           *mesh, TopAbs_FACE );
1910           while ( const TopoDS_Shape* f = fIt->next() )
1911             if (( isAdjFaceMeshed = mesh->GetSubMesh( *f )->IsMeshComputed() ))
1912               break;
1913           if ( isAdjFaceMeshed )
1914             return toSM( error( TCom("Different nb of segment on logically horizontal edges #")
1915                                 << shapeID( botE ) << " and #"
1916                                 << shapeID( topE ) << ": "
1917                                 << tgtSM->GetSubMeshDS()->NbElements() << " != "
1918                                 << srcSM->GetSubMeshDS()->NbElements() ));
1919           tgtSM->ComputeStateEngine( SMESH_subMesh::CLEAN );
1920         }
1921         if ( !tgtSM->IsMeshComputed() )
1922         {
1923           // compute nodes on VERTEXes
1924           SMESH_subMeshIteratorPtr smIt = tgtSM->getDependsOnIterator(/*includeSelf=*/false);
1925           while ( smIt->more() )
1926             smIt->next()->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1927           // project segments
1928           DBGOUT( "COMPUTE H edge (proj) " << tgtSM->GetId());
1929           projector1D->myHyp.SetSourceEdge( TopoDS::Edge( srcSM->GetSubShape() ));
1930           projector1D->InitComputeError();
1931           bool ok = projector1D->Compute( *mesh, tgtSM->GetSubShape() );
1932           if ( !ok )
1933           {
1934             SMESH_ComputeErrorPtr err = projector1D->GetComputeError();
1935             if ( err->IsOK() ) err->myName = COMPERR_ALGO_FAILED;
1936             tgtSM->GetComputeError() = err;
1937             return false;
1938           }
1939         }
1940         tgtSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1941
1942
1943         // Compute quad mesh on wall FACEs
1944         // -------------------------------
1945
1946         // make all EDGES meshed
1947         fSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE );
1948         if ( !fSM->SubMeshesComputed() )
1949           return toSM( error( COMPERR_BAD_INPUT_MESH,
1950                               "Not all edges have valid algorithm and hypothesis"));
1951         // mesh the <face>
1952         quadAlgo->InitComputeError();
1953         DBGOUT( "COMPUTE Quad face " << fSM->GetId());
1954         bool ok = quadAlgo->Compute( *mesh, face );
1955         fSM->GetComputeError() = quadAlgo->GetComputeError();
1956         if ( !ok )
1957           return false;
1958         fSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1959       }
1960       if ( myHelper->GetIsQuadratic() )
1961       {
1962         // fill myHelper with medium nodes built by quadAlgo
1963         for ( SMDS_ElemIteratorPtr fIt = fSM->GetSubMeshDS()->GetElements(); fIt->more(); )
1964           myHelper->AddTLinks( SMDS_Mesh::DownCast<SMDS_MeshFace>( fIt->next() ));
1965       }
1966     }
1967   }
1968
1969   return true;
1970 }
1971
1972 //=======================================================================
1973 //function : findPropagationSource
1974 //purpose  : Returns a source EDGE of propagation to a given EDGE
1975 //=======================================================================
1976
1977 TopoDS_Edge StdMeshers_Prism_3D::findPropagationSource( const TopoDS_Edge& E )
1978 {
1979   if ( myPropagChains )
1980     for ( size_t i = 0; !myPropagChains[i].IsEmpty(); ++i )
1981       if ( myPropagChains[i].Contains( E ))
1982         return TopoDS::Edge( myPropagChains[i].FindKey( 1 ));
1983
1984   return TopoDS_Edge();
1985 }
1986
1987 //=======================================================================
1988 //function : makeQuadsForOutInProjection
1989 //purpose  : Create artificial wall quads for vertical projection between
1990 //           the outer and inner walls
1991 //=======================================================================
1992
1993 void StdMeshers_Prism_3D::makeQuadsForOutInProjection( const Prism_3D::TPrismTopo& thePrism,
1994                                                        multimap< int, int >&       wgt2quad,
1995                                                        map< int, FaceQuadStruct >& iQ2oiQuads)
1996 {
1997   if ( thePrism.NbWires() <= 1 )
1998     return;
1999
2000   std::set< int > doneWires; // processed wires
2001
2002   SMESH_Mesh*      mesh = myHelper->GetMesh();
2003   const bool  isForward = true;
2004   const bool skipMedium = myHelper->GetIsQuadratic();
2005
2006   // make a source side for all projections
2007
2008   multimap< int, int >::reverse_iterator w2q = wgt2quad.rbegin();
2009   const int iQuad = w2q->second;
2010   const int iWire = getWireIndex( thePrism.myWallQuads[ iQuad ].front() );
2011   doneWires.insert( iWire );
2012
2013   UVPtStructVec srcNodes;
2014
2015   Prism_3D::TQuadList::const_iterator quad = thePrism.myWallQuads[ iQuad ].begin();
2016   for ( ; quad != thePrism.myWallQuads[ iQuad ].end(); ++quad )
2017   {
2018     StdMeshers_FaceSidePtr lftSide = (*quad)->side[ QUAD_LEFT_SIDE ];
2019
2020     // assure that all the source (left) EDGEs are meshed
2021     for ( int i = 0; i < lftSide->NbEdges(); ++i )
2022     {
2023       const TopoDS_Edge& srcE = lftSide->Edge(i);
2024       SMESH_subMesh*    srcSM = mesh->GetSubMesh( srcE );
2025       if ( !srcSM->IsMeshComputed() ) {
2026         srcSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE );
2027         srcSM->ComputeStateEngine       ( SMESH_subMesh::COMPUTE );
2028       }
2029       if ( !srcSM->IsMeshComputed() )
2030         return;
2031     }
2032     const UVPtStructVec& subNodes = lftSide->GetUVPtStruct();
2033     UVPtStructVec::const_iterator subBeg = subNodes.begin(), subEnd = subNodes.end();
2034     if ( !srcNodes.empty() ) ++subBeg;
2035     srcNodes.insert( srcNodes.end(), subBeg, subEnd );
2036   }
2037   StdMeshers_FaceSidePtr srcSide = StdMeshers_FaceSide::New( srcNodes );
2038
2039   // make the quads
2040
2041   list< TopoDS_Edge > sideEdges;
2042   TopoDS_Face face;
2043   for ( ++w2q; w2q != wgt2quad.rend(); ++w2q )
2044   {
2045     const int                  iQuad = w2q->second;
2046     const Prism_3D::TQuadList& quads = thePrism.myWallQuads[ iQuad ];
2047     const int                  iWire = getWireIndex( quads.front() );
2048     if ( !doneWires.insert( iWire ).second )
2049       continue;
2050
2051     sideEdges.clear();
2052     for ( quad = quads.begin(); quad != quads.end(); ++quad )
2053     {
2054       StdMeshers_FaceSidePtr lftSide = (*quad)->side[ QUAD_LEFT_SIDE ];
2055       for ( int i = 0; i < lftSide->NbEdges(); ++i )
2056         sideEdges.push_back( lftSide->Edge( i ));
2057       face = lftSide->Face();
2058     }
2059     StdMeshers_FaceSidePtr tgtSide =
2060       StdMeshers_FaceSide::New( face, sideEdges, mesh, isForward, skipMedium, myHelper );
2061
2062     FaceQuadStruct& newQuad = iQ2oiQuads[ iQuad ];
2063     newQuad.side.resize( 4 );
2064     newQuad.side[ QUAD_LEFT_SIDE  ] = srcSide;
2065     newQuad.side[ QUAD_RIGHT_SIDE ] = tgtSide;
2066
2067     wgt2quad.insert( *w2q ); // to process this quad after processing the newQuad
2068   }
2069 }
2070
2071 //=======================================================================
2072 //function : Evaluate
2073 //purpose  :
2074 //=======================================================================
2075
2076 bool StdMeshers_Prism_3D::Evaluate(SMESH_Mesh&         theMesh,
2077                                    const TopoDS_Shape& theShape,
2078                                    MapShapeNbElems&    aResMap)
2079 {
2080   if ( theShape.ShapeType() == TopAbs_COMPOUND )
2081   {
2082     bool ok = true;
2083     for ( TopoDS_Iterator it( theShape ); it.More(); it.Next() )
2084       ok &= Evaluate( theMesh, it.Value(), aResMap );
2085     return ok;
2086   }
2087   SMESH_MesherHelper helper( theMesh );
2088   myHelper = &helper;
2089   myHelper->SetSubShape( theShape );
2090
2091   // find face contains only triangles
2092   vector < SMESH_subMesh * >meshFaces;
2093   TopTools_SequenceOfShape aFaces;
2094   int NumBase = 0, i = 0, NbQFs = 0;
2095   for (TopExp_Explorer exp(theShape, TopAbs_FACE); exp.More(); exp.Next()) {
2096     i++;
2097     aFaces.Append(exp.Current());
2098     SMESH_subMesh *aSubMesh = theMesh.GetSubMesh(exp.Current());
2099     meshFaces.push_back(aSubMesh);
2100     MapShapeNbElemsItr anIt = aResMap.find(meshFaces[i-1]);
2101     if( anIt==aResMap.end() )
2102       return toSM( error( "Submesh can not be evaluated"));
2103
2104     std::vector<int> aVec = (*anIt).second;
2105     int nbtri = Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
2106     int nbqua = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
2107     if( nbtri==0 && nbqua>0 ) {
2108       NbQFs++;
2109     }
2110     if( nbtri>0 ) {
2111       NumBase = i;
2112     }
2113   }
2114
2115   if(NbQFs<4) {
2116     std::vector<int> aResVec(SMDSEntity_Last);
2117     for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
2118     SMESH_subMesh * sm = theMesh.GetSubMesh(theShape);
2119     aResMap.insert(std::make_pair(sm,aResVec));
2120     return toSM( error( "Submesh can not be evaluated" ));
2121   }
2122
2123   if(NumBase==0) NumBase = 1; // only quads => set 1 faces as base
2124
2125   // find number of 1d elems for base face
2126   int nb1d = 0;
2127   TopTools_MapOfShape Edges1;
2128   for (TopExp_Explorer exp(aFaces.Value(NumBase), TopAbs_EDGE); exp.More(); exp.Next()) {
2129     Edges1.Add(exp.Current());
2130     SMESH_subMesh *sm = theMesh.GetSubMesh(exp.Current());
2131     if( sm ) {
2132       MapShapeNbElemsItr anIt = aResMap.find(sm);
2133       if( anIt == aResMap.end() ) continue;
2134       std::vector<int> aVec = (*anIt).second;
2135       nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
2136     }
2137   }
2138   // find face opposite to base face
2139   int OppNum = 0;
2140   for(i=1; i<=6; i++) {
2141     if(i==NumBase) continue;
2142     bool IsOpposite = true;
2143     for(TopExp_Explorer exp(aFaces.Value(i), TopAbs_EDGE); exp.More(); exp.Next()) {
2144       if( Edges1.Contains(exp.Current()) ) {
2145         IsOpposite = false;
2146         break;
2147       }
2148     }
2149     if(IsOpposite) {
2150       OppNum = i;
2151       break;
2152     }
2153   }
2154   // find number of 2d elems on side faces
2155   int nb2d = 0;
2156   for(i=1; i<=6; i++) {
2157     if( i==OppNum || i==NumBase ) continue;
2158     MapShapeNbElemsItr anIt = aResMap.find( meshFaces[i-1] );
2159     if( anIt == aResMap.end() ) continue;
2160     std::vector<int> aVec = (*anIt).second;
2161     nb2d += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
2162   }
2163
2164   MapShapeNbElemsItr anIt = aResMap.find( meshFaces[NumBase-1] );
2165   std::vector<int> aVec = (*anIt).second;
2166   bool IsQuadratic = (aVec[SMDSEntity_Quad_Triangle]>aVec[SMDSEntity_Triangle]) ||
2167                      (aVec[SMDSEntity_Quad_Quadrangle]>aVec[SMDSEntity_Quadrangle]);
2168   int nb2d_face0_3 = Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
2169   int nb2d_face0_4 = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
2170   int nb0d_face0 = aVec[SMDSEntity_Node];
2171   int nb1d_face0_int = ( nb2d_face0_3*3 + nb2d_face0_4*4 - nb1d ) / 2;
2172
2173   std::vector<int> aResVec(SMDSEntity_Last);
2174   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
2175   if(IsQuadratic) {
2176     aResVec[SMDSEntity_Quad_Penta] = nb2d_face0_3 * ( nb2d/nb1d );
2177     aResVec[SMDSEntity_Quad_Hexa] = nb2d_face0_4 * ( nb2d/nb1d );
2178     aResVec[SMDSEntity_Node] = nb0d_face0 * ( 2*nb2d/nb1d - 1 ) - nb1d_face0_int * nb2d/nb1d;
2179   }
2180   else {
2181     aResVec[SMDSEntity_Node] = nb0d_face0 * ( nb2d/nb1d - 1 );
2182     aResVec[SMDSEntity_Penta] = nb2d_face0_3 * ( nb2d/nb1d );
2183     aResVec[SMDSEntity_Hexa] = nb2d_face0_4 * ( nb2d/nb1d );
2184   }
2185   SMESH_subMesh * sm = theMesh.GetSubMesh(theShape);
2186   aResMap.insert(std::make_pair(sm,aResVec));
2187
2188   return true;
2189 }
2190
2191 //================================================================================
2192 /*!
2193  * \brief Create prisms
2194  * \param columns - columns of nodes generated from nodes of a mesh face
2195  * \param helper - helper initialized by mesh and shape to add prisms to
2196  */
2197 //================================================================================
2198
2199 bool StdMeshers_Prism_3D::AddPrisms( vector<const TNodeColumn*> & columns,
2200                                      SMESH_MesherHelper*          helper)
2201 {
2202   size_t nbNodes = columns.size();
2203   size_t nbZ     = columns[0]->size();
2204   if ( nbZ < 2 ) return false;
2205   for ( size_t i = 1; i < nbNodes; ++i )
2206     if ( columns[i]->size() != nbZ )
2207       return false;
2208
2209   // find out orientation
2210   bool isForward = true;
2211   SMDS_VolumeTool vTool;
2212   size_t z = 1;
2213   switch ( nbNodes ) {
2214   case 3: {
2215     SMDS_VolumeOfNodes tmpPenta ( (*columns[0])[z-1], // bottom
2216                                   (*columns[1])[z-1],
2217                                   (*columns[2])[z-1],
2218                                   (*columns[0])[z],   // top
2219                                   (*columns[1])[z],
2220                                   (*columns[2])[z] );
2221     vTool.Set( &tmpPenta );
2222     isForward  = vTool.IsForward();
2223     break;
2224   }
2225   case 4: {
2226     SMDS_VolumeOfNodes tmpHex( (*columns[0])[z-1], (*columns[1])[z-1], // bottom
2227                                (*columns[2])[z-1], (*columns[3])[z-1],
2228                                (*columns[0])[z],   (*columns[1])[z],   // top
2229                                (*columns[2])[z],   (*columns[3])[z] );
2230     vTool.Set( &tmpHex );
2231     isForward  = vTool.IsForward();
2232     break;
2233   }
2234   default:
2235     const int di = (nbNodes+1) / 3;
2236     SMDS_VolumeOfNodes tmpVol ( (*columns[0]   )[z-1],
2237                                 (*columns[di]  )[z-1],
2238                                 (*columns[2*di])[z-1],
2239                                 (*columns[0]   )[z],
2240                                 (*columns[di]  )[z],
2241                                 (*columns[2*di])[z] );
2242     vTool.Set( &tmpVol );
2243     isForward  = vTool.IsForward();
2244   }
2245
2246   // vertical loop on columns
2247
2248   helper->SetElementsOnShape( true );
2249
2250   switch ( nbNodes ) {
2251
2252   case 3: { // ---------- pentahedra
2253     const int i1 = isForward ? 1 : 2;
2254     const int i2 = isForward ? 2 : 1;
2255     for ( z = 1; z < nbZ; ++z )
2256       helper->AddVolume( (*columns[0 ])[z-1], // bottom
2257                          (*columns[i1])[z-1],
2258                          (*columns[i2])[z-1],
2259                          (*columns[0 ])[z],   // top
2260                          (*columns[i1])[z],
2261                          (*columns[i2])[z] );
2262     break;
2263   }
2264   case 4: { // ---------- hexahedra
2265     const int i1 = isForward ? 1 : 3;
2266     const int i3 = isForward ? 3 : 1;
2267     for ( z = 1; z < nbZ; ++z )
2268       helper->AddVolume( (*columns[0])[z-1], (*columns[i1])[z-1], // bottom
2269                          (*columns[2])[z-1], (*columns[i3])[z-1],
2270                          (*columns[0])[z],   (*columns[i1])[z],     // top
2271                          (*columns[2])[z],   (*columns[i3])[z] );
2272     break;
2273   }
2274   case 6: { // ---------- octahedra
2275     const int iBase1 = isForward ? -1 : 0;
2276     const int iBase2 = isForward ?  0 :-1;
2277     for ( z = 1; z < nbZ; ++z )
2278       helper->AddVolume( (*columns[0])[z+iBase1], (*columns[1])[z+iBase1], // bottom or top
2279                          (*columns[2])[z+iBase1], (*columns[3])[z+iBase1],
2280                          (*columns[4])[z+iBase1], (*columns[5])[z+iBase1],
2281                          (*columns[0])[z+iBase2], (*columns[1])[z+iBase2], // top or bottom
2282                          (*columns[2])[z+iBase2], (*columns[3])[z+iBase2],
2283                          (*columns[4])[z+iBase2], (*columns[5])[z+iBase2] );
2284     break;
2285   }
2286   default: // ---------- polyhedra
2287     vector<int> quantities( 2 + nbNodes, 4 );
2288     quantities[0] = quantities[1] = nbNodes;
2289     columns.resize( nbNodes + 1 );
2290     columns[ nbNodes ] = columns[ 0 ];
2291     const int i1 = isForward ? 1 : 3;
2292     const int i3 = isForward ? 3 : 1;
2293     const int iBase1 = isForward ? -1 : 0;
2294     const int iBase2 = isForward ?  0 :-1;
2295     vector<const SMDS_MeshNode*> nodes( 2*nbNodes + 4*nbNodes);
2296     for ( z = 1; z < nbZ; ++z )
2297     {
2298       for ( size_t i = 0; i < nbNodes; ++i ) {
2299         nodes[ i             ] = (*columns[ i ])[z+iBase1]; // bottom or top
2300         nodes[ 2*nbNodes-i-1 ] = (*columns[ i ])[z+iBase2]; // top or bottom
2301         // side
2302         int di = 2*nbNodes + 4*i;
2303         nodes[ di+0 ] = (*columns[i  ])[z  ];
2304         nodes[ di+i1] = (*columns[i+1])[z  ];
2305         nodes[ di+2 ] = (*columns[i+1])[z-1];
2306         nodes[ di+i3] = (*columns[i  ])[z-1];
2307       }
2308       helper->AddPolyhedralVolume( nodes, quantities );
2309     }
2310
2311   } // switch ( nbNodes )
2312
2313   return true;
2314 }
2315
2316 //================================================================================
2317 /*!
2318  * \brief Find correspondence between bottom and top nodes
2319  *  If elements on the bottom and top faces are topologically different,
2320  *  and projection is possible and allowed, perform the projection
2321  *  \retval bool - is a success or not
2322  */
2323 //================================================================================
2324
2325 bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf,
2326                                                  const Prism_3D::TPrismTopo& thePrism)
2327 {
2328   SMESH_subMesh * botSM = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom );
2329   SMESH_subMesh * topSM = myHelper->GetMesh()->GetSubMesh( thePrism.myTop    );
2330
2331   SMESHDS_SubMesh * botSMDS = botSM->GetSubMeshDS();
2332   SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS();
2333
2334   if ( !botSMDS || botSMDS->NbElements() == 0 )
2335   {
2336     _gen->Compute( *myHelper->GetMesh(), botSM->GetSubShape(), /*aShapeOnly=*/true );
2337     botSMDS = botSM->GetSubMeshDS();
2338     if ( !botSMDS || botSMDS->NbElements() == 0 )
2339       return toSM( error(TCom("No elements on face #") << botSM->GetId() ));
2340   }
2341
2342   bool needProject = !topSM->IsMeshComputed();
2343   if ( !needProject &&
2344        (botSMDS->NbElements() != topSMDS->NbElements() ||
2345         botSMDS->NbNodes()    != topSMDS->NbNodes()))
2346   {
2347     MESSAGE("nb elem bot " << botSMDS->NbElements() <<
2348             " top " << ( topSMDS ? topSMDS->NbElements() : 0 ));
2349     MESSAGE("nb node bot " << botSMDS->NbNodes() <<
2350             " top " << ( topSMDS ? topSMDS->NbNodes() : 0 ));
2351     return toSM( error(TCom("Mesh on faces #") << botSM->GetId()
2352                        <<" and #"<< topSM->GetId() << " seems different" ));
2353   }
2354
2355   if ( 0/*needProject && !myProjectTriangles*/ )
2356     return toSM( error(TCom("Mesh on faces #") << botSM->GetId()
2357                        <<" and #"<< topSM->GetId() << " seems different" ));
2358   ///RETURN_BAD_RESULT("Need to project but not allowed");
2359
2360   NSProjUtils::TNodeNodeMap n2nMap;
2361   const NSProjUtils::TNodeNodeMap* n2nMapPtr = & n2nMap;
2362   if ( needProject )
2363   {
2364     if ( !projectBottomToTop( bottomToTopTrsf, thePrism ))
2365       return false;
2366     n2nMapPtr = & TProjction2dAlgo::instance( this )->GetNodesMap();
2367   }
2368
2369   if ( !n2nMapPtr || (int) n2nMapPtr->size() < botSMDS->NbNodes() )
2370   {
2371     // associate top and bottom faces
2372     NSProjUtils::TShapeShapeMap shape2ShapeMap;
2373     const bool sameTopo =
2374       NSProjUtils::FindSubShapeAssociation( thePrism.myBottom, myHelper->GetMesh(),
2375                                             thePrism.myTop,    myHelper->GetMesh(),
2376                                             shape2ShapeMap);
2377     if ( !sameTopo )
2378       for ( size_t iQ = 0; iQ < thePrism.myWallQuads.size(); ++iQ )
2379       {
2380         const Prism_3D::TQuadList& quadList = thePrism.myWallQuads[iQ];
2381         StdMeshers_FaceSidePtr      botSide = quadList.front()->side[ QUAD_BOTTOM_SIDE ];
2382         StdMeshers_FaceSidePtr      topSide = quadList.back ()->side[ QUAD_TOP_SIDE ];
2383         if ( botSide->NbEdges() == topSide->NbEdges() )
2384         {
2385           for ( int iE = 0; iE < botSide->NbEdges(); ++iE )
2386           {
2387             NSProjUtils::InsertAssociation( botSide->Edge( iE ),
2388                                             topSide->Edge( iE ), shape2ShapeMap );
2389             NSProjUtils::InsertAssociation( myHelper->IthVertex( 0, botSide->Edge( iE )),
2390                                             myHelper->IthVertex( 0, topSide->Edge( iE )),
2391                                             shape2ShapeMap );
2392           }
2393         }
2394         else
2395         {
2396           TopoDS_Vertex vb, vt;
2397           StdMeshers_FaceSidePtr sideB, sideT;
2398           vb = myHelper->IthVertex( 0, botSide->Edge( 0 ));
2399           vt = myHelper->IthVertex( 0, topSide->Edge( 0 ));
2400           sideB = quadList.front()->side[ QUAD_LEFT_SIDE ];
2401           sideT = quadList.back ()->side[ QUAD_LEFT_SIDE ];
2402           if ( vb.IsSame( sideB->FirstVertex() ) &&
2403                vt.IsSame( sideT->LastVertex() ))
2404           {
2405             NSProjUtils::InsertAssociation( botSide->Edge( 0 ),
2406                                             topSide->Edge( 0 ), shape2ShapeMap );
2407             NSProjUtils::InsertAssociation( vb, vt, shape2ShapeMap );
2408           }
2409           vb = myHelper->IthVertex( 1, botSide->Edge( botSide->NbEdges()-1 ));
2410           vt = myHelper->IthVertex( 1, topSide->Edge( topSide->NbEdges()-1 ));
2411           sideB = quadList.front()->side[ QUAD_RIGHT_SIDE ];
2412           sideT = quadList.back ()->side[ QUAD_RIGHT_SIDE ];
2413           if ( vb.IsSame( sideB->FirstVertex() ) &&
2414                vt.IsSame( sideT->LastVertex() ))
2415           {
2416             NSProjUtils::InsertAssociation( botSide->Edge( botSide->NbEdges()-1 ),
2417                                             topSide->Edge( topSide->NbEdges()-1 ),
2418                                             shape2ShapeMap );
2419             NSProjUtils::InsertAssociation( vb, vt, shape2ShapeMap );
2420           }
2421         }
2422       }
2423
2424     // Find matching nodes of top and bottom faces
2425     n2nMapPtr = & n2nMap;
2426     if ( ! NSProjUtils::FindMatchingNodesOnFaces( thePrism.myBottom, myHelper->GetMesh(),
2427                                                   thePrism.myTop,    myHelper->GetMesh(),
2428                                                   shape2ShapeMap, n2nMap ))
2429     {
2430       if ( sameTopo )
2431         return toSM( error(TCom("Mesh on faces #") << botSM->GetId()
2432                            <<" and #"<< topSM->GetId() << " seems different" ));
2433       else
2434         return toSM( error(TCom("Topology of faces #") << botSM->GetId()
2435                            <<" and #"<< topSM->GetId() << " seems different" ));
2436     }
2437   }
2438
2439   // Fill myBotToColumnMap
2440
2441   int zSize = myBlock.VerticalSize();
2442   TNodeNodeMap::const_iterator bN_tN = n2nMapPtr->begin();
2443   for ( ; bN_tN != n2nMapPtr->end(); ++bN_tN )
2444   {
2445     const SMDS_MeshNode* botNode = bN_tN->first;
2446     const SMDS_MeshNode* topNode = bN_tN->second;
2447     if ( botNode->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE &&
2448          myBlock.HasNodeColumn( botNode ))
2449       continue; // wall columns are contained in myBlock
2450     // create node column
2451     Prism_3D::TNode bN( botNode );
2452     TNode2ColumnMap::iterator bN_col =
2453       myBotToColumnMap.insert( make_pair ( bN, TNodeColumn() )).first;
2454     TNodeColumn & column = bN_col->second;
2455     column.resize( zSize, 0 );
2456     column.front() = botNode;
2457     column.back()  = topNode;
2458   }
2459   return true;
2460 }
2461
2462 //================================================================================
2463 /*!
2464  * \brief Remove faces from the top face and re-create them by projection from the bottom
2465  * \retval bool - a success or not
2466  */
2467 //================================================================================
2468
2469 bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf &             bottomToTopTrsf,
2470                                               const Prism_3D::TPrismTopo& thePrism )
2471 {
2472   if ( project2dMesh( thePrism.myBottom, thePrism.myTop ))
2473   {
2474     return true;
2475   }
2476   NSProjUtils::TNodeNodeMap& n2nMap =
2477     (NSProjUtils::TNodeNodeMap&) TProjction2dAlgo::instance( this )->GetNodesMap();
2478   n2nMap.clear();
2479
2480   myUseBlock = true;
2481
2482   SMESHDS_Mesh*  meshDS = myHelper->GetMeshDS();
2483   SMESH_subMesh * botSM = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom );
2484   SMESH_subMesh * topSM = myHelper->GetMesh()->GetSubMesh( thePrism.myTop );
2485
2486   SMESHDS_SubMesh * botSMDS = botSM->GetSubMeshDS();
2487   SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS();
2488
2489   if ( topSMDS && topSMDS->NbElements() > 0 )
2490   {
2491     //topSM->ComputeStateEngine( SMESH_subMesh::CLEAN ); -- avoid propagation of events
2492     for ( SMDS_ElemIteratorPtr eIt = topSMDS->GetElements(); eIt->more(); )
2493       meshDS->RemoveFreeElement( eIt->next(), topSMDS, /*fromGroups=*/false );
2494     for ( SMDS_NodeIteratorPtr nIt = topSMDS->GetNodes(); nIt->more(); )
2495       meshDS->RemoveFreeNode( nIt->next(), topSMDS, /*fromGroups=*/false );
2496   }
2497
2498   const TopoDS_Face& botFace = thePrism.myBottom; // oriented within
2499   const TopoDS_Face& topFace = thePrism.myTop;    //    the 3D SHAPE
2500   int topFaceID = meshDS->ShapeToIndex( thePrism.myTop );
2501
2502   SMESH_MesherHelper botHelper( *myHelper->GetMesh() );
2503   botHelper.SetSubShape( botFace );
2504   botHelper.ToFixNodeParameters( true );
2505   bool checkUV;
2506   SMESH_MesherHelper topHelper( *myHelper->GetMesh() );
2507   topHelper.SetSubShape( topFace );
2508   topHelper.ToFixNodeParameters( true );
2509   double distXYZ[4], fixTol = 10 * topHelper.MaxTolerance( topFace );
2510
2511   // Fill myBotToColumnMap
2512
2513   int zSize = myBlock.VerticalSize();
2514   Prism_3D::TNode prevTNode;
2515   SMDS_NodeIteratorPtr nIt = botSMDS->GetNodes();
2516   while ( nIt->more() )
2517   {
2518     const SMDS_MeshNode* botNode = nIt->next();
2519     const SMDS_MeshNode* topNode = 0;
2520     if ( botNode->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE )
2521       continue; // strange
2522
2523     Prism_3D::TNode bN( botNode );
2524     if ( bottomToTopTrsf.Form() == gp_Identity )
2525     {
2526       // compute bottom node params
2527       gp_XYZ paramHint(-1,-1,-1);
2528       if ( prevTNode.IsNeighbor( bN ))
2529       {
2530         paramHint = prevTNode.GetParams();
2531         // double tol = 1e-2 * ( prevTNode.GetCoords() - bN.GetCoords() ).Modulus();
2532         // myBlock.SetTolerance( Min( myBlock.GetTolerance(), tol ));
2533       }
2534       if ( !myBlock.ComputeParameters( bN.GetCoords(), bN.ChangeParams(),
2535                                        ID_BOT_FACE, paramHint ))
2536         return toSM( error(TCom("Can't compute normalized parameters for node ")
2537                            << botNode->GetID() << " on the face #"<< botSM->GetId() ));
2538       prevTNode = bN;
2539       // compute top node coords
2540       gp_XYZ topXYZ; gp_XY topUV;
2541       if ( !myBlock.FacePoint( ID_TOP_FACE, bN.GetParams(), topXYZ ) ||
2542            !myBlock.FaceUV   ( ID_TOP_FACE, bN.GetParams(), topUV ))
2543         return toSM( error(TCom("Can't compute coordinates "
2544                                 "by normalized parameters on the face #")<< topSM->GetId() ));
2545       topNode = meshDS->AddNode( topXYZ.X(),topXYZ.Y(),topXYZ.Z() );
2546       meshDS->SetNodeOnFace( topNode, topFaceID, topUV.X(), topUV.Y() );
2547     }
2548     else // use bottomToTopTrsf
2549     {
2550       gp_XYZ coords = bN.GetCoords();
2551       bottomToTopTrsf.Transforms( coords );
2552       topNode = meshDS->AddNode( coords.X(), coords.Y(), coords.Z() );
2553       gp_XY topUV = botHelper.GetNodeUV( botFace, botNode, 0, &checkUV );
2554       meshDS->SetNodeOnFace( topNode, topFaceID, topUV.X(), topUV.Y() );
2555       distXYZ[0] = -1;
2556       if ( topHelper.CheckNodeUV( topFace, topNode, topUV, fixTol, /*force=*/false, distXYZ ) &&
2557            distXYZ[0] > fixTol && distXYZ[0] < fixTol * 1e+3 )
2558         meshDS->MoveNode( topNode, distXYZ[1], distXYZ[2], distXYZ[3] ); // transform can be inaccurate
2559     }
2560     // create node column
2561     TNode2ColumnMap::iterator bN_col =
2562       myBotToColumnMap.insert( make_pair ( bN, TNodeColumn() )).first;
2563     TNodeColumn & column = bN_col->second;
2564     column.resize( zSize );
2565     column.front() = botNode;
2566     column.back()  = topNode;
2567
2568     n2nMap.insert( n2nMap.end(), make_pair( botNode, topNode ));
2569
2570     if ( _computeCanceled )
2571       return toSM( error( SMESH_ComputeError::New(COMPERR_CANCELED)));
2572   }
2573
2574   // Create top faces
2575
2576   const bool oldSetElemsOnShape = myHelper->SetElementsOnShape( false );
2577
2578   // care of orientation;
2579   // if the bottom faces is orienetd OK then top faces must be reversed
2580   bool reverseTop = true;
2581   if ( myHelper->NbAncestors( botFace, *myBlock.Mesh(), TopAbs_SOLID ) > 1 )
2582     reverseTop = ! myHelper->IsReversedSubMesh( botFace );
2583   int iFrw, iRev, *iPtr = &( reverseTop ? iRev : iFrw );
2584
2585   // loop on bottom mesh faces
2586   SMDS_ElemIteratorPtr faceIt = botSMDS->GetElements();
2587   vector< const SMDS_MeshNode* > nodes;
2588   while ( faceIt->more() )
2589   {
2590     const SMDS_MeshElement* face = faceIt->next();
2591     if ( !face || face->GetType() != SMDSAbs_Face )
2592       continue;
2593
2594     // find top node in columns for each bottom node
2595     int nbNodes = face->NbCornerNodes();
2596     nodes.resize( nbNodes );
2597     for ( iFrw = 0, iRev = nbNodes-1; iFrw < nbNodes; ++iFrw, --iRev )
2598     {
2599       const SMDS_MeshNode* n = face->GetNode( *iPtr );
2600       if ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) {
2601         TNode2ColumnMap::iterator bot_column = myBotToColumnMap.find( n );
2602         if ( bot_column == myBotToColumnMap.end() )
2603           return toSM( error(TCom("No nodes found above node ") << n->GetID() ));
2604         nodes[ iFrw ] = bot_column->second.back();
2605       }
2606       else {
2607         const TNodeColumn* column = myBlock.GetNodeColumn( n );
2608         if ( !column )
2609           return toSM( error(TCom("No side nodes found above node ") << n->GetID() ));
2610         nodes[ iFrw ] = column->back();
2611       }
2612     }
2613     SMDS_MeshElement* newFace = 0;
2614     switch ( nbNodes ) {
2615
2616     case 3: {
2617       newFace = myHelper->AddFace(nodes[0], nodes[1], nodes[2]);
2618       break;
2619     }
2620     case 4: {
2621       newFace = myHelper->AddFace( nodes[0], nodes[1], nodes[2], nodes[3] );
2622       break;
2623     }
2624     default:
2625       newFace = meshDS->AddPolygonalFace( nodes );
2626     }
2627     if ( newFace )
2628       meshDS->SetMeshElementOnShape( newFace, topFaceID );
2629   }
2630
2631   myHelper->SetElementsOnShape( oldSetElemsOnShape );
2632
2633   // Check the projected mesh
2634
2635   if ( thePrism.NbWires() > 1 && // there are holes
2636        topHelper.IsDistorted2D( topSM, /*checkUV=*/false ))
2637   {
2638     SMESH_MeshEditor editor( topHelper.GetMesh() );
2639
2640     // smooth in 2D or 3D?
2641     TopLoc_Location loc;
2642     Handle(Geom_Surface) surface = BRep_Tool::Surface( topFace, loc );
2643     bool isPlanar = GeomLib_IsPlanarSurface( surface ).IsPlanar();
2644
2645     set<const SMDS_MeshNode*> fixedNodes;
2646     TIDSortedElemSet faces;
2647     for ( faceIt = topSMDS->GetElements(); faceIt->more(); )
2648       faces.insert( faces.end(), faceIt->next() );
2649
2650     bool isOk = false;
2651     for ( int isCentroidal = 0; isCentroidal < 2; ++isCentroidal )
2652     {
2653       SMESH_MeshEditor::SmoothMethod algo =
2654         isCentroidal ? SMESH_MeshEditor::CENTROIDAL : SMESH_MeshEditor::LAPLACIAN;
2655
2656       int nbAttempts = isCentroidal ? 1 : 10;
2657       for ( int iAttemp = 0; iAttemp < nbAttempts; ++iAttemp )
2658       {
2659         TIDSortedElemSet workFaces = faces;
2660
2661         // smoothing
2662         editor.Smooth( workFaces, fixedNodes, algo, /*nbIterations=*/ 10,
2663                        /*theTgtAspectRatio=*/1.0, /*the2D=*/!isPlanar);
2664
2665         if (( isOk = !topHelper.IsDistorted2D( topSM, /*checkUV=*/true )) &&
2666             ( !isCentroidal ))
2667           break;
2668       }
2669     }
2670     if ( !isOk )
2671       return toSM( error( TCom("Projection from face #") << botSM->GetId()
2672                           << " to face #" << topSM->GetId()
2673                           << " failed: inverted elements created"));
2674   }
2675
2676   TProjction2dAlgo::instance( this )->SetEventListener( topSM );
2677
2678   return true;
2679 }
2680
2681 //=======================================================================
2682 //function : getSweepTolerance
2683 //purpose  : Compute tolerance to pass to StdMeshers_Sweeper
2684 //=======================================================================
2685
2686 double StdMeshers_Prism_3D::getSweepTolerance( const Prism_3D::TPrismTopo& thePrism )
2687 {
2688   SMESHDS_Mesh*    meshDS = myHelper->GetMeshDS();
2689   SMESHDS_SubMesh * sm[2] = { meshDS->MeshElements( thePrism.myBottom ),
2690                               meshDS->MeshElements( thePrism.myTop )    };
2691   double minDist = 1e100;
2692
2693   vector< SMESH_TNodeXYZ > nodes;
2694   for ( int iSM = 0; iSM < 2; ++iSM )
2695   {
2696     if ( !sm[ iSM ]) continue;
2697
2698     SMDS_ElemIteratorPtr fIt = sm[ iSM ]->GetElements();
2699     while ( fIt->more() )
2700     {
2701       const SMDS_MeshElement* face = fIt->next();
2702       const int            nbNodes = face->NbCornerNodes();
2703       SMDS_ElemIteratorPtr     nIt = face->nodesIterator();
2704
2705       nodes.resize( nbNodes + 1 );
2706       for ( int iN = 0; iN < nbNodes; ++iN )
2707         nodes[ iN ] = nIt->next();
2708       nodes.back() = nodes[0];
2709
2710       // loop on links
2711       double dist2;
2712       for ( int iN = 0; iN < nbNodes; ++iN )
2713       {
2714         if ( nodes[ iN   ]._node->GetPosition()->GetDim() < 2 &&
2715              nodes[ iN+1 ]._node->GetPosition()->GetDim() < 2 )
2716         {
2717           // it's a boundary link; measure distance of other
2718           // nodes to this link
2719           gp_XYZ linkDir = nodes[ iN ] - nodes[ iN+1 ];
2720           double linkLen = linkDir.Modulus();
2721           bool   isDegen = ( linkLen < numeric_limits<double>::min() );
2722           if ( !isDegen ) linkDir /= linkLen;
2723           for ( int iN2 = 0; iN2 < nbNodes; ++iN2 ) // loop on other nodes
2724           {
2725             if ( nodes[ iN2 ] == nodes[ iN ] ||
2726                  nodes[ iN2 ] == nodes[ iN+1 ]) continue;
2727             if ( isDegen )
2728             {
2729               dist2 = ( nodes[ iN ] - nodes[ iN2 ]).SquareModulus();
2730             }
2731             else
2732             {
2733               dist2 = linkDir.CrossSquareMagnitude( nodes[ iN ] - nodes[ iN2 ]);
2734             }
2735             if ( dist2 > numeric_limits<double>::min() )
2736               minDist = Min ( minDist, dist2 );
2737           }
2738         }
2739         // measure length link
2740         else if ( nodes[ iN ]._node < nodes[ iN+1 ]._node ) // not to measure same link twice
2741         {
2742           dist2 = ( nodes[ iN ] - nodes[ iN+1 ]).SquareModulus();
2743           if ( dist2 > numeric_limits<double>::min() )
2744             minDist = Min ( minDist, dist2 );
2745         }
2746       }
2747     }
2748   }
2749   return 0.1 * Sqrt ( minDist );
2750 }
2751
2752 //=======================================================================
2753 //function : isSimpleQuad
2754 //purpose  : check if the bottom FACE is meshable with nice quadrangles,
2755 //           if so the block approach can work rather fast.
2756 //           This is a temporary mean caused by problems in StdMeshers_Sweeper
2757 //=======================================================================
2758
2759 bool StdMeshers_Prism_3D::isSimpleBottom( const Prism_3D::TPrismTopo& thePrism )
2760 {
2761   if ( thePrism.myNbEdgesInWires.front() != 4 )
2762     return false;
2763
2764   // analyse angles between edges
2765   double nbConcaveAng = 0, nbConvexAng = 0;
2766   TopoDS_Face reverseBottom = TopoDS::Face( thePrism.myBottom.Reversed() ); // see initPrism()
2767   TopoDS_Vertex commonV;
2768   const list< TopoDS_Edge >& botEdges = thePrism.myBottomEdges;
2769   list< TopoDS_Edge >::const_iterator edge = botEdges.begin();
2770   while ( edge != botEdges.end() )
2771   {
2772     if ( SMESH_Algo::isDegenerated( *edge ))
2773       return false;
2774     TopoDS_Edge e1 = *edge++;
2775     TopoDS_Edge e2 = ( edge == botEdges.end() ? botEdges.front() : *edge );
2776     if ( ! TopExp::CommonVertex( e1, e2,  commonV ))
2777     {
2778       e2 = botEdges.front();
2779       if ( ! TopExp::CommonVertex( e1, e2,  commonV ))
2780         break;
2781     }
2782     double angle = myHelper->GetAngle( e1, e2, reverseBottom, commonV );
2783     if ( angle < -5 * M_PI/180 )
2784       if ( ++nbConcaveAng > 1 )
2785         return false;
2786     if ( angle > 85 * M_PI/180 )
2787       if ( ++nbConvexAng > 4 )
2788         return false;
2789   }
2790   return true;
2791 }
2792
2793 //=======================================================================
2794 //function : allVerticalEdgesStraight
2795 //purpose  : Defines if all "vertical" EDGEs are straight
2796 //=======================================================================
2797
2798 bool StdMeshers_Prism_3D::allVerticalEdgesStraight( const Prism_3D::TPrismTopo& thePrism )
2799 {
2800   for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
2801   {
2802     const Prism_3D::TQuadList& quads = thePrism.myWallQuads[i];
2803     Prism_3D::TQuadList::const_iterator quadIt = quads.begin();
2804     TopoDS_Edge prevQuadEdge;
2805     for ( ; quadIt != quads.end(); ++quadIt )
2806     {
2807       StdMeshers_FaceSidePtr rightSide = (*quadIt)->side[ QUAD_RIGHT_SIDE ];
2808
2809       if ( !prevQuadEdge.IsNull() &&
2810            !SMESH_Algo::IsContinuous( rightSide->Edge( 0 ), prevQuadEdge ))
2811         return false;
2812
2813       for ( int iE = 0; iE < rightSide->NbEdges(); ++iE )
2814       {
2815         const TopoDS_Edge & rightE = rightSide->Edge( iE );
2816         if ( !SMESH_Algo::IsStraight( rightE, /*degenResult=*/true ))
2817           return false;
2818
2819         if ( iE > 0 &&
2820              !SMESH_Algo::IsContinuous( rightSide->Edge( iE-1 ), rightE ))
2821           return false;
2822
2823         prevQuadEdge = rightE;
2824       }
2825     }
2826   }
2827   return true;
2828 }
2829
2830 //=======================================================================
2831 //function : project2dMesh
2832 //purpose  : Project mesh faces from a source FACE of one prism (theSrcFace)
2833 //           to a source FACE of another prism (theTgtFace)
2834 //=======================================================================
2835
2836 bool StdMeshers_Prism_3D::project2dMesh(const TopoDS_Face& theSrcFace,
2837                                         const TopoDS_Face& theTgtFace)
2838 {
2839   if ( CountEdges( theSrcFace ) != CountEdges( theTgtFace ))
2840     return false;
2841
2842   TProjction2dAlgo* projector2D = TProjction2dAlgo::instance( this );
2843   projector2D->myHyp.SetSourceFace( theSrcFace );
2844   bool ok = projector2D->Compute( *myHelper->GetMesh(), theTgtFace );
2845
2846   SMESH_subMesh* tgtSM = myHelper->GetMesh()->GetSubMesh( theTgtFace );
2847   if ( !ok && tgtSM->GetSubMeshDS() ) {
2848     //tgtSM->ComputeStateEngine( SMESH_subMesh::CLEAN ); -- avoid propagation of events
2849     SMESHDS_Mesh*     meshDS = myHelper->GetMeshDS();
2850     SMESHDS_SubMesh* tgtSMDS = tgtSM->GetSubMeshDS();
2851     for ( SMDS_ElemIteratorPtr eIt = tgtSMDS->GetElements(); eIt->more(); )
2852       meshDS->RemoveFreeElement( eIt->next(), tgtSMDS, /*fromGroups=*/false );
2853     for ( SMDS_NodeIteratorPtr nIt = tgtSMDS->GetNodes(); nIt->more(); )
2854       meshDS->RemoveFreeNode( nIt->next(), tgtSMDS, /*fromGroups=*/false );
2855   }
2856   tgtSM->ComputeStateEngine       ( SMESH_subMesh::CHECK_COMPUTE_STATE );
2857   tgtSM->ComputeSubMeshStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
2858
2859   projector2D->SetEventListener( tgtSM );
2860
2861   return ok;
2862 }
2863
2864 //================================================================================
2865 /*!
2866  * \brief Set projection coordinates of a node to a face and it's sub-shapes
2867  * \param faceID - the face given by in-block ID
2868  * \param params - node normalized parameters
2869  * \retval bool - is a success
2870  */
2871 //================================================================================
2872
2873 bool StdMeshers_Prism_3D::setFaceAndEdgesXYZ( const int faceID, const gp_XYZ& params, int z )
2874 {
2875   // find base and top edges of the face
2876   enum { BASE = 0, TOP, LEFT, RIGHT };
2877   vector< int > edgeVec; // 0-base, 1-top
2878   SMESH_Block::GetFaceEdgesIDs( faceID, edgeVec );
2879
2880   myBlock.EdgePoint( edgeVec[ BASE ], params, myShapeXYZ[ edgeVec[ BASE ]]);
2881   myBlock.EdgePoint( edgeVec[ TOP  ], params, myShapeXYZ[ edgeVec[ TOP ]]);
2882
2883   SHOWYXZ("\nparams ", params);
2884   SHOWYXZ("TOP is " <<edgeVec[ TOP ], myShapeXYZ[ edgeVec[ TOP]]);
2885   SHOWYXZ("BASE is "<<edgeVec[ BASE], myShapeXYZ[ edgeVec[ BASE]]);
2886
2887   if ( faceID == SMESH_Block::ID_Fx0z || faceID == SMESH_Block::ID_Fx1z )
2888   {
2889     myBlock.EdgePoint( edgeVec[ LEFT ], params, myShapeXYZ[ edgeVec[ LEFT ]]);
2890     myBlock.EdgePoint( edgeVec[ RIGHT ], params, myShapeXYZ[ edgeVec[ RIGHT ]]);
2891
2892     SHOWYXZ("VER "<<edgeVec[ LEFT], myShapeXYZ[ edgeVec[ LEFT]]);
2893     SHOWYXZ("VER "<<edgeVec[ RIGHT], myShapeXYZ[ edgeVec[ RIGHT]]);
2894   }
2895   myBlock.FacePoint( faceID, params, myShapeXYZ[ faceID ]);
2896   SHOWYXZ("FacePoint "<<faceID, myShapeXYZ[ faceID]);
2897
2898   return true;
2899 }
2900
2901 //=======================================================================
2902 //function : toSM
2903 //purpose  : If (!isOK), sets the error to a sub-mesh of a current SOLID
2904 //=======================================================================
2905
2906 bool StdMeshers_Prism_3D::toSM( bool isOK )
2907 {
2908   if ( mySetErrorToSM &&
2909        !isOK &&
2910        myHelper &&
2911        !myHelper->GetSubShape().IsNull() &&
2912        myHelper->GetSubShape().ShapeType() == TopAbs_SOLID)
2913   {
2914     SMESH_subMesh* sm = myHelper->GetMesh()->GetSubMesh( myHelper->GetSubShape() );
2915     sm->GetComputeError() = this->GetComputeError();
2916     // clear error in order not to return it twice
2917     _error = COMPERR_OK;
2918     _comment.clear();
2919   }
2920   return isOK;
2921 }
2922
2923 //=======================================================================
2924 //function : shapeID
2925 //purpose  : Return index of a shape
2926 //=======================================================================
2927
2928 int StdMeshers_Prism_3D::shapeID( const TopoDS_Shape& S )
2929 {
2930   if ( S.IsNull() ) return 0;
2931   if ( !myHelper  ) return -3;
2932   return myHelper->GetMeshDS()->ShapeToIndex( S );
2933 }
2934
2935 namespace // utils used by StdMeshers_Prism_3D::IsApplicable()
2936 {
2937   struct EdgeWithNeighbors
2938   {
2939     TopoDS_Edge   _edge;
2940     int           _iBase;     // index in a WIRE with non-base EDGEs excluded
2941     int           _iL, _iR;   // used to connect PrismSide's
2942     int           _iE;        // index in a WIRE
2943     int           _iLE, _iRE; // used to connect EdgeWithNeighbors's
2944     bool          _isBase;    // is used in a base FACE
2945     TopoDS_Vertex _vv[2];     // end VERTEXes
2946     EdgeWithNeighbors(const TopoDS_Edge& E,
2947                       int iE,  int nbE,  int shift,
2948                       int iEE, int nbEE, int shiftE,
2949                       bool isBase, bool setVV ):
2950       _edge( E ),
2951       _iBase( iE + shift ),
2952       _iL ( SMESH_MesherHelper::WrapIndex( iE-1,  Max( 1, nbE  )) + shift ),
2953       _iR ( SMESH_MesherHelper::WrapIndex( iE+1,  Max( 1, nbE  )) + shift ),
2954       _iE ( iEE + shiftE ),
2955       _iLE( SMESH_MesherHelper::WrapIndex( iEE-1, Max( 1, nbEE )) + shiftE ),
2956       _iRE( SMESH_MesherHelper::WrapIndex( iEE+1, Max( 1, nbEE )) + shiftE ),
2957       _isBase( isBase )
2958     {
2959       if ( setVV )
2960       {
2961         Vertex( 0 );
2962         Vertex( 1 );
2963       }
2964     }
2965     EdgeWithNeighbors() {}
2966     bool IsInternal() const { return !_edge.IsNull() && _edge.Orientation() == TopAbs_INTERNAL; }
2967     bool IsConnected( const EdgeWithNeighbors& edge, int iEnd ) const
2968     {
2969       return (( _vv[ iEnd ].IsSame( edge._vv[ 1 - iEnd ])) ||
2970               ( IsInternal() && _vv[ iEnd ].IsSame( edge._vv[ iEnd ])));
2971     }
2972     bool IsConnected( const std::vector< EdgeWithNeighbors > & edges, int iEnd ) const
2973     {
2974       int iEdge = iEnd ? _iRE : _iLE;
2975       return iEdge == _iE ? false : IsConnected( edges[ iEdge ], iEnd );
2976     }
2977     const TopoDS_Vertex& Vertex( int iEnd )
2978     {
2979       if ( _vv[ iEnd ].IsNull() )
2980         _vv[ iEnd ] = SMESH_MesherHelper::IthVertex( iEnd, _edge );
2981       return _vv[ iEnd ];
2982     }
2983   };
2984   // PrismSide contains all FACEs linking a bottom EDGE with a top one.
2985   struct PrismSide
2986   {
2987     TopoDS_Face                 _face;    // a currently treated upper FACE
2988     TopTools_IndexedMapOfShape *_faces;   // all FACEs (pointer because of a private copy constructor)
2989     TopoDS_Edge                 _topEdge; // a current top EDGE
2990     vector< EdgeWithNeighbors >*_edges;   // all EDGEs of _face
2991     int                         _iBotEdge;       // index of _topEdge within _edges
2992     vector< bool >              _isCheckedEdge;  // mark EDGEs whose two owner FACEs found
2993     int                         _nbCheckedEdges; // nb of EDGEs whose location is defined
2994     PrismSide                  *_leftSide;       // neighbor sides
2995     PrismSide                  *_rightSide;
2996     bool                        _isInternal; // whether this side raises from an INTERNAL EDGE
2997     //void SetExcluded() { _leftSide = _rightSide = NULL; }
2998     //bool IsExcluded() const { return !_leftSide; }
2999     const TopoDS_Edge& Edge( int i ) const
3000     {
3001       return (*_edges)[ i ]._edge;
3002     }
3003     int FindEdge( const TopoDS_Edge& E ) const
3004     {
3005       for ( size_t i = 0; i < _edges->size(); ++i )
3006         if ( E.IsSame( Edge( i ))) return i;
3007       return -1;
3008     }
3009     const TopoDS_Vertex& Vertex( int iE, int iEnd ) const
3010     {
3011       return (*_edges)[ iE ].Vertex( iEnd );
3012     }
3013     bool HasVertex( const TopoDS_Vertex& V ) const
3014     {
3015       for ( size_t i = 0; i < _edges->size(); ++i )
3016         if ( V.IsSame( Vertex( i, 0 ))) return true;
3017       return false;
3018     }
3019     bool IsSideFace( const TopTools_ListOfShape& faces,
3020                      const TopoDS_Face&          avoidFace,
3021                      const bool                  checkNeighbors ) const
3022     {
3023       TopTools_ListIteratorOfListOfShape faceIt( faces );
3024       for ( ; faceIt.More(); faceIt.Next() )
3025       {
3026         const TopoDS_Shape& face = faceIt.Value();
3027         if ( !face.IsSame( avoidFace ))
3028         {
3029           if ( _faces->Contains( face )) // avoid returning true for a prism top FACE
3030             return ( !_face.IsNull() || !( face.IsSame( _faces->FindKey( _faces->Extent() ))));
3031         }
3032       }
3033       if ( checkNeighbors )
3034         return (( _leftSide  && _leftSide->IsSideFace ( faces, avoidFace, false )) ||
3035                 ( _rightSide && _rightSide->IsSideFace( faces, avoidFace, false )));
3036
3037       return false;
3038     }
3039   };
3040   //--------------------------------------------------------------------------------
3041   /*!
3042    * \brief Return another faces sharing an edge
3043    */
3044   const TopoDS_Face & getAnotherFace( const TopoDS_Face& face,
3045                                       const TopTools_ListOfShape& faces)
3046   {
3047     TopTools_ListIteratorOfListOfShape faceIt( faces );
3048     for ( ; faceIt.More(); faceIt.Next() )
3049       if ( !face.IsSame( faceIt.Value() ))
3050         return TopoDS::Face( faceIt.Value() );
3051     return face;
3052   }
3053   //--------------------------------------------------------------------------------
3054   /*!
3055    * \brief Return another faces sharing an edge
3056    */
3057   const TopoDS_Face & getAnotherFace( const TopoDS_Face& face,
3058                                       const TopoDS_Edge& edge,
3059                                       TopTools_IndexedDataMapOfShapeListOfShape& facesOfEdge)
3060   {
3061     return getAnotherFace( face, facesOfEdge.FindFromKey( edge ));
3062   }
3063
3064   //--------------------------------------------------------------------------------
3065   /*!
3066    * \brief Return ordered edges of a face
3067    */
3068   //================================================================================
3069   /*!
3070    * \brief Return ordered edges of a face
3071    *  \param [in] face - the face
3072    *  \param [out] edges - return edge (edges from which no vertical faces raise excluded)
3073    *  \param [in] facesOfEdge - faces of each edge
3074    *  \param [in] noHolesAllowed - are multiple wires allowed
3075    */
3076   //================================================================================
3077
3078   bool getEdges( const TopoDS_Face&                         face,
3079                  vector< EdgeWithNeighbors > &              edges,
3080                  TopTools_IndexedDataMapOfShapeListOfShape& facesOfEdge,
3081                  const bool                                 noHolesAllowed)
3082   {
3083     TopoDS_Face f = face;
3084     if ( f.Orientation() != TopAbs_FORWARD &&
3085          f.Orientation() != TopAbs_REVERSED )
3086       f.Orientation( TopAbs_FORWARD );
3087     list< TopoDS_Edge > ee;
3088     list< int >         nbEdgesInWires;
3089     int nbW = SMESH_Block::GetOrderedEdges( f, ee, nbEdgesInWires );
3090     if ( nbW > 1 && noHolesAllowed )
3091       return false;
3092
3093     list< TopoDS_Edge >::iterator   e = ee.begin();
3094     list< int         >::iterator nbE = nbEdgesInWires.begin();
3095     for ( ; nbE != nbEdgesInWires.end(); ++nbE )
3096       for ( int iE = 0; iE < *nbE; ++e, ++iE )
3097         if ( SMESH_Algo::isDegenerated( *e )) // degenerated EDGE is never used
3098         {
3099           e = --ee.erase( e );
3100           --(*nbE);
3101           --iE;
3102         }
3103
3104     int iE, nbTot = 0, iBase, nbBase, nbTotBase = 0;
3105     vector<int> isBase;
3106     edges.clear();
3107     e = ee.begin();
3108     for ( nbE = nbEdgesInWires.begin(); nbE != nbEdgesInWires.end(); ++nbE )
3109     {
3110       nbBase = 0;
3111       isBase.resize( *nbE );
3112       list< TopoDS_Edge >::iterator eIt = e;
3113       for ( iE = 0; iE < *nbE; ++eIt, ++iE )
3114       {
3115         isBase[ iE ] = ( getAnotherFace( face, *eIt, facesOfEdge ) != face );
3116         nbBase += isBase[ iE ];
3117       }
3118       for ( iBase = 0, iE = 0; iE < *nbE; ++e, ++iE )
3119       {
3120         edges.push_back( EdgeWithNeighbors( *e,
3121                                             iBase, nbBase, nbTotBase,
3122                                             iE,    *nbE,   nbTot,
3123                                             isBase[ iE ], nbW > 1 ));
3124         iBase += isBase[ iE ];
3125       }
3126       nbTot     += *nbE;
3127       nbTotBase += nbBase;
3128     }
3129     if ( nbTotBase == 0 )
3130       return false;
3131
3132     // IPAL53099, 54416. Set correct neighbors to INTERNAL EDGEs
3133     if ( nbW > 1 )
3134     {
3135       int iFirst = 0, iLast;
3136       for ( nbE = nbEdgesInWires.begin(); nbE != nbEdgesInWires.end(); ++nbE )
3137       {
3138         iLast = iFirst + *nbE - 1;
3139         bool isConnectOk = ( edges[ iFirst ].IsConnected( edges, 0 ) &&
3140                              edges[ iFirst ].IsConnected( edges, 1 ));
3141         if ( !isConnectOk )
3142         {
3143           for ( iE = iFirst; iE <= iLast; ++iE )
3144           {
3145             if ( !edges[ iE ]._isBase )
3146               continue;
3147             int* iNei[] = { & edges[ iE ]._iL,
3148                             & edges[ iE ]._iR };
3149             for ( int iV = 0; iV < 2; ++iV )
3150             {
3151               if ( edges[ iE ].IsConnected( edges, iV ))
3152                 continue; // Ok - connected to a neighbor EDGE
3153
3154               // look for a connected EDGE
3155               bool found = false;
3156               for ( int iE2 = 0, nbE = edges.size(); iE2 < nbE &&   !found; ++iE2 )
3157                 if (( iE2 != iE ) &&
3158                     ( found = edges[ iE ].IsConnected( edges[ iE2 ], iV )))
3159                 {
3160                   *iNei[ iV ] = edges[ iE2 ]._iBase;
3161                 }
3162               if ( !found )
3163                 *iNei[ iV ] = edges[ iE ]._iBase; // connect to self
3164             }
3165           }
3166         }
3167         iFirst += *nbE;
3168       }
3169     }
3170     return edges.size();
3171   }
3172
3173   //--------------------------------------------------------------------------------
3174   /*!
3175    * \brief Return number of faces sharing given edges
3176    */
3177   // int nbAdjacentFaces( const std::vector< EdgeWithNeighbors >&          edges,
3178   //                      const TopTools_IndexedDataMapOfShapeListOfShape& facesOfEdge )
3179   // {
3180   //   TopTools_MapOfShape adjFaces;
3181
3182   //   for ( size_t i = 0; i < edges.size(); ++i )
3183   //   {
3184   //     TopTools_ListIteratorOfListOfShape faceIt( facesOfEdge.FindFromKey( edges[i]._edge ));
3185   //     for ( ; faceIt.More(); faceIt.Next() )
3186   //       adjFaces.Add( faceIt.Value() );
3187   //   }
3188   //   return adjFaces.Extent();
3189   // }
3190 }
3191
3192 //================================================================================
3193 /*!
3194  * \brief Return true if the algorithm can mesh this shape
3195  *  \param [in] aShape - shape to check
3196  *  \param [in] toCheckAll - if true, this check returns OK if all shapes are OK,
3197  *              else, returns OK if at least one shape is OK
3198  */
3199 //================================================================================
3200
3201 bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckAll)
3202 {
3203   TopExp_Explorer sExp( shape, TopAbs_SOLID );
3204   if ( !sExp.More() )
3205     return false;
3206
3207   for ( ; sExp.More(); sExp.Next() )
3208   {
3209     // check nb shells
3210     TopoDS_Shape shell;
3211     TopExp_Explorer shExp( sExp.Current(), TopAbs_SHELL );
3212     while ( shExp.More() ) {
3213       shell = shExp.Current();
3214       shExp.Next();
3215       if ( shExp.More() && BRep_Tool::IsClosed( shExp.Current() ))
3216         shell.Nullify();
3217     }
3218     if ( shell.IsNull() ) {
3219       if ( toCheckAll ) return false;
3220       continue;
3221     }
3222     // get all faces
3223     TopTools_IndexedMapOfShape allFaces;
3224     TopExp::MapShapes( sExp.Current(), TopAbs_FACE, allFaces );
3225     if ( allFaces.Extent() < 3 ) {
3226       if ( toCheckAll ) return false;
3227       continue;
3228     }
3229     // is a box?
3230     if ( allFaces.Extent() == 6 )
3231     {
3232       TopTools_IndexedMapOfOrientedShape map;
3233       bool isBox = SMESH_Block::FindBlockShapes( TopoDS::Shell( shell ),
3234                                                  TopoDS_Vertex(), TopoDS_Vertex(), map );
3235       if ( isBox ) {
3236         if ( !toCheckAll ) return true;
3237         continue;
3238       }
3239     }
3240 #ifdef _DEBUG_
3241     TopTools_IndexedMapOfShape allShapes; // usage: allShapes.FindIndex( s )
3242     TopExp::MapShapes( shape, allShapes );
3243 #endif
3244
3245     TopTools_IndexedDataMapOfShapeListOfShape facesOfEdge;
3246     TopTools_ListIteratorOfListOfShape faceIt;
3247     TopExp::MapShapesAndAncestors( sExp.Current(), TopAbs_EDGE, TopAbs_FACE , facesOfEdge );
3248     if ( facesOfEdge.IsEmpty() ) {
3249       if ( toCheckAll ) return false;
3250       continue;
3251     }
3252
3253     typedef vector< EdgeWithNeighbors > TEdgeWithNeighborsVec;
3254     vector< TEdgeWithNeighborsVec > faceEdgesVec( allFaces.Extent() + 1 );
3255     const size_t nbEdgesMax = facesOfEdge.Extent() * 2; // there can be seam EDGEs
3256     TopTools_IndexedMapOfShape* facesOfSide = new TopTools_IndexedMapOfShape[ nbEdgesMax ];
3257     SMESHUtils::ArrayDeleter<TopTools_IndexedMapOfShape> delFacesOfSide( facesOfSide );
3258
3259     // try to use each face as a bottom one
3260     bool prismDetected = false;
3261     vector< PrismSide > sides;
3262     for ( int iF = 1; iF < allFaces.Extent() && !prismDetected; ++iF )
3263     {
3264       const TopoDS_Face& botF = TopoDS::Face( allFaces( iF ));
3265
3266       TEdgeWithNeighborsVec& botEdges = faceEdgesVec[ iF ];
3267       if ( botEdges.empty() )
3268         if ( !getEdges( botF, botEdges, facesOfEdge, /*noHoles=*/false ))
3269           break;
3270
3271       int nbBase = 0;
3272       for ( size_t iS = 0; iS < botEdges.size(); ++iS )
3273         nbBase += botEdges[ iS ]._isBase;
3274
3275       if ( allFaces.Extent()-1 <= nbBase )
3276         continue; // all faces are adjacent to botF - no top FACE
3277
3278       // init data of side FACEs
3279       sides.clear();
3280       sides.resize( nbBase );
3281       size_t iS = 0;
3282       for ( size_t iE = 0; iE < botEdges.size(); ++iE )
3283       {
3284         if ( !botEdges[ iE ]._isBase )
3285           continue;
3286         sides[ iS ]._topEdge    = botEdges[ iE ]._edge;
3287         sides[ iS ]._face       = botF;
3288         sides[ iS ]._leftSide   = & sides[ botEdges[ iE ]._iR ];
3289         sides[ iS ]._rightSide  = & sides[ botEdges[ iE ]._iL ];
3290         sides[ iS ]._isInternal = botEdges[ iE ].IsInternal();
3291         sides[ iS ]._faces      = & facesOfSide[ iS ];
3292         sides[ iS ]._faces->Clear();
3293         ++iS;
3294       }
3295
3296       bool isOK = true; // ok for a current botF
3297       bool hasAdvanced = true; // is new data found in a current loop
3298       int  nbFoundSideFaces = 0;
3299       for ( int iLoop = 0; isOK && hasAdvanced; ++iLoop )
3300       {
3301         hasAdvanced = false;
3302         for ( size_t iS = 0; iS < sides.size() && isOK; ++iS )
3303         {
3304           PrismSide& side = sides[ iS ];
3305           if ( side._face.IsNull() )
3306             continue; // probably the prism top face is the last of side._faces
3307
3308           if ( side._topEdge.IsNull() )
3309           {
3310             // find vertical EDGEs --- EDGEs shared with neighbor side FACEs
3311             for ( int is2nd = 0; is2nd < 2 && isOK; ++is2nd ) // 2 adjacent neighbors
3312             {
3313               const PrismSide* adjSide = is2nd ? side._rightSide : side._leftSide;
3314               if ( side._isInternal )
3315               {
3316                 const TopoDS_Vertex& V = side.Vertex( side._iBotEdge, is2nd );
3317                 bool lHasV = side._leftSide ->HasVertex( V );
3318                 bool rHasV = side._rightSide->HasVertex( V );
3319                 if ( lHasV == rHasV )
3320                   adjSide = ( &side == side._leftSide ) ? side._rightSide : side._leftSide;
3321                 else
3322                   adjSide = ( rHasV ) ? side._rightSide : side._leftSide;
3323               }
3324               int di = is2nd ? 1 : -1;
3325               for ( size_t i = 1; i < side._edges->size(); ++i )
3326               {
3327                 int iE = SMESH_MesherHelper::WrapIndex( i*di + side._iBotEdge, side._edges->size());
3328                 if ( side._isCheckedEdge[ iE ] ) continue;
3329                 const TopoDS_Edge&               vertE = side.Edge( iE );
3330                 const TopTools_ListOfShape& neighborFF = facesOfEdge.FindFromKey( vertE );
3331                 bool isEdgeShared = (( adjSide->IsSideFace( neighborFF, side._face,
3332                                                             side._isInternal )) ||
3333                                      ( adjSide == &side &&
3334                                        side._face.IsSame( getAnotherFace( side._face,
3335                                                                           neighborFF ))));
3336                 if ( isEdgeShared ) // vertE is shared with adjSide
3337                 {
3338                   hasAdvanced = true;
3339                   side._isCheckedEdge[ iE ] = true;
3340                   side._nbCheckedEdges++;
3341                   int nbNotCheckedE = side._edges->size() - side._nbCheckedEdges;
3342                   if ( nbNotCheckedE == 1 )
3343                     break;
3344                 }
3345                 else
3346                 {
3347                   if ( i == 1 && iLoop == 0 ) isOK = false;
3348                   break;
3349                 }
3350               }
3351             }
3352             // find a top EDGE
3353             int nbNotCheckedE = side._edges->size() - side._nbCheckedEdges;
3354             if ( nbNotCheckedE == 1 )
3355             {
3356               vector<bool>::iterator ii = std::find( side._isCheckedEdge.begin(),
3357                                                      side._isCheckedEdge.end(), false );
3358               if ( ii != side._isCheckedEdge.end() )
3359               {
3360                 size_t iE = std::distance( side._isCheckedEdge.begin(), ii );
3361                 side._topEdge = side.Edge( iE );
3362               }
3363             }
3364             isOK = ( nbNotCheckedE >= 1 );
3365           }
3366           else //if ( !side._topEdge.IsNull() )
3367           {
3368             // get a next face of a side
3369             const TopoDS_Shape& f = getAnotherFace( side._face, side._topEdge, facesOfEdge );
3370             side._faces->Add( f );
3371             bool stop = false;
3372             if ( f.IsSame( side._face ) || // _topEdge is a seam
3373                  SMESH_MesherHelper::Count( f, TopAbs_WIRE, false ) != 1 )
3374             {
3375               stop = true;
3376             }
3377             else if ( side._leftSide != & side && // not closed side face
3378                       side._leftSide->_faces->Contains( f ))
3379             {
3380               stop = true; // probably f is the prism top face
3381               side._leftSide->_face.Nullify();
3382               side._leftSide->_topEdge.Nullify();
3383             }
3384             else if ( side._rightSide != & side &&
3385                       side._rightSide->_faces->Contains( f ))
3386             {
3387               stop = true; // probably f is the prism top face
3388               side._rightSide->_face.Nullify();
3389               side._rightSide->_topEdge.Nullify();
3390             }
3391             if ( stop )
3392             {
3393               side._face.Nullify();
3394               side._topEdge.Nullify();
3395               continue;
3396             }
3397             side._face  = TopoDS::Face( f );
3398             int faceID  = allFaces.FindIndex( side._face );
3399             side._edges = & faceEdgesVec[ faceID ];
3400             if ( side._edges->empty() )
3401               if ( !getEdges( side._face, * side._edges, facesOfEdge, /*noHoles=*/true ))
3402                 break;
3403             const int nbE = side._edges->size();
3404             if ( nbE >= 4 )
3405             {
3406               hasAdvanced = true;
3407               ++nbFoundSideFaces;
3408               side._iBotEdge = side.FindEdge( side._topEdge );
3409               side._isCheckedEdge.clear();
3410               side._isCheckedEdge.resize( nbE, false );
3411               side._isCheckedEdge[ side._iBotEdge ] = true;
3412               side._nbCheckedEdges = 1; // bottom EDGE is known
3413             }
3414             else // probably a triangular top face found
3415             {
3416               side._face.Nullify();
3417             }
3418             side._topEdge.Nullify();
3419             isOK = ( !side._edges->empty() || side._faces->Extent() > 1 );
3420
3421           } //if ( !side._topEdge.IsNull() )
3422
3423         } // loop on prism sides
3424
3425         if ( nbFoundSideFaces > allFaces.Extent() )
3426         {
3427           isOK = false;
3428         }
3429         if ( iLoop > allFaces.Extent() * 10 )
3430         {
3431           isOK = false;
3432 #ifdef _DEBUG_
3433           cerr << "BUG: infinite loop in StdMeshers_Prism_3D::IsApplicable()" << endl;
3434 #endif
3435         }
3436       } // while hasAdvanced
3437
3438       if ( isOK && sides[0]._faces->Extent() > 1 )
3439       {
3440         const int nbFaces = sides[0]._faces->Extent();
3441         if ( botEdges.size() == 1 ) // cylinder
3442         {
3443           prismDetected = ( nbFaces == allFaces.Extent()-1 );
3444         }
3445         else
3446         {
3447           // check that all face columns end up at the same top face
3448           const TopoDS_Shape& topFace = sides[0]._faces->FindKey( nbFaces );
3449           size_t iS;
3450           for ( iS = 1; iS < sides.size(); ++iS )
3451             if ( ! sides[ iS ]._faces->Contains( topFace ))
3452               break;
3453           if (( prismDetected = ( iS == sides.size() )))
3454           {
3455             // check that bottom and top faces has equal nb of edges
3456             TEdgeWithNeighborsVec& topEdges = faceEdgesVec[ allFaces.FindIndex( topFace )];
3457             if ( topEdges.empty() )
3458               getEdges( TopoDS::Face( topFace ), topEdges, facesOfEdge, /*noHoles=*/false );
3459             prismDetected = ( botEdges.size() == topEdges.size() );
3460           }
3461         }
3462       }
3463     } // loop on allFaces
3464
3465     if ( !prismDetected && toCheckAll ) return false;
3466     if ( prismDetected && !toCheckAll ) return true;
3467
3468   } // loop on solids
3469
3470   return toCheckAll;
3471 }
3472
3473 namespace Prism_3D
3474 {
3475   //================================================================================
3476   /*!
3477    * \brief Return true if this node and other one belong to one face
3478    */
3479   //================================================================================
3480
3481   bool Prism_3D::TNode::IsNeighbor( const Prism_3D::TNode& other ) const
3482   {
3483     if ( !other.myNode || !myNode ) return false;
3484
3485     SMDS_ElemIteratorPtr fIt = other.myNode->GetInverseElementIterator(SMDSAbs_Face);
3486     while ( fIt->more() )
3487       if ( fIt->next()->GetNodeIndex( myNode ) >= 0 )
3488         return true;
3489     return false;
3490   }
3491
3492   //================================================================================
3493   /*!
3494    * \brief Prism initialization
3495    */
3496   //================================================================================
3497
3498   void TPrismTopo::Clear()
3499   {
3500     myShape3D.Nullify();
3501     myTop.Nullify();
3502     myBottom.Nullify();
3503     myWallQuads.clear();
3504     myBottomEdges.clear();
3505     myNbEdgesInWires.clear();
3506     myWallQuads.clear();
3507   }
3508
3509   //================================================================================
3510   /*!
3511    * \brief Set upside-down
3512    */
3513   //================================================================================
3514
3515   void TPrismTopo::SetUpsideDown()
3516   {
3517     std::swap( myBottom, myTop );
3518     myBottomEdges.clear();
3519     std::reverse( myBottomEdges.begin(), myBottomEdges.end() );
3520     for ( size_t i = 0; i < myWallQuads.size(); ++i )
3521     {
3522       myWallQuads[i].reverse();
3523       TQuadList::iterator q = myWallQuads[i].begin();
3524       for ( ; q != myWallQuads[i].end(); ++q )
3525       {
3526         (*q)->shift( 2, /*keepUnitOri=*/true );
3527       }
3528       myBottomEdges.push_back( myWallQuads[i].front()->side[ QUAD_BOTTOM_SIDE ].grid->Edge(0) );
3529     }
3530   }
3531
3532 } // namespace Prism_3D
3533
3534 //================================================================================
3535 /*!
3536  * \brief Constructor. Initialization is needed
3537  */
3538 //================================================================================
3539
3540 StdMeshers_PrismAsBlock::StdMeshers_PrismAsBlock()
3541 {
3542   mySide = 0;
3543 }
3544
3545 StdMeshers_PrismAsBlock::~StdMeshers_PrismAsBlock()
3546 {
3547   Clear();
3548 }
3549 void StdMeshers_PrismAsBlock::Clear()
3550 {
3551   myHelper = 0;
3552   myShapeIDMap.Clear();
3553   myError.reset();
3554
3555   if ( mySide ) {
3556     delete mySide; mySide = 0;
3557   }
3558   myParam2ColumnMaps.clear();
3559   myShapeIndex2ColumnMap.clear();
3560 }
3561
3562 //=======================================================================
3563 //function : initPrism
3564 //purpose  : Analyse shape geometry and mesh.
3565 //           If there are triangles on one of faces, it becomes 'bottom'.
3566 //           thePrism.myBottom can be already set up.
3567 //=======================================================================
3568
3569 bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism,
3570                                     const TopoDS_Shape&   theShape3D,
3571                                     const bool            selectBottom)
3572 {
3573   myHelper->SetSubShape( theShape3D );
3574
3575   SMESH_subMesh* mainSubMesh = myHelper->GetMesh()->GetSubMeshContaining( theShape3D );
3576   if ( !mainSubMesh ) return toSM( error(COMPERR_BAD_INPUT_MESH,"Null submesh of shape3D"));
3577
3578   // detect not-quad FACE sub-meshes of the 3D SHAPE
3579   list< SMESH_subMesh* > notQuadGeomSubMesh;
3580   list< SMESH_subMesh* > notQuadElemSubMesh;
3581   list< SMESH_subMesh* > meshedSubMesh;
3582   int nbFaces = 0;
3583   //
3584   SMESH_subMeshIteratorPtr smIt = mainSubMesh->getDependsOnIterator(false,true);
3585   while ( smIt->more() )
3586   {
3587     SMESH_subMesh* sm = smIt->next();
3588     const TopoDS_Shape& face = sm->GetSubShape();
3589     if      ( face.ShapeType() > TopAbs_FACE ) break;
3590     else if ( face.ShapeType() < TopAbs_FACE ) continue;
3591     nbFaces++;
3592
3593     // is quadrangle FACE?
3594     list< TopoDS_Edge > orderedEdges;
3595     list< int >         nbEdgesInWires;
3596     int nbWires = SMESH_Block::GetOrderedEdges( TopoDS::Face( face ), orderedEdges,
3597                                                 nbEdgesInWires );
3598     if ( nbWires != 1 || nbEdgesInWires.front() != 4 )
3599       notQuadGeomSubMesh.push_back( sm );
3600
3601     // look for a not structured sub-mesh
3602     if ( !sm->IsEmpty() )
3603     {
3604       meshedSubMesh.push_back( sm );
3605       if ( !myHelper->IsSameElemGeometry( sm->GetSubMeshDS(), SMDSGeom_QUADRANGLE ) ||
3606            !myHelper->IsStructured      ( sm ))
3607         notQuadElemSubMesh.push_back( sm );
3608     }
3609   }
3610
3611   int nbNotQuadMeshed = notQuadElemSubMesh.size();
3612   int       nbNotQuad = notQuadGeomSubMesh.size();
3613   bool     hasNotQuad = ( nbNotQuad || nbNotQuadMeshed );
3614
3615   // detect bad cases
3616   if ( nbNotQuadMeshed > 2 )
3617   {
3618     return toSM( error(COMPERR_BAD_INPUT_MESH,
3619                        TCom("More than 2 faces with not quadrangle elements: ")
3620                        <<nbNotQuadMeshed));
3621   }
3622   if ( nbNotQuad > 2 || !thePrism.myBottom.IsNull() )
3623   {
3624     // Issue 0020843 - one of side FACEs is quasi-quadrilateral (not 4 EDGEs).
3625     // Remove from notQuadGeomSubMesh faces meshed with regular grid
3626     int nbQuasiQuads = removeQuasiQuads( notQuadGeomSubMesh, myHelper,
3627                                          TQuadrangleAlgo::instance(this,myHelper) );
3628     nbNotQuad -= nbQuasiQuads;
3629     if ( nbNotQuad > 2 )
3630       return toSM( error(COMPERR_BAD_SHAPE,
3631                          TCom("More than 2 not quadrilateral faces: ") <<nbNotQuad));
3632     hasNotQuad = ( nbNotQuad || nbNotQuadMeshed );
3633   }
3634
3635   // Analyse mesh and topology of FACEs: choose the bottom sub-mesh.
3636   // If there are not quadrangle FACEs, they are top and bottom ones.
3637   // Not quadrangle FACEs must be only on top and bottom.
3638
3639   SMESH_subMesh * botSM = 0;
3640   SMESH_subMesh * topSM = 0;
3641
3642   if ( hasNotQuad ) // can choose a bottom FACE
3643   {
3644     if ( nbNotQuadMeshed > 0 ) botSM = notQuadElemSubMesh.front();
3645     else                       botSM = notQuadGeomSubMesh.front();
3646     if ( nbNotQuadMeshed > 1 ) topSM = notQuadElemSubMesh.back();
3647     else if ( nbNotQuad  > 1 ) topSM = notQuadGeomSubMesh.back();
3648
3649     if ( topSM == botSM ) {
3650       if ( nbNotQuadMeshed > 1 ) topSM = notQuadElemSubMesh.front();
3651       else                       topSM = notQuadGeomSubMesh.front();
3652     }
3653
3654     // detect mesh triangles on wall FACEs
3655     if ( nbNotQuad == 2 && nbNotQuadMeshed > 0 ) {
3656       bool ok = false;
3657       if ( nbNotQuadMeshed == 1 )
3658         ok = ( find( notQuadGeomSubMesh.begin(),
3659                      notQuadGeomSubMesh.end(), botSM ) != notQuadGeomSubMesh.end() );
3660       else
3661         ok = ( notQuadGeomSubMesh == notQuadElemSubMesh );
3662       if ( !ok )
3663         return toSM( error(COMPERR_BAD_INPUT_MESH,
3664                            "Side face meshed with not quadrangle elements"));
3665     }
3666   }
3667
3668   thePrism.myNotQuadOnTop = ( nbNotQuadMeshed > 1 );
3669
3670   // use thePrism.myBottom
3671   if ( !thePrism.myBottom.IsNull() )
3672   {
3673     if ( botSM ) { // <-- not quad geom or mesh on botSM
3674       if ( ! botSM->GetSubShape().IsSame( thePrism.myBottom )) {
3675         std::swap( botSM, topSM );
3676         if ( !botSM || ! botSM->GetSubShape().IsSame( thePrism.myBottom )) {
3677           if ( !selectBottom )
3678             return toSM( error( COMPERR_BAD_INPUT_MESH,
3679                                 "Incompatible non-structured sub-meshes"));
3680           std::swap( botSM, topSM );
3681           thePrism.myBottom = TopoDS::Face( botSM->GetSubShape() );
3682         }
3683       }
3684     }
3685     else if ( !selectBottom ) {
3686       botSM = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom );
3687     }
3688   }
3689   if ( !botSM ) // find a proper bottom
3690   {
3691     bool savedSetErrorToSM = mySetErrorToSM;
3692     mySetErrorToSM = false; // ignore errors in initPrism()
3693
3694     // search among meshed FACEs
3695     list< SMESH_subMesh* >::iterator sm = meshedSubMesh.begin();
3696     for ( ; !botSM && sm != meshedSubMesh.end(); ++sm )
3697     {
3698       thePrism.Clear();
3699       botSM             = *sm;
3700       thePrism.myBottom = TopoDS::Face( botSM->GetSubShape() );
3701       if ( !initPrism( thePrism, theShape3D, /*selectBottom=*/false ))
3702         botSM = NULL;
3703     }
3704     // search among all FACEs
3705     for ( TopExp_Explorer f( theShape3D, TopAbs_FACE ); !botSM && f.More(); f.Next() )
3706     {
3707       int minNbFaces = 2 + myHelper->Count( f.Current(), TopAbs_EDGE, false);
3708       if ( nbFaces < minNbFaces) continue;
3709       thePrism.Clear();
3710       thePrism.myBottom = TopoDS::Face( f.Current() );
3711       botSM             = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom );
3712       if ( !initPrism( thePrism, theShape3D, /*selectBottom=*/false ))
3713         botSM = NULL;
3714     }
3715     mySetErrorToSM = savedSetErrorToSM;
3716     return botSM ? true : toSM( error( COMPERR_BAD_SHAPE ));
3717   }
3718
3719   // find vertex 000 - the one with smallest coordinates (for easy DEBUG :-)
3720   TopoDS_Vertex V000;
3721   double minVal = DBL_MAX, minX = 0, val;
3722   for ( TopExp_Explorer exp( botSM->GetSubShape(), TopAbs_VERTEX );
3723         exp.More(); exp.Next() )
3724   {
3725     const TopoDS_Vertex& v = TopoDS::Vertex( exp.Current() );
3726     gp_Pnt P = BRep_Tool::Pnt( v );
3727     val = P.X() + P.Y() + P.Z();
3728     if ( val < minVal || ( val == minVal && P.X() < minX )) {
3729       V000 = v;
3730       minVal = val;
3731       minX = P.X();
3732     }
3733   }
3734
3735   thePrism.myShape3D = theShape3D;
3736   if ( thePrism.myBottom.IsNull() )
3737     thePrism.myBottom  = TopoDS::Face( botSM->GetSubShape() );
3738   thePrism.myBottom.Orientation( myHelper->GetSubShapeOri( theShape3D, thePrism.myBottom ));
3739   thePrism.myTop.   Orientation( myHelper->GetSubShapeOri( theShape3D, thePrism.myTop ));
3740
3741   // Get ordered bottom edges
3742   TopoDS_Face reverseBottom = // to have order of top EDGEs as in the top FACE
3743     TopoDS::Face( thePrism.myBottom.Reversed() );
3744   SMESH_Block::GetOrderedEdges( reverseBottom,
3745                                 thePrism.myBottomEdges,
3746                                 thePrism.myNbEdgesInWires, V000 );
3747
3748   // Get Wall faces corresponding to the ordered bottom edges and the top FACE
3749   if ( !getWallFaces( thePrism, nbFaces )) // it also sets thePrism.myTop
3750     return false; //toSM( error(COMPERR_BAD_SHAPE, "Can't find side faces"));
3751
3752   if ( topSM )
3753   {
3754     if ( !thePrism.myTop.IsSame( topSM->GetSubShape() ))
3755       return toSM( error
3756                    (notQuadGeomSubMesh.empty() ? COMPERR_BAD_INPUT_MESH : COMPERR_BAD_SHAPE,
3757                     "Non-quadrilateral faces are not opposite"));
3758
3759     // check that the found top and bottom FACEs are opposite
3760     TopTools_IndexedMapOfShape topEdgesMap( thePrism.myBottomEdges.size() );
3761     TopExp::MapShapes( thePrism.myTop, topEdgesMap );
3762     list< TopoDS_Edge >::iterator edge = thePrism.myBottomEdges.begin();
3763     for ( ; edge != thePrism.myBottomEdges.end(); ++edge )
3764       if ( topEdgesMap.Contains( *edge ))
3765         return toSM( error
3766                      (notQuadGeomSubMesh.empty() ? COMPERR_BAD_INPUT_MESH : COMPERR_BAD_SHAPE,
3767                       "Non-quadrilateral faces are not opposite"));
3768   }
3769
3770   if ( thePrism.myBottomEdges.size() > thePrism.myWallQuads.size() )
3771   {
3772     // composite bottom sides => set thePrism upside-down
3773     thePrism.SetUpsideDown();
3774   }
3775
3776   return true;
3777 }
3778
3779 //================================================================================
3780 /*!
3781  * \brief Initialization.
3782  * \param helper - helper loaded with mesh and 3D shape
3783  * \param thePrism - a prism data
3784  * \retval bool - false if a mesh or a shape are KO
3785  */
3786 //================================================================================
3787
3788 bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper*         helper,
3789                                    const Prism_3D::TPrismTopo& thePrism)
3790 {
3791   myHelper = helper;
3792   SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
3793   SMESH_Mesh*     mesh = myHelper->GetMesh();
3794
3795   if ( mySide ) {
3796     delete mySide; mySide = 0;
3797   }
3798   vector< TSideFace* >         sideFaces( NB_WALL_FACES, 0 );
3799   vector< pair< double, double> > params( NB_WALL_FACES );
3800   mySide = new TSideFace( *mesh, sideFaces, params );
3801
3802
3803   SMESH_Block::init();
3804   myShapeIDMap.Clear();
3805   myShapeIndex2ColumnMap.clear();
3806
3807   int wallFaceIds[ NB_WALL_FACES ] = { // to walk around a block
3808     SMESH_Block::ID_Fx0z, SMESH_Block::ID_F1yz,
3809     SMESH_Block::ID_Fx1z, SMESH_Block::ID_F0yz
3810   };
3811
3812   myError = SMESH_ComputeError::New();
3813
3814   myNotQuadOnTop = thePrism.myNotQuadOnTop;
3815
3816   // Find columns of wall nodes and calculate edges' lengths
3817   // --------------------------------------------------------
3818
3819   myParam2ColumnMaps.clear();
3820   myParam2ColumnMaps.resize( thePrism.myBottomEdges.size() ); // total nb edges
3821
3822   size_t iE, nbEdges = thePrism.myNbEdgesInWires.front(); // nb outer edges
3823   vector< double > edgeLength( nbEdges );
3824   multimap< double, int > len2edgeMap;
3825
3826   // for each EDGE: either split into several parts, or join with several next EDGEs
3827   vector<int> nbSplitPerEdge( nbEdges, 0 );
3828   vector<int> nbUnitePerEdge( nbEdges, 0 ); // -1 means "joined to a previous"
3829
3830   // consider continuous straight EDGEs as one side
3831   const int nbSides = countNbSides( thePrism, nbUnitePerEdge, edgeLength );
3832
3833   list< TopoDS_Edge >::const_iterator edgeIt = thePrism.myBottomEdges.begin();
3834   for ( iE = 0; iE < nbEdges; ++iE, ++edgeIt )
3835   {
3836     TParam2ColumnMap & faceColumns = myParam2ColumnMaps[ iE ];
3837
3838     Prism_3D::TQuadList::const_iterator quad = thePrism.myWallQuads[ iE ].begin();
3839     for ( ; quad != thePrism.myWallQuads[ iE ].end(); ++quad )
3840     {
3841       const TopoDS_Edge& quadBot = (*quad)->side[ QUAD_BOTTOM_SIDE ].grid->Edge( 0 );
3842       if ( !myHelper->LoadNodeColumns( faceColumns, (*quad)->face, quadBot, meshDS ))
3843         return error(COMPERR_BAD_INPUT_MESH, TCom("Can't find regular quadrangle mesh ")
3844                      << "on a side face #" << MeshDS()->ShapeToIndex( (*quad)->face ));
3845     }
3846     SHOWYXZ("\np1 F " <<iE, gpXYZ(faceColumns.begin()->second.front() ));
3847     SHOWYXZ("p2 F "   <<iE, gpXYZ(faceColumns.rbegin()->second.front() ));
3848     SHOWYXZ("V First "<<iE, BRep_Tool::Pnt( TopExp::FirstVertex(*edgeIt,true )));
3849
3850     if ( nbSides < NB_WALL_FACES ) // fill map used to split faces
3851       len2edgeMap.insert( make_pair( edgeLength[ iE ], iE )); // sort edges by length
3852   }
3853   // Load columns of internal edges (forming holes)
3854   // and fill map ShapeIndex to TParam2ColumnMap for them
3855   for ( ; edgeIt != thePrism.myBottomEdges.end() ; ++edgeIt, ++iE )
3856   {
3857     TParam2ColumnMap & faceColumns = myParam2ColumnMaps[ iE ];
3858
3859     Prism_3D::TQuadList::const_iterator quad = thePrism.myWallQuads[ iE ].begin();
3860     for ( ; quad != thePrism.myWallQuads[ iE ].end(); ++quad )
3861     {
3862       const TopoDS_Edge& quadBot = (*quad)->side[ QUAD_BOTTOM_SIDE ].grid->Edge( 0 );
3863       if ( !myHelper->LoadNodeColumns( faceColumns, (*quad)->face, quadBot, meshDS ))
3864         return error(COMPERR_BAD_INPUT_MESH, TCom("Can't find regular quadrangle mesh ")
3865                      << "on a side face #" << MeshDS()->ShapeToIndex( (*quad)->face ));
3866     }
3867     if ( !faceColumns.empty() && (int)faceColumns.begin()->second.size() != VerticalSize() )
3868       return error(COMPERR_BAD_INPUT_MESH, "Different 'vertical' discretization");
3869
3870     // edge columns
3871     int id = MeshDS()->ShapeToIndex( *edgeIt );
3872     bool isForward = true; // meaningless for intenal wires
3873     myShapeIndex2ColumnMap[ id ] = make_pair( & faceColumns, isForward );
3874     // columns for vertices
3875     // 1
3876     const SMDS_MeshNode* n0 = faceColumns.begin()->second.front();
3877     id = n0->getshapeId();
3878     myShapeIndex2ColumnMap[ id ] = make_pair( & faceColumns, isForward );
3879     // 2
3880     const SMDS_MeshNode* n1 = faceColumns.rbegin()->second.front();
3881     id = n1->getshapeId();
3882     myShapeIndex2ColumnMap[ id ] = make_pair( & faceColumns, isForward );
3883
3884     // SHOWYXZ("\np1 F " <<iE, gpXYZ(faceColumns.begin()->second.front() ));
3885     // SHOWYXZ("p2 F "   <<iE, gpXYZ(faceColumns.rbegin()->second.front() ));
3886     // SHOWYXZ("V First "<<iE, BRep_Tool::Pnt( TopExp::FirstVertex(*edgeIt,true )));
3887   }
3888
3889   // Create 4 wall faces of a block
3890   // -------------------------------
3891
3892   if ( nbSides <= NB_WALL_FACES ) // ************* Split faces if necessary
3893   {
3894     if ( nbSides != NB_WALL_FACES ) // define how to split
3895     {
3896       if ( len2edgeMap.size() != nbEdges )
3897         RETURN_BAD_RESULT("Uniqueness of edge lengths not assured");
3898
3899       multimap< double, int >::reverse_iterator maxLen_i = len2edgeMap.rbegin();
3900       multimap< double, int >::reverse_iterator midLen_i = ++len2edgeMap.rbegin();
3901
3902       double maxLen = maxLen_i->first;
3903       double midLen = ( len2edgeMap.size() == 1 ) ? 0 : midLen_i->first;
3904       switch ( nbEdges ) {
3905       case 1: // 0-th edge is split into 4 parts
3906         nbSplitPerEdge[ 0 ] = 4;
3907         break;
3908       case 2: // either the longest edge is split into 3 parts, or both edges into halves
3909         if ( maxLen / 3 > midLen / 2 ) {
3910           nbSplitPerEdge[ maxLen_i->second ] = 3;
3911         }
3912         else {
3913           nbSplitPerEdge[ maxLen_i->second ] = 2;
3914           nbSplitPerEdge[ midLen_i->second ] = 2;
3915         }
3916         break;
3917       case 3:
3918         if ( nbSides == 2 )
3919           // split longest into 3 parts
3920           nbSplitPerEdge[ maxLen_i->second ] = 3;
3921         else
3922           // split longest into halves
3923           nbSplitPerEdge[ maxLen_i->second ] = 2;
3924       }
3925     }
3926   }
3927   else // **************************** Unite faces
3928   {
3929     int nbExraFaces = nbSides - 4; // nb of faces to fuse
3930     for ( iE = 0; iE < nbEdges; ++iE )
3931     {
3932       if ( nbUnitePerEdge[ iE ] < 0 )
3933         continue;
3934       // look for already united faces
3935       for ( size_t i = iE; i < iE + nbExraFaces; ++i )
3936       {
3937         if ( nbUnitePerEdge[ i ] > 0 ) // a side including nbUnitePerEdge[i]+1 edge
3938           nbExraFaces += nbUnitePerEdge[ i ];
3939         nbUnitePerEdge[ i ] = -1;
3940       }
3941       nbUnitePerEdge[ iE ] = nbExraFaces;
3942       break;
3943     }
3944   }
3945
3946   // Create TSideFace's
3947   int iSide = 0;
3948   list< TopoDS_Edge >::const_iterator botE = thePrism.myBottomEdges.begin();
3949   for ( iE = 0; iE < nbEdges; ++iE, ++botE )
3950   {
3951     TFaceQuadStructPtr quad = thePrism.myWallQuads[ iE ].front();
3952     const int       nbSplit = nbSplitPerEdge[ iE ];
3953     const int   nbExraFaces = nbUnitePerEdge[ iE ] + 1;
3954     if ( nbSplit > 0 ) // split
3955     {
3956       vector< double > params;
3957       splitParams( nbSplit, &myParam2ColumnMaps[ iE ], params );
3958       const bool isForward =
3959         StdMeshers_PrismAsBlock::IsForwardEdge( myHelper->GetMeshDS(),
3960                                                 myParam2ColumnMaps[iE],
3961                                                 *botE, SMESH_Block::ID_Fx0z );
3962       for ( int i = 0; i < nbSplit; ++i ) {
3963         double f = ( isForward ? params[ i ]   : params[ nbSplit - i-1 ]);
3964         double l = ( isForward ? params[ i+1 ] : params[ nbSplit - i ]);
3965         TSideFace* comp = new TSideFace( *mesh, wallFaceIds[ iSide ],
3966                                          thePrism.myWallQuads[ iE ], *botE,
3967                                          &myParam2ColumnMaps[ iE ], f, l );
3968         mySide->SetComponent( iSide++, comp );
3969       }
3970     }
3971     else if ( nbExraFaces > 1 ) // unite
3972     {
3973       double u0 = 0, sumLen = 0;
3974       for ( size_t i = iE; i < iE + nbExraFaces; ++i )
3975         sumLen += edgeLength[ i ];
3976
3977       vector< TSideFace* >        components( nbExraFaces );
3978       vector< pair< double, double> > params( nbExraFaces );
3979       bool endReached = false;
3980       for ( int i = 0; i < nbExraFaces; ++i, ++botE, ++iE )
3981       {
3982         if ( iE == nbEdges )
3983         {
3984           endReached = true;
3985           botE = thePrism.myBottomEdges.begin();
3986           iE = 0;
3987         }
3988         components[ i ] = new TSideFace( *mesh, wallFaceIds[ iSide ],
3989                                          thePrism.myWallQuads[ iE ], *botE,
3990                                          &myParam2ColumnMaps[ iE ]);
3991         double u1 = u0 + edgeLength[ iE ] / sumLen;
3992         params[ i ] = make_pair( u0 , u1 );
3993         u0 = u1;
3994       }
3995       TSideFace* comp = new TSideFace( *mesh, components, params );
3996       mySide->SetComponent( iSide++, comp );
3997       if ( endReached )
3998         break;
3999       --iE; // for increment in an external loop on iE
4000       --botE;
4001     }
4002     else if ( nbExraFaces < 0 ) // skip already united face
4003     {
4004     }
4005     else // use as is
4006     {
4007       TSideFace* comp = new TSideFace( *mesh, wallFaceIds[ iSide ],
4008                                        thePrism.myWallQuads[ iE ], *botE,
4009                                        &myParam2ColumnMaps[ iE ]);
4010       mySide->SetComponent( iSide++, comp );
4011     }
4012   }
4013
4014
4015   // Fill geometry fields of SMESH_Block
4016   // ------------------------------------
4017
4018   vector< int > botEdgeIdVec;
4019   SMESH_Block::GetFaceEdgesIDs( ID_BOT_FACE, botEdgeIdVec );
4020
4021   bool isForward[NB_WALL_FACES] = { true, true, true, true };
4022   Adaptor2d_Curve2d* botPcurves[NB_WALL_FACES];
4023   Adaptor2d_Curve2d* topPcurves[NB_WALL_FACES];
4024
4025   for ( int iF = 0; iF < NB_WALL_FACES; ++iF )
4026   {
4027     TSideFace * sideFace = mySide->GetComponent( iF );
4028     if ( !sideFace )
4029       RETURN_BAD_RESULT("NULL TSideFace");
4030     int fID = sideFace->FaceID(); // in-block ID
4031
4032     // fill myShapeIDMap
4033     if ( sideFace->InsertSubShapes( myShapeIDMap ) != 8 &&
4034          !sideFace->IsComplex())
4035       MESSAGE( ": Warning : InsertSubShapes() < 8 on side " << iF );
4036
4037     // side faces geometry
4038     Adaptor2d_Curve2d* pcurves[NB_WALL_FACES];
4039     if ( !sideFace->GetPCurves( pcurves ))
4040       RETURN_BAD_RESULT("TSideFace::GetPCurves() failed");
4041
4042     SMESH_Block::TFace& tFace = myFace[ fID - ID_FirstF ];
4043     tFace.Set( fID, sideFace->Surface(), pcurves, isForward );
4044
4045     SHOWYXZ( endl<<"F "<< iF << " id " << fID << " FRW " << sideFace->IsForward(), sideFace->Value(0,0));
4046     // edges 3D geometry
4047     vector< int > edgeIdVec;
4048     SMESH_Block::GetFaceEdgesIDs( fID, edgeIdVec );
4049     for ( int isMax = 0; isMax < 2; ++isMax ) {
4050       {
4051         int eID = edgeIdVec[ isMax ];
4052         SMESH_Block::TEdge& tEdge = myEdge[ eID - ID_FirstE ];
4053         tEdge.Set( eID, sideFace->HorizCurve(isMax), true);
4054         SHOWYXZ(eID<<" HOR"<<isMax<<"(0)", sideFace->HorizCurve(isMax)->Value(0));
4055         SHOWYXZ(eID<<" HOR"<<isMax<<"(1)", sideFace->HorizCurve(isMax)->Value(1));
4056       }
4057       {
4058         int eID = edgeIdVec[ isMax+2 ];
4059         SMESH_Block::TEdge& tEdge = myEdge[ eID - ID_FirstE  ];
4060         tEdge.Set( eID, sideFace->VertiCurve(isMax), true);
4061         SHOWYXZ(eID<<" VER"<<isMax<<"(0)", sideFace->VertiCurve(isMax)->Value(0));
4062         SHOWYXZ(eID<<" VER"<<isMax<<"(1)", sideFace->VertiCurve(isMax)->Value(1));
4063
4064         // corner points
4065         vector< int > vertexIdVec;
4066         SMESH_Block::GetEdgeVertexIDs( eID, vertexIdVec );
4067         myPnt[ vertexIdVec[0] - ID_FirstV ] = tEdge.GetCurve()->Value(0).XYZ();
4068         myPnt[ vertexIdVec[1] - ID_FirstV ] = tEdge.GetCurve()->Value(1).XYZ();
4069       }
4070     }
4071     // pcurves on horizontal faces
4072     for ( iE = 0; iE < NB_WALL_FACES; ++iE ) {
4073       if ( edgeIdVec[ BOTTOM_EDGE ] == botEdgeIdVec[ iE ] ) {
4074         botPcurves[ iE ] = sideFace->HorizPCurve( false, thePrism.myBottom );
4075         topPcurves[ iE ] = sideFace->HorizPCurve( true,  thePrism.myTop );
4076         break;
4077       }
4078     }
4079     //sideFace->dumpNodes( 4 ); // debug
4080   }
4081   // horizontal faces geometry
4082   {
4083     SMESH_Block::TFace& tFace = myFace[ ID_BOT_FACE - ID_FirstF ];
4084     tFace.Set( ID_BOT_FACE, new BRepAdaptor_Surface( thePrism.myBottom ), botPcurves, isForward );
4085     SMESH_Block::Insert( thePrism.myBottom, ID_BOT_FACE, myShapeIDMap );
4086   }
4087   {
4088     SMESH_Block::TFace& tFace = myFace[ ID_TOP_FACE - ID_FirstF ];
4089     tFace.Set( ID_TOP_FACE, new BRepAdaptor_Surface( thePrism.myTop ), topPcurves, isForward );
4090     SMESH_Block::Insert( thePrism.myTop, ID_TOP_FACE, myShapeIDMap );
4091   }
4092   //faceGridToPythonDump( SMESH_Block::ID_Fxy0, 50 );
4093   //faceGridToPythonDump( SMESH_Block::ID_Fxy1 );
4094
4095   // Fill map ShapeIndex to TParam2ColumnMap
4096   // ----------------------------------------
4097
4098   list< TSideFace* > fList;
4099   list< TSideFace* >::iterator fListIt;
4100   fList.push_back( mySide );
4101   for ( fListIt = fList.begin(); fListIt != fList.end(); ++fListIt)
4102   {
4103     int nb = (*fListIt)->NbComponents();
4104     for ( int i = 0; i < nb; ++i ) {
4105       if ( TSideFace* comp = (*fListIt)->GetComponent( i ))
4106         fList.push_back( comp );
4107     }
4108     if ( TParam2ColumnMap* cols = (*fListIt)->GetColumns()) {
4109       // columns for a base edge
4110       int id = MeshDS()->ShapeToIndex( (*fListIt)->BaseEdge() );
4111       bool isForward = (*fListIt)->IsForward();
4112       myShapeIndex2ColumnMap[ id ] = make_pair( cols, isForward );
4113
4114       // columns for vertices
4115       const SMDS_MeshNode* n0 = cols->begin()->second.front();
4116       id = n0->getshapeId();
4117       myShapeIndex2ColumnMap[ id ] = make_pair( cols, isForward );
4118
4119       const SMDS_MeshNode* n1 = cols->rbegin()->second.front();
4120       id = n1->getshapeId();
4121       myShapeIndex2ColumnMap[ id ] = make_pair( cols, !isForward );
4122     }
4123   }
4124
4125 // #define SHOWYXZ(msg, xyz) { gp_Pnt p(xyz); cout << msg << " ("<< p.X() << "; " <<p.Y() << "; " <<p.Z() << ") " <<endl; }
4126
4127 //   double _u[]={ 0.1, 0.1, 0.9, 0.9 };
4128 //   double _v[]={ 0.1, 0.9, 0.1, 0.9 };
4129 //   for ( int z = 0; z < 2; ++z )
4130 //     for ( int i = 0; i < 4; ++i )
4131 //     {
4132 //       //gp_XYZ testPar(0.25, 0.25, 0), testCoord;
4133 //       int iFace = (z ? ID_TOP_FACE : ID_BOT_FACE);
4134 //       gp_XYZ testPar(_u[i], _v[i], z), testCoord;
4135 //       if ( !FacePoint( iFace, testPar, testCoord ))
4136 //         RETURN_BAD_RESULT("TEST FacePoint() FAILED");
4137 //       SHOWYXZ("IN TEST PARAM" , testPar);
4138 //       SHOWYXZ("OUT TEST CORD" , testCoord);
4139 //       if ( !ComputeParameters( testCoord, testPar , iFace))
4140 //         RETURN_BAD_RESULT("TEST ComputeParameters() FAILED");
4141 //       SHOWYXZ("OUT TEST PARAM" , testPar);
4142 //     }
4143   return true;
4144 }
4145
4146 //================================================================================
4147 /*!
4148  * \brief Return pointer to column of nodes
4149  * \param node - bottom node from which the returned column goes up
4150  * \retval const TNodeColumn* - the found column
4151  */
4152 //================================================================================
4153
4154 const TNodeColumn* StdMeshers_PrismAsBlock::GetNodeColumn(const SMDS_MeshNode* node) const
4155 {
4156   int sID = node->getshapeId();
4157
4158   map<int, pair< TParam2ColumnMap*, bool > >::const_iterator col_frw =
4159     myShapeIndex2ColumnMap.find( sID );
4160   if ( col_frw != myShapeIndex2ColumnMap.end() ) {
4161     const TParam2ColumnMap* cols = col_frw->second.first;
4162     TParam2ColumnIt u_col = cols->begin();
4163     for ( ; u_col != cols->end(); ++u_col )
4164       if ( u_col->second[ 0 ] == node )
4165         return & u_col->second;
4166   }
4167   return 0;
4168 }
4169
4170 //=======================================================================
4171 //function : GetLayersTransformation
4172 //purpose  : Return transformations to get coordinates of nodes of each layer
4173 //           by nodes of the bottom. Layer is a set of nodes at a certain step
4174 //           from bottom to top.
4175 //           Transformation to get top node from bottom ones is computed
4176 //           only if the top FACE is not meshed.
4177 //=======================================================================
4178
4179 bool StdMeshers_PrismAsBlock::GetLayersTransformation(vector<gp_Trsf> &           trsf,
4180                                                       const Prism_3D::TPrismTopo& prism) const
4181 {
4182   const bool itTopMeshed = !SubMesh( ID_BOT_FACE )->IsEmpty();
4183   const int zSize = VerticalSize();
4184   if ( zSize < 3 && !itTopMeshed ) return true;
4185   trsf.resize( zSize - 1 );
4186
4187   // Select some node columns by which we will define coordinate system of layers
4188
4189   vector< const TNodeColumn* > columns;
4190   {
4191     bool isReverse;
4192     list< TopoDS_Edge >::const_iterator edgeIt = prism.myBottomEdges.begin();
4193     for ( int iE = 0; iE < prism.myNbEdgesInWires.front(); ++iE, ++edgeIt )
4194     {
4195       if ( SMESH_Algo::isDegenerated( *edgeIt )) continue;
4196       const TParam2ColumnMap* u2colMap =
4197         GetParam2ColumnMap( MeshDS()->ShapeToIndex( *edgeIt ), isReverse );
4198       if ( !u2colMap ) return false;
4199       double f = u2colMap->begin()->first, l = u2colMap->rbegin()->first;
4200       //isReverse = ( edgeIt->Orientation() == TopAbs_REVERSED );
4201       //if ( isReverse ) swap ( f, l ); -- u2colMap takes orientation into account
4202       const int nbCol = 5;
4203       for ( int i = 0; i < nbCol; ++i )
4204       {
4205         double u = f + i/double(nbCol) * ( l - f );
4206         const TNodeColumn* col = & getColumn( u2colMap, u )->second;
4207         if ( columns.empty() || col != columns.back() )
4208           columns.push_back( col );
4209       }
4210     }
4211   }
4212
4213   // Find tolerance to check transformations
4214
4215   double tol2;
4216   {
4217     Bnd_B3d bndBox;
4218     for ( size_t i = 0; i < columns.size(); ++i )
4219       bndBox.Add( gpXYZ( columns[i]->front() ));
4220     tol2 = bndBox.SquareExtent() * 1e-5;
4221   }
4222
4223   // Compute transformations
4224
4225   int xCol = -1;
4226   gp_Trsf fromCsZ, toCs0;
4227   gp_Ax3 cs0 = getLayerCoordSys(0, columns, xCol );
4228   //double dist0 = cs0.Location().Distance( gpXYZ( (*columns[0])[0]));
4229   toCs0.SetTransformation( cs0 );
4230   for ( int z = 1; z < zSize; ++z )
4231   {
4232     gp_Ax3 csZ = getLayerCoordSys(z, columns, xCol );
4233     //double distZ = csZ.Location().Distance( gpXYZ( (*columns[0])[z]));
4234     fromCsZ.SetTransformation( csZ );
4235     fromCsZ.Invert();
4236     gp_Trsf& t = trsf[ z-1 ];
4237     t = fromCsZ * toCs0;
4238     //t.SetScaleFactor( distZ/dist0 ); - it does not work properly, wrong base point
4239
4240     // check a transformation
4241     for ( size_t i = 0; i < columns.size(); ++i )
4242     {
4243       gp_Pnt p0 = gpXYZ( (*columns[i])[0] );
4244       gp_Pnt pz = gpXYZ( (*columns[i])[z] );
4245       t.Transforms( p0.ChangeCoord() );
4246       if ( p0.SquareDistance( pz ) > tol2 )
4247       {
4248         t = gp_Trsf();
4249         return ( z == zSize - 1 ); // OK if fails only bottom->top trsf
4250       }
4251     }
4252   }
4253   return true;
4254 }
4255
4256 //================================================================================
4257 /*!
4258  * \brief Check curve orientation of a bottom edge
4259   * \param meshDS - mesh DS
4260   * \param columnsMap - node columns map of side face
4261   * \param bottomEdge - the bottom edge
4262   * \param sideFaceID - side face in-block ID
4263   * \retval bool - true if orientation coincide with in-block forward orientation
4264  */
4265 //================================================================================
4266
4267 bool StdMeshers_PrismAsBlock::IsForwardEdge(SMESHDS_Mesh*           meshDS,
4268                                             const TParam2ColumnMap& columnsMap,
4269                                             const TopoDS_Edge &     bottomEdge,
4270                                             const int               sideFaceID)
4271 {
4272   bool isForward = false;
4273   if ( SMESH_MesherHelper::IsClosedEdge( bottomEdge ))
4274   {
4275     isForward = ( bottomEdge.Orientation() == TopAbs_FORWARD );
4276   }
4277   else
4278   {
4279     const TNodeColumn&     firstCol = columnsMap.begin()->second;
4280     const SMDS_MeshNode* bottomNode = firstCol[0];
4281     TopoDS_Shape firstVertex = SMESH_MesherHelper::GetSubShapeByNode( bottomNode, meshDS );
4282     isForward = ( firstVertex.IsSame( TopExp::FirstVertex( bottomEdge, true )));
4283   }
4284   // on 2 of 4 sides first vertex is end
4285   if ( sideFaceID == ID_Fx1z || sideFaceID == ID_F0yz )
4286     isForward = !isForward;
4287   return isForward;
4288 }
4289
4290 //=======================================================================
4291 //function : faceGridToPythonDump
4292 //purpose  : Prints a script creating a normal grid on the prism side
4293 //=======================================================================
4294
4295 void StdMeshers_PrismAsBlock::faceGridToPythonDump(const SMESH_Block::TShapeID face,
4296                                                    const int                   nb)
4297 {
4298 #ifdef _DEBUG_
4299   gp_XYZ pOnF[6] = { gp_XYZ(0,0,0), gp_XYZ(0,0,1),
4300                      gp_XYZ(0,0,0), gp_XYZ(0,1,0),
4301                      gp_XYZ(0,0,0), gp_XYZ(1,0,0) };
4302   gp_XYZ p2;
4303   cout << "mesh = smesh.Mesh( 'Face " << face << "')" << endl;
4304   SMESH_Block::TFace& f = myFace[ face - ID_FirstF ];
4305   gp_XYZ params = pOnF[ face - ID_FirstF ];
4306   //const int nb = 10; // nb face rows
4307   for ( int j = 0; j <= nb; ++j )
4308   {
4309     params.SetCoord( f.GetVInd(), double( j )/ nb );
4310     for ( int i = 0; i <= nb; ++i )
4311     {
4312       params.SetCoord( f.GetUInd(), double( i )/ nb );
4313       gp_XYZ p = f.Point( params );
4314       gp_XY uv = f.GetUV( params );
4315       cout << "mesh.AddNode( " << p.X() << ", " << p.Y() << ", " << p.Z() << " )"
4316            << " # " << 1 + i + j * ( nb + 1 )
4317            << " ( " << i << ", " << j << " ) "
4318            << " UV( " << uv.X() << ", " << uv.Y() << " )" << endl;
4319       ShellPoint( params, p2 );
4320       double dist = ( p2 - p ).Modulus();
4321       if ( dist > 1e-4 )
4322         cout << "#### dist from ShellPoint " << dist
4323              << " (" << p2.X() << ", " << p2.Y() << ", " << p2.Z() << " ) " << endl;
4324     }
4325   }
4326   for ( int j = 0; j < nb; ++j )
4327     for ( int i = 0; i < nb; ++i )
4328     {
4329       int n = 1 + i + j * ( nb + 1 );
4330       cout << "mesh.AddFace([ "
4331            << n << ", " << n+1 << ", "
4332            << n+nb+2 << ", " << n+nb+1 << "]) " << endl;
4333     }
4334
4335 #endif
4336 }
4337
4338 //================================================================================
4339 /*!
4340  * \brief Constructor
4341   * \param faceID - in-block ID
4342   * \param face - geom FACE
4343   * \param baseEdge - EDGE proreply oriented in the bottom EDGE !!!
4344   * \param columnsMap - map of node columns
4345   * \param first - first normalized param
4346   * \param last - last normalized param
4347  */
4348 //================================================================================
4349
4350 StdMeshers_PrismAsBlock::TSideFace::TSideFace(SMESH_Mesh&                mesh,
4351                                               const int                  faceID,
4352                                               const Prism_3D::TQuadList& quadList,
4353                                               const TopoDS_Edge&         baseEdge,
4354                                               TParam2ColumnMap*          columnsMap,
4355                                               const double               first,
4356                                               const double               last):
4357   myID( faceID ),
4358   myParamToColumnMap( columnsMap ),
4359   myHelper( mesh )
4360 {
4361   myParams.resize( 1 );
4362   myParams[ 0 ] = make_pair( first, last );
4363   mySurface     = PSurface( new BRepAdaptor_Surface( quadList.front()->face ));
4364   myBaseEdge    = baseEdge;
4365   myIsForward   = StdMeshers_PrismAsBlock::IsForwardEdge( myHelper.GetMeshDS(),
4366                                                           *myParamToColumnMap,
4367                                                           myBaseEdge, myID );
4368   myHelper.SetSubShape( quadList.front()->face );
4369
4370   if ( quadList.size() > 1 ) // side is vertically composite
4371   {
4372     // fill myShapeID2Surf map to enable finding a right surface by any sub-shape ID
4373
4374     SMESHDS_Mesh* meshDS = myHelper.GetMeshDS();
4375
4376     TopTools_IndexedDataMapOfShapeListOfShape subToFaces;
4377     Prism_3D::TQuadList::const_iterator quad = quadList.begin();
4378     for ( ; quad != quadList.end(); ++quad )
4379     {
4380       const TopoDS_Face& face = (*quad)->face;
4381       TopExp::MapShapesAndAncestors( face, TopAbs_VERTEX, TopAbs_FACE, subToFaces );
4382       TopExp::MapShapesAndAncestors( face, TopAbs_EDGE,   TopAbs_FACE, subToFaces );
4383       myShapeID2Surf.insert( make_pair( meshDS->ShapeToIndex( face ),
4384                                         PSurface( new BRepAdaptor_Surface( face ))));
4385     }
4386     for ( int i = 1; i <= subToFaces.Extent(); ++i )
4387     {
4388       const TopoDS_Shape&     sub = subToFaces.FindKey( i );
4389       TopTools_ListOfShape& faces = subToFaces( i );
4390       int subID  = meshDS->ShapeToIndex( sub );
4391       int faceID = meshDS->ShapeToIndex( faces.First() );
4392       myShapeID2Surf.insert ( make_pair( subID, myShapeID2Surf[ faceID ]));
4393     }
4394   }
4395 }
4396
4397 //================================================================================
4398 /*!
4399  * \brief Constructor of a complex side face
4400  */
4401 //================================================================================
4402
4403 StdMeshers_PrismAsBlock::TSideFace::
4404 TSideFace(SMESH_Mesh&                             mesh,
4405           const vector< TSideFace* >&             components,
4406           const vector< pair< double, double> > & params)
4407   :myID( components[0] ? components[0]->myID : 0 ),
4408    myParamToColumnMap( 0 ),
4409    myParams( params ),
4410    myIsForward( true ),
4411    myComponents( components ),
4412    myHelper( mesh )
4413 {
4414   if ( myID == ID_Fx1z || myID == ID_F0yz )
4415   {
4416     // reverse components
4417     std::reverse( myComponents.begin(), myComponents.end() );
4418     std::reverse( myParams.begin(),     myParams.end() );
4419     for ( size_t i = 0; i < myParams.size(); ++i )
4420     {
4421       const double f = myParams[i].first;
4422       const double l = myParams[i].second;
4423       myParams[i] = make_pair( 1. - l, 1. - f );
4424     }
4425   }
4426 }
4427 //================================================================================
4428 /*!
4429  * \brief Copy constructor
4430   * \param other - other side
4431  */
4432 //================================================================================
4433
4434 StdMeshers_PrismAsBlock::TSideFace::TSideFace( const TSideFace& other ):
4435   myID               ( other.myID ),
4436   myParamToColumnMap ( other.myParamToColumnMap ),
4437   mySurface          ( other.mySurface ),
4438   myBaseEdge         ( other.myBaseEdge ),
4439   myShapeID2Surf     ( other.myShapeID2Surf ),
4440   myParams           ( other.myParams ),
4441   myIsForward        ( other.myIsForward ),
4442   myComponents       ( other.myComponents.size() ),
4443   myHelper           ( *other.myHelper.GetMesh() )
4444 {
4445   for ( size_t i = 0 ; i < myComponents.size(); ++i )
4446     myComponents[ i ] = new TSideFace( *other.myComponents[ i ]);
4447 }
4448
4449 //================================================================================
4450 /*!
4451  * \brief Deletes myComponents
4452  */
4453 //================================================================================
4454
4455 StdMeshers_PrismAsBlock::TSideFace::~TSideFace()
4456 {
4457   for ( size_t i = 0 ; i < myComponents.size(); ++i )
4458     if ( myComponents[ i ] )
4459       delete myComponents[ i ];
4460 }
4461
4462 //================================================================================
4463 /*!
4464  * \brief Return geometry of the vertical curve
4465   * \param isMax - true means curve located closer to (1,1,1) block point
4466   * \retval Adaptor3d_Curve* - curve adaptor
4467  */
4468 //================================================================================
4469
4470 Adaptor3d_Curve* StdMeshers_PrismAsBlock::TSideFace::VertiCurve(const bool isMax) const
4471 {
4472   if ( !myComponents.empty() ) {
4473     if ( isMax )
4474       return myComponents.back()->VertiCurve(isMax);
4475     else
4476       return myComponents.front()->VertiCurve(isMax);
4477   }
4478   double f = myParams[0].first, l = myParams[0].second;
4479   if ( !myIsForward ) std::swap( f, l );
4480   return new TVerticalEdgeAdaptor( myParamToColumnMap, isMax ? l : f );
4481 }
4482
4483 //================================================================================
4484 /*!
4485  * \brief Return geometry of the top or bottom curve
4486   * \param isTop -
4487   * \retval Adaptor3d_Curve* -
4488  */
4489 //================================================================================
4490
4491 Adaptor3d_Curve* StdMeshers_PrismAsBlock::TSideFace::HorizCurve(const bool isTop) const
4492 {
4493   return new THorizontalEdgeAdaptor( this, isTop );
4494 }
4495
4496 //================================================================================
4497 /*!
4498  * \brief Return pcurves
4499   * \param pcurv - array of 4 pcurves
4500   * \retval bool - is a success
4501  */
4502 //================================================================================
4503
4504 bool StdMeshers_PrismAsBlock::TSideFace::GetPCurves(Adaptor2d_Curve2d* pcurv[4]) const
4505 {
4506   int iEdge[ 4 ] = { BOTTOM_EDGE, TOP_EDGE, V0_EDGE, V1_EDGE };
4507
4508   for ( int i = 0 ; i < 4 ; ++i ) {
4509     Handle(Geom2d_Line) line;
4510     switch ( iEdge[ i ] ) {
4511     case TOP_EDGE:
4512       line = new Geom2d_Line( gp_Pnt2d( 0, 1 ), gp::DX2d() ); break;
4513     case BOTTOM_EDGE:
4514       line = new Geom2d_Line( gp::Origin2d(), gp::DX2d() ); break;
4515     case V0_EDGE:
4516       line = new Geom2d_Line( gp::Origin2d(), gp::DY2d() ); break;
4517     case V1_EDGE:
4518       line = new Geom2d_Line( gp_Pnt2d( 1, 0 ), gp::DY2d() ); break;
4519     }
4520     pcurv[ i ] = new Geom2dAdaptor_Curve( line, 0, 1 );
4521   }
4522   return true;
4523 }
4524
4525 //================================================================================
4526 /*!
4527  * \brief Returns geometry of pcurve on a horizontal face
4528   * \param isTop - is top or bottom face
4529   * \param horFace - a horizontal face
4530   * \retval Adaptor2d_Curve2d* - curve adaptor
4531  */
4532 //================================================================================
4533
4534 Adaptor2d_Curve2d*
4535 StdMeshers_PrismAsBlock::TSideFace::HorizPCurve(const bool         isTop,
4536                                                 const TopoDS_Face& horFace) const
4537 {
4538   return new TPCurveOnHorFaceAdaptor( this, isTop, horFace );
4539 }
4540
4541 //================================================================================
4542 /*!
4543  * \brief Return a component corresponding to parameter
4544   * \param U - parameter along a horizontal size
4545   * \param localU - parameter along a horizontal size of a component
4546   * \retval TSideFace* - found component
4547  */
4548 //================================================================================
4549
4550 StdMeshers_PrismAsBlock::TSideFace*
4551 StdMeshers_PrismAsBlock::TSideFace::GetComponent(const double U,double & localU) const
4552 {
4553   localU = U;
4554   if ( myComponents.empty() )
4555     return const_cast<TSideFace*>( this );
4556
4557   size_t i;
4558   for ( i = 0; i < myComponents.size(); ++i )
4559     if ( U < myParams[ i ].second )
4560       break;
4561   if ( i >= myComponents.size() )
4562     i = myComponents.size() - 1;
4563
4564   double f = myParams[ i ].first, l = myParams[ i ].second;
4565   localU = ( U - f ) / ( l - f );
4566   return myComponents[ i ];
4567 }
4568
4569 //================================================================================
4570 /*!
4571  * \brief Find node columns for a parameter
4572   * \param U - parameter along a horizontal edge
4573   * \param col1 - the 1st found column
4574   * \param col2 - the 2nd found column
4575   * \retval r - normalized position of U between the found columns
4576  */
4577 //================================================================================
4578
4579 double StdMeshers_PrismAsBlock::TSideFace::GetColumns(const double      U,
4580                                                       TParam2ColumnIt & col1,
4581                                                       TParam2ColumnIt & col2) const
4582 {
4583   double u = U, r = 0;
4584   if ( !myComponents.empty() ) {
4585     TSideFace * comp = GetComponent(U,u);
4586     return comp->GetColumns( u, col1, col2 );
4587   }
4588
4589   if ( !myIsForward )
4590     u = 1 - u;
4591   double f = myParams[0].first, l = myParams[0].second;
4592   u = f + u * ( l - f );
4593
4594   col1 = col2 = getColumn( myParamToColumnMap, u );
4595   if ( ++col2 == myParamToColumnMap->end() ) {
4596     --col2;
4597     r = 0.5;
4598   }
4599   else {
4600     double uf = col1->first;
4601     double ul = col2->first;
4602     r = ( u - uf ) / ( ul - uf );
4603   }
4604   return r;
4605 }
4606
4607 //================================================================================
4608 /*!
4609  * \brief Return all nodes at a given height together with their normalized parameters
4610  *  \param [in] Z - the height of interest
4611  *  \param [out] nodes - map of parameter to node
4612  */
4613 //================================================================================
4614
4615 void StdMeshers_PrismAsBlock::
4616 TSideFace::GetNodesAtZ(const int Z,
4617                        map<double, const SMDS_MeshNode* >& nodes ) const
4618 {
4619   if ( !myComponents.empty() )
4620   {
4621     double u0 = 0.;
4622     for ( size_t i = 0; i < myComponents.size(); ++i )
4623     {
4624       map<double, const SMDS_MeshNode* > nn;
4625       myComponents[i]->GetNodesAtZ( Z, nn );
4626       map<double, const SMDS_MeshNode* >::iterator u2n = nn.begin();
4627       if ( !nodes.empty() && nodes.rbegin()->second == u2n->second )
4628         ++u2n;
4629       const double uRange = myParams[i].second - myParams[i].first;
4630       for ( ; u2n != nn.end(); ++u2n )
4631         nodes.insert( nodes.end(), make_pair( u0 + uRange * u2n->first, u2n->second ));
4632       u0 += uRange;
4633     }
4634   }
4635   else
4636   {
4637     double f = myParams[0].first, l = myParams[0].second;
4638     if ( !myIsForward )
4639       std::swap( f, l );
4640     const double uRange = l - f;
4641     if ( Abs( uRange ) < std::numeric_limits<double>::min() )
4642       return;
4643     TParam2ColumnIt u2col = getColumn( myParamToColumnMap, myParams[0].first + 1e-3 );
4644     for ( ; u2col != myParamToColumnMap->end(); ++u2col )
4645       if ( u2col->first > myParams[0].second + 1e-9 )
4646         break;
4647       else
4648         nodes.insert( nodes.end(),
4649                       make_pair( ( u2col->first - f ) / uRange, u2col->second[ Z ] ));
4650   }
4651 }
4652
4653 //================================================================================
4654 /*!
4655  * \brief Return coordinates by normalized params
4656   * \param U - horizontal param
4657   * \param V - vertical param
4658   * \retval gp_Pnt - result point
4659  */
4660 //================================================================================
4661
4662 gp_Pnt StdMeshers_PrismAsBlock::TSideFace::Value(const Standard_Real U,
4663                                                  const Standard_Real V) const
4664 {
4665   if ( !myComponents.empty() ) {
4666     double u;
4667     TSideFace * comp = GetComponent(U,u);
4668     return comp->Value( u, V );
4669   }
4670
4671   TParam2ColumnIt u_col1, u_col2;
4672   double vR, hR = GetColumns( U, u_col1, u_col2 );
4673
4674   const SMDS_MeshNode* nn[4];
4675
4676   // BEGIN issue 0020680: Bad cell created by Radial prism in center of torus
4677   // Workaround for a wrongly located point returned by mySurface.Value() for
4678   // UV located near boundary of BSpline surface.
4679   // To bypass the problem, we take point from 3D curve of EDGE.
4680   // It solves pb of the bloc_fiss_new.py
4681   const double tol = 1e-3;
4682   if ( V < tol || V+tol >= 1. )
4683   {
4684     nn[0] = V < tol ? u_col1->second.front() : u_col1->second.back();
4685     nn[2] = V < tol ? u_col2->second.front() : u_col2->second.back();
4686     TopoDS_Edge edge;
4687     if ( V < tol )
4688     {
4689       edge = myBaseEdge;
4690     }
4691     else
4692     {
4693       TopoDS_Shape s = myHelper.GetSubShapeByNode( nn[0], myHelper.GetMeshDS() );
4694       if ( s.ShapeType() != TopAbs_EDGE )
4695         s = myHelper.GetSubShapeByNode( nn[2], myHelper.GetMeshDS() );
4696       if ( !s.IsNull() && s.ShapeType() == TopAbs_EDGE )
4697         edge = TopoDS::Edge( s );
4698     }
4699     if ( !edge.IsNull() )
4700     {
4701       double u1 = myHelper.GetNodeU( edge, nn[0], nn[2] );
4702       double u3 = myHelper.GetNodeU( edge, nn[2], nn[0] );
4703       double u = u1 * ( 1 - hR ) + u3 * hR;
4704       TopLoc_Location loc; double f,l;
4705       Handle(Geom_Curve) curve = BRep_Tool::Curve( edge,loc,f,l );
4706       return curve->Value( u ).Transformed( loc );
4707     }
4708   }
4709   // END issue 0020680: Bad cell created by Radial prism in center of torus
4710
4711   vR = getRAndNodes( & u_col1->second, V, nn[0], nn[1] );
4712   vR = getRAndNodes( & u_col2->second, V, nn[2], nn[3] );
4713
4714   if ( !myShapeID2Surf.empty() ) // side is vertically composite
4715   {
4716     // find a FACE on which the 4 nodes lie
4717     TSideFace* me = (TSideFace*) this;
4718     int notFaceID1 = 0, notFaceID2 = 0;
4719     for ( int i = 0; i < 4; ++i )
4720       if ( nn[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) // node on FACE
4721       {
4722         me->mySurface = me->myShapeID2Surf[ nn[i]->getshapeId() ];
4723         notFaceID2 = 0;
4724         break;
4725       }
4726       else if ( notFaceID1 == 0 ) // node on EDGE or VERTEX
4727       {
4728         me->mySurface  = me->myShapeID2Surf[ nn[i]->getshapeId() ];
4729         notFaceID1 = nn[i]->getshapeId();
4730       }
4731       else if ( notFaceID1 != nn[i]->getshapeId() ) // node on other EDGE or VERTEX
4732       {
4733         if ( mySurface != me->myShapeID2Surf[ nn[i]->getshapeId() ])
4734           notFaceID2 = nn[i]->getshapeId();
4735       }
4736     if ( notFaceID2 ) // no nodes of FACE and nodes are on different FACEs
4737     {
4738       SMESHDS_Mesh* meshDS = myHelper.GetMeshDS();
4739       TopoDS_Shape face = myHelper.GetCommonAncestor( meshDS->IndexToShape( notFaceID1 ),
4740                                                        meshDS->IndexToShape( notFaceID2 ),
4741                                                        *myHelper.GetMesh(),
4742                                                        TopAbs_FACE );
4743       if ( face.IsNull() )
4744         throw SALOME_Exception("StdMeshers_PrismAsBlock::TSideFace::Value() face.IsNull()");
4745       int faceID = meshDS->ShapeToIndex( face );
4746       me->mySurface = me->myShapeID2Surf[ faceID ];
4747       if ( !mySurface )
4748         throw SALOME_Exception("StdMeshers_PrismAsBlock::TSideFace::Value() !mySurface");
4749     }
4750   }
4751   ((TSideFace*) this)->myHelper.SetSubShape( mySurface->Face() );
4752
4753   gp_XY uv1 = myHelper.GetNodeUV( mySurface->Face(), nn[0], nn[2]);
4754   gp_XY uv2 = myHelper.GetNodeUV( mySurface->Face(), nn[1], nn[3]);
4755   gp_XY uv12 = uv1 * ( 1 - vR ) + uv2 * vR;
4756
4757   gp_XY uv3 = myHelper.GetNodeUV( mySurface->Face(), nn[2], nn[0]);
4758   gp_XY uv4 = myHelper.GetNodeUV( mySurface->Face(), nn[3], nn[1]);
4759   gp_XY uv34 = uv3 * ( 1 - vR ) + uv4 * vR;
4760
4761   gp_XY uv = uv12 * ( 1 - hR ) + uv34 * hR;
4762
4763   gp_Pnt p = mySurface->Value( uv.X(), uv.Y() );
4764   return p;
4765 }
4766
4767
4768 //================================================================================
4769 /*!
4770  * \brief Return boundary edge
4771   * \param edge - edge index
4772   * \retval TopoDS_Edge - found edge
4773  */
4774 //================================================================================
4775
4776 TopoDS_Edge StdMeshers_PrismAsBlock::TSideFace::GetEdge(const int iEdge) const
4777 {
4778   if ( !myComponents.empty() ) {
4779     switch ( iEdge ) {
4780     case V0_EDGE : return myComponents.front()->GetEdge( iEdge );
4781     case V1_EDGE : return myComponents.back() ->GetEdge( iEdge );
4782     default: return TopoDS_Edge();
4783     }
4784   }
4785   TopoDS_Shape edge;
4786   const SMDS_MeshNode* node = 0;
4787   SMESHDS_Mesh * meshDS = myHelper.GetMesh()->GetMeshDS();
4788   TNodeColumn* column;
4789
4790   switch ( iEdge ) {
4791   case TOP_EDGE:
4792   case BOTTOM_EDGE:
4793     column = & (( ++myParamToColumnMap->begin())->second );
4794     node = ( iEdge == TOP_EDGE ) ? column->back() : column->front();
4795     edge = myHelper.GetSubShapeByNode ( node, meshDS );
4796     if ( edge.ShapeType() == TopAbs_VERTEX ) {
4797       column = & ( myParamToColumnMap->begin()->second );
4798       node = ( iEdge == TOP_EDGE ) ? column->back() : column->front();
4799     }
4800     break;
4801   case V0_EDGE:
4802   case V1_EDGE: {
4803     bool back = ( iEdge == V1_EDGE );
4804     if ( !myIsForward ) back = !back;
4805     if ( back )
4806       column = & ( myParamToColumnMap->rbegin()->second );
4807     else
4808       column = & ( myParamToColumnMap->begin()->second );
4809     if ( column->size() > 0 )
4810       edge = myHelper.GetSubShapeByNode( (*column)[ 1 ], meshDS );
4811     if ( edge.IsNull() || edge.ShapeType() == TopAbs_VERTEX )
4812       node = column->front();
4813     break;
4814   }
4815   default:;
4816   }
4817   if ( !edge.IsNull() && edge.ShapeType() == TopAbs_EDGE )
4818     return TopoDS::Edge( edge );
4819
4820   // find edge by 2 vertices
4821   TopoDS_Shape V1 = edge;
4822   TopoDS_Shape V2 = myHelper.GetSubShapeByNode( node, meshDS );
4823   if ( !V2.IsNull() && V2.ShapeType() == TopAbs_VERTEX && !V2.IsSame( V1 ))
4824   {
4825     TopoDS_Shape ancestor = myHelper.GetCommonAncestor( V1, V2, *myHelper.GetMesh(), TopAbs_EDGE);
4826     if ( !ancestor.IsNull() )
4827       return TopoDS::Edge( ancestor );
4828   }
4829   return TopoDS_Edge();
4830 }
4831
4832 //================================================================================
4833 /*!
4834  * \brief Fill block sub-shapes
4835   * \param shapeMap - map to fill in
4836   * \retval int - nb inserted sub-shapes
4837  */
4838 //================================================================================
4839
4840 int StdMeshers_PrismAsBlock::TSideFace::InsertSubShapes(TBlockShapes& shapeMap) const
4841 {
4842   int nbInserted = 0;
4843
4844   // Insert edges
4845   vector< int > edgeIdVec;
4846   SMESH_Block::GetFaceEdgesIDs( myID, edgeIdVec );
4847
4848   for ( int i = BOTTOM_EDGE; i <=V1_EDGE ; ++i ) {
4849     TopoDS_Edge e = GetEdge( i );
4850     if ( !e.IsNull() ) {
4851       nbInserted += SMESH_Block::Insert( e, edgeIdVec[ i ], shapeMap);
4852     }
4853   }
4854
4855   // Insert corner vertices
4856
4857   TParam2ColumnIt col1, col2 ;
4858   vector< int > vertIdVec;
4859
4860   // from V0 column
4861   SMESH_Block::GetEdgeVertexIDs( edgeIdVec[ V0_EDGE ], vertIdVec);
4862   GetColumns(0, col1, col2 );
4863   const SMDS_MeshNode* node0 = col1->second.front();
4864   const SMDS_MeshNode* node1 = col1->second.back();
4865   TopoDS_Shape v0 = myHelper.GetSubShapeByNode( node0, myHelper.GetMeshDS());
4866   TopoDS_Shape v1 = myHelper.GetSubShapeByNode( node1, myHelper.GetMeshDS());
4867   if ( v0.ShapeType() == TopAbs_VERTEX ) {
4868     nbInserted += SMESH_Block::Insert( v0, vertIdVec[ 0 ], shapeMap);
4869   }
4870   if ( v1.ShapeType() == TopAbs_VERTEX ) {
4871     nbInserted += SMESH_Block::Insert( v1, vertIdVec[ 1 ], shapeMap);
4872   }
4873
4874   // from V1 column
4875   SMESH_Block::GetEdgeVertexIDs( edgeIdVec[ V1_EDGE ], vertIdVec);
4876   GetColumns(1, col1, col2 );
4877   node0 = col2->second.front();
4878   node1 = col2->second.back();
4879   v0 = myHelper.GetSubShapeByNode( node0, myHelper.GetMeshDS());
4880   v1 = myHelper.GetSubShapeByNode( node1, myHelper.GetMeshDS());
4881   if ( v0.ShapeType() == TopAbs_VERTEX ) {
4882     nbInserted += SMESH_Block::Insert( v0, vertIdVec[ 0 ], shapeMap);
4883   }
4884   if ( v1.ShapeType() == TopAbs_VERTEX ) {
4885     nbInserted += SMESH_Block::Insert( v1, vertIdVec[ 1 ], shapeMap);
4886   }
4887
4888 //   TopoDS_Vertex V0, V1, Vcom;
4889 //   TopExp::Vertices( myBaseEdge, V0, V1, true );
4890 //   if ( !myIsForward ) std::swap( V0, V1 );
4891
4892 //   // bottom vertex IDs
4893 //   SMESH_Block::GetEdgeVertexIDs( edgeIdVec[ _u0 ], vertIdVec);
4894 //   SMESH_Block::Insert( V0, vertIdVec[ 0 ], shapeMap);
4895 //   SMESH_Block::Insert( V1, vertIdVec[ 1 ], shapeMap);
4896
4897 //   TopoDS_Edge sideEdge = GetEdge( V0_EDGE );
4898 //   if ( sideEdge.IsNull() || !TopExp::CommonVertex( botEdge, sideEdge, Vcom ))
4899 //     return false;
4900
4901 //   // insert one side edge
4902 //   int edgeID;
4903 //   if ( Vcom.IsSame( V0 )) edgeID = edgeIdVec[ _v0 ];
4904 //   else                    edgeID = edgeIdVec[ _v1 ];
4905 //   SMESH_Block::Insert( sideEdge, edgeID, shapeMap);
4906
4907 //   // top vertex of the side edge
4908 //   SMESH_Block::GetEdgeVertexIDs( edgeID, vertIdVec);
4909 //   TopoDS_Vertex Vtop = TopExp::FirstVertex( sideEdge );
4910 //   if ( Vcom.IsSame( Vtop ))
4911 //     Vtop = TopExp::LastVertex( sideEdge );
4912 //   SMESH_Block::Insert( Vtop, vertIdVec[ 1 ], shapeMap);
4913
4914 //   // other side edge
4915 //   sideEdge = GetEdge( V1_EDGE );
4916 //   if ( sideEdge.IsNull() )
4917 //     return false;
4918 //   if ( edgeID = edgeIdVec[ _v1 ]) edgeID = edgeIdVec[ _v0 ];
4919 //   else                            edgeID = edgeIdVec[ _v1 ];
4920 //   SMESH_Block::Insert( sideEdge, edgeID, shapeMap);
4921
4922 //   // top edge
4923 //   TopoDS_Edge topEdge = GetEdge( TOP_EDGE );
4924 //   SMESH_Block::Insert( topEdge, edgeIdVec[ _u1 ], shapeMap);
4925
4926 //   // top vertex of the other side edge
4927 //   if ( !TopExp::CommonVertex( topEdge, sideEdge, Vcom ))
4928 //     return false;
4929 //   SMESH_Block::GetEdgeVertexIDs( edgeID, vertIdVec );
4930 //   SMESH_Block::Insert( Vcom, vertIdVec[ 1 ], shapeMap);
4931
4932   return nbInserted;
4933 }
4934
4935 //================================================================================
4936 /*!
4937  * \brief Dump ids of nodes of sides
4938  */
4939 //================================================================================
4940
4941 void StdMeshers_PrismAsBlock::TSideFace::dumpNodes(int nbNodes) const
4942 {
4943 #ifdef _DEBUG_
4944   cout << endl << "NODES OF FACE "; SMESH_Block::DumpShapeID( myID, cout ) << endl;
4945   THorizontalEdgeAdaptor* hSize0 = (THorizontalEdgeAdaptor*) HorizCurve(0);
4946   cout << "Horiz side 0: "; hSize0->dumpNodes(nbNodes); cout << endl;
4947   THorizontalEdgeAdaptor* hSize1 = (THorizontalEdgeAdaptor*) HorizCurve(1);
4948   cout << "Horiz side 1: "; hSize1->dumpNodes(nbNodes); cout << endl;
4949   TVerticalEdgeAdaptor* vSide0 = (TVerticalEdgeAdaptor*) VertiCurve(0);
4950   cout << "Verti side 0: "; vSide0->dumpNodes(nbNodes); cout << endl;
4951   TVerticalEdgeAdaptor* vSide1 = (TVerticalEdgeAdaptor*) VertiCurve(1);
4952   cout << "Verti side 1: "; vSide1->dumpNodes(nbNodes); cout << endl;
4953   delete hSize0; delete hSize1; delete vSide0; delete vSide1;
4954 #endif
4955 }
4956
4957 //================================================================================
4958 /*!
4959  * \brief Creates TVerticalEdgeAdaptor
4960   * \param columnsMap - node column map
4961   * \param parameter - normalized parameter
4962  */
4963 //================================================================================
4964
4965 StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::
4966 TVerticalEdgeAdaptor( const TParam2ColumnMap* columnsMap, const double parameter)
4967 {
4968   myNodeColumn = & getColumn( columnsMap, parameter )->second;
4969 }
4970
4971 //================================================================================
4972 /*!
4973  * \brief Return coordinates for the given normalized parameter
4974   * \param U - normalized parameter
4975   * \retval gp_Pnt - coordinates
4976  */
4977 //================================================================================
4978
4979 gp_Pnt StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::Value(const Standard_Real U) const
4980 {
4981   const SMDS_MeshNode* n1;
4982   const SMDS_MeshNode* n2;
4983   double r = getRAndNodes( myNodeColumn, U, n1, n2 );
4984   return gpXYZ(n1) * ( 1 - r ) + gpXYZ(n2) * r;
4985 }
4986
4987 //================================================================================
4988 /*!
4989  * \brief Dump ids of nodes
4990  */
4991 //================================================================================
4992
4993 void StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::dumpNodes(int nbNodes) const
4994 {
4995 #ifdef _DEBUG_
4996   for ( int i = 0; i < nbNodes && i < (int)myNodeColumn->size(); ++i )
4997     cout << (*myNodeColumn)[i]->GetID() << " ";
4998   if ( nbNodes < (int) myNodeColumn->size() )
4999     cout << myNodeColumn->back()->GetID();
5000 #endif
5001 }
5002
5003 //================================================================================
5004 /*!
5005  * \brief Return coordinates for the given normalized parameter
5006   * \param U - normalized parameter
5007   * \retval gp_Pnt - coordinates
5008  */
5009 //================================================================================
5010
5011 gp_Pnt StdMeshers_PrismAsBlock::THorizontalEdgeAdaptor::Value(const Standard_Real U) const
5012 {
5013   return mySide->TSideFace::Value( U, myV );
5014 }
5015
5016 //================================================================================
5017 /*!
5018  * \brief Dump ids of <nbNodes> first nodes and the last one
5019  */
5020 //================================================================================
5021
5022 void StdMeshers_PrismAsBlock::THorizontalEdgeAdaptor::dumpNodes(int nbNodes) const
5023 {
5024 #ifdef _DEBUG_
5025   // Not bedugged code. Last node is sometimes incorrect
5026   const TSideFace* side = mySide;
5027   double u = 0;
5028   if ( mySide->IsComplex() )
5029     side = mySide->GetComponent(0,u);
5030
5031   TParam2ColumnIt col, col2;
5032   TParam2ColumnMap* u2cols = side->GetColumns();
5033   side->GetColumns( u , col, col2 );
5034
5035   int j, i = myV ? mySide->ColumnHeight()-1 : 0;
5036
5037   const SMDS_MeshNode* n = 0;
5038   const SMDS_MeshNode* lastN
5039     = side->IsForward() ? u2cols->rbegin()->second[ i ] : u2cols->begin()->second[ i ];
5040   for ( j = 0; j < nbNodes && n != lastN; ++j )
5041   {
5042     n = col->second[ i ];
5043     cout << n->GetID() << " ";
5044     if ( side->IsForward() )
5045       ++col;
5046     else
5047       --col;
5048   }
5049
5050   // last node
5051   u = 1;
5052   if ( mySide->IsComplex() )
5053     side = mySide->GetComponent(1,u);
5054
5055   side->GetColumns( u , col, col2 );
5056   if ( n != col->second[ i ] )
5057     cout << col->second[ i ]->GetID();
5058 #endif
5059 }
5060
5061 //================================================================================
5062 /*!
5063  * \brief Constructor of TPCurveOnHorFaceAdaptor fills its map of
5064  * normalized parameter to node UV on a horizontal face
5065  *  \param [in] sideFace - lateral prism side
5066  *  \param [in] isTop - is \a horFace top or bottom of the prism
5067  *  \param [in] horFace - top or bottom face of the prism
5068  */
5069 //================================================================================
5070
5071 StdMeshers_PrismAsBlock::
5072 TPCurveOnHorFaceAdaptor::TPCurveOnHorFaceAdaptor( const TSideFace*   sideFace,
5073                                                   const bool         isTop,
5074                                                   const TopoDS_Face& horFace)
5075 {
5076   if ( sideFace && !horFace.IsNull() )
5077   {
5078     //cout << "\n\t FACE " << sideFace->FaceID() << endl;
5079     const int Z = isTop ? sideFace->ColumnHeight() - 1 : 0;
5080     map<double, const SMDS_MeshNode* > u2nodes;
5081     sideFace->GetNodesAtZ( Z, u2nodes );
5082     if ( u2nodes.empty() )
5083       return;
5084
5085     SMESH_MesherHelper helper( *sideFace->GetMesh() );
5086     helper.SetSubShape( horFace );
5087
5088     bool okUV;
5089     gp_XY uv;
5090     double f,l;
5091     Handle(Geom2d_Curve) C2d;
5092     int edgeID = -1;
5093     const double tol = 10 * helper.MaxTolerance( horFace );
5094     const SMDS_MeshNode* prevNode = u2nodes.rbegin()->second;
5095
5096     map<double, const SMDS_MeshNode* >::iterator u2n = u2nodes.begin();
5097     for ( ; u2n != u2nodes.end(); ++u2n )
5098     {
5099       const SMDS_MeshNode* n = u2n->second;
5100       okUV = false;
5101       if ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_EDGE )
5102       {
5103         if ( n->getshapeId() != edgeID )
5104         {
5105           C2d.Nullify();
5106           edgeID = n->getshapeId();
5107           TopoDS_Shape S = helper.GetSubShapeByNode( n, helper.GetMeshDS() );
5108           if ( !S.IsNull() && S.ShapeType() == TopAbs_EDGE )
5109           {
5110             C2d = BRep_Tool::CurveOnSurface( TopoDS::Edge( S ), horFace, f,l );
5111           }
5112         }
5113         if ( !C2d.IsNull() )
5114         {
5115           double u = SMDS_EdgePositionPtr( n->GetPosition() )->GetUParameter();
5116           if ( f <= u && u <= l )
5117           {
5118             uv = C2d->Value( u ).XY();
5119             okUV = helper.CheckNodeUV( horFace, n, uv, tol );
5120           }
5121         }
5122       }
5123       if ( !okUV )
5124         uv = helper.GetNodeUV( horFace, n, prevNode, &okUV );
5125
5126       myUVmap.insert( myUVmap.end(), make_pair( u2n->first, uv ));
5127       // cout << n->getshapeId() << " N " << n->GetID()
5128       //      << " \t" << uv.X() << ", " << uv.Y() << " \t" << u2n->first << endl;
5129
5130       prevNode = n;
5131     }
5132   }
5133 }
5134
5135 //================================================================================
5136 /*!
5137  * \brief Return UV on pcurve for the given normalized parameter
5138   * \param U - normalized parameter
5139   * \retval gp_Pnt - coordinates
5140  */
5141 //================================================================================
5142
5143 gp_Pnt2d StdMeshers_PrismAsBlock::TPCurveOnHorFaceAdaptor::Value(const Standard_Real U) const
5144 {
5145   map< double, gp_XY >::const_iterator i1 = myUVmap.upper_bound( U );
5146
5147   if ( i1 == myUVmap.end() )
5148     return myUVmap.empty() ? gp_XY(0,0) : myUVmap.rbegin()->second;
5149
5150   if ( i1 == myUVmap.begin() )
5151     return (*i1).second;
5152
5153   map< double, gp_XY >::const_iterator i2 = i1--;
5154
5155   double r = ( U - i1->first ) / ( i2->first - i1->first );
5156   return i1->second * ( 1 - r ) + i2->second * r;
5157 }
5158
5159 //================================================================================
5160 /*!
5161  * \brief Projects internal nodes using transformation found by boundary nodes
5162  */
5163 //================================================================================
5164
5165 bool StdMeshers_Sweeper::projectIntPoints(const vector< gp_XYZ >&    fromBndPoints,
5166                                           const vector< gp_XYZ >&    toBndPoints,
5167                                           const vector< gp_XYZ >&    fromIntPoints,
5168                                           vector< gp_XYZ >&          toIntPoints,
5169                                           const double               r,
5170                                           NSProjUtils::TrsfFinder3D& trsf,
5171                                           vector< gp_XYZ > *         bndError)
5172 {
5173   // find transformation
5174   if ( trsf.IsIdentity() && !trsf.Solve( fromBndPoints, toBndPoints ))
5175     return false;
5176
5177   // compute internal points using the found trsf
5178   for ( size_t iP = 0; iP < fromIntPoints.size(); ++iP )
5179   {
5180     toIntPoints[ iP ] = trsf.Transform( fromIntPoints[ iP ]);
5181   }
5182
5183   // compute boundary error
5184   if ( bndError )
5185   {
5186     bndError->resize( fromBndPoints.size() );
5187     gp_XYZ fromTrsf;
5188     for ( size_t iP = 0; iP < fromBndPoints.size(); ++iP )
5189     {
5190       fromTrsf = trsf.Transform( fromBndPoints[ iP ] );
5191       (*bndError)[ iP ]  = toBndPoints[ iP ] - fromTrsf;
5192     }
5193   }
5194
5195   // apply boundary error
5196   if ( bndError && toIntPoints.size() == myTopBotTriangles.size() )
5197   {
5198     for ( size_t iP = 0; iP < toIntPoints.size(); ++iP )
5199     {
5200       const TopBotTriangles& tbTrias = myTopBotTriangles[ iP ];
5201       for ( int i = 0; i < 3; ++i ) // boundary errors at 3 triangle nodes
5202       {
5203         toIntPoints[ iP ] +=
5204           ( (*bndError)[ tbTrias.myBotTriaNodes[i] ] * tbTrias.myBotBC[i] * ( 1 - r ) +
5205             (*bndError)[ tbTrias.myTopTriaNodes[i] ] * tbTrias.myTopBC[i] * ( r     ));
5206       }
5207     }
5208   }
5209
5210   return true;
5211 }
5212
5213 //================================================================================
5214 /*!
5215  * \brief Create internal nodes of the prism by computing an affine transformation
5216  *        from layer to layer
5217  */
5218 //================================================================================
5219
5220 bool StdMeshers_Sweeper::ComputeNodesByTrsf( const double tol,
5221                                              const bool   allowHighBndError)
5222 {
5223   const size_t zSize = myBndColumns[0]->size();
5224   const size_t zSrc = 0, zTgt = zSize-1;
5225   if ( zSize < 3 ) return true;
5226
5227   vector< vector< gp_XYZ > > intPntsOfLayer( zSize ); // node coordinates to compute
5228   // set coordinates of src and tgt nodes
5229   for ( size_t z = 0; z < intPntsOfLayer.size(); ++z )
5230     intPntsOfLayer[ z ].resize( myIntColumns.size() );
5231   for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
5232   {
5233     intPntsOfLayer[ zSrc ][ iP ] = intPoint( iP, zSrc );
5234     intPntsOfLayer[ zTgt ][ iP ] = intPoint( iP, zTgt );
5235   }
5236
5237   // for each internal column find boundary nodes whose error to use for correction
5238   prepareTopBotDelaunay();
5239   bool isErrorCorrectable = findDelaunayTriangles();
5240
5241   // compute coordinates of internal nodes by projecting (transforming) src and tgt
5242   // nodes towards the central layer
5243
5244   vector< NSProjUtils::TrsfFinder3D > trsfOfLayer( zSize );
5245   vector< vector< gp_XYZ > >          bndError( zSize );
5246
5247   // boundary points used to compute an affine transformation from a layer to a next one
5248   vector< gp_XYZ > fromSrcBndPnts( myBndColumns.size() ), fromTgtBndPnts( myBndColumns.size() );
5249   vector< gp_XYZ > toSrcBndPnts  ( myBndColumns.size() ), toTgtBndPnts  ( myBndColumns.size() );
5250   for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
5251   {
5252     fromSrcBndPnts[ iP ] = bndPoint( iP, zSrc );
5253     fromTgtBndPnts[ iP ] = bndPoint( iP, zTgt );
5254   }
5255
5256   size_t zS = zSrc + 1;
5257   size_t zT = zTgt - 1;
5258   for ( ; zS < zT; ++zS, --zT ) // vertical loop on layers
5259   {
5260     for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
5261     {
5262       toSrcBndPnts[ iP ] = bndPoint( iP, zS );
5263       toTgtBndPnts[ iP ] = bndPoint( iP, zT );
5264     }
5265     if (! projectIntPoints( fromSrcBndPnts, toSrcBndPnts,
5266                             intPntsOfLayer[ zS-1 ], intPntsOfLayer[ zS ],
5267                             zS / ( zSize - 1.),
5268                             trsfOfLayer   [ zS-1 ], & bndError[ zS-1 ]))
5269       return false;
5270     if (! projectIntPoints( fromTgtBndPnts, toTgtBndPnts,
5271                             intPntsOfLayer[ zT+1 ], intPntsOfLayer[ zT ],
5272                             zT / ( zSize - 1.),
5273                             trsfOfLayer   [ zT+1 ], & bndError[ zT+1 ]))
5274       return false;
5275
5276     // if ( zT == zTgt - 1 )
5277     // {
5278     //   for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
5279     //   {
5280     //     gp_XYZ fromTrsf = trsfOfLayer   [ zT+1].Transform( fromTgtBndPnts[ iP ] );
5281     //     cout << "mesh.AddNode( "
5282     //          << fromTrsf.X() << ", "
5283     //          << fromTrsf.Y() << ", "
5284     //          << fromTrsf.Z() << ") " << endl;
5285     //   }
5286     //   for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
5287     //     cout << "mesh.AddNode( "
5288     //          << intPntsOfLayer[ zT ][ iP ].X() << ", "
5289     //          << intPntsOfLayer[ zT ][ iP ].Y() << ", "
5290     //          << intPntsOfLayer[ zT ][ iP ].Z() << ") " << endl;
5291     // }
5292
5293     fromTgtBndPnts.swap( toTgtBndPnts );
5294     fromSrcBndPnts.swap( toSrcBndPnts );
5295   }
5296
5297   // Evaluate an error of boundary points
5298
5299   if ( !isErrorCorrectable && !allowHighBndError )
5300   {
5301     for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
5302     {
5303       double sumError = 0;
5304       for ( size_t z = 1; z < zS; ++z ) // loop on layers
5305         sumError += ( bndError[ z-1     ][ iP ].Modulus() +
5306                       bndError[ zSize-z ][ iP ].Modulus() );
5307
5308       if ( sumError > tol )
5309         return false;
5310     }
5311   }
5312
5313   // Compute two projections of internal points to the central layer
5314   // in order to evaluate an error of internal points
5315
5316   bool centerIntErrorIsSmall;
5317   vector< gp_XYZ > centerSrcIntPnts( myIntColumns.size() );
5318   vector< gp_XYZ > centerTgtIntPnts( myIntColumns.size() );
5319
5320   for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
5321   {
5322     toSrcBndPnts[ iP ] = bndPoint( iP, zS );
5323     toTgtBndPnts[ iP ] = bndPoint( iP, zT );
5324   }
5325   if (! projectIntPoints( fromSrcBndPnts, toSrcBndPnts,
5326                           intPntsOfLayer[ zS-1 ], centerSrcIntPnts,
5327                           zS / ( zSize - 1.),
5328                           trsfOfLayer   [ zS-1 ], & bndError[ zS-1 ]))
5329     return false;
5330   if (! projectIntPoints( fromTgtBndPnts, toTgtBndPnts,
5331                           intPntsOfLayer[ zT+1 ], centerTgtIntPnts,
5332                           zT / ( zSize - 1.),
5333                           trsfOfLayer   [ zT+1 ], & bndError[ zT+1 ]))
5334     return false;
5335
5336   // evaluate an error of internal points on the central layer
5337   centerIntErrorIsSmall = true;
5338   if ( zS == zT ) // odd zSize
5339   {
5340     for ( size_t iP = 0; ( iP < myIntColumns.size() && centerIntErrorIsSmall ); ++iP )
5341       centerIntErrorIsSmall =
5342         (centerSrcIntPnts[ iP ] - centerTgtIntPnts[ iP ]).SquareModulus() < tol*tol;
5343   }
5344   else // even zSize
5345   {
5346     for ( size_t iP = 0; ( iP < myIntColumns.size() && centerIntErrorIsSmall ); ++iP )
5347       centerIntErrorIsSmall =
5348         (intPntsOfLayer[ zS-1 ][ iP ] - centerTgtIntPnts[ iP ]).SquareModulus() < tol*tol;
5349   }
5350
5351   // compute final points on the central layer
5352   double r = zS / ( zSize - 1.);
5353   if ( zS == zT )
5354   {
5355     for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
5356     {
5357       intPntsOfLayer[ zS ][ iP ] =
5358         ( 1 - r ) * centerSrcIntPnts[ iP ] + r * centerTgtIntPnts[ iP ];
5359     }
5360   }
5361   else
5362   {
5363     for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
5364     {
5365       intPntsOfLayer[ zS ][ iP ] =
5366         r * intPntsOfLayer[ zS ][ iP ] + ( 1 - r ) * centerSrcIntPnts[ iP ];
5367       intPntsOfLayer[ zT ][ iP ] =
5368         r * intPntsOfLayer[ zT ][ iP ] + ( 1 - r ) * centerTgtIntPnts[ iP ];
5369     }
5370   }
5371
5372   if ( !centerIntErrorIsSmall )
5373   {
5374     // Compensate the central error; continue adding projection
5375     // by going from central layer to the source and target ones
5376
5377     vector< gp_XYZ >& fromSrcIntPnts = centerSrcIntPnts;
5378     vector< gp_XYZ >& fromTgtIntPnts = centerTgtIntPnts;
5379     vector< gp_XYZ >  toSrcIntPnts( myIntColumns.size() );
5380     vector< gp_XYZ >  toTgtIntPnts( myIntColumns.size() );
5381     vector< gp_XYZ >  srcBndError( myBndColumns.size() );
5382     vector< gp_XYZ >  tgtBndError( myBndColumns.size() );
5383
5384     fromTgtBndPnts.swap( toTgtBndPnts );
5385     fromSrcBndPnts.swap( toSrcBndPnts );
5386
5387     for ( ++zS, --zT; zS < zTgt; ++zS, --zT ) // vertical loop on layers
5388     {
5389       // invert transformation
5390       //if ( !trsfOfLayer[ zS+1 ].Invert() )
5391       trsfOfLayer[ zS+1 ] = NSProjUtils::TrsfFinder3D(); // to recompute
5392       //if ( !trsfOfLayer[ zT-1 ].Invert() )
5393       trsfOfLayer[ zT-1 ] = NSProjUtils::TrsfFinder3D();
5394
5395       // project internal nodes and compute bnd error
5396       for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
5397       {
5398         toSrcBndPnts[ iP ] = bndPoint( iP, zS );
5399         toTgtBndPnts[ iP ] = bndPoint( iP, zT );
5400       }
5401       projectIntPoints( fromSrcBndPnts, toSrcBndPnts,
5402                         fromSrcIntPnts, toSrcIntPnts,
5403                         zS / ( zSize - 1.),
5404                         trsfOfLayer[ zS+1 ], & srcBndError );
5405       projectIntPoints( fromTgtBndPnts, toTgtBndPnts,
5406                         fromTgtIntPnts, toTgtIntPnts,
5407                         zT / ( zSize - 1.),
5408                         trsfOfLayer[ zT-1 ], & tgtBndError );
5409
5410       // if ( zS == zTgt - 1 )
5411       // {
5412       //   cout << "mesh2 = smesh.Mesh()" << endl;
5413       //   for ( size_t iP = 0; iP < myBndColumns.size(); ++iP )
5414       //   {
5415       //     gp_XYZ fromTrsf = trsfOfLayer   [ zS+1].Transform( fromSrcBndPnts[ iP ] );
5416       //     cout << "mesh2.AddNode( "
5417       //          << fromTrsf.X() << ", "
5418       //          << fromTrsf.Y() << ", "
5419       //          << fromTrsf.Z() << ") " << endl;
5420       //   }
5421       //   for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
5422       //     cout << "mesh2.AddNode( "
5423       //          << toSrcIntPnts[ iP ].X() << ", "
5424       //          << toSrcIntPnts[ iP ].Y() << ", "
5425       //          << toSrcIntPnts[ iP ].Z() << ") " << endl;
5426       // }
5427
5428       // sum up 2 projections
5429       r = zS / ( zSize - 1.);
5430       vector< gp_XYZ >& zSIntPnts = intPntsOfLayer[ zS ];
5431       vector< gp_XYZ >& zTIntPnts = intPntsOfLayer[ zT ];
5432       for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
5433       {
5434         zSIntPnts[ iP ] = r * zSIntPnts[ iP ]  +  ( 1 - r ) * toSrcIntPnts[ iP ];
5435         zTIntPnts[ iP ] = r * zTIntPnts[ iP ]  +  ( 1 - r ) * toTgtIntPnts[ iP ];
5436       }
5437
5438       fromSrcBndPnts.swap( toSrcBndPnts );
5439       fromSrcIntPnts.swap( toSrcIntPnts );
5440       fromTgtBndPnts.swap( toTgtBndPnts );
5441       fromTgtIntPnts.swap( toTgtIntPnts );
5442     }
5443   }  // if ( !centerIntErrorIsSmall )
5444
5445
5446   //cout << "centerIntErrorIsSmall = " << centerIntErrorIsSmall<< endl;
5447
5448   // Create nodes
5449   for ( size_t iP = 0; iP < myIntColumns.size(); ++iP )
5450   {
5451     vector< const SMDS_MeshNode* > & nodeCol = *myIntColumns[ iP ];
5452     for ( size_t z = zSrc + 1; z < zTgt; ++z ) // vertical loop on layers
5453     {
5454       const gp_XYZ & xyz = intPntsOfLayer[ z ][ iP ];
5455       if ( !( nodeCol[ z ] = myHelper->AddNode( xyz.X(), xyz.Y(), xyz.Z() )))
5456         return false;
5457     }
5458   }
5459
5460   return true;
5461 }
5462
5463 //================================================================================
5464 /*!
5465  * \brief Check if all nodes of each layers have same logical Z
5466  */
5467 //================================================================================
5468
5469 bool StdMeshers_Sweeper::CheckSameZ()
5470 {
5471   myZColumns.resize( myBndColumns.size() );
5472   fillZColumn( myZColumns[0], *myBndColumns[0] );
5473
5474   bool sameZ = true;
5475   const double tol = 0.1 * 1./ myBndColumns[0]->size();
5476
5477   // check columns based on VERTEXes
5478
5479   vector< int > vertexIndex;
5480   vertexIndex.push_back( 0 );
5481   for ( size_t iC = 1; iC < myBndColumns.size() &&  sameZ; ++iC )
5482   {
5483     if ( myBndColumns[iC]->front()->GetPosition()->GetDim() > 0 )
5484       continue; // not on VERTEX
5485
5486     vertexIndex.push_back( iC );
5487     fillZColumn( myZColumns[iC], *myBndColumns[iC] );
5488
5489     for ( size_t iZ = 0; iZ < myZColumns[0].size() &&  sameZ; ++iZ )
5490       sameZ = ( Abs( myZColumns[0][iZ] - myZColumns[iC][iZ]) < tol );
5491   }
5492
5493   // check columns based on EDGEs, one per EDGE
5494
5495   for ( size_t i = 1; i < vertexIndex.size() &&  sameZ; ++i )
5496   {
5497     if ( vertexIndex[i] - vertexIndex[i-1] < 2 )
5498       continue;
5499
5500     int iC = ( vertexIndex[i] + vertexIndex[i-1] ) / 2;
5501     fillZColumn( myZColumns[iC], *myBndColumns[iC] );
5502
5503     for ( size_t iZ = 0; iZ < myZColumns[0].size() &&  sameZ; ++iZ )
5504       sameZ = ( Abs( myZColumns[0][iZ] - myZColumns[iC][iZ]) < tol );
5505   }
5506
5507   if ( sameZ )
5508   {
5509     myZColumns.resize(1);
5510   }
5511   else
5512   {
5513     for ( size_t iC = 1; iC < myBndColumns.size(); ++iC )
5514       fillZColumn( myZColumns[iC], *myBndColumns[iC] );
5515   }
5516
5517   return sameZ;
5518 }
5519
5520 //================================================================================
5521 /*!
5522  * \brief Create internal nodes of the prism all located on straight lines with
5523  *        the same distribution along the lines.
5524  */
5525 //================================================================================
5526
5527 bool StdMeshers_Sweeper::ComputeNodesOnStraightSameZ()
5528 {
5529   TZColumn& z = myZColumns[0];
5530
5531   for ( size_t i = 0; i < myIntColumns.size(); ++i )
5532   {
5533     TNodeColumn& nodes = *myIntColumns[i];
5534     SMESH_NodeXYZ n0( nodes[0] ), n1( nodes.back() );
5535
5536     for ( size_t iZ = 0; iZ < z.size(); ++iZ )
5537     {
5538       gp_XYZ p = n0 * ( 1 - z[iZ] ) + n1 * z[iZ];
5539       nodes[ iZ+1 ] = myHelper->AddNode( p.X(), p.Y(), p.Z() );
5540     }
5541   }
5542
5543   return true;
5544 }
5545
5546 //================================================================================
5547 /*!
5548  * \brief Create internal nodes of the prism all located on straight lines with
5549  *        different distributions along the lines.
5550  */
5551 //================================================================================
5552
5553 bool StdMeshers_Sweeper::ComputeNodesOnStraight()
5554 {
5555   prepareTopBotDelaunay();
5556
5557   const SMDS_MeshNode     *botNode, *topNode;
5558   const BRepMesh_Triangle *topTria;
5559   double botBC[3], topBC[3]; // barycentric coordinates
5560   int    botTriaNodes[3], topTriaNodes[3];
5561   bool   checkUV = true;
5562
5563   int nbInternalNodes = myIntColumns.size();
5564   myBotDelaunay->InitTraversal( nbInternalNodes );
5565
5566   while (( botNode = myBotDelaunay->NextNode( botBC, botTriaNodes )))
5567   {
5568     TNodeColumn* column = myIntColumns[ myNodeID2ColID( botNode->GetID() )];
5569
5570     // find a Delaunay triangle containing the topNode
5571     topNode = column->back();
5572     gp_XY topUV = myHelper->GetNodeUV( myTopFace, topNode, NULL, &checkUV );
5573     // get a starting triangle basing on that top and bot boundary nodes have same index
5574     topTria = myTopDelaunay->GetTriangleNear( botTriaNodes[0] );
5575     topTria = myTopDelaunay->FindTriangle( topUV, topTria, topBC, topTriaNodes );
5576     if ( !topTria )
5577       return false;
5578
5579     // create nodes along a line
5580     SMESH_NodeXYZ botP( botNode ), topP( topNode );
5581     for ( size_t iZ = 0; iZ < myZColumns[0].size(); ++iZ )
5582     {
5583       // use barycentric coordinates as weight of Z of boundary columns
5584       double botZ = 0, topZ = 0;
5585       for ( int i = 0; i < 3; ++i )
5586       {
5587         botZ += botBC[i] * myZColumns[ botTriaNodes[i] ][ iZ ];
5588         topZ += topBC[i] * myZColumns[ topTriaNodes[i] ][ iZ ];
5589       }
5590       double rZ = double( iZ + 1 ) / ( myZColumns[0].size() + 1 );
5591       double z = botZ * ( 1 - rZ ) + topZ * rZ;
5592       gp_XYZ p = botP * ( 1 - z  ) + topP * z;
5593       (*column)[ iZ+1 ] = myHelper->AddNode( p.X(), p.Y(), p.Z() );
5594     }
5595   }
5596
5597   return myBotDelaunay->NbVisitedNodes() == nbInternalNodes;
5598 }
5599
5600 //================================================================================
5601 /*!
5602  * \brief Compute Z of nodes of a straight column
5603  */
5604 //================================================================================
5605
5606 void StdMeshers_Sweeper::fillZColumn( TZColumn&    zColumn,
5607                                       TNodeColumn& nodes )
5608 {
5609   if ( zColumn.size() == nodes.size() - 2 )
5610     return;
5611
5612   gp_Pnt p0 = SMESH_NodeXYZ( nodes[0] );
5613   gp_Vec line( p0, SMESH_NodeXYZ( nodes.back() ));
5614   double len2 = line.SquareMagnitude();
5615
5616   zColumn.resize( nodes.size() - 2 );
5617   for ( size_t i = 0; i < zColumn.size(); ++i )
5618   {
5619     gp_Vec vec( p0, SMESH_NodeXYZ( nodes[ i+1] ));
5620     zColumn[i] = ( line * vec ) / len2; // param [0,1] on the line
5621   }
5622 }
5623
5624 //================================================================================
5625 /*!
5626  * \brief Initialize *Delaunay members
5627  */
5628 //================================================================================
5629
5630 void StdMeshers_Sweeper::prepareTopBotDelaunay()
5631 {
5632   SMESH_MesherHelper* helper[2] = { myHelper, myHelper };
5633   SMESH_MesherHelper botHelper( *myHelper->GetMesh() );
5634   SMESH_MesherHelper topHelper( *myHelper->GetMesh() );
5635   const SMDS_MeshNode* intBotNode = 0;
5636   const SMDS_MeshNode* intTopNode = 0;
5637   if ( myHelper->HasSeam() || myHelper->HasDegeneratedEdges() ) // use individual helpers
5638   {
5639     botHelper.SetSubShape( myBotFace );
5640     topHelper.SetSubShape( myTopFace );
5641     helper[0] = & botHelper;
5642     helper[1] = & topHelper;
5643     if ( !myIntColumns.empty() )
5644     {
5645       TNodeColumn& nodes = *myIntColumns[ myIntColumns.size()/2 ];
5646       intBotNode = nodes[0];
5647       intTopNode = nodes.back();
5648     }
5649   }
5650
5651   UVPtStructVec botUV( myBndColumns.size() );
5652   UVPtStructVec topUV( myBndColumns.size() );
5653   for ( size_t i = 0; i < myBndColumns.size(); ++i )
5654   {
5655     TNodeColumn& nodes = *myBndColumns[i];
5656     botUV[i].node = nodes[0];
5657     botUV[i].SetUV( helper[0]->GetNodeUV( myBotFace, nodes[0], intBotNode ));
5658     topUV[i].node = nodes.back();
5659     topUV[i].SetUV( helper[1]->GetNodeUV( myTopFace, nodes.back(), intTopNode ));
5660     botUV[i].node->setIsMarked( true );
5661   }
5662   TopoDS_Edge dummyE;
5663   SMESH_Mesh* mesh = myHelper->GetMesh();
5664   TSideVector botWires( 1, StdMeshers_FaceSide::New( botUV, myBotFace, dummyE, mesh ));
5665   TSideVector topWires( 1, StdMeshers_FaceSide::New( topUV, myTopFace, dummyE, mesh ));
5666
5667   // Delaunay mesh on the FACEs.
5668   bool checkUV = false;
5669   myBotDelaunay.reset( new NSProjUtils::Delaunay( botWires, checkUV ));
5670   myTopDelaunay.reset( new NSProjUtils::Delaunay( topWires, checkUV ));
5671
5672   if ( myHelper->GetIsQuadratic() )
5673   {
5674     // mark all medium nodes of faces on botFace to avoid their treating
5675     SMESHDS_SubMesh* smDS = myHelper->GetMeshDS()->MeshElements( myBotFace );
5676     SMDS_ElemIteratorPtr eIt = smDS->GetElements();
5677     while ( eIt->more() )
5678     {
5679       const SMDS_MeshElement* e = eIt->next();
5680       for ( int i = e->NbCornerNodes(), nb = e->NbNodes(); i < nb; ++i )
5681         e->GetNode( i )->setIsMarked( true );
5682     }
5683   }
5684
5685   // map to get a node column by a bottom node
5686   myNodeID2ColID.Clear(/*doReleaseMemory=*/false);
5687   myNodeID2ColID.ReSize( myIntColumns.size() );
5688
5689   // un-mark nodes to treat (internal bottom nodes) to be returned by myBotDelaunay
5690   for ( size_t i = 0; i < myIntColumns.size(); ++i )
5691   {
5692     const SMDS_MeshNode* botNode = myIntColumns[i]->front();
5693     botNode->setIsMarked( false );
5694     myNodeID2ColID.Bind( botNode->GetID(), i );
5695   }
5696 }
5697
5698 //================================================================================
5699 /*!
5700  * \brief For each internal node column, find Delaunay triangles including it
5701  *        and Barycentric Coordinates within the triangles. Fill in myTopBotTriangles
5702  */
5703 //================================================================================
5704
5705 bool StdMeshers_Sweeper::findDelaunayTriangles()
5706 {
5707   const SMDS_MeshNode     *botNode, *topNode;
5708   const BRepMesh_Triangle *topTria;
5709   TopBotTriangles          tbTrias;
5710   bool  checkUV = true;
5711
5712   int nbInternalNodes = myIntColumns.size();
5713   myTopBotTriangles.resize( nbInternalNodes );
5714
5715   myBotDelaunay->InitTraversal( nbInternalNodes );
5716
5717   while (( botNode = myBotDelaunay->NextNode( tbTrias.myBotBC, tbTrias.myBotTriaNodes )))
5718   {
5719     int colID = myNodeID2ColID( botNode->GetID() );
5720     TNodeColumn* column = myIntColumns[ colID ];
5721
5722     // find a Delaunay triangle containing the topNode
5723     topNode = column->back();
5724     gp_XY topUV = myHelper->GetNodeUV( myTopFace, topNode, NULL, &checkUV );
5725     // get a starting triangle basing on that top and bot boundary nodes have same index
5726     topTria = myTopDelaunay->GetTriangleNear( tbTrias.myBotTriaNodes[0] );
5727     topTria = myTopDelaunay->FindTriangle( topUV, topTria,
5728                                            tbTrias.myTopBC, tbTrias.myTopTriaNodes );
5729     if ( !topTria )
5730       tbTrias.SetTopByBottom();
5731
5732     myTopBotTriangles[ colID ] = tbTrias;
5733   }
5734
5735   if ( myBotDelaunay->NbVisitedNodes() < nbInternalNodes )
5736   {
5737     myTopBotTriangles.clear();
5738     return false;
5739   }
5740
5741   myBotDelaunay.reset();
5742   myTopDelaunay.reset();
5743   myNodeID2ColID.Clear();
5744
5745   return true;
5746 }
5747
5748 //================================================================================
5749 /*!
5750  * \brief Initialize fields
5751  */
5752 //================================================================================
5753
5754 StdMeshers_Sweeper::TopBotTriangles::TopBotTriangles()
5755 {
5756   myBotBC[0] = myBotBC[1] = myBotBC[2] = myTopBC[0] = myTopBC[1] = myTopBC[2] = 0.;
5757   myBotTriaNodes[0] = myBotTriaNodes[1] = myBotTriaNodes[2] = 0;
5758   myTopTriaNodes[0] = myTopTriaNodes[1] = myTopTriaNodes[2] = 0;
5759 }
5760
5761 //================================================================================
5762 /*!
5763  * \brief Set top data equal to bottom data
5764  */
5765 //================================================================================
5766
5767 void StdMeshers_Sweeper::TopBotTriangles::SetTopByBottom()
5768 {
5769   for ( int i = 0; i < 3; ++i )
5770   {
5771     myTopBC[i]        = myBotBC[i];
5772     myTopTriaNodes[i] = myBotTriaNodes[0];
5773   }
5774 }