src/StdMeshers/StdMeshers_Prism_3D.cxx

   1 // Copyright (C) 2007-2012  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.
  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_MesherHelper.hxx"
  37 #include "StdMeshers_FaceSide.hxx"
  38 #include "StdMeshers_ProjectionSource1D.hxx"
  39 #include "StdMeshers_ProjectionSource2D.hxx"
  40 #include "StdMeshers_ProjectionUtils.hxx"
  41 #include "StdMeshers_Projection_1D.hxx"
  42 #include "StdMeshers_Projection_1D2D.hxx"
  43 #include "StdMeshers_Quadrangle_2D.hxx"
  44
  45 #include "utilities.h"
  46
  47 #include <BRepAdaptor_CompCurve.hxx>
  48 #include <BRep_Tool.hxx>
  49 #include <Bnd_B3d.hxx>
  50 #include <Geom2dAdaptor_Curve.hxx>
  51 #include <Geom2d_Line.hxx>
  52 #include <Geom_Curve.hxx>
  53 #include <TopExp.hxx>
  54 #include <TopExp_Explorer.hxx>
  55 #include <TopTools_ListIteratorOfListOfShape.hxx>
  56 #include <TopTools_ListOfShape.hxx>
  57 #include <TopTools_MapOfShape.hxx>
  58 #include <TopTools_SequenceOfShape.hxx>
  59 #include <TopoDS.hxx>
  60 #include <gp_Ax2.hxx>
  61 #include <gp_Ax3.hxx>
  62
  63 using namespace std;
  64
  65 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
  66 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
  67 #define SHOWYXZ(msg, xyz) // {\
  68 // gp_Pnt p (xyz); \
  69 // cout << msg << " ("<< p.X() << "; " <<p.Y() << "; " <<p.Z() << ") " <<endl;\
  70 // }
  71 #ifdef _DEBUG_
  72 #define DBGOUT(msg) //cout << msg << endl;
  73 #endif
  74
  75 namespace TAssocTool = StdMeshers_ProjectionUtils;
  76
  77 typedef SMESH_Comment TCom;
  78
  79 enum { ID_BOT_FACE = SMESH_Block::ID_Fxy0,
  80        ID_TOP_FACE = SMESH_Block::ID_Fxy1,
  81        BOTTOM_EDGE = 0, TOP_EDGE, V0_EDGE, V1_EDGE, // edge IDs in face
  82        NB_WALL_FACES = 4 }; //
  83
  84 namespace {
  85
  86   //=======================================================================
  87   /*!
  88    * \brief Quadrangle algorithm
  89    */
  90   struct TQuadrangleAlgo : public StdMeshers_Quadrangle_2D
  91   {
  92     TQuadrangleAlgo(int studyId, SMESH_Gen* gen)
  93       : StdMeshers_Quadrangle_2D( gen->GetANewId(), studyId, gen)
  94     {
  95     }
  96     static StdMeshers_Quadrangle_2D* instance( SMESH_Algo*         fatherAlgo,
  97                                                SMESH_MesherHelper* helper=0)
  98     {
  99       static TQuadrangleAlgo* algo = new TQuadrangleAlgo( fatherAlgo->GetStudyId(),
 100                                                           fatherAlgo->GetGen() );
 101       if ( helper &&
 102            algo->myProxyMesh &&
 103            algo->myProxyMesh->GetMesh() != helper->GetMesh() )
 104         algo->myProxyMesh.reset( new SMESH_ProxyMesh( *helper->GetMesh() ));
 105
 106       algo->myQuadStruct.reset();
 107
 108       if ( helper )
 109         algo->_quadraticMesh = helper->GetIsQuadratic();
 110
 111       return algo;
 112     }
 113   };
 114   //=======================================================================
 115   /*!
 116    * \brief Algorithm projecting 1D mesh
 117    */
 118   struct TProjction1dAlgo : public StdMeshers_Projection_1D
 119   {
 120     StdMeshers_ProjectionSource1D myHyp;
 121
 122     TProjction1dAlgo(int studyId, SMESH_Gen* gen)
 123       : StdMeshers_Projection_1D( gen->GetANewId(), studyId, gen),
 124         myHyp( gen->GetANewId(), studyId, gen)
 125     {
 126       StdMeshers_Projection_1D::_sourceHypo = & myHyp;
 127     }
 128     static TProjction1dAlgo* instance( SMESH_Algo* fatherAlgo )
 129     {
 130       static TProjction1dAlgo* algo = new TProjction1dAlgo( fatherAlgo->GetStudyId(),
 131                                                             fatherAlgo->GetGen() );
 132       return algo;
 133     }
 134   };
 135   //=======================================================================
 136   /*!
 137    * \brief Algorithm projecting 2D mesh
 138    */
 139   struct TProjction2dAlgo : public StdMeshers_Projection_1D2D
 140   {
 141     StdMeshers_ProjectionSource2D myHyp;
 142
 143     TProjction2dAlgo(int studyId, SMESH_Gen* gen)
 144       : StdMeshers_Projection_1D2D( gen->GetANewId(), studyId, gen),
 145         myHyp( gen->GetANewId(), studyId, gen)
 146     {
 147       StdMeshers_Projection_2D::_sourceHypo = & myHyp;
 148     }
 149     static TProjction2dAlgo* instance( SMESH_Algo* fatherAlgo )
 150     {
 151       static TProjction2dAlgo* algo = new TProjction2dAlgo( fatherAlgo->GetStudyId(),
 152                                                             fatherAlgo->GetGen() );
 153       return algo;
 154     }
 155   };
 156
 157   //================================================================================
 158   /*!
 159    * \brief Make \a botE be the BOTTOM_SIDE of \a quad.
 160    *        Return false if the BOTTOM_SIDE is composite
 161    */
 162   //================================================================================
 163
 164   bool setBottomEdge( const TopoDS_Edge&   botE,
 165                       faceQuadStruct::Ptr& quad,
 166                       const TopoDS_Shape&  face)
 167   {
 168     quad->side[ QUAD_TOP_SIDE  ]->Reverse();
 169     quad->side[ QUAD_LEFT_SIDE ]->Reverse();
 170     int edgeIndex = 0;
 171     for ( size_t i = 0; i < quad->side.size(); ++i )
 172     {
 173       StdMeshers_FaceSide* quadSide = quad->side[i];
 174       for ( int iE = 0; iE < quadSide->NbEdges(); ++iE )
 175         if ( botE.IsSame( quadSide->Edge( iE )))
 176         {
 177           if ( quadSide->NbEdges() > 1 )
 178             return false;
 179           edgeIndex = i;
 180           i = quad->side.size(); // to quit from the outer loop
 181           break;
 182         }
 183     }
 184     if ( edgeIndex != QUAD_BOTTOM_SIDE )
 185       quad->shift( quad->side.size() - edgeIndex, /*keepUnitOri=*/false );
 186
 187     quad->face = TopoDS::Face( face );
 188
 189     return true;
 190   }
 191
 192   //================================================================================
 193   /*!
 194    * \brief Return iterator pointing to node column for the given parameter
 195    * \param columnsMap - node column map
 196    * \param parameter - parameter
 197    * \retval TParam2ColumnMap::iterator - result
 198    *
 199    * it returns closest left column
 200    */
 201   //================================================================================
 202
 203   TParam2ColumnIt getColumn( const TParam2ColumnMap* columnsMap,
 204                              const double            parameter )
 205   {
 206     TParam2ColumnIt u_col = columnsMap->upper_bound( parameter );
 207     if ( u_col != columnsMap->begin() )
 208       --u_col;
 209     return u_col; // return left column
 210   }
 211
 212   //================================================================================
 213   /*!
 214    * \brief Return nodes around given parameter and a ratio
 215    * \param column - node column
 216    * \param param - parameter
 217    * \param node1 - lower node
 218    * \param node2 - upper node
 219    * \retval double - ratio
 220    */
 221   //================================================================================
 222
 223   double getRAndNodes( const TNodeColumn*     column,
 224                        const double           param,
 225                        const SMDS_MeshNode* & node1,
 226                        const SMDS_MeshNode* & node2)
 227   {
 228     if ( param >= 1.0 || column->size() == 1) {
 229       node1 = node2 = column->back();
 230       return 0;
 231     }
 232
 233     int i = int( param * ( column->size() - 1 ));
 234     double u0 = double( i )/ double( column->size() - 1 );
 235     double r = ( param - u0 ) * ( column->size() - 1 );
 236
 237     node1 = (*column)[ i ];
 238     node2 = (*column)[ i + 1];
 239     return r;
 240   }
 241
 242   //================================================================================
 243   /*!
 244    * \brief Compute boundary parameters of face parts
 245     * \param nbParts - nb of parts to split columns into
 246     * \param columnsMap - node columns of the face to split
 247     * \param params - computed parameters
 248    */
 249   //================================================================================
 250
 251   void splitParams( const int               nbParts,
 252                     const TParam2ColumnMap* columnsMap,
 253                     vector< double > &      params)
 254   {
 255     params.clear();
 256     params.reserve( nbParts + 1 );
 257     TParam2ColumnIt last_par_col = --columnsMap->end();
 258     double par = columnsMap->begin()->first; // 0.
 259     double parLast = last_par_col->first;
 260     params.push_back( par );
 261     for ( int i = 0; i < nbParts - 1; ++ i )
 262     {
 263       double partSize = ( parLast - par ) / double ( nbParts - i );
 264       TParam2ColumnIt par_col = getColumn( columnsMap, par + partSize );
 265       if ( par_col->first == par ) {
 266         ++par_col;
 267         if ( par_col == last_par_col ) {
 268           while ( i < nbParts - 1 )
 269             params.push_back( par + partSize * i++ );
 270           break;
 271         }
 272       }
 273       par = par_col->first;
 274       params.push_back( par );
 275     }
 276     params.push_back( parLast ); // 1.
 277   }
 278
 279   //================================================================================
 280   /*!
 281    * \brief Return coordinate system for z-th layer of nodes
 282    */
 283   //================================================================================
 284
 285   gp_Ax2 getLayerCoordSys(const int                           z,
 286                           const vector< const TNodeColumn* >& columns,
 287                           int&                                xColumn)
 288   {
 289     // gravity center of a layer
 290     gp_XYZ O(0,0,0);
 291     int vertexCol = -1;
 292     for ( int i = 0; i < columns.size(); ++i )
 293     {
 294       O += gpXYZ( (*columns[ i ])[ z ]);
 295       if ( vertexCol < 0 &&
 296            columns[ i ]->front()->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX )
 297         vertexCol = i;
 298     }
 299     O /= columns.size();
 300
 301     // Z axis
 302     gp_Vec Z(0,0,0);
 303     int iPrev = columns.size()-1;
 304     for ( int i = 0; i < columns.size(); ++i )
 305     {
 306       gp_Vec v1( O, gpXYZ( (*columns[ iPrev ])[ z ]));
 307       gp_Vec v2( O, gpXYZ( (*columns[ i ]    )[ z ]));
 308       Z += v1 ^ v2;
 309       iPrev = i;
 310     }
 311
 312     if ( vertexCol >= 0 )
 313     {
 314       O = gpXYZ( (*columns[ vertexCol ])[ z ]);
 315     }
 316     if ( xColumn < 0 || xColumn >= columns.size() )
 317     {
 318       // select a column for X dir
 319       double maxDist = 0;
 320       for ( int i = 0; i < columns.size(); ++i )
 321       {
 322         double dist = ( O - gpXYZ((*columns[ i ])[ z ])).SquareModulus();
 323         if ( dist > maxDist )
 324         {
 325           xColumn = i;
 326           maxDist = dist;
 327         }
 328       }
 329     }
 330
 331     // X axis
 332     gp_Vec X( O, gpXYZ( (*columns[ xColumn ])[ z ]));
 333
 334     return gp_Ax2( O, Z, X);
 335   }
 336
 337   //================================================================================
 338   /*!
 339    * \brief Removes submeshes that are or can be meshed with regular grid from given list
 340    *  \retval int - nb of removed submeshes
 341    */
 342   //================================================================================
 343
 344   int removeQuasiQuads(list< SMESH_subMesh* >&   notQuadSubMesh,
 345                        SMESH_MesherHelper*       helper,
 346                        StdMeshers_Quadrangle_2D* quadAlgo)
 347   {
 348     int nbRemoved = 0;
 349     //SMESHDS_Mesh* mesh = notQuadSubMesh.front()->GetFather()->GetMeshDS();
 350     list< SMESH_subMesh* >::iterator smIt = notQuadSubMesh.begin();
 351     while ( smIt != notQuadSubMesh.end() )
 352     {
 353       SMESH_subMesh* faceSm = *smIt;
 354       SMESHDS_SubMesh* faceSmDS = faceSm->GetSubMeshDS();
 355       int nbQuads = faceSmDS ? faceSmDS->NbElements() : 0;
 356       bool toRemove;
 357       if ( nbQuads > 0 )
 358         toRemove = helper->IsStructured( faceSm );
 359       else
 360         toRemove = quadAlgo->CheckNbEdges( *helper->GetMesh(),
 361                                            faceSm->GetSubShape() );
 362       nbRemoved += toRemove;
 363       if ( toRemove )
 364         smIt = notQuadSubMesh.erase( smIt );
 365       else
 366         ++smIt;
 367     }
 368
 369     return nbRemoved;
 370   }
 371
 372 } // namespace
 373
 374 //=======================================================================
 375 //function : StdMeshers_Prism_3D
 376 //purpose  :
 377 //=======================================================================
 378
 379 StdMeshers_Prism_3D::StdMeshers_Prism_3D(int hypId, int studyId, SMESH_Gen* gen)
 380   :SMESH_3D_Algo(hypId, studyId, gen)
 381 {
 382   _name                    = "Prism_3D";
 383   _shapeType               = (1 << TopAbs_SOLID); // 1 bit per shape type
 384   _onlyUnaryInput          = false; // accept all SOLIDs at once
 385   _requireDiscreteBoundary = false; // mesh FACEs and EDGEs by myself
 386   _supportSubmeshes        = true;  // "source" FACE must be meshed by other algo
 387   _neededLowerHyps[ 1 ]    = true;  // suppress warning on hiding a global 1D algo
 388   _neededLowerHyps[ 2 ]    = true;  // suppress warning on hiding a global 2D algo
 389
 390   //myProjectTriangles       = false;
 391   mySetErrorToSM           = true;  // to pass an error to a sub-mesh of a current solid or not
 392 }
 393
 394 //================================================================================
 395 /*!
 396  * \brief Destructor
 397  */
 398 //================================================================================
 399
 400 StdMeshers_Prism_3D::~StdMeshers_Prism_3D()
 401 {}
 402
 403 //=======================================================================
 404 //function : CheckHypothesis
 405 //purpose  :
 406 //=======================================================================
 407
 408 bool StdMeshers_Prism_3D::CheckHypothesis(SMESH_Mesh&                          aMesh,
 409                                           const TopoDS_Shape&                  aShape,
 410                                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 411 {
 412   // Check shape geometry
 413 /*  PAL16229
 414   aStatus = SMESH_Hypothesis::HYP_BAD_GEOMETRY;
 415
 416   // find not quadrangle faces
 417   list< TopoDS_Shape > notQuadFaces;
 418   int nbEdge, nbWire, nbFace = 0;
 419   TopExp_Explorer exp( aShape, TopAbs_FACE );
 420   for ( ; exp.More(); exp.Next() ) {
 421     ++nbFace;
 422     const TopoDS_Shape& face = exp.Current();
 423     nbEdge = TAssocTool::Count( face, TopAbs_EDGE, 0 );
 424     nbWire = TAssocTool::Count( face, TopAbs_WIRE, 0 );
 425     if (  nbEdge!= 4 || nbWire!= 1 ) {
 426       if ( !notQuadFaces.empty() ) {
 427         if ( TAssocTool::Count( notQuadFaces.back(), TopAbs_EDGE, 0 ) != nbEdge ||
 428              TAssocTool::Count( notQuadFaces.back(), TopAbs_WIRE, 0 ) != nbWire )
 429           RETURN_BAD_RESULT("Different not quad faces");
 430       }
 431       notQuadFaces.push_back( face );
 432     }
 433   }
 434   if ( !notQuadFaces.empty() )
 435   {
 436     if ( notQuadFaces.size() != 2 )
 437       RETURN_BAD_RESULT("Bad nb not quad faces: " << notQuadFaces.size());
 438
 439     // check total nb faces
 440     nbEdge = TAssocTool::Count( notQuadFaces.back(), TopAbs_EDGE, 0 );
 441     if ( nbFace != nbEdge + 2 )
 442       RETURN_BAD_RESULT("Bad nb of faces: " << nbFace << " but must be " << nbEdge + 2);
 443   }
 444 */
 445   // no hypothesis
 446   aStatus = SMESH_Hypothesis::HYP_OK;
 447   return true;
 448 }
 449
 450 //=======================================================================
 451 //function : Compute
 452 //purpose  : Compute mesh on a COMPOUND of SOLIDs
 453 //=======================================================================
 454
 455 bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theShape)
 456 {
 457   SMESH_MesherHelper helper( theMesh );
 458   myHelper = &helper;
 459
 460   int nbSolids = helper.Count( theShape, TopAbs_SOLID, /*skipSame=*/false );
 461   if ( nbSolids < 1 )
 462     return true;
 463
 464   TopTools_IndexedDataMapOfShapeListOfShape faceToSolids;
 465   TopExp::MapShapesAndAncestors( theShape, TopAbs_FACE, TopAbs_SOLID, faceToSolids );
 466
 467   // look for meshed FACEs ("source" FACEs) that must be prism bottoms
 468   list< TopoDS_Face > meshedFaces, notQuadMeshedFaces, notQuadFaces;
 469   const bool meshHasQuads = ( theMesh.NbQuadrangles() > 0 );
 470   //StdMeshers_Quadrangle_2D* quadAlgo = TQuadrangleAlgo::instance( this );
 471   for ( int iF = 1; iF < faceToSolids.Extent(); ++iF )
 472   {
 473     const TopoDS_Face& face = TopoDS::Face( faceToSolids.FindKey( iF ));
 474     SMESH_subMesh*   faceSM = theMesh.GetSubMesh( face );
 475     if ( !faceSM->IsEmpty() )
 476     {
 477       if ( !meshHasQuads ||
 478            !helper.IsSameElemGeometry( faceSM->GetSubMeshDS(), SMDSGeom_QUADRANGLE ) ||
 479            !helper.IsStructured( faceSM )
 480            )
 481         notQuadMeshedFaces.push_front( face );
 482       else if ( myHelper->Count( face, TopAbs_EDGE, /*ignoreSame=*/false ) != 4 )
 483         meshedFaces.push_front( face );
 484       else
 485         meshedFaces.push_back( face );
 486     }
 487     // not add not quadrilateral FACE as we can't compute it
 488     // else if ( !quadAlgo->CheckNbEdges( theMesh, face ))
 489     // // not add not quadrilateral FACE as it can be a prism side
 490     // // else if ( myHelper->Count( face, TopAbs_EDGE, /*ignoreSame=*/false ) != 4 )
 491     // {
 492     //   notQuadFaces.push_back( face );
 493     // }
 494   }
 495   // notQuadFaces are of medium priority, put them before ordinary meshed faces
 496   meshedFaces.splice( meshedFaces.begin(), notQuadFaces );
 497   // notQuadMeshedFaces are of highest priority, put them before notQuadFaces
 498   meshedFaces.splice( meshedFaces.begin(), notQuadMeshedFaces );
 499
 500   Prism_3D::TPrismTopo prism;
 501
 502   if ( nbSolids == 1 )
 503   {
 504     if ( !meshedFaces.empty() )
 505       prism.myBottom = meshedFaces.front();
 506     return ( initPrism( prism, TopExp_Explorer( theShape, TopAbs_SOLID ).Current() ) &&
 507              compute( prism ));
 508   }
 509
 510   TopTools_MapOfShape meshedSolids;
 511   list< Prism_3D::TPrismTopo > meshedPrism;
 512   TopTools_ListIteratorOfListOfShape solidIt;
 513
 514   while ( meshedSolids.Extent() < nbSolids )
 515   {
 516     if ( _computeCanceled )
 517       return toSM( error( SMESH_ComputeError::New(COMPERR_CANCELED)));
 518
 519     // compute prisms having avident computed source FACE
 520     while ( !meshedFaces.empty() )
 521     {
 522       TopoDS_Face face = meshedFaces.front();
 523       meshedFaces.pop_front();
 524       TopTools_ListOfShape& solidList = faceToSolids.ChangeFromKey( face );
 525       while ( !solidList.IsEmpty() )
 526       {
 527         TopoDS_Shape solid = solidList.First();
 528         solidList.RemoveFirst();
 529         if ( meshedSolids.Add( solid ))
 530         {
 531           prism.Clear();
 532           prism.myBottom = face;
 533           if ( !initPrism( prism, solid ) ||
 534                !compute( prism ))
 535             return false;
 536
 537           meshedFaces.push_front( prism.myTop );
 538           meshedPrism.push_back( prism );
 539         }
 540       }
 541     }
 542     if ( meshedSolids.Extent() == nbSolids )
 543       break;
 544
 545     // below in the loop we try to find source FACEs somehow
 546
 547     // project mesh from source FACEs of computed prisms to
 548     // prisms sharing wall FACEs
 549     list< Prism_3D::TPrismTopo >::iterator prismIt = meshedPrism.begin();
 550     for ( ; prismIt != meshedPrism.end(); ++prismIt )
 551     {
 552       for ( size_t iW = 0; iW < prismIt->myWallQuads.size(); ++iW )
 553       {
 554         Prism_3D::TQuadList::iterator wQuad = prismIt->myWallQuads[iW].begin();
 555         for ( ; wQuad != prismIt->myWallQuads[iW].end(); ++ wQuad )
 556         {
 557           const TopoDS_Face& wFace = (*wQuad)->face;
 558           TopTools_ListOfShape& solidList = faceToSolids.ChangeFromKey( wFace );
 559           solidIt.Initialize( solidList );
 560           while ( solidIt.More() )
 561           {
 562             const TopoDS_Shape& solid = solidIt.Value();
 563             if ( meshedSolids.Contains( solid )) {
 564               solidList.Remove( solidIt );
 565               continue; // already computed prism
 566             }
 567             // find a source FACE of the SOLID: it's a FACE sharing a bottom EDGE with wFace
 568             const TopoDS_Edge& wEdge = (*wQuad)->side[ QUAD_TOP_SIDE ]->Edge(0);
 569             PShapeIteratorPtr faceIt = myHelper->GetAncestors( wEdge, *myHelper->GetMesh(),
 570                                                                TopAbs_FACE);
 571             while ( const TopoDS_Shape* f = faceIt->next() )
 572             {
 573               const TopoDS_Face& candidateF = TopoDS::Face( *f );
 574               prism.Clear();
 575               prism.myBottom  = candidateF;
 576               mySetErrorToSM = false;
 577               if ( !myHelper->IsSubShape( candidateF, prismIt->myShape3D ) &&
 578                    !myHelper->GetMesh()->GetSubMesh( candidateF )->IsMeshComputed() &&
 579                    initPrism( prism, solid ) &&
 580                    project2dMesh( prismIt->myBottom, candidateF))
 581               {
 582                 mySetErrorToSM = true;
 583                 if ( !compute( prism ))
 584                   return false;
 585                 meshedFaces.push_front( prism.myTop );
 586                 meshedFaces.push_front( prism.myBottom );
 587                 meshedPrism.push_back( prism );
 588                 meshedSolids.Add( solid );
 589               }
 590               InitComputeError();
 591             }
 592             mySetErrorToSM = true;
 593             InitComputeError();
 594             if ( meshedSolids.Contains( solid ))
 595               solidList.Remove( solidIt );
 596             else
 597               solidIt.Next();
 598           }
 599         }
 600       }
 601       if ( !meshedFaces.empty() )
 602         break; // to compute prisms with avident sources
 603     }
 604
 605     // find FACEs with local 1D hyps, which has to be computed by now,
 606     // or at least any computed FACEs
 607     for ( int iF = 1; ( meshedFaces.empty() && iF < faceToSolids.Extent() ); ++iF )
 608     {
 609       const TopoDS_Face&               face = TopoDS::Face( faceToSolids.FindKey( iF ));
 610       const TopTools_ListOfShape& solidList = faceToSolids.FindFromKey( face );
 611       if ( solidList.IsEmpty() ) continue;
 612       SMESH_subMesh*                 faceSM = theMesh.GetSubMesh( face );
 613       if ( !faceSM->IsEmpty() )
 614       {
 615         meshedFaces.push_back( face ); // lower priority
 616       }
 617       else
 618       {
 619         bool allSubMeComputed = true;
 620         SMESH_subMeshIteratorPtr smIt = faceSM->getDependsOnIterator(false,true);
 621         while ( smIt->more() && allSubMeComputed )
 622           allSubMeComputed = smIt->next()->IsMeshComputed();
 623         if ( allSubMeComputed )
 624         {
 625           faceSM->ComputeStateEngine( SMESH_subMesh::COMPUTE );
 626           if ( !faceSM->IsEmpty() )
 627             meshedFaces.push_front( face ); // higher priority
 628           else
 629             faceSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
 630         }
 631       }
 632     }
 633
 634
 635     // TODO. there are other ways to find out the source FACE:
 636     // propagation, topological similarity, ect.
 637
 638     // simply try to mesh all not meshed SOLIDs
 639     if ( meshedFaces.empty() )
 640     {
 641       for ( TopExp_Explorer solid( theShape, TopAbs_SOLID ); solid.More(); solid.Next() )
 642       {
 643         mySetErrorToSM = false;
 644         prism.Clear();
 645         if ( !meshedSolids.Contains( solid.Current() ) &&
 646              initPrism( prism, solid.Current() ))
 647         {
 648           mySetErrorToSM = true;
 649           if ( !compute( prism ))
 650             return false;
 651           meshedFaces.push_front( prism.myTop );
 652           meshedFaces.push_front( prism.myBottom );
 653           meshedPrism.push_back( prism );
 654           meshedSolids.Add( solid.Current() );
 655         }
 656         mySetErrorToSM = true;
 657       }
 658     }
 659
 660     if ( meshedFaces.empty() ) // set same error to 10 not-computed solids
 661     {
 662       SMESH_ComputeErrorPtr err = SMESH_ComputeError::New
 663         ( COMPERR_BAD_INPUT_MESH, "No meshed source face found", this );
 664
 665       const int maxNbErrors = 10; // limit nb errors not to overload the Compute dialog
 666       TopExp_Explorer solid( theShape, TopAbs_SOLID );
 667       for ( int i = 0; ( i < maxNbErrors && solid.More() ); ++i, solid.Next() )
 668         if ( !meshedSolids.Contains( solid.Current() ))
 669         {
 670           SMESH_subMesh* sm = theMesh.GetSubMesh( solid.Current() );
 671           sm->GetComputeError() = err;
 672         }
 673       return error( err );
 674     }
 675   }
 676   return true;
 677 }
 678
 679 //================================================================================
 680 /*!
 681  * \brief Find wall faces by bottom edges
 682  */
 683 //================================================================================
 684
 685 bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism,
 686                                         const int              totalNbFaces)
 687 {
 688   thePrism.myWallQuads.clear();
 689
 690   SMESH_Mesh* mesh = myHelper->GetMesh();
 691
 692   StdMeshers_Quadrangle_2D* quadAlgo = TQuadrangleAlgo::instance( this, myHelper );
 693
 694   TopTools_MapOfShape faceMap;
 695   TopTools_IndexedDataMapOfShapeListOfShape edgeToFaces;
 696   TopExp::MapShapesAndAncestors( thePrism.myShape3D,
 697                                  TopAbs_EDGE, TopAbs_FACE, edgeToFaces );
 698
 699   // ------------------------------
 700   // Get the 1st row of wall FACEs
 701   // ------------------------------
 702
 703   list< TopoDS_Edge >::iterator edge = thePrism.myBottomEdges.begin();
 704   std::list< int >::iterator     nbE = thePrism.myNbEdgesInWires.begin();
 705   int iE = 0;
 706   while ( edge != thePrism.myBottomEdges.end() )
 707   {
 708     ++iE;
 709     if ( BRep_Tool::Degenerated( *edge ))
 710     {
 711       edge = thePrism.myBottomEdges.erase( edge );
 712       --iE;
 713       --(*nbE);
 714     }
 715     else
 716     {
 717       TopTools_ListIteratorOfListOfShape faceIt( edgeToFaces.FindFromKey( *edge ));
 718       for ( ; faceIt.More(); faceIt.Next() )
 719       {
 720         const TopoDS_Face& face = TopoDS::Face( faceIt.Value() );
 721         if ( !thePrism.myBottom.IsSame( face ))
 722         {
 723           Prism_3D::TQuadList quadList( 1, quadAlgo->CheckNbEdges( *mesh, face ));
 724           if ( !quadList.back() )
 725             return toSM( error(TCom("Side face #") << shapeID( face )
 726                                << " not meshable with quadrangles"));
 727           if ( ! setBottomEdge( *edge, quadList.back(), face ))
 728             return toSM( error(TCom("Composite 'horizontal' edges are not supported")));
 729           thePrism.myWallQuads.push_back( quadList );
 730           faceMap.Add( face );
 731           break;
 732         }
 733       }
 734       ++edge;
 735     }
 736     if ( iE == *nbE )
 737     {
 738       iE = 0;
 739       ++nbE;
 740     }
 741   }
 742
 743   // -------------------------
 744   // Find the rest wall FACEs
 745   // -------------------------
 746
 747   // Compose a vector of indixes of right neighbour FACE for each wall FACE
 748   // that is not so evident in case of several WIREs in the bottom FACE
 749   thePrism.myRightQuadIndex.clear();
 750   for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
 751     thePrism.myRightQuadIndex.push_back( i+1 );
 752   list< int >::iterator nbEinW = thePrism.myNbEdgesInWires.begin();
 753   for ( int iLeft = 0; nbEinW != thePrism.myNbEdgesInWires.end(); ++nbEinW )
 754   {
 755     thePrism.myRightQuadIndex[ iLeft + *nbEinW - 1 ] = iLeft; // 1st EDGE index of a current WIRE
 756     iLeft += *nbEinW;
 757   }
 758
 759   while ( totalNbFaces - faceMap.Extent() > 2 )
 760   {
 761     // find wall FACEs adjacent to each of wallQuads by the right side EDGE
 762     int nbKnownFaces;
 763     do {
 764       nbKnownFaces = faceMap.Extent();
 765       StdMeshers_FaceSide *rightSide, *topSide; // sides of the quad
 766       for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
 767       {
 768         rightSide = thePrism.myWallQuads[i].back()->side[ QUAD_RIGHT_SIDE ];
 769         for ( int iE = 0; iE < rightSide->NbEdges(); ++iE ) // rightSide can be composite
 770         {
 771           const TopoDS_Edge & rightE = rightSide->Edge( iE );
 772           TopTools_ListIteratorOfListOfShape face( edgeToFaces.FindFromKey( rightE ));
 773           for ( ; face.More(); face.Next() )
 774             if ( faceMap.Add( face.Value() ))
 775             {
 776               // a new wall FACE encountered, store it in thePrism.myWallQuads
 777               const int iRight = thePrism.myRightQuadIndex[i];
 778               topSide = thePrism.myWallQuads[ iRight ].back()->side[ QUAD_TOP_SIDE ];
 779               const TopoDS_Edge&   newBotE = topSide->Edge(0);
 780               const TopoDS_Shape& newWallF = face.Value();
 781               thePrism.myWallQuads[ iRight ].push_back( quadAlgo->CheckNbEdges( *mesh, newWallF ));
 782               if ( !thePrism.myWallQuads[ iRight ].back() )
 783                 return toSM( error(TCom("Side face #") << shapeID( newWallF ) <<
 784                                    " not meshable with quadrangles"));
 785               if ( ! setBottomEdge( newBotE, thePrism.myWallQuads[ iRight ].back(), newWallF ))
 786                 return toSM( error(TCom("Composite 'horizontal' edges are not supported")));
 787             }
 788         }
 789       }
 790     } while ( nbKnownFaces != faceMap.Extent() );
 791
 792     // find wall FACEs adjacent to each of thePrism.myWallQuads by the top side EDGE
 793     if ( totalNbFaces - faceMap.Extent() > 2 )
 794     {
 795       for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
 796       {
 797         StdMeshers_FaceSide* topSide = thePrism.myWallQuads[i].back()->side[ QUAD_TOP_SIDE ];
 798         const TopoDS_Edge &     topE = topSide->Edge( 0 );
 799         if ( topSide->NbEdges() > 1 )
 800           return toSM( error(COMPERR_BAD_SHAPE, TCom("Side face #") <<
 801                              shapeID( thePrism.myWallQuads[i].back()->face )
 802                              << " has a composite top edge"));
 803         TopTools_ListIteratorOfListOfShape faceIt( edgeToFaces.FindFromKey( topE ));
 804         for ( ; faceIt.More(); faceIt.Next() )
 805           if ( faceMap.Add( faceIt.Value() ))
 806           {
 807             // a new wall FACE encountered, store it in wallQuads
 808             thePrism.myWallQuads[ i ].push_back( quadAlgo->CheckNbEdges( *mesh, faceIt.Value() ));
 809             if ( !thePrism.myWallQuads[ i ].back() )
 810               return toSM( error(TCom("Side face #") << shapeID( faceIt.Value() ) <<
 811                                  " not meshable with quadrangles"));
 812             if ( ! setBottomEdge( topE, thePrism.myWallQuads[ i ].back(), faceIt.Value() ))
 813               return toSM( error(TCom("Composite 'horizontal' edges are not supported")));
 814             if ( totalNbFaces - faceMap.Extent() == 2 )
 815             {
 816               i = thePrism.myWallQuads.size(); // to quit from the outer loop
 817               break;
 818             }
 819           }
 820       }
 821     }
 822   } // while ( totalNbFaces - faceMap.Extent() > 2 )
 823
 824   // ------------------
 825   // Find the top FACE
 826   // ------------------
 827
 828   if ( thePrism.myTop.IsNull() )
 829   {
 830     // now only top and bottom FACEs are not in the faceMap
 831     faceMap.Add( thePrism.myBottom );
 832     for ( TopExp_Explorer f( thePrism.myShape3D, TopAbs_FACE );f.More(); f.Next() )
 833       if ( !faceMap.Contains( f.Current() )) {
 834         thePrism.myTop = TopoDS::Face( f.Current() );
 835         break;
 836       }
 837     if ( thePrism.myTop.IsNull() )
 838       return toSM( error("Top face not found"));
 839   }
 840
 841   // Check that the top FACE shares all the top EDGEs
 842   for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i )
 843   {
 844     StdMeshers_FaceSide* topSide = thePrism.myWallQuads[i].back()->side[ QUAD_TOP_SIDE ];
 845     const TopoDS_Edge &     topE = topSide->Edge( 0 );
 846     if ( !myHelper->IsSubShape( topE, thePrism.myTop ))
 847       return toSM( error( TCom("Wrong source face (#") << shapeID( thePrism.myBottom )));
 848   }
 849
 850   return true;
 851 }
 852
 853 //=======================================================================
 854 //function : compute
 855 //purpose  : Compute mesh on a SOLID
 856 //=======================================================================
 857
 858 bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism)
 859 {
 860   myHelper->IsQuadraticSubMesh( thePrism.myShape3D );
 861   if ( _computeCanceled )
 862     return toSM( error( SMESH_ComputeError::New(COMPERR_CANCELED)));
 863
 864   // Make all side FACEs of thePrism meshed with quads
 865   if ( !computeWalls( thePrism ))
 866     return false;
 867
 868   // Analyse mesh and geometry to find block sub-shapes and submeshes
 869   if ( !myBlock.Init( myHelper, thePrism ))
 870     return toSM( error( myBlock.GetError()));
 871
 872   SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
 873
 874   int volumeID = meshDS->ShapeToIndex( thePrism.myShape3D );
 875
 876
 877   // To compute coordinates of a node inside a block, it is necessary to know
 878   // 1. normalized parameters of the node by which
 879   // 2. coordinates of node projections on all block sub-shapes are computed
 880
 881   // So we fill projections on vertices at once as they are same for all nodes
 882   myShapeXYZ.resize( myBlock.NbSubShapes() );
 883   for ( int iV = SMESH_Block::ID_FirstV; iV < SMESH_Block::ID_FirstE; ++iV ) {
 884     myBlock.VertexPoint( iV, myShapeXYZ[ iV ]);
 885     SHOWYXZ("V point " <<iV << " ", myShapeXYZ[ iV ]);
 886   }
 887
 888   // Projections on the top and bottom faces are taken from nodes existing
 889   // on these faces; find correspondence between bottom and top nodes
 890   myBotToColumnMap.clear();
 891   if ( !assocOrProjBottom2Top() ) // it also fills myBotToColumnMap
 892     return false;
 893
 894
 895   // Create nodes inside the block
 896
 897   // try to use transformation (issue 0020680)
 898   vector<gp_Trsf> trsf;
 899   if ( myBlock.GetLayersTransformation( trsf, thePrism ))
 900   {
 901     // loop on nodes inside the bottom face
 902     TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin();
 903     for ( ; bot_column != myBotToColumnMap.end(); ++bot_column )
 904     {
 905       const Prism_3D::TNode& tBotNode = bot_column->first; // bottom TNode
 906       if ( tBotNode.GetPositionType() != SMDS_TOP_FACE )
 907         continue; // node is not inside face
 908
 909       // column nodes; middle part of the column are zero pointers
 910       TNodeColumn& column = bot_column->second;
 911       TNodeColumn::iterator columnNodes = column.begin();
 912       for ( int z = 0; columnNodes != column.end(); ++columnNodes, ++z)
 913       {
 914         const SMDS_MeshNode* & node = *columnNodes;
 915         if ( node ) continue; // skip bottom or top node
 916
 917         gp_XYZ coords = tBotNode.GetCoords();
 918         trsf[z-1].Transforms( coords );
 919         node = meshDS->AddNode( coords.X(), coords.Y(), coords.Z() );
 920         meshDS->SetNodeInVolume( node, volumeID );
 921       }
 922     } // loop on bottom nodes
 923   }
 924   else // use block approach
 925   {
 926     // loop on nodes inside the bottom face
 927     Prism_3D::TNode prevBNode;
 928     TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin();
 929     for ( ; bot_column != myBotToColumnMap.end(); ++bot_column )
 930     {
 931       const Prism_3D::TNode& tBotNode = bot_column->first; // bottom TNode
 932       if ( tBotNode.GetPositionType() != SMDS_TOP_FACE )
 933         continue; // node is not inside face
 934
 935       // column nodes; middle part of the column are zero pointers
 936       TNodeColumn& column = bot_column->second;
 937
 938       // compute bottom node parameters
 939       gp_XYZ paramHint(-1,-1,-1);
 940       if ( prevBNode.IsNeighbor( tBotNode ))
 941         paramHint = prevBNode.GetParams();
 942       if ( !myBlock.ComputeParameters( tBotNode.GetCoords(), tBotNode.ChangeParams(),
 943                                        ID_BOT_FACE, paramHint ))
 944         return toSM( error(TCom("Can't compute normalized parameters for node ")
 945                            << tBotNode.myNode->GetID() << " on the face #"
 946                            << myBlock.SubMesh( ID_BOT_FACE )->GetId() ));
 947       prevBNode = tBotNode;
 948
 949       myShapeXYZ[ ID_BOT_FACE ] = tBotNode.GetCoords();
 950       gp_XYZ botParams          = tBotNode.GetParams();
 951
 952       // compute top node parameters
 953       myShapeXYZ[ ID_TOP_FACE ] = gpXYZ( column.back() );
 954       gp_XYZ topParams = botParams;
 955       topParams.SetZ( 1 );
 956       if ( column.size() > 2 ) {
 957         gp_Pnt topCoords = myShapeXYZ[ ID_TOP_FACE ];
 958         if ( !myBlock.ComputeParameters( topCoords, topParams, ID_TOP_FACE, topParams ))
 959           return toSM( error(TCom("Can't compute normalized parameters ")
 960                              << "for node " << column.back()->GetID()
 961                              << " on the face #"<< column.back()->getshapeId() ));
 962       }
 963
 964       // vertical loop
 965       TNodeColumn::iterator columnNodes = column.begin();
 966       for ( int z = 0; columnNodes != column.end(); ++columnNodes, ++z)
 967       {
 968         const SMDS_MeshNode* & node = *columnNodes;
 969         if ( node ) continue; // skip bottom or top node
 970
 971         // params of a node to create
 972         double rz = (double) z / (double) ( column.size() - 1 );
 973         gp_XYZ params = botParams * ( 1 - rz ) + topParams * rz;
 974
 975         // set coords on all faces and nodes
 976         const int nbSideFaces = 4;
 977         int sideFaceIDs[nbSideFaces] = { SMESH_Block::ID_Fx0z,
 978                                          SMESH_Block::ID_Fx1z,
 979                                          SMESH_Block::ID_F0yz,
 980                                          SMESH_Block::ID_F1yz };
 981         for ( int iF = 0; iF < nbSideFaces; ++iF )
 982           if ( !setFaceAndEdgesXYZ( sideFaceIDs[ iF ], params, z ))
 983             return false;
 984
 985         // compute coords for a new node
 986         gp_XYZ coords;
 987         if ( !SMESH_Block::ShellPoint( params, myShapeXYZ, coords ))
 988           return toSM( error("Can't compute coordinates by normalized parameters"));
 989
 990         SHOWYXZ("TOPFacePoint ",myShapeXYZ[ ID_TOP_FACE]);
 991         SHOWYXZ("BOT Node "<< tBotNode.myNode->GetID(),gpXYZ(tBotNode.myNode));
 992         SHOWYXZ("ShellPoint ",coords);
 993
 994         // create a node
 995         node = meshDS->AddNode( coords.X(), coords.Y(), coords.Z() );
 996         meshDS->SetNodeInVolume( node, volumeID );
 997       }
 998     } // loop on bottom nodes
 999   }
1000
1001   // Create volumes
1002
1003   SMESHDS_SubMesh* smDS = myBlock.SubMeshDS( ID_BOT_FACE );
1004   if ( !smDS ) return toSM( error(COMPERR_BAD_INPUT_MESH, "Null submesh"));
1005
1006   // loop on bottom mesh faces
1007   SMDS_ElemIteratorPtr faceIt = smDS->GetElements();
1008   while ( faceIt->more() )
1009   {
1010     const SMDS_MeshElement* face = faceIt->next();
1011     if ( !face || face->GetType() != SMDSAbs_Face )
1012       continue;
1013
1014     // find node columns for each node
1015     int nbNodes = face->NbCornerNodes();
1016     vector< const TNodeColumn* > columns( nbNodes );
1017     for ( int i = 0; i < nbNodes; ++i )
1018     {
1019       const SMDS_MeshNode* n = face->GetNode( i );
1020       if ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) {
1021         TNode2ColumnMap::iterator bot_column = myBotToColumnMap.find( n );
1022         if ( bot_column == myBotToColumnMap.end() )
1023           return toSM( error(TCom("No nodes found above node ") << n->GetID() ));
1024         columns[ i ] = & bot_column->second;
1025       }
1026       else {
1027         columns[ i ] = myBlock.GetNodeColumn( n );
1028         if ( !columns[ i ] )
1029           return toSM( error(TCom("No side nodes found above node ") << n->GetID() ));
1030       }
1031     }
1032     // create prisms
1033     AddPrisms( columns, myHelper );
1034
1035   } // loop on bottom mesh faces
1036
1037   // clear data
1038   myBotToColumnMap.clear();
1039   myBlock.Clear();
1040
1041   return true;
1042 }
1043
1044 //=======================================================================
1045 //function : computeWalls
1046 //purpose  : Compute 2D mesh on walls FACEs of a prism
1047 //=======================================================================
1048
1049 bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism)
1050 {
1051   SMESH_Mesh*     mesh = myHelper->GetMesh();
1052   SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
1053   DBGOUT( endl << "COMPUTE Prism " << meshDS->ShapeToIndex( thePrism.myShape3D ));
1054
1055   TProjction1dAlgo* projector1D = TProjction1dAlgo::instance( this );
1056   StdMeshers_Quadrangle_2D* quadAlgo = TQuadrangleAlgo::instance( this, myHelper );
1057
1058   SMESH_HypoFilter hyp1dFilter( SMESH_HypoFilter::IsAlgo(),/*not=*/true);
1059   hyp1dFilter.And( SMESH_HypoFilter::HasDim( 1 ));
1060   hyp1dFilter.And( SMESH_HypoFilter::IsMoreLocalThan( thePrism.myShape3D, *mesh ));
1061
1062   // Discretize equally 'vertical' EDGEs
1063   // -----------------------------------
1064   // find source FACE sides for projection: either already computed ones or
1065   // the 'most composite' ones
1066   multimap< int, int > wgt2quad;
1067   for ( size_t iW = 0; iW != thePrism.myWallQuads.size(); ++iW )
1068   {
1069     Prism_3D::TQuadList::const_iterator quad = thePrism.myWallQuads[iW].begin();
1070     int wgt = 0; // "weight"
1071     for ( ; quad != thePrism.myWallQuads[iW].end(); ++quad )
1072     {
1073       StdMeshers_FaceSide* lftSide = (*quad)->side[ QUAD_LEFT_SIDE ];
1074       for ( int i = 0; i < lftSide->NbEdges(); ++i )
1075       {
1076         ++wgt;
1077         const TopoDS_Edge& E = lftSide->Edge(i);
1078         if ( mesh->GetSubMesh( E )->IsMeshComputed() )
1079           wgt += 10;
1080         else if ( mesh->GetHypothesis( E, hyp1dFilter, true )) // local hypothesis!
1081           wgt += 100;
1082       }
1083     }
1084     wgt2quad.insert( make_pair( wgt, iW ));
1085
1086     // in quadratic mesh, pass ignoreMediumNodes to quad sides
1087     if ( myHelper->GetIsQuadratic() )
1088     {
1089       quad = thePrism.myWallQuads[iW].begin();
1090       for ( ; quad != thePrism.myWallQuads[iW].end(); ++quad )
1091         for ( int i = 0; i < NB_QUAD_SIDES; ++i )
1092           (*quad)->side[ i ]->SetIgnoreMediumNodes( true );
1093     }
1094   }
1095
1096   // Project 'vertical' EDGEs, from left to right
1097   multimap< int, int >::reverse_iterator w2q = wgt2quad.rbegin();
1098   for ( ; w2q != wgt2quad.rend(); ++w2q )
1099   {
1100     const int iW = w2q->second;
1101     const Prism_3D::TQuadList&         quads = thePrism.myWallQuads[ iW ];
1102     Prism_3D::TQuadList::const_iterator quad = quads.begin();
1103     for ( ; quad != quads.end(); ++quad )
1104     {
1105       StdMeshers_FaceSide* rgtSide = (*quad)->side[ QUAD_RIGHT_SIDE ]; // tgt
1106       StdMeshers_FaceSide* lftSide = (*quad)->side[ QUAD_LEFT_SIDE ];  // src
1107       bool swapLeftRight = ( lftSide->NbSegments( /*update=*/true ) == 0 &&
1108                              rgtSide->NbSegments( /*update=*/true )  > 0 );
1109       if ( swapLeftRight )
1110         std::swap( lftSide, rgtSide );
1111
1112       // assure that all the source (left) EDGEs are meshed
1113       int nbSrcSegments = 0;
1114       for ( int i = 0; i < lftSide->NbEdges(); ++i )
1115       {
1116         const TopoDS_Edge& srcE = lftSide->Edge(i);
1117         SMESH_subMesh*    srcSM = mesh->GetSubMesh( srcE );
1118         if ( !srcSM->IsMeshComputed() ) {
1119           DBGOUT( "COMPUTE V edge " << srcSM->GetId() );
1120           srcSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE );
1121           srcSM->ComputeStateEngine       ( SMESH_subMesh::COMPUTE );
1122           if ( !srcSM->IsMeshComputed() )
1123             return false;
1124         }
1125         nbSrcSegments += srcSM->GetSubMeshDS()->NbElements();
1126       }
1127       // check target EDGEs
1128       int nbTgtMeshed = 0, nbTgtSegments = 0;
1129       vector< bool > isTgtEdgeComputed( rgtSide->NbEdges() );
1130       for ( int i = 0; i < rgtSide->NbEdges(); ++i )
1131       {
1132         const TopoDS_Edge& tgtE = rgtSide->Edge(i);
1133         SMESH_subMesh*    tgtSM = mesh->GetSubMesh( tgtE );
1134         if (( isTgtEdgeComputed[ i ] = tgtSM->IsMeshComputed() )) {
1135           ++nbTgtMeshed;
1136           nbTgtSegments += tgtSM->GetSubMeshDS()->NbElements();
1137         }
1138       }
1139       if ( rgtSide->NbEdges() == nbTgtMeshed ) // all tgt EDGEs meshed
1140       {
1141         if ( nbTgtSegments != nbSrcSegments )
1142         {
1143           for ( int i = 0; i < lftSide->NbEdges(); ++i )
1144             addBadInputElements( meshDS->MeshElements( lftSide->Edge( i )));
1145           for ( int i = 0; i < rgtSide->NbEdges(); ++i )
1146             addBadInputElements( meshDS->MeshElements( rgtSide->Edge( i )));
1147           return toSM( error( TCom("Different nb of segment on logically vertical edges #")
1148                               << shapeID( lftSide->Edge(0) ) << " and #"
1149                               << shapeID( rgtSide->Edge(0) ) << ": "
1150                               << nbSrcSegments << " != " << nbTgtSegments ));
1151         }
1152         continue;
1153       }
1154       // Compute 'vertical projection'
1155       if ( nbTgtMeshed == 0 )
1156       {
1157         // compute nodes on target VERTEXes
1158         const UVPtStructVec&  srcNodeStr = lftSide->GetUVPtStruct();
1159         if ( srcNodeStr.size() == 0 )
1160           return toSM( error( TCom("Invalid node positions on edge #") <<
1161                               shapeID( lftSide->Edge(0) )));
1162         vector< SMDS_MeshNode* > newNodes( srcNodeStr.size() );
1163         for ( int is2ndV = 0; is2ndV < 2; ++is2ndV )
1164         {
1165           const TopoDS_Edge& E = rgtSide->Edge( is2ndV ? rgtSide->NbEdges()-1 : 0 );
1166           TopoDS_Vertex      v = myHelper->IthVertex( is2ndV, E );
1167           mesh->GetSubMesh( v )->ComputeStateEngine( SMESH_subMesh::COMPUTE );
1168           const SMDS_MeshNode* n = SMESH_Algo::VertexNode( v, meshDS );
1169           newNodes[ is2ndV ? 0 : newNodes.size()-1 ] = (SMDS_MeshNode*) n;
1170         }
1171
1172         // compute nodes on target EDGEs
1173         DBGOUT( "COMPUTE V edge (proj) " << shapeID( lftSide->Edge(0)));
1174         rgtSide->Reverse(); // direct it same as the lftSide
1175         myHelper->SetElementsOnShape( false ); // myHelper holds the prism shape
1176         TopoDS_Edge tgtEdge;
1177         for ( size_t iN = 1; iN < srcNodeStr.size()-1; ++iN ) // add nodes
1178         {
1179           gp_Pnt       p = rgtSide->Value3d  ( srcNodeStr[ iN ].normParam );
1180           double       u = rgtSide->Parameter( srcNodeStr[ iN ].normParam, tgtEdge );
1181           newNodes[ iN ] = meshDS->AddNode( p.X(), p.Y(), p.Z() );
1182           meshDS->SetNodeOnEdge( newNodes[ iN ], tgtEdge, u );
1183         }
1184         for ( size_t iN = 1; iN < srcNodeStr.size(); ++iN ) // add segments
1185         {
1186           // find an EDGE to set a new segment
1187           std::pair<int, TopAbs_ShapeEnum> id2type =
1188             myHelper->GetMediumPos( newNodes[ iN-1 ], newNodes[ iN ] );
1189           if ( id2type.second != TopAbs_EDGE )
1190           {
1191             // new nodes are on different EDGEs; put one of them on VERTEX
1192             const int      edgeIndex = rgtSide->EdgeIndex( srcNodeStr[ iN-1 ].normParam );
1193             const double vertexParam = rgtSide->LastParameter( edgeIndex );
1194             const gp_Pnt           p = BRep_Tool::Pnt( rgtSide->LastVertex( edgeIndex ));
1195             const int         isPrev = ( Abs( srcNodeStr[ iN-1 ].normParam - vertexParam ) <
1196                                          Abs( srcNodeStr[ iN   ].normParam - vertexParam ));
1197             meshDS->UnSetNodeOnShape( newNodes[ iN-isPrev ] );
1198             meshDS->SetNodeOnVertex ( newNodes[ iN-isPrev ], rgtSide->LastVertex( edgeIndex ));
1199             meshDS->MoveNode        ( newNodes[ iN-isPrev ], p.X(), p.Y(), p.Z() );
1200             id2type.first = newNodes[ iN-(1-isPrev) ]->getshapeId();
1201           }
1202           SMDS_MeshElement* newEdge = myHelper->AddEdge( newNodes[ iN-1 ], newNodes[ iN ] );
1203           meshDS->SetMeshElementOnShape( newEdge, id2type.first );
1204         }
1205         myHelper->SetElementsOnShape( true );
1206         for ( int i = 0; i < rgtSide->NbEdges(); ++i ) // update state of sub-meshes
1207         {
1208           const TopoDS_Edge& E = rgtSide->Edge( i );
1209           SMESH_subMesh* tgtSM = mesh->GetSubMesh( E );
1210           tgtSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1211         }
1212
1213         // to continue projection from the just computed side as a source
1214         if ( !swapLeftRight && rgtSide->NbEdges() > 1 && w2q->second == iW )
1215         {
1216           std::pair<int,int> wgt2quadKeyVal( w2q->first + 1, thePrism.myRightQuadIndex[ iW ]);
1217           wgt2quad.insert( wgt2quadKeyVal ); // it will be skipped by ++w2q
1218           wgt2quad.insert( wgt2quadKeyVal );
1219           w2q = wgt2quad.rbegin();
1220         }
1221       }
1222       else
1223       {
1224         // HOPE assigned hypotheses are OK, so that equal nb of segments will be generated
1225         //return toSM( error("Partial projection not implemented"));
1226       }
1227     } // loop on quads of a composite wall side
1228   } // loop on the ordered wall sides
1229
1230
1231
1232   for ( size_t iW = 0; iW != thePrism.myWallQuads.size(); ++iW )
1233   {
1234     Prism_3D::TQuadList::const_iterator quad = thePrism.myWallQuads[iW].begin();
1235     for ( ; quad != thePrism.myWallQuads[iW].end(); ++quad )
1236     {
1237       // Top EDGEs must be projections from the bottom ones
1238       // to compute stuctured quad mesh on wall FACEs
1239       // ---------------------------------------------------
1240       {
1241         const TopoDS_Edge& botE = (*quad)->side[ QUAD_BOTTOM_SIDE ]->Edge(0);
1242         const TopoDS_Edge& topE = (*quad)->side[ QUAD_TOP_SIDE    ]->Edge(0);
1243         SMESH_subMesh*    botSM = mesh->GetSubMesh( botE );
1244         SMESH_subMesh*    topSM = mesh->GetSubMesh( topE );
1245         SMESH_subMesh*    srcSM = botSM;
1246         SMESH_subMesh*    tgtSM = topSM;
1247         if ( !srcSM->IsMeshComputed() && topSM->IsMeshComputed() )
1248           std::swap( srcSM, tgtSM );
1249
1250         if ( !srcSM->IsMeshComputed() )
1251         {
1252           DBGOUT( "COMPUTE H edge " << srcSM->GetId());
1253           srcSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE ); // nodes on VERTEXes
1254           srcSM->ComputeStateEngine( SMESH_subMesh::COMPUTE );        // segments on the EDGE
1255         }
1256         srcSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1257
1258         if ( !tgtSM->IsMeshComputed() )
1259         {
1260           // compute nodes on VERTEXes
1261           tgtSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE );
1262           // project segments
1263           DBGOUT( "COMPUTE H edge (proj) " << tgtSM->GetId());
1264           projector1D->myHyp.SetSourceEdge( TopoDS::Edge( srcSM->GetSubShape() ));
1265           projector1D->InitComputeError();
1266           bool ok = projector1D->Compute( *mesh, tgtSM->GetSubShape() );
1267           if ( !ok )
1268           {
1269             SMESH_ComputeErrorPtr err = projector1D->GetComputeError();
1270             if ( err->IsOK() ) err->myName = COMPERR_ALGO_FAILED;
1271             tgtSM->GetComputeError() = err;
1272             return false;
1273           }
1274         }
1275         tgtSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1276       }
1277
1278       // Compute quad mesh on wall FACEs
1279       // -------------------------------
1280       const TopoDS_Face& face = (*quad)->face;
1281       SMESH_subMesh* fSM = mesh->GetSubMesh( face );
1282       if ( ! fSM->IsMeshComputed() )
1283       {
1284         // make all EDGES meshed
1285         fSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE );
1286         if ( !fSM->SubMeshesComputed() )
1287           return toSM( error( COMPERR_BAD_INPUT_MESH,
1288                               "Not all edges have valid algorithm and hypothesis"));
1289         // mesh the <face>
1290         quadAlgo->InitComputeError();
1291         DBGOUT( "COMPUTE Quad face " << fSM->GetId());
1292         bool ok = quadAlgo->Compute( *mesh, face );
1293         fSM->GetComputeError() = quadAlgo->GetComputeError();
1294         if ( !ok )
1295           return false;
1296         fSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1297       }
1298       if ( myHelper->GetIsQuadratic() )
1299       {
1300         // fill myHelper with medium nodes built by quadAlgo
1301         SMDS_ElemIteratorPtr fIt = fSM->GetSubMeshDS()->GetElements();
1302         while ( fIt->more() )
1303           myHelper->AddTLinks( dynamic_cast<const SMDS_MeshFace*>( fIt->next() ));
1304       }
1305     }
1306   }
1307
1308   return true;
1309 }
1310
1311 //=======================================================================
1312 //function : Evaluate
1313 //purpose  :
1314 //=======================================================================
1315
1316 bool StdMeshers_Prism_3D::Evaluate(SMESH_Mesh&         theMesh,
1317                                    const TopoDS_Shape& theShape,
1318                                    MapShapeNbElems&    aResMap)
1319 {
1320   if ( theShape.ShapeType() == TopAbs_COMPOUND )
1321   {
1322     bool ok = true;
1323     for ( TopoDS_Iterator it( theShape ); it.More(); it.Next() )
1324       ok &= Evaluate( theMesh, it.Value(), aResMap );
1325     return ok;
1326   }
1327   SMESH_MesherHelper helper( theMesh );
1328   myHelper = &helper;
1329   myHelper->SetSubShape( theShape );
1330
1331   // find face contains only triangles
1332   vector < SMESH_subMesh * >meshFaces;
1333   TopTools_SequenceOfShape aFaces;
1334   int NumBase = 0, i = 0, NbQFs = 0;
1335   for (TopExp_Explorer exp(theShape, TopAbs_FACE); exp.More(); exp.Next()) {
1336     i++;
1337     aFaces.Append(exp.Current());
1338     SMESH_subMesh *aSubMesh = theMesh.GetSubMesh(exp.Current());
1339     meshFaces.push_back(aSubMesh);
1340     MapShapeNbElemsItr anIt = aResMap.find(meshFaces[i-1]);
1341     if( anIt==aResMap.end() )
1342       return toSM( error( "Submesh can not be evaluated"));
1343
1344     std::vector<int> aVec = (*anIt).second;
1345     int nbtri = Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
1346     int nbqua = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
1347     if( nbtri==0 && nbqua>0 ) {
1348       NbQFs++;
1349     }
1350     if( nbtri>0 ) {
1351       NumBase = i;
1352     }
1353   }
1354
1355   if(NbQFs<4) {
1356     std::vector<int> aResVec(SMDSEntity_Last);
1357     for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1358     SMESH_subMesh * sm = theMesh.GetSubMesh(theShape);
1359     aResMap.insert(std::make_pair(sm,aResVec));
1360     return toSM( error( "Submesh can not be evaluated" ));
1361   }
1362
1363   if(NumBase==0) NumBase = 1; // only quads => set 1 faces as base
1364
1365   // find number of 1d elems for base face
1366   int nb1d = 0;
1367   TopTools_MapOfShape Edges1;
1368   for (TopExp_Explorer exp(aFaces.Value(NumBase), TopAbs_EDGE); exp.More(); exp.Next()) {
1369     Edges1.Add(exp.Current());
1370     SMESH_subMesh *sm = theMesh.GetSubMesh(exp.Current());
1371     if( sm ) {
1372       MapShapeNbElemsItr anIt = aResMap.find(sm);
1373       if( anIt == aResMap.end() ) continue;
1374       std::vector<int> aVec = (*anIt).second;
1375       nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
1376     }
1377   }
1378   // find face opposite to base face
1379   int OppNum = 0;
1380   for(i=1; i<=6; i++) {
1381     if(i==NumBase) continue;
1382     bool IsOpposite = true;
1383     for(TopExp_Explorer exp(aFaces.Value(i), TopAbs_EDGE); exp.More(); exp.Next()) {
1384       if( Edges1.Contains(exp.Current()) ) {
1385         IsOpposite = false;
1386         break;
1387       }
1388     }
1389     if(IsOpposite) {
1390       OppNum = i;
1391       break;
1392     }
1393   }
1394   // find number of 2d elems on side faces
1395   int nb2d = 0;
1396   for(i=1; i<=6; i++) {
1397     if( i==OppNum || i==NumBase ) continue;
1398     MapShapeNbElemsItr anIt = aResMap.find( meshFaces[i-1] );
1399     if( anIt == aResMap.end() ) continue;
1400     std::vector<int> aVec = (*anIt).second;
1401     nb2d += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
1402   }
1403
1404   MapShapeNbElemsItr anIt = aResMap.find( meshFaces[NumBase-1] );
1405   std::vector<int> aVec = (*anIt).second;
1406   bool IsQuadratic = (aVec[SMDSEntity_Quad_Triangle]>aVec[SMDSEntity_Triangle]) ||
1407                      (aVec[SMDSEntity_Quad_Quadrangle]>aVec[SMDSEntity_Quadrangle]);
1408   int nb2d_face0_3 = Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
1409   int nb2d_face0_4 = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
1410   int nb0d_face0 = aVec[SMDSEntity_Node];
1411   int nb1d_face0_int = ( nb2d_face0_3*3 + nb2d_face0_4*4 - nb1d ) / 2;
1412
1413   std::vector<int> aResVec(SMDSEntity_Last);
1414   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
1415   if(IsQuadratic) {
1416     aResVec[SMDSEntity_Quad_Penta] = nb2d_face0_3 * ( nb2d/nb1d );
1417     aResVec[SMDSEntity_Quad_Hexa] = nb2d_face0_4 * ( nb2d/nb1d );
1418     aResVec[SMDSEntity_Node] = nb0d_face0 * ( 2*nb2d/nb1d - 1 ) - nb1d_face0_int * nb2d/nb1d;
1419   }
1420   else {
1421     aResVec[SMDSEntity_Node] = nb0d_face0 * ( nb2d/nb1d - 1 );
1422     aResVec[SMDSEntity_Penta] = nb2d_face0_3 * ( nb2d/nb1d );
1423     aResVec[SMDSEntity_Hexa] = nb2d_face0_4 * ( nb2d/nb1d );
1424   }
1425   SMESH_subMesh * sm = theMesh.GetSubMesh(theShape);
1426   aResMap.insert(std::make_pair(sm,aResVec));
1427
1428   return true;
1429 }
1430
1431 //================================================================================
1432 /*!
1433  * \brief Create prisms
1434  * \param columns - columns of nodes generated from nodes of a mesh face
1435  * \param helper - helper initialized by mesh and shape to add prisms to
1436  */
1437 //================================================================================
1438
1439 void StdMeshers_Prism_3D::AddPrisms( vector<const TNodeColumn*> & columns,
1440                                      SMESH_MesherHelper*          helper)
1441 {
1442   int nbNodes = columns.size();
1443   int nbZ     = columns[0]->size();
1444   if ( nbZ < 2 ) return;
1445
1446   // find out orientation
1447   bool isForward = true;
1448   SMDS_VolumeTool vTool;
1449   int z = 1;
1450   switch ( nbNodes ) {
1451   case 3: {
1452     SMDS_VolumeOfNodes tmpPenta ( (*columns[0])[z-1], // bottom
1453                                   (*columns[1])[z-1],
1454                                   (*columns[2])[z-1],
1455                                   (*columns[0])[z],   // top
1456                                   (*columns[1])[z],
1457                                   (*columns[2])[z] );
1458     vTool.Set( &tmpPenta );
1459     isForward  = vTool.IsForward();
1460     break;
1461   }
1462   case 4: {
1463     SMDS_VolumeOfNodes tmpHex( (*columns[0])[z-1], (*columns[1])[z-1], // bottom
1464                                (*columns[2])[z-1], (*columns[3])[z-1],
1465                                (*columns[0])[z],   (*columns[1])[z],   // top
1466                                (*columns[2])[z],   (*columns[3])[z] );
1467     vTool.Set( &tmpHex );
1468     isForward  = vTool.IsForward();
1469     break;
1470   }
1471   default:
1472     const int di = (nbNodes+1) / 3;
1473     SMDS_VolumeOfNodes tmpVol ( (*columns[0]   )[z-1],
1474                                 (*columns[di]  )[z-1],
1475                                 (*columns[2*di])[z-1],
1476                                 (*columns[0]   )[z],
1477                                 (*columns[di]  )[z],
1478                                 (*columns[2*di])[z] );
1479     vTool.Set( &tmpVol );
1480     isForward  = vTool.IsForward();
1481   }
1482
1483   // vertical loop on columns
1484
1485   helper->SetElementsOnShape( true );
1486
1487   switch ( nbNodes ) {
1488
1489   case 3: { // ---------- pentahedra
1490     const int i1 = isForward ? 1 : 2;
1491     const int i2 = isForward ? 2 : 1;
1492     for ( z = 1; z < nbZ; ++z )
1493       helper->AddVolume( (*columns[0 ])[z-1], // bottom
1494                          (*columns[i1])[z-1],
1495                          (*columns[i2])[z-1],
1496                          (*columns[0 ])[z],   // top
1497                          (*columns[i1])[z],
1498                          (*columns[i2])[z] );
1499     break;
1500   }
1501   case 4: { // ---------- hexahedra
1502     const int i1 = isForward ? 1 : 3;
1503     const int i3 = isForward ? 3 : 1;
1504     for ( z = 1; z < nbZ; ++z )
1505       helper->AddVolume( (*columns[0])[z-1], (*columns[i1])[z-1], // bottom
1506                          (*columns[2])[z-1], (*columns[i3])[z-1],
1507                          (*columns[0])[z],   (*columns[i1])[z],     // top
1508                          (*columns[2])[z],   (*columns[i3])[z] );
1509     break;
1510   }
1511   case 6: { // ---------- octahedra
1512     const int iBase1 = isForward ? -1 : 0;
1513     const int iBase2 = isForward ?  0 :-1;
1514     for ( z = 1; z < nbZ; ++z )
1515       helper->AddVolume( (*columns[0])[z+iBase1], (*columns[1])[z+iBase1], // bottom or top
1516                          (*columns[2])[z+iBase1], (*columns[3])[z+iBase1],
1517                          (*columns[4])[z+iBase1], (*columns[5])[z+iBase1],
1518                          (*columns[0])[z+iBase2], (*columns[1])[z+iBase2], // top or bottom
1519                          (*columns[2])[z+iBase2], (*columns[3])[z+iBase2],
1520                          (*columns[4])[z+iBase2], (*columns[5])[z+iBase2] );
1521     break;
1522   }
1523   default: // ---------- polyhedra
1524     vector<int> quantities( 2 + nbNodes, 4 );
1525     quantities[0] = quantities[1] = nbNodes;
1526     columns.resize( nbNodes + 1 );
1527     columns[ nbNodes ] = columns[ 0 ];
1528     const int i1 = isForward ? 1 : 3;
1529     const int i3 = isForward ? 3 : 1;
1530     const int iBase1 = isForward ? -1 : 0;
1531     const int iBase2 = isForward ?  0 :-1;
1532     vector<const SMDS_MeshNode*> nodes( 2*nbNodes + 4*nbNodes);
1533     for ( z = 1; z < nbZ; ++z )
1534     {
1535       for ( int i = 0; i < nbNodes; ++i ) {
1536         nodes[ i             ] = (*columns[ i ])[z+iBase1]; // bottom or top
1537         nodes[ 2*nbNodes-i-1 ] = (*columns[ i ])[z+iBase2]; // top or bottom
1538         // side
1539         int di = 2*nbNodes + 4*i;
1540         nodes[ di+0 ] = (*columns[i  ])[z  ];
1541         nodes[ di+i1] = (*columns[i+1])[z  ];
1542         nodes[ di+2 ] = (*columns[i+1])[z-1];
1543         nodes[ di+i3] = (*columns[i  ])[z-1];
1544       }
1545       helper->AddPolyhedralVolume( nodes, quantities );
1546     }
1547
1548   } // switch ( nbNodes )
1549 }
1550
1551 //================================================================================
1552 /*!
1553  * \brief Find correspondence between bottom and top nodes
1554  *  If elements on the bottom and top faces are topologically different,
1555  *  and projection is possible and allowed, perform the projection
1556  *  \retval bool - is a success or not
1557  */
1558 //================================================================================
1559
1560 bool StdMeshers_Prism_3D::assocOrProjBottom2Top()
1561 {
1562   SMESH_subMesh * botSM = myBlock.SubMesh( ID_BOT_FACE );
1563   SMESH_subMesh * topSM = myBlock.SubMesh( ID_TOP_FACE );
1564
1565   SMESHDS_SubMesh * botSMDS = botSM->GetSubMeshDS();
1566   SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS();
1567
1568   if ( !botSMDS || botSMDS->NbElements() == 0 )
1569   {
1570     _gen->Compute( *myHelper->GetMesh(), botSM->GetSubShape() );
1571     botSMDS = botSM->GetSubMeshDS();
1572     if ( !botSMDS || botSMDS->NbElements() == 0 )
1573       return toSM( error(TCom("No elements on face #") << botSM->GetId() ));
1574   }
1575
1576   bool needProject = !topSM->IsMeshComputed();
1577   if ( !needProject &&
1578        (botSMDS->NbElements() != topSMDS->NbElements() ||
1579         botSMDS->NbNodes()    != topSMDS->NbNodes()))
1580   {
1581     MESSAGE("nb elem bot " << botSMDS->NbElements() <<
1582             " top " << ( topSMDS ? topSMDS->NbElements() : 0 ));
1583     MESSAGE("nb node bot " << botSMDS->NbNodes() <<
1584             " top " << ( topSMDS ? topSMDS->NbNodes() : 0 ));
1585     return toSM( error(TCom("Mesh on faces #") << botSM->GetId()
1586                        <<" and #"<< topSM->GetId() << " seems different" ));
1587   }
1588
1589   if ( 0/*needProject && !myProjectTriangles*/ )
1590     return toSM( error(TCom("Mesh on faces #") << botSM->GetId()
1591                        <<" and #"<< topSM->GetId() << " seems different" ));
1592   ///RETURN_BAD_RESULT("Need to project but not allowed");
1593
1594   if ( needProject )
1595   {
1596     return projectBottomToTop();
1597   }
1598
1599   TopoDS_Face botFace = TopoDS::Face( myBlock.Shape( ID_BOT_FACE ));
1600   TopoDS_Face topFace = TopoDS::Face( myBlock.Shape( ID_TOP_FACE ));
1601   // associate top and bottom faces
1602   TAssocTool::TShapeShapeMap shape2ShapeMap;
1603   if ( !TAssocTool::FindSubShapeAssociation( botFace, myBlock.Mesh(),
1604                                              topFace, myBlock.Mesh(),
1605                                              shape2ShapeMap) )
1606     return toSM( error(TCom("Topology of faces #") << botSM->GetId()
1607                        <<" and #"<< topSM->GetId() << " seems different" ));
1608
1609   // Find matching nodes of top and bottom faces
1610   TNodeNodeMap n2nMap;
1611   if ( ! TAssocTool::FindMatchingNodesOnFaces( botFace, myBlock.Mesh(),
1612                                                topFace, myBlock.Mesh(),
1613                                                shape2ShapeMap, n2nMap ))
1614     return toSM( error(TCom("Mesh on faces #") << botSM->GetId()
1615                        <<" and #"<< topSM->GetId() << " seems different" ));
1616
1617   // Fill myBotToColumnMap
1618
1619   int zSize = myBlock.VerticalSize();
1620   //TNode prevTNode;
1621   TNodeNodeMap::iterator bN_tN = n2nMap.begin();
1622   for ( ; bN_tN != n2nMap.end(); ++bN_tN )
1623   {
1624     const SMDS_MeshNode* botNode = bN_tN->first;
1625     const SMDS_MeshNode* topNode = bN_tN->second;
1626     if ( botNode->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE )
1627       continue; // wall columns are contained in myBlock
1628     // create node column
1629     Prism_3D::TNode bN( botNode );
1630     TNode2ColumnMap::iterator bN_col =
1631       myBotToColumnMap.insert( make_pair ( bN, TNodeColumn() )).first;
1632     TNodeColumn & column = bN_col->second;
1633     column.resize( zSize );
1634     column.front() = botNode;
1635     column.back()  = topNode;
1636   }
1637   return true;
1638 }
1639
1640 //================================================================================
1641 /*!
1642  * \brief Remove quadrangles from the top face and
1643  * create triangles there by projection from the bottom
1644  * \retval bool - a success or not
1645  */
1646 //================================================================================
1647
1648 bool StdMeshers_Prism_3D::projectBottomToTop()
1649 {
1650   SMESHDS_Mesh*  meshDS = myBlock.MeshDS();
1651   SMESH_subMesh * botSM = myBlock.SubMesh( ID_BOT_FACE );
1652   SMESH_subMesh * topSM = myBlock.SubMesh( ID_TOP_FACE );
1653
1654   SMESHDS_SubMesh * botSMDS = botSM->GetSubMeshDS();
1655   SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS();
1656
1657   if ( topSMDS && topSMDS->NbElements() > 0 )
1658     topSM->ComputeStateEngine( SMESH_subMesh::CLEAN );
1659
1660   const TopoDS_Shape& botFace = myBlock.Shape( ID_BOT_FACE ); // oriented within the 3D SHAPE
1661   const TopoDS_Shape& topFace = myBlock.Shape( ID_TOP_FACE);
1662   int topFaceID = meshDS->ShapeToIndex( topFace );
1663
1664   // Fill myBotToColumnMap
1665
1666   int zSize = myBlock.VerticalSize();
1667   Prism_3D::TNode prevTNode;
1668   SMDS_NodeIteratorPtr nIt = botSMDS->GetNodes();
1669   while ( nIt->more() )
1670   {
1671     const SMDS_MeshNode* botNode = nIt->next();
1672     if ( botNode->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE )
1673       continue; // strange
1674     // compute bottom node params
1675     Prism_3D::TNode bN( botNode );
1676     gp_XYZ paramHint(-1,-1,-1);
1677     if ( prevTNode.IsNeighbor( bN ))
1678       paramHint = prevTNode.GetParams();
1679     if ( !myBlock.ComputeParameters( bN.GetCoords(), bN.ChangeParams(),
1680                                      ID_BOT_FACE, paramHint ))
1681       return toSM( error(TCom("Can't compute normalized parameters for node ")
1682                          << botNode->GetID() << " on the face #"<< botSM->GetId() ));
1683     prevTNode = bN;
1684     // compute top node coords
1685     gp_XYZ topXYZ; gp_XY topUV;
1686     if ( !myBlock.FacePoint( ID_TOP_FACE, bN.GetParams(), topXYZ ) ||
1687          !myBlock.FaceUV   ( ID_TOP_FACE, bN.GetParams(), topUV ))
1688       return toSM( error(TCom("Can't compute coordinates "
1689                               "by normalized parameters on the face #")<< topSM->GetId() ));
1690     SMDS_MeshNode * topNode = meshDS->AddNode( topXYZ.X(),topXYZ.Y(),topXYZ.Z() );
1691     meshDS->SetNodeOnFace( topNode, topFaceID, topUV.X(), topUV.Y() );
1692     // create node column
1693     TNode2ColumnMap::iterator bN_col =
1694       myBotToColumnMap.insert( make_pair ( bN, TNodeColumn() )).first;
1695     TNodeColumn & column = bN_col->second;
1696     column.resize( zSize );
1697     column.front() = botNode;
1698     column.back()  = topNode;
1699   }
1700
1701   // Create top faces
1702
1703   const bool oldSetElemsOnShape = myHelper->SetElementsOnShape( false );
1704
1705   // care of orientation;
1706   // if the bottom faces is orienetd OK then top faces must be reversed
1707   bool reverseTop = true;
1708   if ( myHelper->NbAncestors( botFace, *myBlock.Mesh(), TopAbs_SOLID ) > 1 )
1709     reverseTop = ! myHelper->IsReversedSubMesh( TopoDS::Face( botFace ));
1710   int iFrw, iRev, *iPtr = &( reverseTop ? iRev : iFrw );
1711
1712   // loop on bottom mesh faces
1713   SMDS_ElemIteratorPtr faceIt = botSMDS->GetElements();
1714   vector< const SMDS_MeshNode* > nodes;
1715   while ( faceIt->more() )
1716   {
1717     const SMDS_MeshElement* face = faceIt->next();
1718     if ( !face || face->GetType() != SMDSAbs_Face )
1719       continue;
1720
1721     // find top node in columns for each bottom node
1722     int nbNodes = face->NbCornerNodes();
1723     nodes.resize( nbNodes );
1724     for ( iFrw = 0, iRev = nbNodes-1; iFrw < nbNodes; ++iFrw, --iRev )
1725     {
1726       const SMDS_MeshNode* n = face->GetNode( *iPtr );
1727       if ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) {
1728         TNode2ColumnMap::iterator bot_column = myBotToColumnMap.find( n );
1729         if ( bot_column == myBotToColumnMap.end() )
1730           return toSM( error(TCom("No nodes found above node ") << n->GetID() ));
1731         nodes[ iFrw ] = bot_column->second.back();
1732       }
1733       else {
1734         const TNodeColumn* column = myBlock.GetNodeColumn( n );
1735         if ( !column )
1736           return toSM( error(TCom("No side nodes found above node ") << n->GetID() ));
1737         nodes[ iFrw ] = column->back();
1738       }
1739     }
1740     SMDS_MeshElement* newFace = 0;
1741     switch ( nbNodes ) {
1742
1743     case 3: {
1744       newFace = myHelper->AddFace(nodes[0], nodes[1], nodes[2]);
1745       break;
1746       }
1747     case 4: {
1748       newFace = myHelper->AddFace( nodes[0], nodes[1], nodes[2], nodes[3] );
1749       break;
1750       }
1751     default:
1752       newFace = meshDS->AddPolygonalFace( nodes );
1753     }
1754     if ( newFace )
1755       meshDS->SetMeshElementOnShape( newFace, topFaceID );
1756   }
1757
1758   myHelper->SetElementsOnShape( oldSetElemsOnShape );
1759
1760   return true;
1761 }
1762
1763 //=======================================================================
1764 //function : project2dMesh
1765 //purpose  : Project mesh faces from a source FACE of one prism (theSrcFace)
1766 //           to a source FACE of another prism (theTgtFace)
1767 //=======================================================================
1768
1769 bool StdMeshers_Prism_3D::project2dMesh(const TopoDS_Face& theSrcFace,
1770                                         const TopoDS_Face& theTgtFace)
1771 {
1772   TProjction2dAlgo* projector2D = TProjction2dAlgo::instance( this );
1773   projector2D->myHyp.SetSourceFace( theSrcFace );
1774   bool ok = projector2D->Compute( *myHelper->GetMesh(), theTgtFace );
1775
1776   SMESH_subMesh* tgtSM = myHelper->GetMesh()->GetSubMesh( theTgtFace );
1777   tgtSM->ComputeStateEngine       ( SMESH_subMesh::CHECK_COMPUTE_STATE );
1778   tgtSM->ComputeSubMeshStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
1779
1780   return ok;
1781 }
1782
1783 //================================================================================
1784 /*!
1785  * \brief Set projection coordinates of a node to a face and it's sub-shapes
1786  * \param faceID - the face given by in-block ID
1787  * \param params - node normalized parameters
1788  * \retval bool - is a success
1789  */
1790 //================================================================================
1791
1792 bool StdMeshers_Prism_3D::setFaceAndEdgesXYZ( const int faceID, const gp_XYZ& params, int z )
1793 {
1794   // find base and top edges of the face
1795   enum { BASE = 0, TOP, LEFT, RIGHT };
1796   vector< int > edgeVec; // 0-base, 1-top
1797   SMESH_Block::GetFaceEdgesIDs( faceID, edgeVec );
1798
1799   myBlock.EdgePoint( edgeVec[ BASE ], params, myShapeXYZ[ edgeVec[ BASE ]]);
1800   myBlock.EdgePoint( edgeVec[ TOP  ], params, myShapeXYZ[ edgeVec[ TOP ]]);
1801
1802   SHOWYXZ("\nparams ", params);
1803   SHOWYXZ("TOP is " <<edgeVec[ TOP ], myShapeXYZ[ edgeVec[ TOP]]);
1804   SHOWYXZ("BASE is "<<edgeVec[ BASE], myShapeXYZ[ edgeVec[ BASE]]);
1805
1806   if ( faceID == SMESH_Block::ID_Fx0z || faceID == SMESH_Block::ID_Fx1z )
1807   {
1808     myBlock.EdgePoint( edgeVec[ LEFT ], params, myShapeXYZ[ edgeVec[ LEFT ]]);
1809     myBlock.EdgePoint( edgeVec[ RIGHT ], params, myShapeXYZ[ edgeVec[ RIGHT ]]);
1810
1811     SHOWYXZ("VER "<<edgeVec[ LEFT], myShapeXYZ[ edgeVec[ LEFT]]);
1812     SHOWYXZ("VER "<<edgeVec[ RIGHT], myShapeXYZ[ edgeVec[ RIGHT]]);
1813   }
1814   myBlock.FacePoint( faceID, params, myShapeXYZ[ faceID ]);
1815   SHOWYXZ("FacePoint "<<faceID, myShapeXYZ[ faceID]);
1816
1817   return true;
1818 }
1819
1820 //=======================================================================
1821 //function : toSM
1822 //purpose  : If (!isOK), sets the error to a sub-mesh of a current SOLID
1823 //=======================================================================
1824
1825 bool StdMeshers_Prism_3D::toSM( bool isOK )
1826 {
1827   if ( mySetErrorToSM &&
1828        !isOK &&
1829        myHelper &&
1830        !myHelper->GetSubShape().IsNull() &&
1831        myHelper->GetSubShape().ShapeType() == TopAbs_SOLID)
1832   {
1833     SMESH_subMesh* sm = myHelper->GetMesh()->GetSubMesh( myHelper->GetSubShape() );
1834     sm->GetComputeError() = this->GetComputeError();
1835     // clear error in order not to return it twice
1836     _error = COMPERR_OK;
1837     _comment.clear();
1838   }
1839   return isOK;
1840 }
1841
1842 //=======================================================================
1843 //function : shapeID
1844 //purpose  : Return index of a shape
1845 //=======================================================================
1846
1847 int StdMeshers_Prism_3D::shapeID( const TopoDS_Shape& S )
1848 {
1849   if ( S.IsNull() ) return 0;
1850   if ( !myHelper  ) return -3;
1851   return myHelper->GetMeshDS()->ShapeToIndex( S );
1852 }
1853
1854 namespace Prism_3D
1855 {
1856   //================================================================================
1857   /*!
1858    * \brief Return true if this node and other one belong to one face
1859    */
1860   //================================================================================
1861
1862   bool Prism_3D::TNode::IsNeighbor( const Prism_3D::TNode& other ) const
1863   {
1864     if ( !other.myNode || !myNode ) return false;
1865
1866     SMDS_ElemIteratorPtr fIt = other.myNode->GetInverseElementIterator(SMDSAbs_Face);
1867     while ( fIt->more() )
1868       if ( fIt->next()->GetNodeIndex( myNode ) >= 0 )
1869         return true;
1870     return false;
1871   }
1872
1873   //================================================================================
1874   /*!
1875    * \brief Prism initialization
1876    */
1877   //================================================================================
1878
1879   void TPrismTopo::Clear()
1880   {
1881     myShape3D.Nullify();
1882     myTop.Nullify();
1883     myBottom.Nullify();
1884     myWallQuads.clear();
1885     myBottomEdges.clear();
1886     myNbEdgesInWires.clear();
1887     myWallQuads.clear();
1888   }
1889
1890 } // namespace Prism_3D
1891
1892 //================================================================================
1893 /*!
1894  * \brief Constructor. Initialization is needed
1895  */
1896 //================================================================================
1897
1898 StdMeshers_PrismAsBlock::StdMeshers_PrismAsBlock()
1899 {
1900   mySide = 0;
1901 }
1902
1903 StdMeshers_PrismAsBlock::~StdMeshers_PrismAsBlock()
1904 {
1905   Clear();
1906 }
1907 void StdMeshers_PrismAsBlock::Clear()
1908 {
1909   myHelper = 0;
1910   myShapeIDMap.Clear();
1911   myError.reset();
1912
1913   if ( mySide ) {
1914     delete mySide; mySide = 0;
1915   }
1916   myParam2ColumnMaps.clear();
1917   myShapeIndex2ColumnMap.clear();
1918 }
1919
1920 //=======================================================================
1921 //function : initPrism
1922 //purpose  : Analyse shape geometry and mesh.
1923 //           If there are triangles on one of faces, it becomes 'bottom'.
1924 //           thePrism.myBottom can be already set up.
1925 //=======================================================================
1926
1927 bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism,
1928                                     const TopoDS_Shape&   shape3D)
1929 {
1930   myHelper->SetSubShape( shape3D );
1931
1932   SMESH_subMesh* mainSubMesh = myHelper->GetMesh()->GetSubMeshContaining( shape3D );
1933   if ( !mainSubMesh ) return toSM( error(COMPERR_BAD_INPUT_MESH,"Null submesh of shape3D"));
1934
1935   // detect not-quad FACE sub-meshes of the 3D SHAPE
1936   list< SMESH_subMesh* > notQuadGeomSubMesh;
1937   list< SMESH_subMesh* > notQuadElemSubMesh;
1938   int nbFaces = 0;
1939   //
1940   SMESH_subMesh* anyFaceSM = 0;
1941   SMESH_subMeshIteratorPtr smIt = mainSubMesh->getDependsOnIterator(false,true);
1942   while ( smIt->more() )
1943   {
1944     SMESH_subMesh* sm = smIt->next();
1945     const TopoDS_Shape& face = sm->GetSubShape();
1946     if      ( face.ShapeType() > TopAbs_FACE ) break;
1947     else if ( face.ShapeType() < TopAbs_FACE ) continue;
1948     nbFaces++;
1949     anyFaceSM = sm;
1950
1951     // is quadrangle FACE?
1952     list< TopoDS_Edge > orderedEdges;
1953     list< int >         nbEdgesInWires;
1954     int nbWires = SMESH_Block::GetOrderedEdges( TopoDS::Face( face ), orderedEdges,
1955                                                 nbEdgesInWires );
1956     if ( nbWires != 1 || nbEdgesInWires.front() != 4 )
1957       notQuadGeomSubMesh.push_back( sm );
1958
1959     // look for not quadrangle mesh elements
1960     if ( SMESHDS_SubMesh* smDS = sm->GetSubMeshDS() )
1961       if ( !myHelper->IsSameElemGeometry( smDS, SMDSGeom_QUADRANGLE ))
1962         notQuadElemSubMesh.push_back( sm );
1963   }
1964
1965   int nbNotQuadMeshed = notQuadElemSubMesh.size();
1966   int       nbNotQuad = notQuadGeomSubMesh.size();
1967   bool     hasNotQuad = ( nbNotQuad || nbNotQuadMeshed );
1968
1969   // detect bad cases
1970   if ( nbNotQuadMeshed > 2 )
1971   {
1972     return toSM( error(COMPERR_BAD_INPUT_MESH,
1973                        TCom("More than 2 faces with not quadrangle elements: ")
1974                        <<nbNotQuadMeshed));
1975   }
1976   if ( nbNotQuad > 2 || !thePrism.myBottom.IsNull() )
1977   {
1978     // Issue 0020843 - one of side FACEs is quasi-quadrilateral (not 4 EDGEs).
1979     // Remove from notQuadGeomSubMesh faces meshed with regular grid
1980     int nbQuasiQuads = removeQuasiQuads( notQuadGeomSubMesh, myHelper,
1981                                          TQuadrangleAlgo::instance(this,myHelper) );
1982     nbNotQuad -= nbQuasiQuads;
1983     if ( nbNotQuad > 2 )
1984       return toSM( error(COMPERR_BAD_SHAPE,
1985                          TCom("More than 2 not quadrilateral faces: ") <<nbNotQuad));
1986     hasNotQuad = ( nbNotQuad || nbNotQuadMeshed );
1987   }
1988
1989   // Analyse mesh and topology of FACEs: choose the bottom sub-mesh.
1990   // If there are not quadrangle FACEs, they are top and bottom ones.
1991   // Not quadrangle FACEs must be only on top and bottom.
1992
1993   SMESH_subMesh * botSM = 0;
1994   SMESH_subMesh * topSM = 0;
1995
1996   if ( hasNotQuad ) // can chose a bottom FACE
1997   {
1998     if ( nbNotQuadMeshed > 0 ) botSM = notQuadElemSubMesh.front();
1999     else                       botSM = notQuadGeomSubMesh.front();
2000     if ( nbNotQuadMeshed > 1 ) topSM = notQuadElemSubMesh.back();
2001     else if ( nbNotQuad  > 1 ) topSM = notQuadGeomSubMesh.back();
2002
2003     if ( topSM == botSM ) {
2004       if ( nbNotQuadMeshed > 1 ) topSM = notQuadElemSubMesh.front();
2005       else                       topSM = notQuadGeomSubMesh.front();
2006     }
2007
2008     // detect mesh triangles on wall FACEs
2009     if ( nbNotQuad == 2 && nbNotQuadMeshed > 0 ) {
2010       bool ok = false;
2011       if ( nbNotQuadMeshed == 1 )
2012         ok = ( find( notQuadGeomSubMesh.begin(),
2013                      notQuadGeomSubMesh.end(), botSM ) != notQuadGeomSubMesh.end() );
2014       else
2015         ok = ( notQuadGeomSubMesh == notQuadElemSubMesh );
2016       if ( !ok )
2017         return toSM( error(COMPERR_BAD_INPUT_MESH,
2018                            "Side face meshed with not quadrangle elements"));
2019     }
2020   }
2021
2022   thePrism.myNotQuadOnTop = ( nbNotQuadMeshed > 1 );
2023
2024   // use thePrism.myBottom
2025   if ( !thePrism.myBottom.IsNull() )
2026   {
2027     if ( botSM ) {
2028       if ( ! botSM->GetSubShape().IsSame( thePrism.myBottom )) {
2029         std::swap( botSM, topSM );
2030         if ( ! botSM->GetSubShape().IsSame( thePrism.myBottom ))
2031           return toSM( error( COMPERR_BAD_INPUT_MESH,
2032                               "Incompatible non-structured sub-meshes"));
2033       }
2034     }
2035     else {
2036       botSM = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom );
2037     }
2038   }
2039   else if ( !botSM ) // find a proper bottom
2040   {
2041     // composite walls or not prism shape
2042     for ( TopExp_Explorer f( shape3D, TopAbs_FACE ); f.More(); f.Next() )
2043     {
2044       int minNbFaces = 2 + myHelper->Count( f.Current(), TopAbs_EDGE, false);
2045       if ( nbFaces >= minNbFaces)
2046       {
2047         thePrism.Clear();
2048         thePrism.myBottom = TopoDS::Face( f.Current() );
2049         if ( initPrism( thePrism, shape3D ))
2050           return true;
2051       }
2052       return toSM( error( COMPERR_BAD_SHAPE ));
2053     }
2054   }
2055
2056   // find vertex 000 - the one with smallest coordinates (for easy DEBUG :-)
2057   TopoDS_Vertex V000;
2058   double minVal = DBL_MAX, minX, val;
2059   for ( TopExp_Explorer exp( botSM->GetSubShape(), TopAbs_VERTEX );
2060         exp.More(); exp.Next() )
2061   {
2062     const TopoDS_Vertex& v = TopoDS::Vertex( exp.Current() );
2063     gp_Pnt P = BRep_Tool::Pnt( v );
2064     val = P.X() + P.Y() + P.Z();
2065     if ( val < minVal || ( val == minVal && P.X() < minX )) {
2066       V000 = v;
2067       minVal = val;
2068       minX = P.X();
2069     }
2070   }
2071
2072   thePrism.myShape3D = shape3D;
2073   if ( thePrism.myBottom.IsNull() )
2074     thePrism.myBottom  = TopoDS::Face( botSM->GetSubShape() );
2075   thePrism.myBottom.Orientation( myHelper->GetSubShapeOri( shape3D,
2076                                                            thePrism.myBottom ));
2077   // Get ordered bottom edges
2078   TopoDS_Face reverseBottom = // to have order of top EDGEs as in the top FACE
2079     TopoDS::Face( thePrism.myBottom.Reversed() );
2080   SMESH_Block::GetOrderedEdges( reverseBottom,
2081                                 thePrism.myBottomEdges,
2082                                 thePrism.myNbEdgesInWires, V000 );
2083
2084   // Get Wall faces corresponding to the ordered bottom edges and the top FACE
2085   if ( !getWallFaces( thePrism, nbFaces ))
2086     return false; //toSM( error(COMPERR_BAD_SHAPE, "Can't find side faces"));
2087
2088   if ( topSM )
2089   {
2090     if ( !thePrism.myTop.IsSame( topSM->GetSubShape() ))
2091       return toSM( error
2092                    (notQuadGeomSubMesh.empty() ? COMPERR_BAD_INPUT_MESH : COMPERR_BAD_SHAPE,
2093                     "Non-quadrilateral faces are not opposite"));
2094
2095     // check that the found top and bottom FACEs are opposite
2096     list< TopoDS_Edge >::iterator edge = thePrism.myBottomEdges.begin();
2097     for ( ; edge != thePrism.myBottomEdges.end(); ++edge )
2098       if ( myHelper->IsSubShape( *edge, thePrism.myTop ))
2099         return toSM( error
2100                      (notQuadGeomSubMesh.empty() ? COMPERR_BAD_INPUT_MESH : COMPERR_BAD_SHAPE,
2101                       "Non-quadrilateral faces are not opposite"));
2102   }
2103
2104   return true;
2105 }
2106
2107 //================================================================================
2108 /*!
2109  * \brief Initialization.
2110  * \param helper - helper loaded with mesh and 3D shape
2111  * \param thePrism - a prosm data
2112  * \retval bool - false if a mesh or a shape are KO
2113  */
2114 //================================================================================
2115
2116 bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper*         helper,
2117                                    const Prism_3D::TPrismTopo& thePrism)
2118 {
2119   if ( mySide ) {
2120     delete mySide; mySide = 0;
2121   }
2122   vector< TSideFace* >         sideFaces( NB_WALL_FACES, 0 );
2123   vector< pair< double, double> > params( NB_WALL_FACES );
2124   mySide = new TSideFace( sideFaces, params );
2125
2126   myHelper = helper;
2127   SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
2128
2129   SMESH_Block::init();
2130   myShapeIDMap.Clear();
2131   myShapeIndex2ColumnMap.clear();
2132
2133   int wallFaceIds[ NB_WALL_FACES ] = { // to walk around a block
2134     SMESH_Block::ID_Fx0z, SMESH_Block::ID_F1yz,
2135     SMESH_Block::ID_Fx1z, SMESH_Block::ID_F0yz
2136   };
2137
2138   myError = SMESH_ComputeError::New();
2139
2140   myNotQuadOnTop = thePrism.myNotQuadOnTop;
2141
2142   // Find columns of wall nodes and calculate edges' lengths
2143   // --------------------------------------------------------
2144
2145   myParam2ColumnMaps.clear();
2146   myParam2ColumnMaps.resize( thePrism.myBottomEdges.size() ); // total nb edges
2147
2148   size_t iE, nbEdges = thePrism.myNbEdgesInWires.front(); // nb outer edges
2149   vector< double >    edgeLength( nbEdges );
2150   multimap< double, int > len2edgeMap;
2151
2152   list< TopoDS_Edge >::const_iterator edgeIt = thePrism.myBottomEdges.begin();
2153   for ( iE = 0; iE < nbEdges; ++iE, ++edgeIt )
2154   {
2155     TParam2ColumnMap & faceColumns = myParam2ColumnMaps[ iE ];
2156
2157     Prism_3D::TQuadList::const_iterator quad = thePrism.myWallQuads[ iE ].begin();
2158     for ( ; quad != thePrism.myWallQuads[ iE ].end(); ++quad )
2159     {
2160       const TopoDS_Edge& quadBot = (*quad)->side[ QUAD_BOTTOM_SIDE ]->Edge( 0 );
2161       if ( !myHelper->LoadNodeColumns( faceColumns, (*quad)->face, quadBot, meshDS ))
2162         return error(COMPERR_BAD_INPUT_MESH, TCom("Can't find regular quadrangle mesh ")
2163                      << "on a side face #" << MeshDS()->ShapeToIndex( (*quad)->face ));
2164     }
2165     SHOWYXZ("\np1 F " <<iE, gpXYZ(faceColumns.begin()->second.front() ));
2166     SHOWYXZ("p2 F "   <<iE, gpXYZ(faceColumns.rbegin()->second.front() ));
2167     SHOWYXZ("V First "<<iE, BRep_Tool::Pnt( TopExp::FirstVertex(*edgeIt,true )));
2168
2169     edgeLength[ iE ] = SMESH_Algo::EdgeLength( *edgeIt );
2170
2171     if ( nbEdges < NB_WALL_FACES ) // fill map used to split faces
2172     {
2173       SMESHDS_SubMesh* smDS = meshDS->MeshElements( *edgeIt);
2174       if ( !smDS )
2175         return error(COMPERR_BAD_INPUT_MESH, TCom("Null submesh on the edge #")
2176                      << MeshDS()->ShapeToIndex( *edgeIt ));
2177       len2edgeMap.insert( make_pair( edgeLength[ iE ], iE ));
2178     }
2179   }
2180   // Load columns of internal edges (forming holes)
2181   // and fill map ShapeIndex to TParam2ColumnMap for them
2182   for ( ; edgeIt != thePrism.myBottomEdges.end() ; ++edgeIt, ++iE )
2183   {
2184     TParam2ColumnMap & faceColumns = myParam2ColumnMaps[ iE ];
2185
2186     Prism_3D::TQuadList::const_iterator quad = thePrism.myWallQuads[ iE ].begin();
2187     for ( ; quad != thePrism.myWallQuads[ iE ].end(); ++quad )
2188     {
2189       const TopoDS_Edge& quadBot = (*quad)->side[ QUAD_BOTTOM_SIDE ]->Edge( 0 );
2190       if ( !myHelper->LoadNodeColumns( faceColumns, (*quad)->face, quadBot, meshDS ))
2191         return error(COMPERR_BAD_INPUT_MESH, TCom("Can't find regular quadrangle mesh ")
2192                      << "on a side face #" << MeshDS()->ShapeToIndex( (*quad)->face ));
2193     }
2194     // edge columns
2195     int id = MeshDS()->ShapeToIndex( *edgeIt );
2196     bool isForward = true; // meaningless for intenal wires
2197     myShapeIndex2ColumnMap[ id ] = make_pair( & faceColumns, isForward );
2198     // columns for vertices
2199     // 1
2200     const SMDS_MeshNode* n0 = faceColumns.begin()->second.front();
2201     id = n0->getshapeId();
2202     myShapeIndex2ColumnMap[ id ] = make_pair( & faceColumns, isForward );
2203     // 2
2204     const SMDS_MeshNode* n1 = faceColumns.rbegin()->second.front();
2205     id = n1->getshapeId();
2206     myShapeIndex2ColumnMap[ id ] = make_pair( & faceColumns, isForward );
2207
2208     // SHOWYXZ("\np1 F " <<iE, gpXYZ(faceColumns.begin()->second.front() ));
2209     // SHOWYXZ("p2 F "   <<iE, gpXYZ(faceColumns.rbegin()->second.front() ));
2210     // SHOWYXZ("V First "<<iE, BRep_Tool::Pnt( TopExp::FirstVertex(*edgeIt,true )));
2211   }
2212
2213   // Create 4 wall faces of a block
2214   // -------------------------------
2215
2216   if ( nbEdges <= NB_WALL_FACES ) // ************* Split faces if necessary
2217   {
2218     map< int, int > iE2nbSplit;
2219     if ( nbEdges != NB_WALL_FACES ) // define how to split
2220     {
2221       if ( len2edgeMap.size() != nbEdges )
2222         RETURN_BAD_RESULT("Uniqueness of edge lengths not assured");
2223
2224       multimap< double, int >::reverse_iterator maxLen_i = len2edgeMap.rbegin();
2225       multimap< double, int >::reverse_iterator midLen_i = ++len2edgeMap.rbegin();
2226
2227       double maxLen = maxLen_i->first;
2228       double midLen = ( len2edgeMap.size() == 1 ) ? 0 : midLen_i->first;
2229       switch ( nbEdges ) {
2230       case 1: // 0-th edge is split into 4 parts
2231         iE2nbSplit.insert( make_pair( 0, 4 )); break;
2232       case 2: // either the longest edge is split into 3 parts, or both edges into halves
2233         if ( maxLen / 3 > midLen / 2 ) {
2234           iE2nbSplit.insert( make_pair( maxLen_i->second, 3 ));
2235         }
2236         else {
2237           iE2nbSplit.insert( make_pair( maxLen_i->second, 2 ));
2238           iE2nbSplit.insert( make_pair( midLen_i->second, 2 ));
2239         }
2240         break;
2241       case 3:
2242         // split longest into halves
2243         iE2nbSplit.insert( make_pair( maxLen_i->second, 2 ));
2244       }
2245     }
2246     // Create TSideFace's
2247     int iSide = 0;
2248     list< TopoDS_Edge >::const_iterator botE = thePrism.myBottomEdges.begin();
2249     for ( iE = 0; iE < nbEdges; ++iE, ++botE )
2250     {
2251       TFaceQuadStructPtr quad = thePrism.myWallQuads[ iE ].front();
2252       // split?
2253       map< int, int >::iterator i_nb = iE2nbSplit.find( iE );
2254       if ( i_nb != iE2nbSplit.end() ) {
2255         // split!
2256         int nbSplit = i_nb->second;
2257         vector< double > params;
2258         splitParams( nbSplit, &myParam2ColumnMaps[ iE ], params );
2259         const bool isForward =
2260           StdMeshers_PrismAsBlock::IsForwardEdge( myHelper->GetMeshDS(),
2261                                                   myParam2ColumnMaps[iE],
2262                                                   *botE, SMESH_Block::ID_Fx0z );
2263         for ( int i = 0; i < nbSplit; ++i ) {
2264           double f = ( isForward ? params[ i ]   : params[ nbSplit - i-1 ]);
2265           double l = ( isForward ? params[ i+1 ] : params[ nbSplit - i ]);
2266           TSideFace* comp = new TSideFace( myHelper, wallFaceIds[ iSide ],
2267                                            thePrism.myWallQuads[ iE ], *botE,
2268                                            &myParam2ColumnMaps[ iE ], f, l );
2269           mySide->SetComponent( iSide++, comp );
2270         }
2271       }
2272       else {
2273         TSideFace* comp = new TSideFace( myHelper, wallFaceIds[ iSide ],
2274                                          thePrism.myWallQuads[ iE ], *botE,
2275                                          &myParam2ColumnMaps[ iE ]);
2276         mySide->SetComponent( iSide++, comp );
2277       }
2278     }
2279   }
2280   else { // **************************** Unite faces
2281
2282     // unite first faces
2283     int nbExraFaces = nbEdges - 3;
2284     int iSide = 0, iE;
2285     double u0 = 0, sumLen = 0;
2286     for ( iE = 0; iE < nbExraFaces; ++iE )
2287       sumLen += edgeLength[ iE ];
2288
2289     vector< TSideFace* >        components( nbExraFaces );
2290     vector< pair< double, double> > params( nbExraFaces );
2291     list< TopoDS_Edge >::const_iterator botE = thePrism.myBottomEdges.begin();
2292     for ( iE = 0; iE < nbExraFaces; ++iE, ++botE )
2293     {
2294       components[ iE ] = new TSideFace( myHelper, wallFaceIds[ iSide ],
2295                                         thePrism.myWallQuads[ iE ], *botE,
2296                                         &myParam2ColumnMaps[ iE ]);
2297       double u1 = u0 + edgeLength[ iE ] / sumLen;
2298       params[ iE ] = make_pair( u0 , u1 );
2299       u0 = u1;
2300     }
2301     mySide->SetComponent( iSide++, new TSideFace( components, params ));
2302
2303     // fill the rest faces
2304     for ( ; iE < nbEdges; ++iE, ++botE )
2305     {
2306       TSideFace* comp = new TSideFace( myHelper, wallFaceIds[ iSide ],
2307                                        thePrism.myWallQuads[ iE ], *botE,
2308                                        &myParam2ColumnMaps[ iE ]);
2309       mySide->SetComponent( iSide++, comp );
2310     }
2311   }
2312
2313
2314   // Fill geometry fields of SMESH_Block
2315   // ------------------------------------
2316
2317   vector< int > botEdgeIdVec;
2318   SMESH_Block::GetFaceEdgesIDs( ID_BOT_FACE, botEdgeIdVec );
2319
2320   bool isForward[NB_WALL_FACES] = { true, true, true, true };
2321   Adaptor2d_Curve2d* botPcurves[NB_WALL_FACES];
2322   Adaptor2d_Curve2d* topPcurves[NB_WALL_FACES];
2323
2324   for ( int iF = 0; iF < NB_WALL_FACES; ++iF )
2325   {
2326     TSideFace * sideFace = mySide->GetComponent( iF );
2327     if ( !sideFace )
2328       RETURN_BAD_RESULT("NULL TSideFace");
2329     int fID = sideFace->FaceID(); // in-block ID
2330
2331     // fill myShapeIDMap
2332     if ( sideFace->InsertSubShapes( myShapeIDMap ) != 8 &&
2333          !sideFace->IsComplex())
2334       MESSAGE( ": Warning : InsertSubShapes() < 8 on side " << iF );
2335
2336     // side faces geometry
2337     Adaptor2d_Curve2d* pcurves[NB_WALL_FACES];
2338     if ( !sideFace->GetPCurves( pcurves ))
2339       RETURN_BAD_RESULT("TSideFace::GetPCurves() failed");
2340
2341     SMESH_Block::TFace& tFace = myFace[ fID - ID_FirstF ];
2342     tFace.Set( fID, sideFace->Surface(), pcurves, isForward );
2343
2344     SHOWYXZ( endl<<"F "<< iF << " id " << fID << " FRW " << sideFace->IsForward(), sideFace->Value(0,0));
2345     // edges 3D geometry
2346     vector< int > edgeIdVec;
2347     SMESH_Block::GetFaceEdgesIDs( fID, edgeIdVec );
2348     for ( int isMax = 0; isMax < 2; ++isMax ) {
2349       {
2350         int eID = edgeIdVec[ isMax ];
2351         SMESH_Block::TEdge& tEdge = myEdge[ eID - ID_FirstE ];
2352         tEdge.Set( eID, sideFace->HorizCurve(isMax), true);
2353         SHOWYXZ(eID<<" HOR"<<isMax<<"(0)", sideFace->HorizCurve(isMax)->Value(0));
2354         SHOWYXZ(eID<<" HOR"<<isMax<<"(1)", sideFace->HorizCurve(isMax)->Value(1));
2355       }
2356       {
2357         int eID = edgeIdVec[ isMax+2 ];
2358         SMESH_Block::TEdge& tEdge = myEdge[ eID - ID_FirstE  ];
2359         tEdge.Set( eID, sideFace->VertiCurve(isMax), true);
2360         SHOWYXZ(eID<<" VER"<<isMax<<"(0)", sideFace->VertiCurve(isMax)->Value(0));
2361         SHOWYXZ(eID<<" VER"<<isMax<<"(1)", sideFace->VertiCurve(isMax)->Value(1));
2362
2363         // corner points
2364         vector< int > vertexIdVec;
2365         SMESH_Block::GetEdgeVertexIDs( eID, vertexIdVec );
2366         myPnt[ vertexIdVec[0] - ID_FirstV ] = tEdge.GetCurve()->Value(0).XYZ();
2367         myPnt[ vertexIdVec[1] - ID_FirstV ] = tEdge.GetCurve()->Value(1).XYZ();
2368       }
2369     }
2370     // pcurves on horizontal faces
2371     for ( iE = 0; iE < NB_WALL_FACES; ++iE ) {
2372       if ( edgeIdVec[ BOTTOM_EDGE ] == botEdgeIdVec[ iE ] ) {
2373         botPcurves[ iE ] = sideFace->HorizPCurve( false, thePrism.myBottom );
2374         topPcurves[ iE ] = sideFace->HorizPCurve( true,  thePrism.myTop );
2375         break;
2376       }
2377     }
2378     //sideFace->dumpNodes( 4 ); // debug
2379   }
2380   // horizontal faces geometry
2381   {
2382     SMESH_Block::TFace& tFace = myFace[ ID_BOT_FACE - ID_FirstF ];
2383     tFace.Set( ID_BOT_FACE, new BRepAdaptor_Surface( thePrism.myBottom ), botPcurves, isForward );
2384     SMESH_Block::Insert( thePrism.myBottom, ID_BOT_FACE, myShapeIDMap );
2385   }
2386   {
2387     SMESH_Block::TFace& tFace = myFace[ ID_TOP_FACE - ID_FirstF ];
2388     tFace.Set( ID_TOP_FACE, new BRepAdaptor_Surface( thePrism.myTop ), topPcurves, isForward );
2389     SMESH_Block::Insert( thePrism.myTop, ID_TOP_FACE, myShapeIDMap );
2390   }
2391
2392   // Fill map ShapeIndex to TParam2ColumnMap
2393   // ----------------------------------------
2394
2395   list< TSideFace* > fList;
2396   list< TSideFace* >::iterator fListIt;
2397   fList.push_back( mySide );
2398   for ( fListIt = fList.begin(); fListIt != fList.end(); ++fListIt)
2399   {
2400     int nb = (*fListIt)->NbComponents();
2401     for ( int i = 0; i < nb; ++i ) {
2402       if ( TSideFace* comp = (*fListIt)->GetComponent( i ))
2403         fList.push_back( comp );
2404     }
2405     if ( TParam2ColumnMap* cols = (*fListIt)->GetColumns()) {
2406       // columns for a base edge
2407       int id = MeshDS()->ShapeToIndex( (*fListIt)->BaseEdge() );
2408       bool isForward = (*fListIt)->IsForward();
2409       myShapeIndex2ColumnMap[ id ] = make_pair( cols, isForward );
2410
2411       // columns for vertices
2412       const SMDS_MeshNode* n0 = cols->begin()->second.front();
2413       id = n0->getshapeId();
2414       myShapeIndex2ColumnMap[ id ] = make_pair( cols, isForward );
2415
2416       const SMDS_MeshNode* n1 = cols->rbegin()->second.front();
2417       id = n1->getshapeId();
2418       myShapeIndex2ColumnMap[ id ] = make_pair( cols, !isForward );
2419     }
2420   }
2421
2422 //   gp_XYZ testPar(0.25, 0.25, 0), testCoord;
2423 //   if ( !FacePoint( ID_BOT_FACE, testPar, testCoord ))
2424 //     RETURN_BAD_RESULT("TEST FacePoint() FAILED");
2425 //   SHOWYXZ("IN TEST PARAM" , testPar);
2426 //   SHOWYXZ("OUT TEST CORD" , testCoord);
2427 //   if ( !ComputeParameters( testCoord, testPar , ID_BOT_FACE))
2428 //     RETURN_BAD_RESULT("TEST ComputeParameters() FAILED");
2429 //   SHOWYXZ("OUT TEST PARAM" , testPar);
2430
2431   return true;
2432 }
2433
2434 //================================================================================
2435 /*!
2436  * \brief Return pointer to column of nodes
2437  * \param node - bottom node from which the returned column goes up
2438  * \retval const TNodeColumn* - the found column
2439  */
2440 //================================================================================
2441
2442 const TNodeColumn* StdMeshers_PrismAsBlock::GetNodeColumn(const SMDS_MeshNode* node) const
2443 {
2444   int sID = node->getshapeId();
2445
2446   map<int, pair< TParam2ColumnMap*, bool > >::const_iterator col_frw =
2447     myShapeIndex2ColumnMap.find( sID );
2448   if ( col_frw != myShapeIndex2ColumnMap.end() ) {
2449     const TParam2ColumnMap* cols = col_frw->second.first;
2450     TParam2ColumnIt u_col = cols->begin();
2451     for ( ; u_col != cols->end(); ++u_col )
2452       if ( u_col->second[ 0 ] == node )
2453         return & u_col->second;
2454   }
2455   return 0;
2456 }
2457
2458 //=======================================================================
2459 //function : GetLayersTransformation
2460 //purpose  : Return transformations to get coordinates of nodes of each layer
2461 //           by nodes of the bottom. Layer is a set of nodes at a certain step
2462 //           from bottom to top.
2463 //=======================================================================
2464
2465 bool StdMeshers_PrismAsBlock::GetLayersTransformation(vector<gp_Trsf> &           trsf,
2466                                                       const Prism_3D::TPrismTopo& prism) const
2467 {
2468   const int zSize = VerticalSize();
2469   if ( zSize < 3 ) return true;
2470   trsf.resize( zSize - 2 );
2471
2472   // Select some node columns by which we will define coordinate system of layers
2473
2474   vector< const TNodeColumn* > columns;
2475   {
2476     bool isReverse;
2477     list< TopoDS_Edge >::const_iterator edgeIt = prism.myBottomEdges.begin();
2478     for ( int iE = 0; iE < prism.myNbEdgesInWires.front(); ++iE, ++edgeIt )
2479     {
2480       if ( BRep_Tool::Degenerated( *edgeIt )) continue;
2481       const TParam2ColumnMap* u2colMap =
2482         GetParam2ColumnMap( MeshDS()->ShapeToIndex( *edgeIt ), isReverse );
2483       if ( !u2colMap ) return false;
2484       double f = u2colMap->begin()->first, l = u2colMap->rbegin()->first;
2485       //isReverse = ( edgeIt->Orientation() == TopAbs_REVERSED );
2486       //if ( isReverse ) swap ( f, l ); -- u2colMap takes orientation into account
2487       const int nbCol = 5;
2488       for ( int i = 0; i < nbCol; ++i )
2489       {
2490         double u = f + i/double(nbCol) * ( l - f );
2491         const TNodeColumn* col = & getColumn( u2colMap, u )->second;
2492         if ( columns.empty() || col != columns.back() )
2493           columns.push_back( col );
2494       }
2495     }
2496   }
2497
2498   // Find tolerance to check transformations
2499
2500   double tol2;
2501   {
2502     Bnd_B3d bndBox;
2503     for ( int i = 0; i < columns.size(); ++i )
2504       bndBox.Add( gpXYZ( columns[i]->front() ));
2505     tol2 = bndBox.SquareExtent() * 1e-5;
2506   }
2507
2508   // Compute transformations
2509
2510   int xCol = -1;
2511   gp_Trsf fromCsZ, toCs0;
2512   gp_Ax3 cs0 = getLayerCoordSys(0, columns, xCol );
2513   //double dist0 = cs0.Location().Distance( gpXYZ( (*columns[0])[0]));
2514   toCs0.SetTransformation( cs0 );
2515   for ( int z = 1; z < zSize-1; ++z )
2516   {
2517     gp_Ax3 csZ = getLayerCoordSys(z, columns, xCol );
2518     //double distZ = csZ.Location().Distance( gpXYZ( (*columns[0])[z]));
2519     fromCsZ.SetTransformation( csZ );
2520     fromCsZ.Invert();
2521     gp_Trsf& t = trsf[ z-1 ];
2522     t = fromCsZ * toCs0;
2523     //t.SetScaleFactor( distZ/dist0 ); - it does not work properly, wrong base point
2524
2525     // check a transformation
2526     for ( int i = 0; i < columns.size(); ++i )
2527     {
2528       gp_Pnt p0 = gpXYZ( (*columns[i])[0] );
2529       gp_Pnt pz = gpXYZ( (*columns[i])[z] );
2530       t.Transforms( p0.ChangeCoord() );
2531       if ( p0.SquareDistance( pz ) > tol2 )
2532         return false;
2533     }
2534   }
2535   return true;
2536 }
2537
2538 //================================================================================
2539 /*!
2540  * \brief Check curve orientation of a bootom edge
2541   * \param meshDS - mesh DS
2542   * \param columnsMap - node columns map of side face
2543   * \param bottomEdge - the bootom edge
2544   * \param sideFaceID - side face in-block ID
2545   * \retval bool - true if orientation coinside with in-block forward orientation
2546  */
2547 //================================================================================
2548
2549 bool StdMeshers_PrismAsBlock::IsForwardEdge(SMESHDS_Mesh*           meshDS,
2550                                             const TParam2ColumnMap& columnsMap,
2551                                             const TopoDS_Edge &     bottomEdge,
2552                                             const int               sideFaceID)
2553 {
2554   bool isForward = false;
2555   if ( SMESH_MesherHelper::IsClosedEdge( bottomEdge ))
2556   {
2557     isForward = ( bottomEdge.Orientation() == TopAbs_FORWARD );
2558   }
2559   else
2560   {
2561     const TNodeColumn&     firstCol = columnsMap.begin()->second;
2562     const SMDS_MeshNode* bottomNode = firstCol[0];
2563     TopoDS_Shape firstVertex = SMESH_MesherHelper::GetSubShapeByNode( bottomNode, meshDS );
2564     isForward = ( firstVertex.IsSame( TopExp::FirstVertex( bottomEdge, true )));
2565   }
2566   // on 2 of 4 sides first vertex is end
2567   if ( sideFaceID == ID_Fx1z || sideFaceID == ID_F0yz )
2568     isForward = !isForward;
2569   return isForward;
2570 }
2571
2572 //================================================================================
2573 /*!
2574  * \brief Constructor
2575   * \param faceID - in-block ID
2576   * \param face - geom FACE
2577   * \param baseEdge - EDGE proreply oriented in the bottom EDGE !!!
2578   * \param columnsMap - map of node columns
2579   * \param first - first normalized param
2580   * \param last - last normalized param
2581  */
2582 //================================================================================
2583
2584 StdMeshers_PrismAsBlock::TSideFace::TSideFace(SMESH_MesherHelper*        helper,
2585                                               const int                  faceID,
2586                                               const Prism_3D::TQuadList& quadList,
2587                                               const TopoDS_Edge&         baseEdge,
2588                                               TParam2ColumnMap*          columnsMap,
2589                                               const double               first,
2590                                               const double               last):
2591   myID( faceID ),
2592   myParamToColumnMap( columnsMap ),
2593   myHelper( helper )
2594 {
2595   myParams.resize( 1 );
2596   myParams[ 0 ] = make_pair( first, last );
2597   mySurface     = PSurface( new BRepAdaptor_Surface( quadList.front()->face ));
2598   myBaseEdge    = baseEdge;
2599   myIsForward   = StdMeshers_PrismAsBlock::IsForwardEdge( myHelper->GetMeshDS(),
2600                                                           *myParamToColumnMap,
2601                                                           myBaseEdge, myID );
2602   if ( quadList.size() > 1 ) // side is vertically composite
2603   {
2604     // fill myShapeID2Surf map to enable finding a right surface by any sub-shape ID
2605
2606     SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
2607
2608     TopTools_IndexedDataMapOfShapeListOfShape subToFaces;
2609     Prism_3D::TQuadList::const_iterator quad = quadList.begin();
2610     for ( ; quad != quadList.end(); ++quad )
2611     {
2612       const TopoDS_Face& face = (*quad)->face;
2613       TopExp::MapShapesAndAncestors( face, TopAbs_VERTEX, TopAbs_FACE, subToFaces );
2614       TopExp::MapShapesAndAncestors( face, TopAbs_EDGE,   TopAbs_FACE, subToFaces );
2615       myShapeID2Surf.insert( make_pair( meshDS->ShapeToIndex( face ),
2616                                         PSurface( new BRepAdaptor_Surface( face ))));
2617     }
2618     for ( int i = 1; i <= subToFaces.Extent(); ++i )
2619     {
2620       const TopoDS_Shape&     sub = subToFaces.FindKey( i );
2621       TopTools_ListOfShape& faces = subToFaces( i );
2622       int subID  = meshDS->ShapeToIndex( sub );
2623       int faceID = meshDS->ShapeToIndex( faces.First() );
2624       myShapeID2Surf.insert ( make_pair( subID, myShapeID2Surf[ faceID ]));
2625     }
2626   }
2627 }
2628
2629 //================================================================================
2630 /*!
2631  * \brief Constructor of complex side face
2632  */
2633 //================================================================================
2634
2635 StdMeshers_PrismAsBlock::TSideFace::
2636 TSideFace(const vector< TSideFace* >&             components,
2637           const vector< pair< double, double> > & params)
2638   :myID( components[0] ? components[0]->myID : 0 ),
2639    myParamToColumnMap( 0 ),
2640    myParams( params ),
2641    myIsForward( true ),
2642    myComponents( components ),
2643    myHelper( components[0] ? components[0]->myHelper : 0 )
2644 {}
2645 //================================================================================
2646 /*!
2647  * \brief Copy constructor
2648   * \param other - other side
2649  */
2650 //================================================================================
2651
2652 StdMeshers_PrismAsBlock::TSideFace::TSideFace( const TSideFace& other )
2653 {
2654   myID               = other.myID;
2655   mySurface          = other.mySurface;
2656   myBaseEdge         = other.myBaseEdge;
2657   myParams           = other.myParams;
2658   myIsForward        = other.myIsForward;
2659   myHelper           = other.myHelper;
2660   myParamToColumnMap = other.myParamToColumnMap;
2661
2662   myComponents.resize( other.myComponents.size());
2663   for (int i = 0 ; i < myComponents.size(); ++i )
2664     myComponents[ i ] = new TSideFace( *other.myComponents[ i ]);
2665 }
2666
2667 //================================================================================
2668 /*!
2669  * \brief Deletes myComponents
2670  */
2671 //================================================================================
2672
2673 StdMeshers_PrismAsBlock::TSideFace::~TSideFace()
2674 {
2675   for (int i = 0 ; i < myComponents.size(); ++i )
2676     if ( myComponents[ i ] )
2677       delete myComponents[ i ];
2678 }
2679
2680 //================================================================================
2681 /*!
2682  * \brief Return geometry of the vertical curve
2683   * \param isMax - true means curve located closer to (1,1,1) block point
2684   * \retval Adaptor3d_Curve* - curve adaptor
2685  */
2686 //================================================================================
2687
2688 Adaptor3d_Curve* StdMeshers_PrismAsBlock::TSideFace::VertiCurve(const bool isMax) const
2689 {
2690   if ( !myComponents.empty() ) {
2691     if ( isMax )
2692       return myComponents.back()->VertiCurve(isMax);
2693     else
2694       return myComponents.front()->VertiCurve(isMax);
2695   }
2696   double f = myParams[0].first, l = myParams[0].second;
2697   if ( !myIsForward ) std::swap( f, l );
2698   return new TVerticalEdgeAdaptor( myParamToColumnMap, isMax ? l : f );
2699 }
2700
2701 //================================================================================
2702 /*!
2703  * \brief Return geometry of the top or bottom curve
2704   * \param isTop -
2705   * \retval Adaptor3d_Curve* -
2706  */
2707 //================================================================================
2708
2709 Adaptor3d_Curve* StdMeshers_PrismAsBlock::TSideFace::HorizCurve(const bool isTop) const
2710 {
2711   return new THorizontalEdgeAdaptor( this, isTop );
2712 }
2713
2714 //================================================================================
2715 /*!
2716  * \brief Return pcurves
2717   * \param pcurv - array of 4 pcurves
2718   * \retval bool - is a success
2719  */
2720 //================================================================================
2721
2722 bool StdMeshers_PrismAsBlock::TSideFace::GetPCurves(Adaptor2d_Curve2d* pcurv[4]) const
2723 {
2724   int iEdge[ 4 ] = { BOTTOM_EDGE, TOP_EDGE, V0_EDGE, V1_EDGE };
2725
2726   for ( int i = 0 ; i < 4 ; ++i ) {
2727     Handle(Geom2d_Line) line;
2728     switch ( iEdge[ i ] ) {
2729     case TOP_EDGE:
2730       line = new Geom2d_Line( gp_Pnt2d( 0, 1 ), gp::DX2d() ); break;
2731     case BOTTOM_EDGE:
2732       line = new Geom2d_Line( gp::Origin2d(), gp::DX2d() ); break;
2733     case V0_EDGE:
2734       line = new Geom2d_Line( gp::Origin2d(), gp::DY2d() ); break;
2735     case V1_EDGE:
2736       line = new Geom2d_Line( gp_Pnt2d( 1, 0 ), gp::DY2d() ); break;
2737     }
2738     pcurv[ i ] = new Geom2dAdaptor_Curve( line, 0, 1 );
2739   }
2740   return true;
2741 }
2742
2743 //================================================================================
2744 /*!
2745  * \brief Returns geometry of pcurve on a horizontal face
2746   * \param isTop - is top or bottom face
2747   * \param horFace - a horizontal face
2748   * \retval Adaptor2d_Curve2d* - curve adaptor
2749  */
2750 //================================================================================
2751
2752 Adaptor2d_Curve2d*
2753 StdMeshers_PrismAsBlock::TSideFace::HorizPCurve(const bool         isTop,
2754                                                 const TopoDS_Face& horFace) const
2755 {
2756   return new TPCurveOnHorFaceAdaptor( this, isTop, horFace );
2757 }
2758
2759 //================================================================================
2760 /*!
2761  * \brief Return a component corresponding to parameter
2762   * \param U - parameter along a horizontal size
2763   * \param localU - parameter along a horizontal size of a component
2764   * \retval TSideFace* - found component
2765  */
2766 //================================================================================
2767
2768 StdMeshers_PrismAsBlock::TSideFace*
2769 StdMeshers_PrismAsBlock::TSideFace::GetComponent(const double U,double & localU) const
2770 {
2771   localU = U;
2772   if ( myComponents.empty() )
2773     return const_cast<TSideFace*>( this );
2774
2775   int i;
2776   for ( i = 0; i < myComponents.size(); ++i )
2777     if ( U < myParams[ i ].second )
2778       break;
2779   if ( i >= myComponents.size() )
2780     i = myComponents.size() - 1;
2781
2782   double f = myParams[ i ].first, l = myParams[ i ].second;
2783   localU = ( U - f ) / ( l - f );
2784   return myComponents[ i ];
2785 }
2786
2787 //================================================================================
2788 /*!
2789  * \brief Find node columns for a parameter
2790   * \param U - parameter along a horizontal edge
2791   * \param col1 - the 1st found column
2792   * \param col2 - the 2nd found column
2793   * \retval r - normalized position of U between the found columns
2794  */
2795 //================================================================================
2796
2797 double StdMeshers_PrismAsBlock::TSideFace::GetColumns(const double      U,
2798                                                       TParam2ColumnIt & col1,
2799                                                       TParam2ColumnIt & col2) const
2800 {
2801   double u = U, r = 0;
2802   if ( !myComponents.empty() ) {
2803     TSideFace * comp = GetComponent(U,u);
2804     return comp->GetColumns( u, col1, col2 );
2805   }
2806
2807   if ( !myIsForward )
2808     u = 1 - u;
2809   double f = myParams[0].first, l = myParams[0].second;
2810   u = f + u * ( l - f );
2811
2812   col1 = col2 = getColumn( myParamToColumnMap, u );
2813   if ( ++col2 == myParamToColumnMap->end() ) {
2814     --col2;
2815     r = 0.5;
2816   }
2817   else {
2818     double uf = col1->first;
2819     double ul = col2->first;
2820     r = ( u - uf ) / ( ul - uf );
2821   }
2822   return r;
2823 }
2824
2825 //================================================================================
2826 /*!
2827  * \brief Return coordinates by normalized params
2828   * \param U - horizontal param
2829   * \param V - vertical param
2830   * \retval gp_Pnt - result point
2831  */
2832 //================================================================================
2833
2834 gp_Pnt StdMeshers_PrismAsBlock::TSideFace::Value(const Standard_Real U,
2835                                                  const Standard_Real V) const
2836 {
2837   if ( !myComponents.empty() ) {
2838     double u;
2839     TSideFace * comp = GetComponent(U,u);
2840     return comp->Value( u, V );
2841   }
2842
2843   TParam2ColumnIt u_col1, u_col2;
2844   double vR, hR = GetColumns( U, u_col1, u_col2 );
2845
2846   const SMDS_MeshNode* nn[4];
2847
2848   // BEGIN issue 0020680: Bad cell created by Radial prism in center of torus
2849   // Workaround for a wrongly located point returned by mySurface.Value() for
2850   // UV located near boundary of BSpline surface.
2851   // To bypass the problem, we take point from 3D curve of EDGE.
2852   // It solves pb of the bloc_fiss_new.py
2853   const double tol = 1e-3;
2854   if ( V < tol || V+tol >= 1. )
2855   {
2856     nn[0] = V < tol ? u_col1->second.front() : u_col1->second.back();
2857     nn[2] = V < tol ? u_col2->second.front() : u_col2->second.back();
2858     TopoDS_Edge edge;
2859     if ( V < tol )
2860     {
2861       edge = myBaseEdge;
2862     }
2863     else
2864     {
2865       TopoDS_Shape s = myHelper->GetSubShapeByNode( nn[0], myHelper->GetMeshDS() );
2866       if ( s.ShapeType() != TopAbs_EDGE )
2867         s = myHelper->GetSubShapeByNode( nn[2], myHelper->GetMeshDS() );
2868       if ( s.ShapeType() == TopAbs_EDGE )
2869         edge = TopoDS::Edge( s );
2870     }
2871     if ( !edge.IsNull() )
2872     {
2873       double u1 = myHelper->GetNodeU( edge, nn[0] );
2874       double u3 = myHelper->GetNodeU( edge, nn[2] );
2875       double u = u1 * ( 1 - hR ) + u3 * hR;
2876       TopLoc_Location loc; double f,l;
2877       Handle(Geom_Curve) curve = BRep_Tool::Curve( edge,loc,f,l );
2878       return curve->Value( u ).Transformed( loc );
2879     }
2880   }
2881   // END issue 0020680: Bad cell created by Radial prism in center of torus
2882
2883   vR = getRAndNodes( & u_col1->second, V, nn[0], nn[1] );
2884   vR = getRAndNodes( & u_col2->second, V, nn[2], nn[3] );
2885
2886   if ( !myShapeID2Surf.empty() ) // side is vertically composite
2887   {
2888     // find a FACE on which the 4 nodes lie
2889     TSideFace* me = (TSideFace*) this;
2890     int notFaceID1 = 0, notFaceID2 = 0;
2891     for ( int i = 0; i < 4; ++i )
2892       if ( nn[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) // node on FACE
2893       {
2894         me->mySurface = me->myShapeID2Surf[ nn[i]->getshapeId() ];
2895         notFaceID2 = 0;
2896         break;
2897       }
2898       else if ( notFaceID1 == 0 ) // node on EDGE or VERTEX
2899       {
2900         me->mySurface  = me->myShapeID2Surf[ nn[i]->getshapeId() ];
2901         notFaceID1 = nn[i]->getshapeId();
2902       }
2903       else if ( notFaceID1 != nn[i]->getshapeId() ) // node on other EDGE or VERTEX
2904       {
2905         if ( mySurface != me->myShapeID2Surf[ nn[i]->getshapeId() ])
2906           notFaceID2 = nn[i]->getshapeId();
2907       }
2908     if ( notFaceID2 ) // no nodes of FACE and nodes are on different FACEs
2909     {
2910       SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
2911       TopoDS_Shape face = myHelper->GetCommonAncestor( meshDS->IndexToShape( notFaceID1 ),
2912                                                        meshDS->IndexToShape( notFaceID2 ),
2913                                                        *myHelper->GetMesh(),
2914                                                        TopAbs_FACE );
2915       if ( face.IsNull() )
2916         throw SALOME_Exception("StdMeshers_PrismAsBlock::TSideFace::Value() face.IsNull()");
2917       int faceID = meshDS->ShapeToIndex( face );
2918       me->mySurface = me->myShapeID2Surf[ faceID ];
2919       if ( !mySurface )
2920         throw SALOME_Exception("StdMeshers_PrismAsBlock::TSideFace::Value() !mySurface");
2921     }
2922   }
2923
2924   gp_XY uv1 = myHelper->GetNodeUV( mySurface->Face(), nn[0], nn[2]);
2925   gp_XY uv2 = myHelper->GetNodeUV( mySurface->Face(), nn[1], nn[3]);
2926   gp_XY uv12 = uv1 * ( 1 - vR ) + uv2 * vR;
2927
2928   gp_XY uv3 = myHelper->GetNodeUV( mySurface->Face(), nn[2], nn[0]);
2929   gp_XY uv4 = myHelper->GetNodeUV( mySurface->Face(), nn[3], nn[1]);
2930   gp_XY uv34 = uv3 * ( 1 - vR ) + uv4 * vR;
2931
2932   gp_XY uv = uv12 * ( 1 - hR ) + uv34 * hR;
2933
2934   gp_Pnt p = mySurface->Value( uv.X(), uv.Y() );
2935   return p;
2936 }
2937
2938
2939 //================================================================================
2940 /*!
2941  * \brief Return boundary edge
2942   * \param edge - edge index
2943   * \retval TopoDS_Edge - found edge
2944  */
2945 //================================================================================
2946
2947 TopoDS_Edge StdMeshers_PrismAsBlock::TSideFace::GetEdge(const int iEdge) const
2948 {
2949   if ( !myComponents.empty() ) {
2950     switch ( iEdge ) {
2951     case V0_EDGE : return myComponents.front()->GetEdge( iEdge );
2952     case V1_EDGE : return myComponents.back() ->GetEdge( iEdge );
2953     default: return TopoDS_Edge();
2954     }
2955   }
2956   TopoDS_Shape edge;
2957   const SMDS_MeshNode* node = 0;
2958   SMESHDS_Mesh * meshDS = myHelper->GetMesh()->GetMeshDS();
2959   TNodeColumn* column;
2960
2961   switch ( iEdge ) {
2962   case TOP_EDGE:
2963   case BOTTOM_EDGE:
2964     column = & (( ++myParamToColumnMap->begin())->second );
2965     node = ( iEdge == TOP_EDGE ) ? column->back() : column->front();
2966     edge = myHelper->GetSubShapeByNode ( node, meshDS );
2967     if ( edge.ShapeType() == TopAbs_VERTEX ) {
2968       column = & ( myParamToColumnMap->begin()->second );
2969       node = ( iEdge == TOP_EDGE ) ? column->back() : column->front();
2970     }
2971     break;
2972   case V0_EDGE:
2973   case V1_EDGE: {
2974     bool back = ( iEdge == V1_EDGE );
2975     if ( !myIsForward ) back = !back;
2976     if ( back )
2977       column = & ( myParamToColumnMap->rbegin()->second );
2978     else
2979       column = & ( myParamToColumnMap->begin()->second );
2980     if ( column->size() > 0 )
2981       edge = myHelper->GetSubShapeByNode( (*column)[ 1 ], meshDS );
2982     if ( edge.IsNull() || edge.ShapeType() == TopAbs_VERTEX )
2983       node = column->front();
2984     break;
2985   }
2986   default:;
2987   }
2988   if ( !edge.IsNull() && edge.ShapeType() == TopAbs_EDGE )
2989     return TopoDS::Edge( edge );
2990
2991   // find edge by 2 vertices
2992   TopoDS_Shape V1 = edge;
2993   TopoDS_Shape V2 = myHelper->GetSubShapeByNode( node, meshDS );
2994   if ( !V2.IsNull() && V2.ShapeType() == TopAbs_VERTEX && !V2.IsSame( V1 ))
2995   {
2996     TopoDS_Shape ancestor = myHelper->GetCommonAncestor( V1, V2, *myHelper->GetMesh(), TopAbs_EDGE);
2997     if ( !ancestor.IsNull() )
2998       return TopoDS::Edge( ancestor );
2999   }
3000   return TopoDS_Edge();
3001 }
3002
3003 //================================================================================
3004 /*!
3005  * \brief Fill block sub-shapes
3006   * \param shapeMap - map to fill in
3007   * \retval int - nb inserted sub-shapes
3008  */
3009 //================================================================================
3010
3011 int StdMeshers_PrismAsBlock::TSideFace::InsertSubShapes(TBlockShapes& shapeMap) const
3012 {
3013   int nbInserted = 0;
3014
3015   // Insert edges
3016   vector< int > edgeIdVec;
3017   SMESH_Block::GetFaceEdgesIDs( myID, edgeIdVec );
3018
3019   for ( int i = BOTTOM_EDGE; i <=V1_EDGE ; ++i ) {
3020     TopoDS_Edge e = GetEdge( i );
3021     if ( !e.IsNull() ) {
3022       nbInserted += SMESH_Block::Insert( e, edgeIdVec[ i ], shapeMap);
3023     }
3024   }
3025
3026   // Insert corner vertices
3027
3028   TParam2ColumnIt col1, col2 ;
3029   vector< int > vertIdVec;
3030
3031   // from V0 column
3032   SMESH_Block::GetEdgeVertexIDs( edgeIdVec[ V0_EDGE ], vertIdVec);
3033   GetColumns(0, col1, col2 );
3034   const SMDS_MeshNode* node0 = col1->second.front();
3035   const SMDS_MeshNode* node1 = col1->second.back();
3036   TopoDS_Shape v0 = myHelper->GetSubShapeByNode( node0, myHelper->GetMeshDS());
3037   TopoDS_Shape v1 = myHelper->GetSubShapeByNode( node1, myHelper->GetMeshDS());
3038   if ( v0.ShapeType() == TopAbs_VERTEX ) {
3039     nbInserted += SMESH_Block::Insert( v0, vertIdVec[ 0 ], shapeMap);
3040   }
3041   if ( v1.ShapeType() == TopAbs_VERTEX ) {
3042     nbInserted += SMESH_Block::Insert( v1, vertIdVec[ 1 ], shapeMap);
3043   }
3044
3045   // from V1 column
3046   SMESH_Block::GetEdgeVertexIDs( edgeIdVec[ V1_EDGE ], vertIdVec);
3047   GetColumns(1, col1, col2 );
3048   node0 = col2->second.front();
3049   node1 = col2->second.back();
3050   v0 = myHelper->GetSubShapeByNode( node0, myHelper->GetMeshDS());
3051   v1 = myHelper->GetSubShapeByNode( node1, myHelper->GetMeshDS());
3052   if ( v0.ShapeType() == TopAbs_VERTEX ) {
3053     nbInserted += SMESH_Block::Insert( v0, vertIdVec[ 0 ], shapeMap);
3054   }
3055   if ( v1.ShapeType() == TopAbs_VERTEX ) {
3056     nbInserted += SMESH_Block::Insert( v1, vertIdVec[ 1 ], shapeMap);
3057   }
3058
3059 //   TopoDS_Vertex V0, V1, Vcom;
3060 //   TopExp::Vertices( myBaseEdge, V0, V1, true );
3061 //   if ( !myIsForward ) std::swap( V0, V1 );
3062
3063 //   // bottom vertex IDs
3064 //   SMESH_Block::GetEdgeVertexIDs( edgeIdVec[ _u0 ], vertIdVec);
3065 //   SMESH_Block::Insert( V0, vertIdVec[ 0 ], shapeMap);
3066 //   SMESH_Block::Insert( V1, vertIdVec[ 1 ], shapeMap);
3067
3068 //   TopoDS_Edge sideEdge = GetEdge( V0_EDGE );
3069 //   if ( sideEdge.IsNull() || !TopExp::CommonVertex( botEdge, sideEdge, Vcom ))
3070 //     return false;
3071
3072 //   // insert one side edge
3073 //   int edgeID;
3074 //   if ( Vcom.IsSame( V0 )) edgeID = edgeIdVec[ _v0 ];
3075 //   else                    edgeID = edgeIdVec[ _v1 ];
3076 //   SMESH_Block::Insert( sideEdge, edgeID, shapeMap);
3077
3078 //   // top vertex of the side edge
3079 //   SMESH_Block::GetEdgeVertexIDs( edgeID, vertIdVec);
3080 //   TopoDS_Vertex Vtop = TopExp::FirstVertex( sideEdge );
3081 //   if ( Vcom.IsSame( Vtop ))
3082 //     Vtop = TopExp::LastVertex( sideEdge );
3083 //   SMESH_Block::Insert( Vtop, vertIdVec[ 1 ], shapeMap);
3084
3085 //   // other side edge
3086 //   sideEdge = GetEdge( V1_EDGE );
3087 //   if ( sideEdge.IsNull() )
3088 //     return false;
3089 //   if ( edgeID = edgeIdVec[ _v1 ]) edgeID = edgeIdVec[ _v0 ];
3090 //   else                            edgeID = edgeIdVec[ _v1 ];
3091 //   SMESH_Block::Insert( sideEdge, edgeID, shapeMap);
3092
3093 //   // top edge
3094 //   TopoDS_Edge topEdge = GetEdge( TOP_EDGE );
3095 //   SMESH_Block::Insert( topEdge, edgeIdVec[ _u1 ], shapeMap);
3096
3097 //   // top vertex of the other side edge
3098 //   if ( !TopExp::CommonVertex( topEdge, sideEdge, Vcom ))
3099 //     return false;
3100 //   SMESH_Block::GetEdgeVertexIDs( edgeID, vertIdVec );
3101 //   SMESH_Block::Insert( Vcom, vertIdVec[ 1 ], shapeMap);
3102
3103   return nbInserted;
3104 }
3105
3106 //================================================================================
3107 /*!
3108  * \brief Dump ids of nodes of sides
3109  */
3110 //================================================================================
3111
3112 void StdMeshers_PrismAsBlock::TSideFace::dumpNodes(int nbNodes) const
3113 {
3114 #ifdef _DEBUG_
3115   cout << endl << "NODES OF FACE "; SMESH_Block::DumpShapeID( myID, cout ) << endl;
3116   THorizontalEdgeAdaptor* hSize0 = (THorizontalEdgeAdaptor*) HorizCurve(0);
3117   cout << "Horiz side 0: "; hSize0->dumpNodes(nbNodes); cout << endl;
3118   THorizontalEdgeAdaptor* hSize1 = (THorizontalEdgeAdaptor*) HorizCurve(1);
3119   cout << "Horiz side 1: "; hSize1->dumpNodes(nbNodes); cout << endl;
3120   TVerticalEdgeAdaptor* vSide0 = (TVerticalEdgeAdaptor*) VertiCurve(0);
3121   cout << "Verti side 0: "; vSide0->dumpNodes(nbNodes); cout << endl;
3122   TVerticalEdgeAdaptor* vSide1 = (TVerticalEdgeAdaptor*) VertiCurve(1);
3123   cout << "Verti side 1: "; vSide1->dumpNodes(nbNodes); cout << endl;
3124   delete hSize0; delete hSize1; delete vSide0; delete vSide1;
3125 #endif
3126 }
3127
3128 //================================================================================
3129 /*!
3130  * \brief Creates TVerticalEdgeAdaptor
3131   * \param columnsMap - node column map
3132   * \param parameter - normalized parameter
3133  */
3134 //================================================================================
3135
3136 StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::
3137 TVerticalEdgeAdaptor( const TParam2ColumnMap* columnsMap, const double parameter)
3138 {
3139   myNodeColumn = & getColumn( columnsMap, parameter )->second;
3140 }
3141
3142 //================================================================================
3143 /*!
3144  * \brief Return coordinates for the given normalized parameter
3145   * \param U - normalized parameter
3146   * \retval gp_Pnt - coordinates
3147  */
3148 //================================================================================
3149
3150 gp_Pnt StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::Value(const Standard_Real U) const
3151 {
3152   const SMDS_MeshNode* n1;
3153   const SMDS_MeshNode* n2;
3154   double r = getRAndNodes( myNodeColumn, U, n1, n2 );
3155   return gpXYZ(n1) * ( 1 - r ) + gpXYZ(n2) * r;
3156 }
3157
3158 //================================================================================
3159 /*!
3160  * \brief Dump ids of nodes
3161  */
3162 //================================================================================
3163
3164 void StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::dumpNodes(int nbNodes) const
3165 {
3166 #ifdef _DEBUG_
3167   for ( int i = 0; i < nbNodes && i < myNodeColumn->size(); ++i )
3168     cout << (*myNodeColumn)[i]->GetID() << " ";
3169   if ( nbNodes < myNodeColumn->size() )
3170     cout << myNodeColumn->back()->GetID();
3171 #endif
3172 }
3173
3174 //================================================================================
3175 /*!
3176  * \brief Return coordinates for the given normalized parameter
3177   * \param U - normalized parameter
3178   * \retval gp_Pnt - coordinates
3179  */
3180 //================================================================================
3181
3182 gp_Pnt StdMeshers_PrismAsBlock::THorizontalEdgeAdaptor::Value(const Standard_Real U) const
3183 {
3184   return mySide->TSideFace::Value( U, myV );
3185 }
3186
3187 //================================================================================
3188 /*!
3189  * \brief Dump ids of <nbNodes> first nodes and the last one
3190  */
3191 //================================================================================
3192
3193 void StdMeshers_PrismAsBlock::THorizontalEdgeAdaptor::dumpNodes(int nbNodes) const
3194 {
3195 #ifdef _DEBUG_
3196   // Not bedugged code. Last node is sometimes incorrect
3197   const TSideFace* side = mySide;
3198   double u = 0;
3199   if ( mySide->IsComplex() )
3200     side = mySide->GetComponent(0,u);
3201
3202   TParam2ColumnIt col, col2;
3203   TParam2ColumnMap* u2cols = side->GetColumns();
3204   side->GetColumns( u , col, col2 );
3205
3206   int j, i = myV ? mySide->ColumnHeight()-1 : 0;
3207
3208   const SMDS_MeshNode* n = 0;
3209   const SMDS_MeshNode* lastN
3210     = side->IsForward() ? u2cols->rbegin()->second[ i ] : u2cols->begin()->second[ i ];
3211   for ( j = 0; j < nbNodes && n != lastN; ++j )
3212   {
3213     n = col->second[ i ];
3214     cout << n->GetID() << " ";
3215     if ( side->IsForward() )
3216       ++col;
3217     else
3218       --col;
3219   }
3220
3221   // last node
3222   u = 1;
3223   if ( mySide->IsComplex() )
3224     side = mySide->GetComponent(1,u);
3225
3226   side->GetColumns( u , col, col2 );
3227   if ( n != col->second[ i ] )
3228     cout << col->second[ i ]->GetID();
3229 #endif
3230 }
3231 //================================================================================
3232 /*!
3233  * \brief Return UV on pcurve for the given normalized parameter
3234   * \param U - normalized parameter
3235   * \retval gp_Pnt - coordinates
3236  */
3237 //================================================================================
3238
3239 gp_Pnt2d StdMeshers_PrismAsBlock::TPCurveOnHorFaceAdaptor::Value(const Standard_Real U) const
3240 {
3241   TParam2ColumnIt u_col1, u_col2;
3242   double r = mySide->GetColumns( U, u_col1, u_col2 );
3243   gp_XY uv1 = mySide->GetNodeUV( myFace, u_col1->second[ myZ ]);
3244   gp_XY uv2 = mySide->GetNodeUV( myFace, u_col2->second[ myZ ]);
3245   return uv1 * ( 1 - r ) + uv2 * r;
3246 }