src/StdMeshers/StdMeshers_Prism_3D.cxx

   1 //  Copyright (C) 2007-2010  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 //  SMESH SMESH : implementaion of SMESH idl descriptions
  24 // File      : StdMeshers_Prism_3D.cxx
  25 // Module    : SMESH
  26 // Created   : Fri Oct 20 11:37:07 2006
  27 // Author    : Edward AGAPOV (eap)
  28 //
  29 #include "StdMeshers_Prism_3D.hxx"
  30
  31 #include "StdMeshers_ProjectionUtils.hxx"
  32 #include "SMESH_MesherHelper.hxx"
  33 #include "SMDS_VolumeTool.hxx"
  34 #include "SMDS_VolumeOfNodes.hxx"
  35 #include "SMDS_EdgePosition.hxx"
  36 #include "SMESH_Comment.hxx"
  37
  38 #include "utilities.h"
  39
  40 #include <BRep_Tool.hxx>
  41 #include <Bnd_B3d.hxx>
  42 #include <Geom2dAdaptor_Curve.hxx>
  43 #include <Geom2d_Line.hxx>
  44 #include <Geom_Curve.hxx>
  45 #include <TopExp.hxx>
  46 #include <TopExp_Explorer.hxx>
  47 #include <TopTools_ListIteratorOfListOfShape.hxx>
  48 #include <TopTools_MapOfShape.hxx>
  49 #include <TopTools_SequenceOfShape.hxx>
  50 #include <TopoDS.hxx>
  51 #include <gp_Ax2.hxx>
  52 #include <gp_Ax3.hxx>
  53
  54 using namespace std;
  55
  56 #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; }
  57 #define gpXYZ(n) gp_XYZ(n->X(),n->Y(),n->Z())
  58 #define SHOWYXZ(msg, xyz) // {\
  59 // gp_Pnt p (xyz); \
  60 // cout << msg << " ("<< p.X() << "; " <<p.Y() << "; " <<p.Z() << ") " <<endl;\
  61 // }
  62
  63 typedef StdMeshers_ProjectionUtils TAssocTool;
  64 typedef SMESH_Comment              TCom;
  65
  66 enum { ID_BOT_FACE = SMESH_Block::ID_Fxy0,
  67        ID_TOP_FACE = SMESH_Block::ID_Fxy1,
  68        BOTTOM_EDGE = 0, TOP_EDGE, V0_EDGE, V1_EDGE, // edge IDs in face
  69        NB_WALL_FACES = 4 }; //
  70
  71 namespace {
  72
  73   //================================================================================
  74   /*!
  75    * \brief Return iterator pointing to node column for the given parameter
  76    * \param columnsMap - node column map
  77    * \param parameter - parameter
  78    * \retval TParam2ColumnMap::iterator - result
  79    *
  80    * it returns closest left column
  81    */
  82   //================================================================================
  83
  84   TParam2ColumnIt getColumn( const TParam2ColumnMap* columnsMap,
  85                              const double            parameter )
  86   {
  87     TParam2ColumnIt u_col = columnsMap->upper_bound( parameter );
  88     if ( u_col != columnsMap->begin() )
  89       --u_col;
  90     return u_col; // return left column
  91   }
  92
  93   //================================================================================
  94   /*!
  95    * \brief Return nodes around given parameter and a ratio
  96    * \param column - node column
  97    * \param param - parameter
  98    * \param node1 - lower node
  99    * \param node2 - upper node
 100    * \retval double - ratio
 101    */
 102   //================================================================================
 103
 104   double getRAndNodes( const TNodeColumn*     column,
 105                        const double           param,
 106                        const SMDS_MeshNode* & node1,
 107                        const SMDS_MeshNode* & node2)
 108   {
 109     if ( param >= 1.0 || column->size() == 1) {
 110       node1 = node2 = column->back();
 111       return 0;
 112     }
 113
 114     int i = int( param * ( column->size() - 1 ));
 115     double u0 = double( i )/ double( column->size() - 1 );
 116     double r = ( param - u0 ) * ( column->size() - 1 );
 117
 118     node1 = (*column)[ i ];
 119     node2 = (*column)[ i + 1];
 120     return r;
 121   }
 122
 123   //================================================================================
 124   /*!
 125    * \brief Compute boundary parameters of face parts
 126     * \param nbParts - nb of parts to split columns into
 127     * \param columnsMap - node columns of the face to split
 128     * \param params - computed parameters
 129    */
 130   //================================================================================
 131
 132   void splitParams( const int               nbParts,
 133                     const TParam2ColumnMap* columnsMap,
 134                     vector< double > &      params)
 135   {
 136     params.clear();
 137     params.reserve( nbParts + 1 );
 138     TParam2ColumnIt last_par_col = --columnsMap->end();
 139     double par = columnsMap->begin()->first; // 0.
 140     double parLast = last_par_col->first;
 141     params.push_back( par );
 142     for ( int i = 0; i < nbParts - 1; ++ i )
 143     {
 144       double partSize = ( parLast - par ) / double ( nbParts - i );
 145       TParam2ColumnIt par_col = getColumn( columnsMap, par + partSize );
 146       if ( par_col->first == par ) {
 147         ++par_col;
 148         if ( par_col == last_par_col ) {
 149           while ( i < nbParts - 1 )
 150             params.push_back( par + partSize * i++ );
 151           break;
 152         }
 153       }
 154       par = par_col->first;
 155       params.push_back( par );
 156     }
 157     params.push_back( parLast ); // 1.
 158   }
 159
 160   //================================================================================
 161   /*!
 162    * \brief Return coordinate system for z-th layer of nodes
 163    */
 164   //================================================================================
 165
 166   gp_Ax2 getLayerCoordSys(const int                           z,
 167                           const vector< const TNodeColumn* >& columns,
 168                           int&                                xColumn)
 169   {
 170     // gravity center of a layer
 171     gp_XYZ O(0,0,0);
 172     int vertexCol = -1;
 173     for ( int i = 0; i < columns.size(); ++i )
 174     {
 175       O += gpXYZ( (*columns[ i ])[ z ]);
 176       if ( vertexCol < 0 &&
 177            columns[ i ]->front()->GetPosition()->GetTypeOfPosition() == SMDS_TOP_VERTEX )
 178         vertexCol = i;
 179     }
 180     O /= columns.size();
 181
 182     // Z axis
 183     gp_Vec Z(0,0,0);
 184     int iPrev = columns.size()-1;
 185     for ( int i = 0; i < columns.size(); ++i )
 186     {
 187       gp_Vec v1( O, gpXYZ( (*columns[ iPrev ])[ z ]));
 188       gp_Vec v2( O, gpXYZ( (*columns[ i ]    )[ z ]));
 189       Z += v1 ^ v2;
 190       iPrev = i;
 191     }
 192
 193     if ( vertexCol >= 0 )
 194     {
 195       O = gpXYZ( (*columns[ vertexCol ])[ z ]);
 196     }
 197     if ( xColumn < 0 || xColumn >= columns.size() )
 198     {
 199       // select a column for X dir
 200       double maxDist = 0;
 201       for ( int i = 0; i < columns.size(); ++i )
 202       {
 203         double dist = ( O - gpXYZ((*columns[ i ])[ z ])).SquareModulus();
 204         if ( dist > maxDist )
 205         {
 206           xColumn = i;
 207           maxDist = dist;
 208         }
 209       }
 210     }
 211
 212     // X axis
 213     gp_Vec X( O, gpXYZ( (*columns[ xColumn ])[ z ]));
 214
 215     return gp_Ax2( O, Z, X);
 216   }
 217 }
 218
 219 //=======================================================================
 220 //function : StdMeshers_Prism_3D
 221 //purpose  :
 222 //=======================================================================
 223
 224 StdMeshers_Prism_3D::StdMeshers_Prism_3D(int hypId, int studyId, SMESH_Gen* gen)
 225   :SMESH_3D_Algo(hypId, studyId, gen)
 226 {
 227   _name = "Prism_3D";
 228   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);       // 1 bit per shape type
 229   myProjectTriangles = false;
 230 }
 231
 232 //================================================================================
 233 /*!
 234  * \brief Destructor
 235  */
 236 //================================================================================
 237
 238 StdMeshers_Prism_3D::~StdMeshers_Prism_3D()
 239 {}
 240
 241 //=======================================================================
 242 //function : CheckHypothesis
 243 //purpose  :
 244 //=======================================================================
 245
 246 bool StdMeshers_Prism_3D::CheckHypothesis(SMESH_Mesh&                          aMesh,
 247                                           const TopoDS_Shape&                  aShape,
 248                                           SMESH_Hypothesis::Hypothesis_Status& aStatus)
 249 {
 250   // Check shape geometry
 251 /*  PAL16229
 252   aStatus = SMESH_Hypothesis::HYP_BAD_GEOMETRY;
 253
 254   // find not quadrangle faces
 255   list< TopoDS_Shape > notQuadFaces;
 256   int nbEdge, nbWire, nbFace = 0;
 257   TopExp_Explorer exp( aShape, TopAbs_FACE );
 258   for ( ; exp.More(); exp.Next() ) {
 259     ++nbFace;
 260     const TopoDS_Shape& face = exp.Current();
 261     nbEdge = TAssocTool::Count( face, TopAbs_EDGE, 0 );
 262     nbWire = TAssocTool::Count( face, TopAbs_WIRE, 0 );
 263     if (  nbEdge!= 4 || nbWire!= 1 ) {
 264       if ( !notQuadFaces.empty() ) {
 265         if ( TAssocTool::Count( notQuadFaces.back(), TopAbs_EDGE, 0 ) != nbEdge ||
 266              TAssocTool::Count( notQuadFaces.back(), TopAbs_WIRE, 0 ) != nbWire )
 267           RETURN_BAD_RESULT("Different not quad faces");
 268       }
 269       notQuadFaces.push_back( face );
 270     }
 271   }
 272   if ( !notQuadFaces.empty() )
 273   {
 274     if ( notQuadFaces.size() != 2 )
 275       RETURN_BAD_RESULT("Bad nb not quad faces: " << notQuadFaces.size());
 276
 277     // check total nb faces
 278     nbEdge = TAssocTool::Count( notQuadFaces.back(), TopAbs_EDGE, 0 );
 279     if ( nbFace != nbEdge + 2 )
 280       RETURN_BAD_RESULT("Bad nb of faces: " << nbFace << " but must be " << nbEdge + 2);
 281   }
 282 */
 283   // no hypothesis
 284   aStatus = SMESH_Hypothesis::HYP_OK;
 285   return true;
 286 }
 287
 288 //=======================================================================
 289 //function : Compute
 290 //purpose  :
 291 //=======================================================================
 292
 293 bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theShape)
 294 {
 295   SMESH_MesherHelper helper( theMesh );
 296   myHelper = &helper;
 297
 298   myHelper->IsQuadraticSubMesh( theShape );
 299
 300   // Analyse mesh and geomerty to find block subshapes and submeshes
 301   if ( !myBlock.Init( myHelper, theShape ))
 302     return error( myBlock.GetError());
 303
 304   SMESHDS_Mesh* meshDS = theMesh.GetMeshDS();
 305
 306   int volumeID = meshDS->ShapeToIndex( theShape );
 307
 308
 309   // To compute coordinates of a node inside a block, it is necessary to know
 310   // 1. normalized parameters of the node by which
 311   // 2. coordinates of node projections on all block sub-shapes are computed
 312
 313   // So we fill projections on vertices at once as they are same for all nodes
 314   myShapeXYZ.resize( myBlock.NbSubShapes() );
 315   for ( int iV = SMESH_Block::ID_FirstV; iV < SMESH_Block::ID_FirstE; ++iV ) {
 316     myBlock.VertexPoint( iV, myShapeXYZ[ iV ]);
 317     SHOWYXZ("V point " <<iV << " ", myShapeXYZ[ iV ]);
 318   }
 319
 320   // Projections on the top and bottom faces are taken from nodes existing
 321   // on these faces; find correspondence between bottom and top nodes
 322   myBotToColumnMap.clear();
 323   if ( !assocOrProjBottom2Top() ) // it also fill myBotToColumnMap
 324     return false;
 325
 326
 327   // Create nodes inside the block
 328
 329   // try to use transformation (issue 0020680)
 330   vector<gp_Trsf> trsf;
 331   if ( myBlock.GetLayersTransformation(trsf))
 332   {
 333     // loop on nodes inside the bottom face
 334     TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin();
 335     for ( ; bot_column != myBotToColumnMap.end(); ++bot_column )
 336     {
 337       const TNode& tBotNode = bot_column->first; // bottom TNode
 338       if ( tBotNode.GetPositionType() != SMDS_TOP_FACE )
 339         continue; // node is not inside face
 340
 341       // column nodes; middle part of the column are zero pointers
 342       TNodeColumn& column = bot_column->second;
 343       TNodeColumn::iterator columnNodes = column.begin();
 344       for ( int z = 0; columnNodes != column.end(); ++columnNodes, ++z)
 345       {
 346         const SMDS_MeshNode* & node = *columnNodes;
 347         if ( node ) continue; // skip bottom or top node
 348
 349         gp_XYZ coords = tBotNode.GetCoords();
 350         trsf[z-1].Transforms( coords );
 351         node = meshDS->AddNode( coords.X(), coords.Y(), coords.Z() );
 352         meshDS->SetNodeInVolume( node, volumeID );
 353       }
 354     } // loop on bottom nodes
 355   }
 356   else // use block approach
 357   {
 358     // loop on nodes inside the bottom face
 359     TNode prevBNode;
 360     TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin();
 361     for ( ; bot_column != myBotToColumnMap.end(); ++bot_column )
 362     {
 363       const TNode& tBotNode = bot_column->first; // bottom TNode
 364       if ( tBotNode.GetPositionType() != SMDS_TOP_FACE )
 365         continue; // node is not inside face
 366
 367       // column nodes; middle part of the column are zero pointers
 368       TNodeColumn& column = bot_column->second;
 369
 370       // compute bottom node parameters
 371       gp_XYZ paramHint(-1,-1,-1);
 372       if ( prevBNode.IsNeighbor( tBotNode ))
 373         paramHint = prevBNode.GetParams();
 374       if ( !myBlock.ComputeParameters( tBotNode.GetCoords(), tBotNode.ChangeParams(),
 375                                        ID_BOT_FACE, paramHint ))
 376         return error(TCom("Can't compute normalized parameters for node ")
 377                      << tBotNode.myNode->GetID() << " on the face #"
 378                      << myBlock.SubMesh( ID_BOT_FACE )->GetId() );
 379       prevBNode = tBotNode;
 380
 381       myShapeXYZ[ ID_BOT_FACE ] = tBotNode.GetCoords();
 382       gp_XYZ botParams          = tBotNode.GetParams();
 383
 384       // compute top node parameters
 385       myShapeXYZ[ ID_TOP_FACE ] = gpXYZ( column.back() );
 386       gp_XYZ topParams = botParams;
 387       topParams.SetZ( 1 );
 388       if ( column.size() > 2 ) {
 389         gp_Pnt topCoords = myShapeXYZ[ ID_TOP_FACE ];
 390         if ( !myBlock.ComputeParameters( topCoords, topParams, ID_TOP_FACE, topParams ))
 391           return error(TCom("Can't compute normalized parameters ")
 392                        << "for node " << column.back()->GetID()
 393                        << " on the face #"<< column.back()->GetPosition()->GetShapeId() );
 394       }
 395
 396       // vertical loop
 397       TNodeColumn::iterator columnNodes = column.begin();
 398       for ( int z = 0; columnNodes != column.end(); ++columnNodes, ++z)
 399       {
 400         const SMDS_MeshNode* & node = *columnNodes;
 401         if ( node ) continue; // skip bottom or top node
 402
 403         // params of a node to create
 404         double rz = (double) z / (double) ( column.size() - 1 );
 405         gp_XYZ params = botParams * ( 1 - rz ) + topParams * rz;
 406
 407         // set coords on all faces and nodes
 408         const int nbSideFaces = 4;
 409         int sideFaceIDs[nbSideFaces] = { SMESH_Block::ID_Fx0z,
 410                                          SMESH_Block::ID_Fx1z,
 411                                          SMESH_Block::ID_F0yz,
 412                                          SMESH_Block::ID_F1yz };
 413         for ( int iF = 0; iF < nbSideFaces; ++iF )
 414           if ( !setFaceAndEdgesXYZ( sideFaceIDs[ iF ], params, z ))
 415             return false;
 416
 417         // compute coords for a new node
 418         gp_XYZ coords;
 419         if ( !SMESH_Block::ShellPoint( params, myShapeXYZ, coords ))
 420           return error("Can't compute coordinates by normalized parameters");
 421
 422         SHOWYXZ("TOPFacePoint ",myShapeXYZ[ ID_TOP_FACE]);
 423         SHOWYXZ("BOT Node "<< tBotNode.myNode->GetID(),gpXYZ(tBotNode.myNode));
 424         SHOWYXZ("ShellPoint ",coords);
 425
 426         // create a node
 427         node = meshDS->AddNode( coords.X(), coords.Y(), coords.Z() );
 428         meshDS->SetNodeInVolume( node, volumeID );
 429       }
 430     } // loop on bottom nodes
 431   }
 432
 433   // Create volumes
 434
 435   SMESHDS_SubMesh* smDS = myBlock.SubMeshDS( ID_BOT_FACE );
 436   if ( !smDS ) return error(COMPERR_BAD_INPUT_MESH, "Null submesh");
 437
 438   // loop on bottom mesh faces
 439   SMDS_ElemIteratorPtr faceIt = smDS->GetElements();
 440   while ( faceIt->more() )
 441   {
 442     const SMDS_MeshElement* face = faceIt->next();
 443     if ( !face || face->GetType() != SMDSAbs_Face )
 444       continue;
 445     int nbNodes = face->NbNodes();
 446     if ( face->IsQuadratic() )
 447       nbNodes /= 2;
 448
 449     // find node columns for each node
 450     vector< const TNodeColumn* > columns( nbNodes );
 451     for ( int i = 0; i < nbNodes; ++i )
 452     {
 453       const SMDS_MeshNode* n = face->GetNode( i );
 454       if ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) {
 455         TNode2ColumnMap::iterator bot_column = myBotToColumnMap.find( n );
 456         if ( bot_column == myBotToColumnMap.end() )
 457           return error(TCom("No nodes found above node ") << n->GetID() );
 458         columns[ i ] = & bot_column->second;
 459       }
 460       else {
 461         columns[ i ] = myBlock.GetNodeColumn( n );
 462         if ( !columns[ i ] )
 463           return error(TCom("No side nodes found above node ") << n->GetID() );
 464       }
 465     }
 466     // create prisms
 467     AddPrisms( columns, myHelper );
 468
 469   } // loop on bottom mesh faces
 470
 471   return true;
 472 }
 473
 474
 475 //=======================================================================
 476 //function : Evaluate
 477 //purpose  :
 478 //=======================================================================
 479
 480 bool StdMeshers_Prism_3D::Evaluate(SMESH_Mesh& theMesh,
 481                                    const TopoDS_Shape& theShape,
 482                                    MapShapeNbElems& aResMap)
 483 {
 484   // find face contains only triangles
 485   vector < SMESH_subMesh * >meshFaces;
 486   TopTools_SequenceOfShape aFaces;
 487   int NumBase = 0, i = 0, NbQFs = 0;
 488   for (TopExp_Explorer exp(theShape, TopAbs_FACE); exp.More(); exp.Next()) {
 489     i++;
 490     aFaces.Append(exp.Current());
 491     SMESH_subMesh *aSubMesh = theMesh.GetSubMesh(exp.Current());
 492     meshFaces.push_back(aSubMesh);
 493     MapShapeNbElemsItr anIt = aResMap.find(meshFaces[i-1]);
 494     if( anIt==aResMap.end() ) {
 495       SMESH_ComputeErrorPtr& smError = aSubMesh->GetComputeError();
 496       smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
 497       return false;
 498     }
 499     std::vector<int> aVec = (*anIt).second;
 500     int nbtri = Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
 501     int nbqua = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
 502     if( nbtri==0 && nbqua>0 ) {
 503       NbQFs++;
 504     }
 505     if( nbtri>0 ) {
 506       NumBase = i;
 507     }
 508   }
 509
 510   if(NbQFs<4) {
 511     std::vector<int> aResVec(SMDSEntity_Last);
 512     for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
 513     SMESH_subMesh * sm = theMesh.GetSubMesh(theShape);
 514     aResMap.insert(std::make_pair(sm,aResVec));
 515     SMESH_ComputeErrorPtr& smError = sm->GetComputeError();
 516     smError.reset( new SMESH_ComputeError(COMPERR_ALGO_FAILED,"Submesh can not be evaluated",this));
 517     return false;
 518   }
 519
 520   if(NumBase==0) NumBase = 1; // only quads => set 1 faces as base
 521
 522   // find number of 1d elems for base face
 523   int nb1d = 0;
 524   TopTools_MapOfShape Edges1;
 525   for (TopExp_Explorer exp(aFaces.Value(NumBase), TopAbs_EDGE); exp.More(); exp.Next()) {
 526     Edges1.Add(exp.Current());
 527     SMESH_subMesh *sm = theMesh.GetSubMesh(exp.Current());
 528     if( sm ) {
 529       MapShapeNbElemsItr anIt = aResMap.find(sm);
 530       if( anIt == aResMap.end() ) continue;
 531       std::vector<int> aVec = (*anIt).second;
 532       nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
 533     }
 534   }
 535   // find face opposite to base face
 536   int OppNum = 0;
 537   for(i=1; i<=6; i++) {
 538     if(i==NumBase) continue;
 539     bool IsOpposite = true;
 540     for(TopExp_Explorer exp(aFaces.Value(i), TopAbs_EDGE); exp.More(); exp.Next()) {
 541       if( Edges1.Contains(exp.Current()) ) {
 542         IsOpposite = false;
 543         break;
 544       }
 545     }
 546     if(IsOpposite) {
 547       OppNum = i;
 548       break;
 549     }
 550   }
 551   // find number of 2d elems on side faces
 552   int nb2d = 0;
 553   for(i=1; i<=6; i++) {
 554     if( i==OppNum || i==NumBase ) continue;
 555     MapShapeNbElemsItr anIt = aResMap.find( meshFaces[i-1] );
 556     if( anIt == aResMap.end() ) continue;
 557     std::vector<int> aVec = (*anIt).second;
 558     nb2d += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
 559   }
 560
 561   MapShapeNbElemsItr anIt = aResMap.find( meshFaces[NumBase-1] );
 562   std::vector<int> aVec = (*anIt).second;
 563   bool IsQuadratic = (aVec[SMDSEntity_Quad_Triangle]>aVec[SMDSEntity_Triangle]) ||
 564                      (aVec[SMDSEntity_Quad_Quadrangle]>aVec[SMDSEntity_Quadrangle]);
 565   int nb2d_face0_3 = Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
 566   int nb2d_face0_4 = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
 567   int nb0d_face0 = aVec[SMDSEntity_Node];
 568   int nb1d_face0_int = ( nb2d_face0_3*3 + nb2d_face0_4*4 - nb1d ) / 2;
 569
 570   std::vector<int> aResVec(SMDSEntity_Last);
 571   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
 572   if(IsQuadratic) {
 573     aResVec[SMDSEntity_Quad_Penta] = nb2d_face0_3 * ( nb2d/nb1d );
 574     aResVec[SMDSEntity_Quad_Hexa] = nb2d_face0_4 * ( nb2d/nb1d );
 575     aResVec[SMDSEntity_Node] = nb0d_face0 * ( 2*nb2d/nb1d - 1 ) - nb1d_face0_int * nb2d/nb1d;
 576   }
 577   else {
 578     aResVec[SMDSEntity_Node] = nb0d_face0 * ( nb2d/nb1d - 1 );
 579     aResVec[SMDSEntity_Penta] = nb2d_face0_3 * ( nb2d/nb1d );
 580     aResVec[SMDSEntity_Hexa] = nb2d_face0_4 * ( nb2d/nb1d );
 581   }
 582   SMESH_subMesh * sm = theMesh.GetSubMesh(theShape);
 583   aResMap.insert(std::make_pair(sm,aResVec));
 584
 585   return true;
 586 }
 587
 588
 589 //================================================================================
 590 /*!
 591  * \brief Create prisms
 592  * \param columns - columns of nodes generated from nodes of a mesh face
 593  * \param helper - helper initialized by mesh and shape to add prisms to
 594  */
 595 //================================================================================
 596
 597 void StdMeshers_Prism_3D::AddPrisms( vector<const TNodeColumn*> & columns,
 598                                      SMESH_MesherHelper*          helper)
 599 {
 600   SMESHDS_Mesh * meshDS = helper->GetMeshDS();
 601   int shapeID = helper->GetSubShapeID();
 602
 603   int nbNodes = columns.size();
 604   int nbZ     = columns[0]->size();
 605   if ( nbZ < 2 ) return;
 606
 607   // find out orientation
 608   bool isForward = true;
 609   SMDS_VolumeTool vTool;
 610   int z = 1;
 611   switch ( nbNodes ) {
 612   case 3: {
 613     const SMDS_MeshNode* botNodes[3] = { (*columns[0])[z-1],
 614                                          (*columns[1])[z-1],
 615                                          (*columns[2])[z-1] };
 616     const SMDS_MeshNode* topNodes[3] = { (*columns[0])[z],
 617                                          (*columns[1])[z],
 618                                          (*columns[2])[z] };
 619     SMDS_VolumeOfNodes tmpVol ( botNodes[0], botNodes[1], botNodes[2],
 620                                 topNodes[0], topNodes[1], topNodes[2]);
 621     vTool.Set( &tmpVol );
 622     isForward  = vTool.IsForward();
 623     break;
 624   }
 625   case 4: {
 626     const SMDS_MeshNode* botNodes[4] = { (*columns[0])[z-1], (*columns[1])[z-1],
 627                                          (*columns[2])[z-1], (*columns[3])[z-1] };
 628     const SMDS_MeshNode* topNodes[4] = { (*columns[0])[z], (*columns[1])[z],
 629                                          (*columns[2])[z], (*columns[3])[z] };
 630     SMDS_VolumeOfNodes tmpVol ( botNodes[0], botNodes[1], botNodes[2], botNodes[3],
 631                                 topNodes[0], topNodes[1], topNodes[2], topNodes[3]);
 632     vTool.Set( &tmpVol );
 633     isForward  = vTool.IsForward();
 634     break;
 635   }
 636   }
 637
 638   // vertical loop on columns
 639   for ( z = 1; z < nbZ; ++z )
 640   {
 641     SMDS_MeshElement* vol = 0;
 642     switch ( nbNodes ) {
 643
 644     case 3: {
 645       const SMDS_MeshNode* botNodes[3] = { (*columns[0])[z-1],
 646                                            (*columns[1])[z-1],
 647                                            (*columns[2])[z-1] };
 648       const SMDS_MeshNode* topNodes[3] = { (*columns[0])[z],
 649                                            (*columns[1])[z],
 650                                            (*columns[2])[z] };
 651       if ( isForward )
 652         vol = helper->AddVolume( botNodes[0], botNodes[1], botNodes[2],
 653                                  topNodes[0], topNodes[1], topNodes[2]);
 654       else
 655         vol = helper->AddVolume( topNodes[0], topNodes[1], topNodes[2],
 656                                  botNodes[0], botNodes[1], botNodes[2]);
 657       break;
 658       }
 659     case 4: {
 660       const SMDS_MeshNode* botNodes[4] = { (*columns[0])[z-1], (*columns[1])[z-1],
 661                                            (*columns[2])[z-1], (*columns[3])[z-1] };
 662       const SMDS_MeshNode* topNodes[4] = { (*columns[0])[z], (*columns[1])[z],
 663                                            (*columns[2])[z], (*columns[3])[z] };
 664       if ( isForward )
 665         vol = helper->AddVolume( botNodes[0], botNodes[1], botNodes[2], botNodes[3],
 666                                  topNodes[0], topNodes[1], topNodes[2], topNodes[3]);
 667       else
 668         vol = helper->AddVolume( topNodes[0], topNodes[1], topNodes[2], topNodes[3],
 669                                  botNodes[0], botNodes[1], botNodes[2], botNodes[3]);
 670       break;
 671       }
 672     default:
 673       // polyhedron
 674       vector<const SMDS_MeshNode*> nodes( 2*nbNodes + 4*nbNodes);
 675       vector<int> quantities( 2 + nbNodes, 4 );
 676       quantities[0] = quantities[1] = nbNodes;
 677       columns.resize( nbNodes + 1 );
 678       columns[ nbNodes ] = columns[ 0 ];
 679       for ( int i = 0; i < nbNodes; ++i ) {
 680         nodes[ i         ] = (*columns[ i ])[z-1]; // bottom
 681         nodes[ i+nbNodes ] = (*columns[ i ])[z  ]; // top
 682         // side
 683         int di = 2*nbNodes + 4*i - 1;
 684         nodes[ di   ] = (*columns[i  ])[z-1];
 685         nodes[ di+1 ] = (*columns[i+1])[z-1];
 686         nodes[ di+2 ] = (*columns[i+1])[z  ];
 687         nodes[ di+3 ] = (*columns[i  ])[z  ];
 688       }
 689       vol = meshDS->AddPolyhedralVolume( nodes, quantities );
 690     }
 691     if ( vol && shapeID > 0 )
 692       meshDS->SetMeshElementOnShape( vol, shapeID );
 693   }
 694 }
 695
 696 //================================================================================
 697 /*!
 698  * \brief Find correspondence between bottom and top nodes
 699  *  If elements on the bottom and top faces are topologically different,
 700  *  and projection is possible and allowed, perform the projection
 701  *  \retval bool - is a success or not
 702  */
 703 //================================================================================
 704
 705 bool StdMeshers_Prism_3D::assocOrProjBottom2Top()
 706 {
 707   SMESH_subMesh * botSM = myBlock.SubMesh( ID_BOT_FACE );
 708   SMESH_subMesh * topSM = myBlock.SubMesh( ID_TOP_FACE );
 709
 710   SMESHDS_SubMesh * botSMDS = botSM->GetSubMeshDS();
 711   SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS();
 712
 713   if ( !botSMDS || botSMDS->NbElements() == 0 )
 714     return error(TCom("No elememts on face #") << botSM->GetId());
 715
 716   bool needProject = false;
 717   if ( !topSMDS ||
 718        botSMDS->NbElements() != topSMDS->NbElements() ||
 719        botSMDS->NbNodes()    != topSMDS->NbNodes())
 720   {
 721     if ( myBlock.HasNotQuadElemOnTop() )
 722       return error(TCom("Mesh on faces #") << botSM->GetId()
 723                    <<" and #"<< topSM->GetId() << " seems different" );
 724     needProject = true;
 725   }
 726
 727   if ( 0/*needProject && !myProjectTriangles*/ )
 728     return error(TCom("Mesh on faces #") << botSM->GetId()
 729                  <<" and #"<< topSM->GetId() << " seems different" );
 730   ///RETURN_BAD_RESULT("Need to project but not allowed");
 731
 732   if ( needProject )
 733   {
 734     return projectBottomToTop();
 735   }
 736
 737   TopoDS_Face botFace = TopoDS::Face( myBlock.Shape( ID_BOT_FACE ));
 738   TopoDS_Face topFace = TopoDS::Face( myBlock.Shape( ID_TOP_FACE ));
 739   // associate top and bottom faces
 740   TAssocTool::TShapeShapeMap shape2ShapeMap;
 741   if ( !TAssocTool::FindSubShapeAssociation( botFace, myBlock.Mesh(),
 742                                              topFace, myBlock.Mesh(),
 743                                              shape2ShapeMap) )
 744     return error(TCom("Topology of faces #") << botSM->GetId()
 745                  <<" and #"<< topSM->GetId() << " seems different" );
 746
 747   // Find matching nodes of top and bottom faces
 748   TNodeNodeMap n2nMap;
 749   if ( ! TAssocTool::FindMatchingNodesOnFaces( botFace, myBlock.Mesh(),
 750                                                topFace, myBlock.Mesh(),
 751                                                shape2ShapeMap, n2nMap ))
 752     return error(TCom("Mesh on faces #") << botSM->GetId()
 753                  <<" and #"<< topSM->GetId() << " seems different" );
 754
 755   // Fill myBotToColumnMap
 756
 757   int zSize = myBlock.VerticalSize();
 758   //TNode prevTNode;
 759   TNodeNodeMap::iterator bN_tN = n2nMap.begin();
 760   for ( ; bN_tN != n2nMap.end(); ++bN_tN )
 761   {
 762     const SMDS_MeshNode* botNode = bN_tN->first;
 763     const SMDS_MeshNode* topNode = bN_tN->second;
 764     if ( botNode->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE )
 765       continue; // wall columns are contained in myBlock
 766     // compute bottom node params
 767     TNode bN( botNode );
 768 //     if ( zSize > 2 ) {
 769 //       gp_XYZ paramHint(-1,-1,-1);
 770 //       if ( prevTNode.IsNeighbor( bN ))
 771 //         paramHint = prevTNode.GetParams();
 772 //       if ( !myBlock.ComputeParameters( bN.GetCoords(), bN.ChangeParams(),
 773 //                                        ID_BOT_FACE, paramHint ))
 774 //         return error(TCom("Can't compute normalized parameters for node ")
 775 //                      << botNode->GetID() << " on the face #"<< botSM->GetId() );
 776 //       prevTNode = bN;
 777 //     }
 778     // create node column
 779     TNode2ColumnMap::iterator bN_col =
 780       myBotToColumnMap.insert( make_pair ( bN, TNodeColumn() )).first;
 781     TNodeColumn & column = bN_col->second;
 782     column.resize( zSize );
 783     column.front() = botNode;
 784     column.back()  = topNode;
 785   }
 786   return true;
 787 }
 788
 789 //================================================================================
 790 /*!
 791  * \brief Remove quadrangles from the top face and
 792  * create triangles there by projection from the bottom
 793  * \retval bool - a success or not
 794  */
 795 //================================================================================
 796
 797 bool StdMeshers_Prism_3D::projectBottomToTop()
 798 {
 799   SMESH_subMesh * botSM = myBlock.SubMesh( ID_BOT_FACE );
 800   SMESH_subMesh * topSM = myBlock.SubMesh( ID_TOP_FACE );
 801
 802   SMESHDS_SubMesh * botSMDS = botSM->GetSubMeshDS();
 803   SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS();
 804
 805   if ( topSMDS )
 806     topSM->ComputeStateEngine( SMESH_subMesh::CLEAN );
 807
 808   SMESHDS_Mesh* meshDS = myBlock.MeshDS();
 809   int shapeID = myHelper->GetSubShapeID();
 810   int topFaceID = meshDS->ShapeToIndex( topSM->GetSubShape() );
 811
 812   // Fill myBotToColumnMap
 813
 814   int zSize = myBlock.VerticalSize();
 815   TNode prevTNode;
 816   SMDS_NodeIteratorPtr nIt = botSMDS->GetNodes();
 817   while ( nIt->more() )
 818   {
 819     const SMDS_MeshNode* botNode = nIt->next();
 820     if ( botNode->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE )
 821       continue; // strange
 822     // compute bottom node params
 823     TNode bN( botNode );
 824     gp_XYZ paramHint(-1,-1,-1);
 825     if ( prevTNode.IsNeighbor( bN ))
 826       paramHint = prevTNode.GetParams();
 827     if ( !myBlock.ComputeParameters( bN.GetCoords(), bN.ChangeParams(),
 828                                      ID_BOT_FACE, paramHint ))
 829       return error(TCom("Can't compute normalized parameters for node ")
 830                    << botNode->GetID() << " on the face #"<< botSM->GetId() );
 831     prevTNode = bN;
 832     // compute top node coords
 833     gp_XYZ topXYZ; gp_XY topUV;
 834     if ( !myBlock.FacePoint( ID_TOP_FACE, bN.GetParams(), topXYZ ) ||
 835          !myBlock.FaceUV   ( ID_TOP_FACE, bN.GetParams(), topUV ))
 836       return error(TCom("Can't compute coordinates "
 837                         "by normalized parameters on the face #")<< topSM->GetId() );
 838     SMDS_MeshNode * topNode = meshDS->AddNode( topXYZ.X(),topXYZ.Y(),topXYZ.Z() );
 839     meshDS->SetNodeOnFace( topNode, topFaceID, topUV.X(), topUV.Y() );
 840     // create node column
 841     TNode2ColumnMap::iterator bN_col =
 842       myBotToColumnMap.insert( make_pair ( bN, TNodeColumn() )).first;
 843     TNodeColumn & column = bN_col->second;
 844     column.resize( zSize );
 845     column.front() = botNode;
 846     column.back()  = topNode;
 847   }
 848
 849   // Create top faces
 850
 851   // loop on bottom mesh faces
 852   SMDS_ElemIteratorPtr faceIt = botSMDS->GetElements();
 853   while ( faceIt->more() )
 854   {
 855     const SMDS_MeshElement* face = faceIt->next();
 856     if ( !face || face->GetType() != SMDSAbs_Face )
 857       continue;
 858     int nbNodes = face->NbNodes();
 859     if ( face->IsQuadratic() )
 860       nbNodes /= 2;
 861
 862     // find top node in columns for each bottom node
 863     vector< const SMDS_MeshNode* > nodes( nbNodes );
 864     for ( int i = 0; i < nbNodes; ++i )
 865     {
 866       const SMDS_MeshNode* n = face->GetNode( nbNodes - i - 1 );
 867       if ( n->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) {
 868         TNode2ColumnMap::iterator bot_column = myBotToColumnMap.find( n );
 869         if ( bot_column == myBotToColumnMap.end() )
 870           return error(TCom("No nodes found above node ") << n->GetID() );
 871         nodes[ i ] = bot_column->second.back();
 872       }
 873       else {
 874         const TNodeColumn* column = myBlock.GetNodeColumn( n );
 875         if ( !column )
 876           return error(TCom("No side nodes found above node ") << n->GetID() );
 877         nodes[ i ] = column->back();
 878       }
 879     }
 880     // create a face, with reversed orientation
 881     SMDS_MeshElement* newFace = 0;
 882     switch ( nbNodes ) {
 883
 884     case 3: {
 885       newFace = myHelper->AddFace(nodes[0], nodes[1], nodes[2]);
 886       break;
 887       }
 888     case 4: {
 889       newFace = myHelper->AddFace( nodes[0], nodes[1], nodes[2], nodes[3] );
 890       break;
 891       }
 892     default:
 893       newFace = meshDS->AddPolygonalFace( nodes );
 894     }
 895     if ( newFace && shapeID > 0 )
 896       meshDS->SetMeshElementOnShape( newFace, shapeID );
 897   }
 898
 899   return true;
 900 }
 901
 902 //================================================================================
 903 /*!
 904  * \brief Set projection coordinates of a node to a face and it's subshapes
 905  * \param faceID - the face given by in-block ID
 906  * \param params - node normalized parameters
 907  * \retval bool - is a success
 908  */
 909 //================================================================================
 910
 911 bool StdMeshers_Prism_3D::setFaceAndEdgesXYZ( const int faceID, const gp_XYZ& params, int z )
 912 {
 913   // find base and top edges of the face
 914   enum { BASE = 0, TOP, LEFT, RIGHT };
 915   vector< int > edgeVec; // 0-base, 1-top
 916   SMESH_Block::GetFaceEdgesIDs( faceID, edgeVec );
 917
 918   myBlock.EdgePoint( edgeVec[ BASE ], params, myShapeXYZ[ edgeVec[ BASE ]]);
 919   myBlock.EdgePoint( edgeVec[ TOP  ], params, myShapeXYZ[ edgeVec[ TOP ]]);
 920
 921   SHOWYXZ("\nparams ", params);
 922   SHOWYXZ("TOP is " <<edgeVec[ TOP ], myShapeXYZ[ edgeVec[ TOP]]);
 923   SHOWYXZ("BASE is "<<edgeVec[ BASE], myShapeXYZ[ edgeVec[ BASE]]);
 924
 925   if ( faceID == SMESH_Block::ID_Fx0z || faceID == SMESH_Block::ID_Fx1z )
 926   {
 927     myBlock.EdgePoint( edgeVec[ LEFT ], params, myShapeXYZ[ edgeVec[ LEFT ]]);
 928     myBlock.EdgePoint( edgeVec[ RIGHT ], params, myShapeXYZ[ edgeVec[ RIGHT ]]);
 929
 930     SHOWYXZ("VER "<<edgeVec[ LEFT], myShapeXYZ[ edgeVec[ LEFT]]);
 931     SHOWYXZ("VER "<<edgeVec[ RIGHT], myShapeXYZ[ edgeVec[ RIGHT]]);
 932   }
 933   myBlock.FacePoint( faceID, params, myShapeXYZ[ faceID ]);
 934   SHOWYXZ("FacePoint "<<faceID, myShapeXYZ[ faceID]);
 935
 936   return true;
 937 }
 938
 939 //================================================================================
 940 /*!
 941  * \brief Return true if this node and other one belong to one face
 942  */
 943 //================================================================================
 944
 945 bool TNode::IsNeighbor( const TNode& other ) const
 946 {
 947   if ( !other.myNode || !myNode ) return false;
 948
 949   SMDS_ElemIteratorPtr fIt = other.myNode->GetInverseElementIterator(SMDSAbs_Face);
 950   while ( fIt->more() )
 951     if ( fIt->next()->GetNodeIndex( myNode ) >= 0 )
 952       return true;
 953   return false;
 954 }
 955
 956 //================================================================================
 957 /*!
 958  * \brief Constructor. Initialization is needed
 959  */
 960 //================================================================================
 961
 962 StdMeshers_PrismAsBlock::StdMeshers_PrismAsBlock()
 963 {
 964   mySide = 0;
 965 }
 966
 967 StdMeshers_PrismAsBlock::~StdMeshers_PrismAsBlock()
 968 {
 969   if ( mySide ) {
 970     delete mySide; mySide = 0;
 971   }
 972 }
 973
 974 //================================================================================
 975 /*!
 976  * \brief Initialization.
 977  * \param helper - helper loaded with mesh and 3D shape
 978  * \param shape3D - a closed shell or solid
 979  * \retval bool - false if a mesh or a shape are KO
 980  */
 981 //================================================================================
 982
 983 bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper,
 984                                    const TopoDS_Shape& shape3D)
 985 {
 986   if ( mySide ) {
 987     delete mySide; mySide = 0;
 988   }
 989   vector< TSideFace* > sideFaces( NB_WALL_FACES, 0 );
 990   vector< pair< double, double> > params ( NB_WALL_FACES );
 991   mySide = new TSideFace( sideFaces, params );
 992
 993   myHelper = helper;
 994   SMESHDS_Mesh* meshDS = myHelper->GetMeshDS();
 995
 996   SMESH_Block::init();
 997   myShapeIDMap.Clear();
 998   myShapeIndex2ColumnMap.clear();
 999
1000   int wallFaceIds[ NB_WALL_FACES ] = { // to walk around a block
1001     SMESH_Block::ID_Fx0z, SMESH_Block::ID_F1yz,
1002     SMESH_Block::ID_Fx1z, SMESH_Block::ID_F0yz
1003   };
1004
1005   myError = SMESH_ComputeError::New();
1006
1007   // -------------------------------------------------------------
1008   // Look for top and bottom faces: not quadrangle ones or meshed
1009   // with not quadrangle elements
1010   // -------------------------------------------------------------
1011
1012   list< SMESH_subMesh* > notQuadGeomSubMesh;
1013   list< SMESH_subMesh* > notQuadElemSubMesh;
1014   int nbFaces = 0;
1015   //
1016   SMESH_subMesh* mainSubMesh = myHelper->GetMesh()->GetSubMeshContaining( shape3D );
1017   if ( !mainSubMesh ) return error(COMPERR_BAD_INPUT_MESH,"Null submesh of shape3D");
1018
1019   // analyse face submeshes
1020   SMESH_subMeshIteratorPtr smIt = mainSubMesh->getDependsOnIterator(false,false);
1021   while ( smIt->more() )
1022   {
1023     SMESH_subMesh* sm = smIt->next();
1024     const TopoDS_Shape& face = sm->GetSubShape();
1025     if ( face.ShapeType() != TopAbs_FACE )
1026       continue;
1027     nbFaces++;
1028
1029     // is quadrangle face?
1030     list< TopoDS_Edge > orderedEdges;
1031     list< int >         nbEdgesInWires;
1032     TopoDS_Vertex       V000;
1033     int nbWires = GetOrderedEdges( TopoDS::Face( face ),
1034                                    V000, orderedEdges, nbEdgesInWires );
1035     if ( nbWires != 1 || nbEdgesInWires.front() != 4 )
1036       notQuadGeomSubMesh.push_back( sm );
1037
1038     // look for not quadrangle mesh elements
1039     if ( SMESHDS_SubMesh* smDS = sm->GetSubMeshDS() ) {
1040       bool hasNotQuad = false;
1041       SMDS_ElemIteratorPtr eIt = smDS->GetElements();
1042       while ( eIt->more() && !hasNotQuad ) {
1043         const SMDS_MeshElement* elem = eIt->next();
1044         if ( elem->GetType() == SMDSAbs_Face ) {
1045           int nbNodes = elem->NbNodes();
1046           if ( elem->IsQuadratic() )
1047             nbNodes /= 2;
1048           hasNotQuad = ( nbNodes != 4 );
1049         }
1050       }
1051       if ( hasNotQuad )
1052         notQuadElemSubMesh.push_back( sm );
1053     }
1054     else {
1055       return error(COMPERR_BAD_INPUT_MESH,TCom("Not meshed face #")<<sm->GetId());
1056     }
1057     // check if a quadrangle face is meshed with a quadranglar grid
1058     if ( notQuadGeomSubMesh.back() != sm &&
1059          notQuadElemSubMesh.back() != sm )
1060     {
1061       // count nb edges on face sides
1062       vector< int > nbEdges;
1063       nbEdges.reserve( nbEdgesInWires.front() );
1064       for ( list< TopoDS_Edge >::iterator edge = orderedEdges.begin();
1065             edge != orderedEdges.end(); ++edge )
1066       {
1067         if ( SMESHDS_SubMesh* smDS = meshDS->MeshElements( *edge ))
1068           nbEdges.push_back ( smDS->NbElements() );
1069         else
1070           nbEdges.push_back ( 0 );
1071       }
1072       int nbQuads = sm->GetSubMeshDS()->NbElements();
1073       if ( nbEdges[0] *  nbEdges[1] != nbQuads ||
1074            nbEdges[0] != nbEdges[2] ||
1075            nbEdges[1] != nbEdges[3] )
1076         notQuadElemSubMesh.push_back( sm );
1077     }
1078   }
1079
1080   // ----------------------------------------------------------------------
1081   // Analyse faces mesh and topology: choose the bottom submesh.
1082   // If there are not quadrangle geom faces, they are top and bottom ones.
1083   // Not quadrangle geom faces must be only on top and bottom.
1084   // ----------------------------------------------------------------------
1085
1086   SMESH_subMesh * botSM = 0;
1087   SMESH_subMesh * topSM = 0;
1088
1089   int nbNotQuad       = notQuadGeomSubMesh.size();
1090   int nbNotQuadMeshed = notQuadElemSubMesh.size();
1091   bool hasNotQuad = ( nbNotQuad || nbNotQuadMeshed );
1092
1093   // detect bad cases
1094   if ( nbNotQuad > 0 && nbNotQuad != 2 )
1095     return error(COMPERR_BAD_SHAPE,
1096                  TCom("More than 2 not quadrilateral faces: ")
1097                  <<nbNotQuad);
1098   if ( nbNotQuadMeshed > 2 )
1099     return error(COMPERR_BAD_INPUT_MESH,
1100                  TCom("More than 2 faces with not quadrangle elements: ")
1101                  <<nbNotQuadMeshed);
1102
1103   // get found submeshes
1104   if ( hasNotQuad )
1105   {
1106     if ( nbNotQuadMeshed > 0 ) botSM = notQuadElemSubMesh.front();
1107     else                       botSM = notQuadGeomSubMesh.front();
1108     if ( nbNotQuadMeshed > 1 ) topSM = notQuadElemSubMesh.back();
1109     else if ( nbNotQuad  > 1 ) topSM = notQuadGeomSubMesh.back();
1110   }
1111   // detect other bad cases
1112   if ( nbNotQuad == 2 && nbNotQuadMeshed > 0 ) {
1113     bool ok = false;
1114     if ( nbNotQuadMeshed == 1 )
1115       ok = ( find( notQuadGeomSubMesh.begin(),
1116                    notQuadGeomSubMesh.end(), botSM ) != notQuadGeomSubMesh.end() );
1117     else
1118       ok = ( notQuadGeomSubMesh == notQuadElemSubMesh );
1119     if ( !ok )
1120       return error(COMPERR_BAD_INPUT_MESH, "Side face meshed with not quadrangle elements");
1121   }
1122
1123   myNotQuadOnTop = ( nbNotQuadMeshed > 1 );
1124
1125   // ----------------------------------------------------------
1126
1127   if ( nbNotQuad == 0 ) // Standard block of 6 quadrangle faces ?
1128   {
1129     // SMESH_Block will perform geometry analysis, we need just to find 2
1130     // connected vertices on top and bottom
1131
1132     TopoDS_Vertex Vbot, Vtop;
1133     if ( nbNotQuadMeshed > 0 ) // Look for vertices
1134     {
1135       TopTools_IndexedMapOfShape edgeMap;
1136       TopExp::MapShapes( botSM->GetSubShape(), TopAbs_EDGE, edgeMap );
1137       // vertex 1 is any vertex of the bottom face
1138       Vbot = TopExp::FirstVertex( TopoDS::Edge( edgeMap( 1 )));
1139       // vertex 2 is end vertex of edge sharing Vbot and not belonging to the bottom face
1140       TopTools_ListIteratorOfListOfShape ancestIt = Mesh()->GetAncestors( Vbot );
1141       for ( ; Vtop.IsNull() && ancestIt.More(); ancestIt.Next() )
1142       {
1143         const TopoDS_Shape & ancestor = ancestIt.Value();
1144         if ( ancestor.ShapeType() == TopAbs_EDGE && !edgeMap.FindIndex( ancestor ))
1145         {
1146           TopoDS_Vertex V1, V2;
1147           TopExp::Vertices( TopoDS::Edge( ancestor ), V1, V2);
1148           if      ( Vbot.IsSame ( V1 )) Vtop = V2;
1149           else if ( Vbot.IsSame ( V2 )) Vtop = V1;
1150           // check that Vtop belongs to shape3D
1151           TopExp_Explorer exp( shape3D, TopAbs_VERTEX );
1152           for ( ; exp.More(); exp.Next() )
1153             if ( Vtop.IsSame( exp.Current() ))
1154               break;
1155           if ( !exp.More() )
1156             Vtop.Nullify();
1157         }
1158       }
1159     }
1160     // get shell from shape3D
1161     TopoDS_Shell shell;
1162     TopExp_Explorer exp( shape3D, TopAbs_SHELL );
1163     int nbShell = 0;
1164     for ( ; exp.More(); exp.Next(), ++nbShell )
1165       shell = TopoDS::Shell( exp.Current() );
1166 //     if ( nbShell != 1 )
1167 //       RETURN_BAD_RESULT("There must be 1 shell in the block");
1168
1169     // Load geometry in SMESH_Block
1170     if ( !SMESH_Block::FindBlockShapes( shell, Vbot, Vtop, myShapeIDMap )) {
1171       if ( !hasNotQuad )
1172         return error(COMPERR_BAD_SHAPE, "Can't detect top and bottom of a prism");
1173     }
1174     else {
1175       if ( !botSM ) botSM = Mesh()->GetSubMeshContaining( myShapeIDMap( ID_BOT_FACE ));
1176       if ( !topSM ) topSM = Mesh()->GetSubMeshContaining( myShapeIDMap( ID_TOP_FACE ));
1177     }
1178
1179   } // end  Standard block of 6 quadrangle faces
1180   // --------------------------------------------------------
1181
1182   // Here the top and bottom faces are found
1183   if ( nbNotQuadMeshed == 2 ) // roughly check correspondence of horiz meshes
1184   {
1185 //     SMESHDS_SubMesh* topSMDS = topSM->GetSubMeshDS();
1186 //     SMESHDS_SubMesh* botSMDS = botSM->GetSubMeshDS();
1187 //     if ( topSMDS->NbNodes() != botSMDS->NbNodes() ||
1188 //          topSMDS->NbElements() != botSMDS->NbElements() )
1189 //       RETURN_BAD_RESULT("Top mesh doesn't correspond to bottom one");
1190   }
1191
1192   // ---------------------------------------------------------
1193   // If there are not quadrangle geom faces, we emulate
1194   // a block of 6 quadrangle faces.
1195   // Load SMESH_Block with faces and edges geometry
1196   // ---------------------------------------------------------
1197
1198
1199   // find vertex 000 - the one with smallest coordinates (for easy DEBUG :-)
1200   TopoDS_Vertex V000;
1201   double minVal = DBL_MAX, minX, val;
1202   for ( TopExp_Explorer exp( botSM->GetSubShape(), TopAbs_VERTEX );
1203         exp.More(); exp.Next() )
1204   {
1205     const TopoDS_Vertex& v = TopoDS::Vertex( exp.Current() );
1206     gp_Pnt P = BRep_Tool::Pnt( v );
1207     val = P.X() + P.Y() + P.Z();
1208     if ( val < minVal || ( val == minVal && P.X() < minX )) {
1209       V000 = v;
1210       minVal = val;
1211       minX = P.X();
1212     }
1213   }
1214
1215   // Get ordered bottom edges
1216   list< TopoDS_Edge > orderedEdges;
1217   list< int >         nbVertexInWires;
1218   SMESH_Block::GetOrderedEdges( TopoDS::Face( botSM->GetSubShape().Reversed() ),
1219                                 V000, orderedEdges, nbVertexInWires );
1220 //   if ( nbVertexInWires.size() != 1 )
1221 //     RETURN_BAD_RESULT("Wrong prism geometry");
1222
1223   // Get Wall faces corresponding to the ordered bottom edges
1224   list< TopoDS_Face > wallFaces;
1225   if ( !GetWallFaces( Mesh(), shape3D, botSM->GetSubShape(), orderedEdges, wallFaces))
1226     return error(COMPERR_BAD_SHAPE, "Can't find side faces");
1227
1228   // Find columns of wall nodes and calculate edges' lengths
1229   // --------------------------------------------------------
1230
1231   myParam2ColumnMaps.clear();
1232   myParam2ColumnMaps.resize( orderedEdges.size() ); // total nb edges
1233
1234   int iE, nbEdges = nbVertexInWires.front(); // nb outer edges
1235   vector< double > edgeLength( nbEdges );
1236   map< double, int > len2edgeMap;
1237
1238   list< TopoDS_Edge >::iterator edgeIt = orderedEdges.begin();
1239   list< TopoDS_Face >::iterator faceIt = wallFaces.begin();
1240   for ( iE = 0; iE < nbEdges; ++edgeIt, ++faceIt )
1241   {
1242     TParam2ColumnMap & faceColumns = myParam2ColumnMaps[ iE ];
1243     if ( !myHelper->LoadNodeColumns( faceColumns, *faceIt, *edgeIt, meshDS ))
1244       return error(COMPERR_BAD_INPUT_MESH, TCom("Can't find regular quadrangle mesh ")
1245                    << "on a side face #" << MeshDS()->ShapeToIndex( *faceIt ));
1246
1247     SHOWYXZ("\np1 F "<<iE, gpXYZ(faceColumns.begin()->second.front() ));
1248     SHOWYXZ("p2 F "<<iE, gpXYZ(faceColumns.rbegin()->second.front() ));
1249     SHOWYXZ("V First "<<iE, BRep_Tool::Pnt( TopExp::FirstVertex(*edgeIt,true )));
1250
1251     edgeLength[ iE ] = SMESH_Algo::EdgeLength( *edgeIt );
1252
1253     if ( nbEdges < NB_WALL_FACES ) // fill map used to split faces
1254     {
1255       SMESHDS_SubMesh* smDS = meshDS->MeshElements( *edgeIt);
1256       if ( !smDS )
1257         return error(COMPERR_BAD_INPUT_MESH, TCom("Null submesh on the edge #")
1258                      << MeshDS()->ShapeToIndex( *edgeIt ));
1259       // assure length uniqueness
1260       edgeLength[ iE ] *= smDS->NbNodes() + edgeLength[ iE ] / ( 1000 + iE );
1261       len2edgeMap[ edgeLength[ iE ]] = iE;
1262     }
1263     ++iE;
1264   }
1265   // Load columns of internal edges (forming holes)
1266   // and fill map ShapeIndex to TParam2ColumnMap for them
1267   for ( ; edgeIt != orderedEdges.end() ; ++edgeIt, ++faceIt )
1268   {
1269     TParam2ColumnMap & faceColumns = myParam2ColumnMaps[ iE ];
1270     if ( !myHelper->LoadNodeColumns( faceColumns, *faceIt, *edgeIt, meshDS ))
1271       return error(COMPERR_BAD_INPUT_MESH, TCom("Can't find regular quadrangle mesh ")
1272                    << "on a side face #" << MeshDS()->ShapeToIndex( *faceIt ));
1273     // edge columns
1274     int id = MeshDS()->ShapeToIndex( *edgeIt );
1275     bool isForward = true; // meaningless for intenal wires
1276     myShapeIndex2ColumnMap[ id ] = make_pair( & faceColumns, isForward );
1277     // columns for vertices
1278     // 1
1279     const SMDS_MeshNode* n0 = faceColumns.begin()->second.front();
1280     id = n0->GetPosition()->GetShapeId();
1281     myShapeIndex2ColumnMap[ id ] = make_pair( & faceColumns, isForward );
1282     // 2
1283     const SMDS_MeshNode* n1 = faceColumns.rbegin()->second.front();
1284     id = n1->GetPosition()->GetShapeId();
1285     myShapeIndex2ColumnMap[ id ] = make_pair( & faceColumns, isForward );
1286 //     SHOWYXZ("\np1 F "<<iE, gpXYZ(faceColumns.begin()->second.front() ));
1287 //     SHOWYXZ("p2 F "<<iE, gpXYZ(faceColumns.rbegin()->second.front() ));
1288 //     SHOWYXZ("V First "<<iE, BRep_Tool::Pnt( TopExp::FirstVertex(*edgeIt,true )));
1289     ++iE;
1290   }
1291
1292   // Create 4 wall faces of a block
1293   // -------------------------------
1294
1295   if ( nbEdges <= NB_WALL_FACES ) // ************* Split faces if necessary
1296   {
1297     map< int, int > iE2nbSplit;
1298     if ( nbEdges != NB_WALL_FACES ) // define how to split
1299     {
1300       if ( len2edgeMap.size() != nbEdges )
1301         RETURN_BAD_RESULT("Uniqueness of edge lengths not assured");
1302       map< double, int >::reverse_iterator maxLen_i = len2edgeMap.rbegin();
1303       map< double, int >::reverse_iterator midLen_i = ++len2edgeMap.rbegin();
1304       double maxLen = maxLen_i->first;
1305       double midLen = ( len2edgeMap.size() == 1 ) ? 0 : midLen_i->first;
1306       switch ( nbEdges ) {
1307       case 1: // 0-th edge is split into 4 parts
1308         iE2nbSplit.insert( make_pair( 0, 4 )); break;
1309       case 2: // either the longest edge is split into 3 parts, or both edges into halves
1310         if ( maxLen / 3 > midLen / 2 ) {
1311           iE2nbSplit.insert( make_pair( maxLen_i->second, 3 ));
1312         }
1313         else {
1314           iE2nbSplit.insert( make_pair( maxLen_i->second, 2 ));
1315           iE2nbSplit.insert( make_pair( midLen_i->second, 2 ));
1316         }
1317         break;
1318       case 3:
1319         // split longest into halves
1320         iE2nbSplit.insert( make_pair( maxLen_i->second, 2 ));
1321       }
1322     }
1323     // Create TSideFace's
1324     faceIt = wallFaces.begin();
1325     edgeIt = orderedEdges.begin();
1326     int iSide = 0;
1327     for ( iE = 0; iE < nbEdges; ++edgeIt, ++faceIt )
1328     {
1329      // split?
1330       map< int, int >::iterator i_nb = iE2nbSplit.find( iE );
1331       if ( i_nb != iE2nbSplit.end() ) {
1332         // split!
1333         int nbSplit = i_nb->second;
1334         vector< double > params;
1335         splitParams( nbSplit, &myParam2ColumnMaps[ iE ], params );
1336         bool isForward = ( edgeIt->Orientation() == TopAbs_FORWARD );
1337         for ( int i = 0; i < nbSplit; ++i ) {
1338           double f = ( isForward ? params[ i ] : params[ nbSplit - i-1 ]);
1339           double l = ( isForward ? params[ i+1 ] : params[ nbSplit - i ]);
1340           TSideFace* comp = new TSideFace( myHelper, wallFaceIds[ iSide ],
1341                                            *faceIt, *edgeIt,
1342                                            &myParam2ColumnMaps[ iE ], f, l );
1343           mySide->SetComponent( iSide++, comp );
1344         }
1345       }
1346       else {
1347         TSideFace* comp = new TSideFace( myHelper, wallFaceIds[ iSide ],
1348                                          *faceIt, *edgeIt,
1349                                          &myParam2ColumnMaps[ iE ]);
1350         mySide->SetComponent( iSide++, comp );
1351       }
1352       ++iE;
1353     }
1354   }
1355   else { // **************************** Unite faces
1356
1357     // unite first faces
1358     int nbExraFaces = nbEdges - 3;
1359     int iSide = 0, iE;
1360     double u0 = 0, sumLen = 0;
1361     for ( iE = 0; iE < nbExraFaces; ++iE )
1362       sumLen += edgeLength[ iE ];
1363
1364     vector< TSideFace* > components( nbExraFaces );
1365     vector< pair< double, double> > params( nbExraFaces );
1366     faceIt = wallFaces.begin();
1367     edgeIt = orderedEdges.begin();
1368     for ( iE = 0; iE < nbExraFaces; ++edgeIt, ++faceIt )
1369     {
1370       components[ iE ] = new TSideFace( myHelper, wallFaceIds[ iSide ],
1371                                         *faceIt, *edgeIt,
1372                                         &myParam2ColumnMaps[ iE ]);
1373       double u1 = u0 + edgeLength[ iE ] / sumLen;
1374       params[ iE ] = make_pair( u0 , u1 );
1375       u0 = u1;
1376       ++iE;
1377     }
1378     mySide->SetComponent( iSide++, new TSideFace( components, params ));
1379
1380     // fill the rest faces
1381     for ( ; iE < nbEdges; ++faceIt, ++edgeIt )
1382     {
1383       TSideFace* comp = new TSideFace( myHelper, wallFaceIds[ iSide ],
1384                                        *faceIt, *edgeIt,
1385                                        &myParam2ColumnMaps[ iE ]);
1386       mySide->SetComponent( iSide++, comp );
1387       ++iE;
1388     }
1389   }
1390
1391
1392   // Fill geometry fields of SMESH_Block
1393   // ------------------------------------
1394
1395   TopoDS_Face botF = TopoDS::Face( botSM->GetSubShape() );
1396   TopoDS_Face topF = TopoDS::Face( topSM->GetSubShape() );
1397
1398   vector< int > botEdgeIdVec;
1399   SMESH_Block::GetFaceEdgesIDs( ID_BOT_FACE, botEdgeIdVec );
1400
1401   bool isForward[NB_WALL_FACES] = { true, true, true, true };
1402   Adaptor2d_Curve2d* botPcurves[NB_WALL_FACES];
1403   Adaptor2d_Curve2d* topPcurves[NB_WALL_FACES];
1404
1405   for ( int iF = 0; iF < NB_WALL_FACES; ++iF )
1406   {
1407     TSideFace * sideFace = mySide->GetComponent( iF );
1408     if ( !sideFace )
1409       RETURN_BAD_RESULT("NULL TSideFace");
1410     int fID = sideFace->FaceID();
1411
1412     // fill myShapeIDMap
1413     if ( sideFace->InsertSubShapes( myShapeIDMap ) != 8 &&
1414          !sideFace->IsComplex())
1415       MESSAGE( ": Warning : InsertSubShapes() < 8 on side " << iF );
1416
1417     // side faces geometry
1418     Adaptor2d_Curve2d* pcurves[NB_WALL_FACES];
1419     if ( !sideFace->GetPCurves( pcurves ))
1420       RETURN_BAD_RESULT("TSideFace::GetPCurves() failed");
1421
1422     SMESH_Block::TFace& tFace = myFace[ fID - ID_FirstF ];
1423     tFace.Set( fID, sideFace->Surface(), pcurves, isForward );
1424
1425     SHOWYXZ( endl<<"F "<< iF << " id " << fID << " FRW " << sideFace->IsForward(), sideFace->Value(0,0));
1426     // edges 3D geometry
1427     vector< int > edgeIdVec;
1428     SMESH_Block::GetFaceEdgesIDs( fID, edgeIdVec );
1429     for ( int isMax = 0; isMax < 2; ++isMax ) {
1430       {
1431         int eID = edgeIdVec[ isMax ];
1432         SMESH_Block::TEdge& tEdge = myEdge[ eID - ID_FirstE ];
1433         tEdge.Set( eID, sideFace->HorizCurve(isMax), true);
1434         SHOWYXZ(eID<<" HOR"<<isMax<<"(0)", sideFace->HorizCurve(isMax)->Value(0));
1435         SHOWYXZ(eID<<" HOR"<<isMax<<"(1)", sideFace->HorizCurve(isMax)->Value(1));
1436       }
1437       {
1438         int eID = edgeIdVec[ isMax+2 ];
1439         SMESH_Block::TEdge& tEdge = myEdge[ eID - ID_FirstE  ];
1440         tEdge.Set( eID, sideFace->VertiCurve(isMax), true);
1441         SHOWYXZ(eID<<" VER"<<isMax<<"(0)", sideFace->VertiCurve(isMax)->Value(0));
1442         SHOWYXZ(eID<<" VER"<<isMax<<"(1)", sideFace->VertiCurve(isMax)->Value(1));
1443
1444         // corner points
1445         vector< int > vertexIdVec;
1446         SMESH_Block::GetEdgeVertexIDs( eID, vertexIdVec );
1447         myPnt[ vertexIdVec[0] - ID_FirstV ] = tEdge.GetCurve()->Value(0).XYZ();
1448         myPnt[ vertexIdVec[1] - ID_FirstV ] = tEdge.GetCurve()->Value(1).XYZ();
1449       }
1450     }
1451     // pcurves on horizontal faces
1452     for ( iE = 0; iE < NB_WALL_FACES; ++iE ) {
1453       if ( edgeIdVec[ BOTTOM_EDGE ] == botEdgeIdVec[ iE ] ) {
1454         botPcurves[ iE ] = sideFace->HorizPCurve( false, botF );
1455         topPcurves[ iE ] = sideFace->HorizPCurve( true,  topF );
1456         break;
1457       }
1458     }
1459     //sideFace->dumpNodes( 4 ); // debug
1460   }
1461   // horizontal faces geometry
1462   {
1463     SMESH_Block::TFace& tFace = myFace[ ID_BOT_FACE - ID_FirstF ];
1464     tFace.Set( ID_BOT_FACE, new BRepAdaptor_Surface( botF ), botPcurves, isForward );
1465     SMESH_Block::Insert( botF, ID_BOT_FACE, myShapeIDMap );
1466   }
1467   {
1468     SMESH_Block::TFace& tFace = myFace[ ID_TOP_FACE - ID_FirstF ];
1469     tFace.Set( ID_TOP_FACE, new BRepAdaptor_Surface( topF ), topPcurves, isForward );
1470     SMESH_Block::Insert( topF, ID_TOP_FACE, myShapeIDMap );
1471   }
1472
1473   // Fill map ShapeIndex to TParam2ColumnMap
1474   // ----------------------------------------
1475
1476   list< TSideFace* > fList;
1477   list< TSideFace* >::iterator fListIt;
1478   fList.push_back( mySide );
1479   for ( fListIt = fList.begin(); fListIt != fList.end(); ++fListIt)
1480   {
1481     int nb = (*fListIt)->NbComponents();
1482     for ( int i = 0; i < nb; ++i ) {
1483       if ( TSideFace* comp = (*fListIt)->GetComponent( i ))
1484         fList.push_back( comp );
1485     }
1486     if ( TParam2ColumnMap* cols = (*fListIt)->GetColumns()) {
1487       // columns for a base edge
1488       int id = MeshDS()->ShapeToIndex( (*fListIt)->BaseEdge() );
1489       bool isForward = (*fListIt)->IsForward();
1490       myShapeIndex2ColumnMap[ id ] = make_pair( cols, isForward );
1491
1492       // columns for vertices
1493       const SMDS_MeshNode* n0 = cols->begin()->second.front();
1494       id = n0->GetPosition()->GetShapeId();
1495       myShapeIndex2ColumnMap[ id ] = make_pair( cols, isForward );
1496
1497       const SMDS_MeshNode* n1 = cols->rbegin()->second.front();
1498       id = n1->GetPosition()->GetShapeId();
1499       myShapeIndex2ColumnMap[ id ] = make_pair( cols, !isForward );
1500     }
1501   }
1502
1503 //   gp_XYZ testPar(0.25, 0.25, 0), testCoord;
1504 //   if ( !FacePoint( ID_BOT_FACE, testPar, testCoord ))
1505 //     RETURN_BAD_RESULT("TEST FacePoint() FAILED");
1506 //   SHOWYXZ("IN TEST PARAM" , testPar);
1507 //   SHOWYXZ("OUT TEST CORD" , testCoord);
1508 //   if ( !ComputeParameters( testCoord, testPar , ID_BOT_FACE))
1509 //     RETURN_BAD_RESULT("TEST ComputeParameters() FAILED");
1510 //   SHOWYXZ("OUT TEST PARAM" , testPar);
1511
1512   return true;
1513 }
1514
1515 //================================================================================
1516 /*!
1517  * \brief Return pointer to column of nodes
1518  * \param node - bottom node from which the returned column goes up
1519  * \retval const TNodeColumn* - the found column
1520  */
1521 //================================================================================
1522
1523 const TNodeColumn* StdMeshers_PrismAsBlock::GetNodeColumn(const SMDS_MeshNode* node) const
1524 {
1525   int sID = node->GetPosition()->GetShapeId();
1526
1527   map<int, pair< TParam2ColumnMap*, bool > >::const_iterator col_frw =
1528     myShapeIndex2ColumnMap.find( sID );
1529   if ( col_frw != myShapeIndex2ColumnMap.end() ) {
1530     const TParam2ColumnMap* cols = col_frw->second.first;
1531     TParam2ColumnIt u_col = cols->begin();
1532     for ( ; u_col != cols->end(); ++u_col )
1533       if ( u_col->second[ 0 ] == node )
1534         return & u_col->second;
1535   }
1536   return 0;
1537 }
1538
1539 //=======================================================================
1540 //function : GetLayersTransformation
1541 //purpose  : Return transformations to get coordinates of nodes of each layer
1542 //           by nodes of the bottom. Layer is a set of nodes at a certain step
1543 //           from bottom to top.
1544 //=======================================================================
1545
1546 bool StdMeshers_PrismAsBlock::GetLayersTransformation(vector<gp_Trsf> & trsf) const
1547 {
1548   const int zSize = VerticalSize();
1549   if ( zSize < 3 ) return true;
1550   trsf.resize( zSize - 2 );
1551
1552   // Select some node columns by which we will define coordinate system of layers
1553
1554   vector< const TNodeColumn* > columns;
1555   {
1556     const TopoDS_Shape& baseFace = Shape(ID_BOT_FACE);
1557     list< TopoDS_Edge > orderedEdges;
1558     list< int >         nbEdgesInWires;
1559     GetOrderedEdges( TopoDS::Face( baseFace ), TopoDS_Vertex(), orderedEdges, nbEdgesInWires );
1560     bool isReverse;
1561     list< TopoDS_Edge >::iterator edgeIt = orderedEdges.begin();
1562     for ( int iE = 0; iE < nbEdgesInWires.front(); ++iE, ++edgeIt )
1563     {
1564       const TParam2ColumnMap& u2colMap =
1565         GetParam2ColumnMap( myHelper->GetMeshDS()->ShapeToIndex( *edgeIt ), isReverse );
1566       isReverse = ( edgeIt->Orientation() == TopAbs_REVERSED );
1567       double f = u2colMap.begin()->first, l = u2colMap.rbegin()->first;
1568       if ( isReverse ) swap ( f, l );
1569       const int nbCol = 5;
1570       for ( int i = 0; i < nbCol; ++i )
1571       {
1572         double u = f + i/double(nbCol) * ( l - f );
1573         const TNodeColumn* col = & getColumn( & u2colMap, u )->second;
1574         if ( columns.empty() || col != columns.back() )
1575           columns.push_back( col );
1576       }
1577     }
1578   }
1579
1580   // Find tolerance to check transformations
1581
1582   double tol2;
1583   {
1584     Bnd_B3d bndBox;
1585     for ( int i = 0; i < columns.size(); ++i )
1586       bndBox.Add( gpXYZ( columns[i]->front() ));
1587     tol2 = bndBox.SquareExtent() * 4 * 1e-4;
1588   }
1589
1590   // Compute transformations
1591
1592   int xCol = -1;
1593   gp_Trsf fromCsZ, toCs0;
1594   gp_Ax3 cs0 = getLayerCoordSys(0, columns, xCol );
1595   //double dist0 = cs0.Location().Distance( gpXYZ( (*columns[0])[0]));
1596   toCs0.SetTransformation( cs0 );
1597   for ( int z = 1; z < zSize-1; ++z )
1598   {
1599     gp_Ax3 csZ = getLayerCoordSys(z, columns, xCol );
1600     //double distZ = csZ.Location().Distance( gpXYZ( (*columns[0])[z]));
1601     fromCsZ.SetTransformation( csZ );
1602     fromCsZ.Invert();
1603     gp_Trsf& t = trsf[ z-1 ];
1604     t = fromCsZ * toCs0;
1605     //t.SetScaleFactor( distZ/dist0 ); - it does not work properly, wrong base point
1606
1607     // check a transformation
1608     for ( int i = 0; i < columns.size(); ++i )
1609     {
1610       gp_Pnt p0 = gpXYZ( (*columns[i])[0] );
1611       gp_Pnt pz = gpXYZ( (*columns[i])[z] );
1612       t.Transforms( p0.ChangeCoord() );
1613       if ( p0.SquareDistance( pz ) > tol2 )
1614         return false;
1615     }
1616   }
1617   return true;
1618 }
1619
1620 //================================================================================
1621 /*!
1622  * \brief Check curve orientation of a bootom edge
1623   * \param meshDS - mesh DS
1624   * \param columnsMap - node columns map of side face
1625   * \param bottomEdge - the bootom edge
1626   * \param sideFaceID - side face in-block ID
1627   * \retval bool - true if orientation coinside with in-block froward orientation
1628  */
1629 //================================================================================
1630
1631 bool StdMeshers_PrismAsBlock::IsForwardEdge(SMESHDS_Mesh*           meshDS,
1632                                             const TParam2ColumnMap& columnsMap,
1633                                             const TopoDS_Edge &     bottomEdge,
1634                                             const int               sideFaceID)
1635 {
1636   bool isForward = false;
1637   if ( TAssocTool::IsClosedEdge( bottomEdge ))
1638   {
1639     isForward = ( bottomEdge.Orientation() == TopAbs_FORWARD );
1640   }
1641   else
1642   {
1643     const TNodeColumn& firstCol = columnsMap.begin()->second;
1644     const SMDS_MeshNode* bottomNode = firstCol[0];
1645     TopoDS_Shape firstVertex = SMESH_MesherHelper::GetSubShapeByNode( bottomNode, meshDS );
1646     isForward = ( firstVertex.IsSame( TopExp::FirstVertex( bottomEdge, true )));
1647   }
1648   // on 2 of 4 sides first vertex is end
1649   if ( sideFaceID == ID_Fx1z || sideFaceID == ID_F0yz )
1650     isForward = !isForward;
1651   return isForward;
1652 }
1653
1654 //================================================================================
1655 /*!
1656  * \brief Find wall faces by bottom edges
1657  * \param mesh - the mesh
1658  * \param mainShape - the prism
1659  * \param bottomFace - the bottom face
1660  * \param bottomEdges - edges bounding the bottom face
1661  * \param wallFaces - faces list to fill in
1662  */
1663 //================================================================================
1664
1665 bool StdMeshers_PrismAsBlock::GetWallFaces( SMESH_Mesh*                     mesh,
1666                                             const TopoDS_Shape &            mainShape,
1667                                             const TopoDS_Shape &            bottomFace,
1668                                             const std::list< TopoDS_Edge >& bottomEdges,
1669                                             std::list< TopoDS_Face >&       wallFaces)
1670 {
1671   wallFaces.clear();
1672
1673   TopTools_IndexedMapOfShape faceMap;
1674   TopExp::MapShapes( mainShape, TopAbs_FACE, faceMap );
1675
1676   list< TopoDS_Edge >::const_iterator edge = bottomEdges.begin();
1677   for ( ; edge != bottomEdges.end(); ++edge )
1678   {
1679     TopTools_ListIteratorOfListOfShape ancestIt = mesh->GetAncestors( *edge );
1680     for ( ; ancestIt.More(); ancestIt.Next() )
1681     {
1682       const TopoDS_Shape& ancestor = ancestIt.Value();
1683       if ( ancestor.ShapeType() == TopAbs_FACE && // face
1684            !bottomFace.IsSame( ancestor ) &&      // not bottom
1685            faceMap.FindIndex( ancestor ))         // belongs to the prism
1686       {
1687         wallFaces.push_back( TopoDS::Face( ancestor ));
1688         break;
1689       }
1690     }
1691   }
1692   return ( wallFaces.size() == bottomEdges.size() );
1693 }
1694
1695 //================================================================================
1696 /*!
1697  * \brief Constructor
1698   * \param faceID - in-block ID
1699   * \param face - geom face
1700   * \param columnsMap - map of node columns
1701   * \param first - first normalized param
1702   * \param last - last normalized param
1703  */
1704 //================================================================================
1705
1706 StdMeshers_PrismAsBlock::TSideFace::TSideFace(SMESH_MesherHelper* helper,
1707                                               const int           faceID,
1708                                               const TopoDS_Face&  face,
1709                                               const TopoDS_Edge&  baseEdge,
1710                                               TParam2ColumnMap*   columnsMap,
1711                                               const double        first,
1712                                               const double        last):
1713   myID( faceID ),
1714   myParamToColumnMap( columnsMap ),
1715   myBaseEdge( baseEdge ),
1716   myHelper( helper )
1717 {
1718   mySurface.Initialize( face );
1719   myParams.resize( 1 );
1720   myParams[ 0 ] = make_pair( first, last );
1721   myIsForward = StdMeshers_PrismAsBlock::IsForwardEdge( myHelper->GetMeshDS(),
1722                                                         *myParamToColumnMap,
1723                                                         myBaseEdge, myID );
1724 }
1725
1726 //================================================================================
1727 /*!
1728  * \brief Constructor of complex side face
1729  */
1730 //================================================================================
1731
1732 StdMeshers_PrismAsBlock::TSideFace::
1733 TSideFace(const vector< TSideFace* >&             components,
1734           const vector< pair< double, double> > & params)
1735   :myID( components[0] ? components[0]->myID : 0 ),
1736    myParamToColumnMap( 0 ),
1737    myParams( params ),
1738    myIsForward( true ),
1739    myComponents( components ),
1740    myHelper( components[0] ? components[0]->myHelper : 0 )
1741 {}
1742 //================================================================================
1743 /*!
1744  * \brief Copy constructor
1745   * \param other - other side
1746  */
1747 //================================================================================
1748
1749 StdMeshers_PrismAsBlock::TSideFace::TSideFace( const TSideFace& other )
1750 {
1751   myID               = other.myID;
1752   mySurface          = other.mySurface;
1753   myBaseEdge         = other.myBaseEdge;
1754   myParams           = other.myParams;
1755   myIsForward        = other.myIsForward;
1756   myHelper           = other.myHelper;
1757   myParamToColumnMap = other.myParamToColumnMap;
1758
1759   myComponents.resize( other.myComponents.size());
1760   for (int i = 0 ; i < myComponents.size(); ++i )
1761     myComponents[ i ] = new TSideFace( *other.myComponents[ i ]);
1762 }
1763
1764 //================================================================================
1765 /*!
1766  * \brief Deletes myComponents
1767  */
1768 //================================================================================
1769
1770 StdMeshers_PrismAsBlock::TSideFace::~TSideFace()
1771 {
1772   for (int i = 0 ; i < myComponents.size(); ++i )
1773     if ( myComponents[ i ] )
1774       delete myComponents[ i ];
1775 }
1776
1777 //================================================================================
1778 /*!
1779  * \brief Return geometry of the vertical curve
1780   * \param isMax - true means curve located closer to (1,1,1) block point
1781   * \retval Adaptor3d_Curve* - curve adaptor
1782  */
1783 //================================================================================
1784
1785 Adaptor3d_Curve* StdMeshers_PrismAsBlock::TSideFace::VertiCurve(const bool isMax) const
1786 {
1787   if ( !myComponents.empty() ) {
1788     if ( isMax )
1789       return myComponents.back()->VertiCurve(isMax);
1790     else
1791       return myComponents.front()->VertiCurve(isMax);
1792   }
1793   double f = myParams[0].first, l = myParams[0].second;
1794   if ( !myIsForward ) std::swap( f, l );
1795   return new TVerticalEdgeAdaptor( myParamToColumnMap, isMax ? l : f );
1796 }
1797
1798 //================================================================================
1799 /*!
1800  * \brief Return geometry of the top or bottom curve
1801   * \param isTop -
1802   * \retval Adaptor3d_Curve* -
1803  */
1804 //================================================================================
1805
1806 Adaptor3d_Curve* StdMeshers_PrismAsBlock::TSideFace::HorizCurve(const bool isTop) const
1807 {
1808   return new THorizontalEdgeAdaptor( this, isTop );
1809 }
1810
1811 //================================================================================
1812 /*!
1813  * \brief Return pcurves
1814   * \param pcurv - array of 4 pcurves
1815   * \retval bool - is a success
1816  */
1817 //================================================================================
1818
1819 bool StdMeshers_PrismAsBlock::TSideFace::GetPCurves(Adaptor2d_Curve2d* pcurv[4]) const
1820 {
1821   int iEdge[ 4 ] = { BOTTOM_EDGE, TOP_EDGE, V0_EDGE, V1_EDGE };
1822
1823   for ( int i = 0 ; i < 4 ; ++i ) {
1824     Handle(Geom2d_Line) line;
1825     switch ( iEdge[ i ] ) {
1826     case TOP_EDGE:
1827       line = new Geom2d_Line( gp_Pnt2d( 0, 1 ), gp::DX2d() ); break;
1828     case BOTTOM_EDGE:
1829       line = new Geom2d_Line( gp::Origin2d(), gp::DX2d() ); break;
1830     case V0_EDGE:
1831       line = new Geom2d_Line( gp::Origin2d(), gp::DY2d() ); break;
1832     case V1_EDGE:
1833       line = new Geom2d_Line( gp_Pnt2d( 1, 0 ), gp::DY2d() ); break;
1834     }
1835     pcurv[ i ] = new Geom2dAdaptor_Curve( line, 0, 1 );
1836   }
1837   return true;
1838 }
1839
1840 //================================================================================
1841 /*!
1842  * \brief Returns geometry of pcurve on a horizontal face
1843   * \param isTop - is top or bottom face
1844   * \param horFace - a horizontal face
1845   * \retval Adaptor2d_Curve2d* - curve adaptor
1846  */
1847 //================================================================================
1848
1849 Adaptor2d_Curve2d*
1850 StdMeshers_PrismAsBlock::TSideFace::HorizPCurve(const bool         isTop,
1851                                                 const TopoDS_Face& horFace) const
1852 {
1853   return new TPCurveOnHorFaceAdaptor( this, isTop, horFace );
1854 }
1855
1856 //================================================================================
1857 /*!
1858  * \brief Return a component corresponding to parameter
1859   * \param U - parameter along a horizontal size
1860   * \param localU - parameter along a horizontal size of a component
1861   * \retval TSideFace* - found component
1862  */
1863 //================================================================================
1864
1865 StdMeshers_PrismAsBlock::TSideFace*
1866 StdMeshers_PrismAsBlock::TSideFace::GetComponent(const double U,double & localU) const
1867 {
1868   localU = U;
1869   if ( myComponents.empty() )
1870     return const_cast<TSideFace*>( this );
1871
1872   int i;
1873   for ( i = 0; i < myComponents.size(); ++i )
1874     if ( U < myParams[ i ].second )
1875       break;
1876   if ( i >= myComponents.size() )
1877     i = myComponents.size() - 1;
1878
1879   double f = myParams[ i ].first, l = myParams[ i ].second;
1880   localU = ( U - f ) / ( l - f );
1881   return myComponents[ i ];
1882 }
1883
1884 //================================================================================
1885 /*!
1886  * \brief Find node columns for a parameter
1887   * \param U - parameter along a horizontal edge
1888   * \param col1 - the 1st found column
1889   * \param col2 - the 2nd found column
1890   * \retval r - normalized position of U between the found columns
1891  */
1892 //================================================================================
1893
1894 double StdMeshers_PrismAsBlock::TSideFace::GetColumns(const double      U,
1895                                                       TParam2ColumnIt & col1,
1896                                                       TParam2ColumnIt & col2) const
1897 {
1898   double u = U, r = 0;
1899   if ( !myComponents.empty() ) {
1900     TSideFace * comp = GetComponent(U,u);
1901     return comp->GetColumns( u, col1, col2 );
1902   }
1903
1904   if ( !myIsForward )
1905     u = 1 - u;
1906   double f = myParams[0].first, l = myParams[0].second;
1907   u = f + u * ( l - f );
1908
1909   col1 = col2 = getColumn( myParamToColumnMap, u );
1910   if ( ++col2 == myParamToColumnMap->end() ) {
1911     --col2;
1912     r = 0.5;
1913   }
1914   else {
1915     double uf = col1->first;
1916     double ul = col2->first;
1917     r = ( u - uf ) / ( ul - uf );
1918   }
1919   return r;
1920 }
1921
1922 //================================================================================
1923 /*!
1924  * \brief Return coordinates by normalized params
1925   * \param U - horizontal param
1926   * \param V - vertical param
1927   * \retval gp_Pnt - result point
1928  */
1929 //================================================================================
1930
1931 gp_Pnt StdMeshers_PrismAsBlock::TSideFace::Value(const Standard_Real U,
1932                                                  const Standard_Real V) const
1933 {
1934   if ( !myComponents.empty() ) {
1935     double u;
1936     TSideFace * comp = GetComponent(U,u);
1937     return comp->Value( u, V );
1938   }
1939
1940   TParam2ColumnIt u_col1, u_col2;
1941   double vR, hR = GetColumns( U, u_col1, u_col2 );
1942
1943   const SMDS_MeshNode* n1 = 0;
1944   const SMDS_MeshNode* n2 = 0;
1945   const SMDS_MeshNode* n3 = 0;
1946   const SMDS_MeshNode* n4 = 0;
1947
1948   // BEGIN issue 0020680: EDF 1252 SMESH: Bad cell created by Radial prism in center of torus
1949   // Workaround for a wrongly located point returned by mySurface.Value() for
1950   // UV located near boundary of BSpline surface.
1951   // To bypass the problem, we take point from 3D curve of edge.
1952   // It solves pb of the bloc_fiss_new.py
1953   const double tol = 1e-3;
1954   if ( V < tol || V+tol >= 1. )
1955   {
1956     n1 = V < tol ? u_col1->second.front() : u_col1->second.back();
1957     n3 = V < tol ? u_col2->second.front() : u_col2->second.back();
1958     TopoDS_Edge edge;
1959     if ( V < tol )
1960     {
1961       edge = myBaseEdge;
1962     }
1963     else
1964     {
1965       TopoDS_Shape s = myHelper->GetSubShapeByNode( n1, myHelper->GetMeshDS() );
1966       if ( s.ShapeType() != TopAbs_EDGE )
1967         s = myHelper->GetSubShapeByNode( n3, myHelper->GetMeshDS() );
1968       if ( s.ShapeType() == TopAbs_EDGE )
1969         edge = TopoDS::Edge( s );
1970     }
1971     if ( !edge.IsNull() )
1972     {
1973       double u1 = myHelper->GetNodeU( edge, n1 );
1974       double u3 = myHelper->GetNodeU( edge, n3 );
1975       double u = u1 * ( 1 - hR ) + u3 * hR;
1976       TopLoc_Location loc; double f,l;
1977       Handle(Geom_Curve) curve = BRep_Tool::Curve( edge,loc,f,l );
1978       return curve->Value( u ).Transformed( loc );
1979     }
1980   }
1981   // END issue 0020680: EDF 1252 SMESH: Bad cell created by Radial prism in center of torus
1982
1983   vR = getRAndNodes( & u_col1->second, V, n1, n2 );
1984   vR = getRAndNodes( & u_col2->second, V, n3, n4 );
1985
1986   gp_XY uv1 = myHelper->GetNodeUV( mySurface.Face(), n1, n4);
1987   gp_XY uv2 = myHelper->GetNodeUV( mySurface.Face(), n2, n3);
1988   gp_XY uv12 = uv1 * ( 1 - vR ) + uv2 * vR;
1989
1990   gp_XY uv3 = myHelper->GetNodeUV( mySurface.Face(), n3, n2);
1991   gp_XY uv4 = myHelper->GetNodeUV( mySurface.Face(), n4, n1);
1992   gp_XY uv34 = uv3 * ( 1 - vR ) + uv4 * vR;
1993
1994   gp_XY uv = uv12 * ( 1 - hR ) + uv34 * hR;
1995
1996   gp_Pnt p = mySurface.Value( uv.X(), uv.Y() );
1997   return p;
1998 }
1999
2000
2001 //================================================================================
2002 /*!
2003  * \brief Return boundary edge
2004   * \param edge - edge index
2005   * \retval TopoDS_Edge - found edge
2006  */
2007 //================================================================================
2008
2009 TopoDS_Edge StdMeshers_PrismAsBlock::TSideFace::GetEdge(const int iEdge) const
2010 {
2011   if ( !myComponents.empty() ) {
2012     switch ( iEdge ) {
2013     case V0_EDGE : return myComponents.front()->GetEdge( iEdge );
2014     case V1_EDGE : return myComponents.back() ->GetEdge( iEdge );
2015     default: return TopoDS_Edge();
2016     }
2017   }
2018   TopoDS_Shape edge;
2019   const SMDS_MeshNode* node = 0;
2020   SMESHDS_Mesh * meshDS = myHelper->GetMesh()->GetMeshDS();
2021   TNodeColumn* column;
2022
2023   switch ( iEdge ) {
2024   case TOP_EDGE:
2025   case BOTTOM_EDGE:
2026     column = & (( ++myParamToColumnMap->begin())->second );
2027     node = ( iEdge == TOP_EDGE ) ? column->back() : column->front();
2028     edge = myHelper->GetSubShapeByNode ( node, meshDS );
2029     if ( edge.ShapeType() == TopAbs_VERTEX ) {
2030       column = & ( myParamToColumnMap->begin()->second );
2031       node = ( iEdge == TOP_EDGE ) ? column->back() : column->front();
2032     }
2033     break;
2034   case V0_EDGE:
2035   case V1_EDGE: {
2036     bool back = ( iEdge == V1_EDGE );
2037     if ( !myIsForward ) back = !back;
2038     if ( back )
2039       column = & ( myParamToColumnMap->rbegin()->second );
2040     else
2041       column = & ( myParamToColumnMap->begin()->second );
2042     if ( column->size() > 0 )
2043       edge = myHelper->GetSubShapeByNode( (*column)[ 1 ], meshDS );
2044     if ( edge.IsNull() || edge.ShapeType() == TopAbs_VERTEX )
2045       node = column->front();
2046     break;
2047   }
2048   default:;
2049   }
2050   if ( !edge.IsNull() && edge.ShapeType() == TopAbs_EDGE )
2051     return TopoDS::Edge( edge );
2052
2053   // find edge by 2 vertices
2054   TopoDS_Shape V1 = edge;
2055   TopoDS_Shape V2 = myHelper->GetSubShapeByNode( node, meshDS );
2056   if ( V2.ShapeType() == TopAbs_VERTEX && !V2.IsSame( V1 ))
2057   {
2058     TopTools_ListIteratorOfListOfShape ancestIt =
2059       myHelper->GetMesh()->GetAncestors( V1 );
2060     for ( ; ancestIt.More(); ancestIt.Next() )
2061     {
2062       const TopoDS_Shape & ancestor = ancestIt.Value();
2063       if ( ancestor.ShapeType() == TopAbs_EDGE )
2064         for ( TopExp_Explorer e( ancestor, TopAbs_VERTEX ); e.More(); e.Next() )
2065           if ( V2.IsSame( e.Current() ))
2066             return TopoDS::Edge( ancestor );
2067     }
2068   }
2069   return TopoDS_Edge();
2070 }
2071
2072 //================================================================================
2073 /*!
2074  * \brief Fill block subshapes
2075   * \param shapeMap - map to fill in
2076   * \retval int - nb inserted subshapes
2077  */
2078 //================================================================================
2079
2080 int StdMeshers_PrismAsBlock::TSideFace::InsertSubShapes(TBlockShapes& shapeMap) const
2081 {
2082   int nbInserted = 0;
2083
2084   // Insert edges
2085   vector< int > edgeIdVec;
2086   SMESH_Block::GetFaceEdgesIDs( myID, edgeIdVec );
2087
2088   for ( int i = BOTTOM_EDGE; i <=V1_EDGE ; ++i ) {
2089     TopoDS_Edge e = GetEdge( i );
2090     if ( !e.IsNull() ) {
2091       nbInserted += SMESH_Block::Insert( e, edgeIdVec[ i ], shapeMap);
2092     }
2093   }
2094
2095   // Insert corner vertices
2096
2097   TParam2ColumnIt col1, col2 ;
2098   vector< int > vertIdVec;
2099
2100   // from V0 column
2101   SMESH_Block::GetEdgeVertexIDs( edgeIdVec[ V0_EDGE ], vertIdVec);
2102   GetColumns(0, col1, col2 );
2103   const SMDS_MeshNode* node0 = col1->second.front();
2104   const SMDS_MeshNode* node1 = col1->second.back();
2105   TopoDS_Shape v0 = myHelper->GetSubShapeByNode( node0, myHelper->GetMeshDS());
2106   TopoDS_Shape v1 = myHelper->GetSubShapeByNode( node1, myHelper->GetMeshDS());
2107   if ( v0.ShapeType() == TopAbs_VERTEX ) {
2108     nbInserted += SMESH_Block::Insert( v0, vertIdVec[ 0 ], shapeMap);
2109   }
2110   if ( v1.ShapeType() == TopAbs_VERTEX ) {
2111     nbInserted += SMESH_Block::Insert( v1, vertIdVec[ 1 ], shapeMap);
2112   }
2113
2114   // from V1 column
2115   SMESH_Block::GetEdgeVertexIDs( edgeIdVec[ V1_EDGE ], vertIdVec);
2116   GetColumns(1, col1, col2 );
2117   node0 = col2->second.front();
2118   node1 = col2->second.back();
2119   v0 = myHelper->GetSubShapeByNode( node0, myHelper->GetMeshDS());
2120   v1 = myHelper->GetSubShapeByNode( node1, myHelper->GetMeshDS());
2121   if ( v0.ShapeType() == TopAbs_VERTEX ) {
2122     nbInserted += SMESH_Block::Insert( v0, vertIdVec[ 0 ], shapeMap);
2123   }
2124   if ( v1.ShapeType() == TopAbs_VERTEX ) {
2125     nbInserted += SMESH_Block::Insert( v1, vertIdVec[ 1 ], shapeMap);
2126   }
2127
2128 //   TopoDS_Vertex V0, V1, Vcom;
2129 //   TopExp::Vertices( myBaseEdge, V0, V1, true );
2130 //   if ( !myIsForward ) std::swap( V0, V1 );
2131
2132 //   // bottom vertex IDs
2133 //   SMESH_Block::GetEdgeVertexIDs( edgeIdVec[ _u0 ], vertIdVec);
2134 //   SMESH_Block::Insert( V0, vertIdVec[ 0 ], shapeMap);
2135 //   SMESH_Block::Insert( V1, vertIdVec[ 1 ], shapeMap);
2136
2137 //   TopoDS_Edge sideEdge = GetEdge( V0_EDGE );
2138 //   if ( sideEdge.IsNull() || !TopExp::CommonVertex( botEdge, sideEdge, Vcom ))
2139 //     return false;
2140
2141 //   // insert one side edge
2142 //   int edgeID;
2143 //   if ( Vcom.IsSame( V0 )) edgeID = edgeIdVec[ _v0 ];
2144 //   else                    edgeID = edgeIdVec[ _v1 ];
2145 //   SMESH_Block::Insert( sideEdge, edgeID, shapeMap);
2146
2147 //   // top vertex of the side edge
2148 //   SMESH_Block::GetEdgeVertexIDs( edgeID, vertIdVec);
2149 //   TopoDS_Vertex Vtop = TopExp::FirstVertex( sideEdge );
2150 //   if ( Vcom.IsSame( Vtop ))
2151 //     Vtop = TopExp::LastVertex( sideEdge );
2152 //   SMESH_Block::Insert( Vtop, vertIdVec[ 1 ], shapeMap);
2153
2154 //   // other side edge
2155 //   sideEdge = GetEdge( V1_EDGE );
2156 //   if ( sideEdge.IsNull() )
2157 //     return false;
2158 //   if ( edgeID = edgeIdVec[ _v1 ]) edgeID = edgeIdVec[ _v0 ];
2159 //   else                            edgeID = edgeIdVec[ _v1 ];
2160 //   SMESH_Block::Insert( sideEdge, edgeID, shapeMap);
2161
2162 //   // top edge
2163 //   TopoDS_Edge topEdge = GetEdge( TOP_EDGE );
2164 //   SMESH_Block::Insert( topEdge, edgeIdVec[ _u1 ], shapeMap);
2165
2166 //   // top vertex of the other side edge
2167 //   if ( !TopExp::CommonVertex( topEdge, sideEdge, Vcom ))
2168 //     return false;
2169 //   SMESH_Block::GetEdgeVertexIDs( edgeID, vertIdVec );
2170 //   SMESH_Block::Insert( Vcom, vertIdVec[ 1 ], shapeMap);
2171
2172   return nbInserted;
2173 }
2174
2175 //================================================================================
2176 /*!
2177  * \brief Dump ids of nodes of sides
2178  */
2179 //================================================================================
2180
2181 void StdMeshers_PrismAsBlock::TSideFace::dumpNodes(int nbNodes) const
2182 {
2183 #ifdef _DEBUG_
2184   cout << endl << "NODES OF FACE "; SMESH_Block::DumpShapeID( myID, cout ) << endl;
2185   THorizontalEdgeAdaptor* hSize0 = (THorizontalEdgeAdaptor*) HorizCurve(0);
2186   cout << "Horiz side 0: "; hSize0->dumpNodes(nbNodes); cout << endl;
2187   THorizontalEdgeAdaptor* hSize1 = (THorizontalEdgeAdaptor*) HorizCurve(1);
2188   cout << "Horiz side 1: "; hSize1->dumpNodes(nbNodes); cout << endl;
2189   TVerticalEdgeAdaptor* vSide0 = (TVerticalEdgeAdaptor*) VertiCurve(0);
2190   cout << "Verti side 0: "; vSide0->dumpNodes(nbNodes); cout << endl;
2191   TVerticalEdgeAdaptor* vSide1 = (TVerticalEdgeAdaptor*) VertiCurve(1);
2192   cout << "Verti side 1: "; vSide1->dumpNodes(nbNodes); cout << endl;
2193   delete hSize0; delete hSize1; delete vSide0; delete vSide1;
2194 #endif
2195 }
2196
2197 //================================================================================
2198 /*!
2199  * \brief Creates TVerticalEdgeAdaptor
2200   * \param columnsMap - node column map
2201   * \param parameter - normalized parameter
2202  */
2203 //================================================================================
2204
2205 StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::
2206 TVerticalEdgeAdaptor( const TParam2ColumnMap* columnsMap, const double parameter)
2207 {
2208   myNodeColumn = & getColumn( columnsMap, parameter )->second;
2209 }
2210
2211 //================================================================================
2212 /*!
2213  * \brief Return coordinates for the given normalized parameter
2214   * \param U - normalized parameter
2215   * \retval gp_Pnt - coordinates
2216  */
2217 //================================================================================
2218
2219 gp_Pnt StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::Value(const Standard_Real U) const
2220 {
2221   const SMDS_MeshNode* n1;
2222   const SMDS_MeshNode* n2;
2223   double r = getRAndNodes( myNodeColumn, U, n1, n2 );
2224   return gpXYZ(n1) * ( 1 - r ) + gpXYZ(n2) * r;
2225 }
2226
2227 //================================================================================
2228 /*!
2229  * \brief Dump ids of nodes
2230  */
2231 //================================================================================
2232
2233 void StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::dumpNodes(int nbNodes) const
2234 {
2235 #ifdef _DEBUG_
2236   for ( int i = 0; i < nbNodes && i < myNodeColumn->size(); ++i )
2237     cout << (*myNodeColumn)[i]->GetID() << " ";
2238   if ( nbNodes < myNodeColumn->size() )
2239     cout << myNodeColumn->back()->GetID();
2240 #endif
2241 }
2242
2243 //================================================================================
2244 /*!
2245  * \brief Return coordinates for the given normalized parameter
2246   * \param U - normalized parameter
2247   * \retval gp_Pnt - coordinates
2248  */
2249 //================================================================================
2250
2251 gp_Pnt StdMeshers_PrismAsBlock::THorizontalEdgeAdaptor::Value(const Standard_Real U) const
2252 {
2253   return mySide->TSideFace::Value( U, myV );
2254 }
2255
2256 //================================================================================
2257 /*!
2258  * \brief Dump ids of <nbNodes> first nodes and the last one
2259  */
2260 //================================================================================
2261
2262 void StdMeshers_PrismAsBlock::THorizontalEdgeAdaptor::dumpNodes(int nbNodes) const
2263 {
2264 #ifdef _DEBUG_
2265   // Not bedugged code. Last node is sometimes incorrect
2266   const TSideFace* side = mySide;
2267   double u = 0;
2268   if ( mySide->IsComplex() )
2269     side = mySide->GetComponent(0,u);
2270
2271   TParam2ColumnIt col, col2;
2272   TParam2ColumnMap* u2cols = side->GetColumns();
2273   side->GetColumns( u , col, col2 );
2274
2275   int j, i = myV ? mySide->ColumnHeight()-1 : 0;
2276
2277   const SMDS_MeshNode* n = 0;
2278   const SMDS_MeshNode* lastN
2279     = side->IsForward() ? u2cols->rbegin()->second[ i ] : u2cols->begin()->second[ i ];
2280   for ( j = 0; j < nbNodes && n != lastN; ++j )
2281   {
2282     n = col->second[ i ];
2283     cout << n->GetID() << " ";
2284     if ( side->IsForward() )
2285       ++col;
2286     else
2287       --col;
2288   }
2289
2290   // last node
2291   u = 1;
2292   if ( mySide->IsComplex() )
2293     side = mySide->GetComponent(1,u);
2294
2295   side->GetColumns( u , col, col2 );
2296   if ( n != col->second[ i ] )
2297     cout << col->second[ i ]->GetID();
2298 #endif
2299 }
2300 //================================================================================
2301 /*!
2302  * \brief Return UV on pcurve for the given normalized parameter
2303   * \param U - normalized parameter
2304   * \retval gp_Pnt - coordinates
2305  */
2306 //================================================================================
2307
2308 gp_Pnt2d StdMeshers_PrismAsBlock::TPCurveOnHorFaceAdaptor::Value(const Standard_Real U) const
2309 {
2310   TParam2ColumnIt u_col1, u_col2;
2311   double r = mySide->GetColumns( U, u_col1, u_col2 );
2312   gp_XY uv1 = mySide->GetNodeUV( myFace, u_col1->second[ myZ ]);
2313   gp_XY uv2 = mySide->GetNodeUV( myFace, u_col2->second[ myZ ]);
2314   return uv1 * ( 1 - r ) + uv2 * r;
2315 }