src/StdMeshers/StdMeshers_Penta_3D.cxx

   1 //  SMESH StdMeshers_Penta_3D implementaion of SMESH idl descriptions
   2 //
   3 //  Copyright (C) 2003  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 //
  24 //  File   : StdMeshers_Penta_3D.cxx
  25 //  Module : SMESH
  26
  27 using namespace std;
  28
  29 #include "StdMeshers_Penta_3D.hxx"
  30
  31 #include "utilities.h"
  32 #include "Utils_ExceptHandlers.hxx"
  33
  34 #include "SMDS_EdgePosition.hxx"
  35 #include "SMDS_MeshElement.hxx"
  36 #include "SMDS_VolumeOfNodes.hxx"
  37 #include "SMDS_VolumeTool.hxx"
  38 #include "SMESHDS_SubMesh.hxx"
  39 #include "SMESH_Mesh.hxx"
  40 #include "SMESH_MeshEditor.hxx"
  41 #include "SMESH_subMesh.hxx"
  42 #include "SMESH_subMeshEventListener.hxx"
  43 #include "SMESH_Comment.hxx"
  44
  45 #include <BRep_Tool.hxx>
  46 #include <TopExp.hxx>
  47 #include <TopExp_Explorer.hxx>
  48 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
  49 #include <TopTools_IndexedMapOfShape.hxx>
  50 #include <TopTools_ListIteratorOfListOfShape.hxx>
  51 #include <TopTools_ListOfShape.hxx>
  52 #include <TopoDS.hxx>
  53 #include <TopoDS_Edge.hxx>
  54 #include <TopoDS_Shell.hxx>
  55 #include <TopoDS_Vertex.hxx>
  56 #include <gp_Pnt.hxx>
  57
  58 #include <stdio.h>
  59 #include <algorithm>
  60
  61 using namespace std;
  62
  63 typedef map < int, int, less<int> >::iterator   \
  64   StdMeshers_IteratorOfDataMapOfIntegerInteger;
  65
  66 enum { NB_WALL_FACES = 4 };
  67
  68 //=======================================================================
  69 //function : StdMeshers_Penta_3D
  70 //purpose  :
  71 //=======================================================================
  72 StdMeshers_Penta_3D::StdMeshers_Penta_3D()
  73 : myErrorStatus(SMESH_ComputeError::New())
  74 {
  75   myTol3D=0.1;
  76   myWallNodesMaps.resize( SMESH_Block::NbFaces() );
  77   myShapeXYZ.resize( SMESH_Block::NbSubShapes() );
  78   myTool = 0;
  79 }
  80
  81 //=======================================================================
  82 //function : ~StdMeshers_Penta_3D
  83 //purpose  :
  84 //=======================================================================
  85
  86 StdMeshers_Penta_3D::~StdMeshers_Penta_3D()
  87 {
  88 }
  89
  90 //=======================================================================
  91 //function : Compute
  92 //purpose  :
  93 //=======================================================================
  94 bool StdMeshers_Penta_3D::Compute(SMESH_Mesh& aMesh,
  95                                   const TopoDS_Shape& aShape)
  96 {
  97   MESSAGE("StdMeshers_Penta_3D::Compute()");
  98   //
  99   bool bOK=false;
 100   //
 101   myShape=aShape;
 102   SetMesh(aMesh);
 103   //
 104   CheckData();
 105   if (!myErrorStatus->IsOK()) {
 106     return bOK;
 107   }
 108
 109   SMESH_MesherHelper helper(aMesh);
 110   myTool = &helper;
 111   myCreateQuadratic = myTool->IsQuadraticSubMesh(aShape);
 112
 113   //
 114   MakeBlock();
 115   if (!myErrorStatus->IsOK()) {
 116     return bOK;
 117   }
 118   //
 119   ClearMeshOnFxy1();
 120   if (!myErrorStatus->IsOK()) {
 121     return bOK;
 122   }
 123   //
 124   MakeNodes();
 125   if (!myErrorStatus->IsOK()) {
 126     return bOK;
 127   }
 128   //
 129   MakeConnectingMap();
 130   //
 131   MakeMeshOnFxy1();
 132   if (!myErrorStatus->IsOK()) {
 133     return bOK;
 134   }
 135   //
 136   MakeVolumeMesh();
 137   //
 138   return !bOK;
 139 }
 140
 141 //=======================================================================
 142 //function : MakeNodes
 143 //purpose  :
 144 //=======================================================================
 145 void StdMeshers_Penta_3D::MakeNodes()
 146 {
 147   const int aNbSIDs=9;
 148   int i, j, k, ij, iNbN, aNodeID, aSize, iErr;
 149   double aX, aY, aZ;
 150   SMESH_Block::TShapeID aSID, aSIDs[aNbSIDs]={
 151     SMESH_Block::ID_V000, SMESH_Block::ID_V100,
 152     SMESH_Block::ID_V110, SMESH_Block::ID_V010,
 153     SMESH_Block::ID_Ex00, SMESH_Block::ID_E1y0,
 154     SMESH_Block::ID_Ex10, SMESH_Block::ID_E0y0,
 155     SMESH_Block::ID_Fxy0
 156   };
 157   //
 158   SMESH_Mesh* pMesh=GetMesh();
 159   //
 160   // 1. Define the sizes of mesh
 161   //
 162   // 1.1 Horizontal size
 163   myJSize=0;
 164   for (i=0; i<aNbSIDs; ++i) {
 165     const TopoDS_Shape& aS = myBlock.Shape(aSIDs[i]);
 166     SMESH_subMesh *aSubMesh = pMesh->GetSubMeshContaining(aS);
 167     ASSERT(aSubMesh);
 168     SMESHDS_SubMesh *aSM = aSubMesh->GetSubMeshDS();
 169     if(!myCreateQuadratic) {
 170       iNbN = aSM->NbNodes();
 171     }
 172     else {
 173       iNbN = 0;
 174       SMDS_NodeIteratorPtr itn = aSM->GetNodes();
 175       while(itn->more()) {
 176         const SMDS_MeshNode* aNode = itn->next();
 177         if(myTool->IsMedium(aNode))
 178           continue;
 179         iNbN++;
 180       }
 181     }
 182     myJSize += iNbN;
 183   }
 184   //printf("***  Horizontal: number of nodes summary=%d\n", myJSize);
 185   //
 186   // 1.2 Vertical size
 187   myISize=2;
 188   {
 189     const TopoDS_Shape& aS=myBlock.Shape(SMESH_Block::ID_E00z);
 190     SMESH_subMesh *aSubMesh = pMesh->GetSubMeshContaining(aS);
 191     ASSERT(aSubMesh);
 192     SMESHDS_SubMesh *aSM = aSubMesh->GetSubMeshDS();
 193     if(!myCreateQuadratic) {
 194       iNbN = aSM->NbNodes();
 195     }
 196     else {
 197       iNbN = 0;
 198       SMDS_NodeIteratorPtr itn = aSM->GetNodes();
 199       while(itn->more()) {
 200         const SMDS_MeshNode* aNode = itn->next();
 201         if(myTool->IsMedium(aNode))
 202           continue;
 203         iNbN++;
 204       }
 205     }
 206     myISize += iNbN;
 207   }
 208   //printf("***  Vertical: number of nodes on edges and vertices=%d\n", myISize);
 209   //
 210   aSize=myISize*myJSize;
 211   myTNodes.resize(aSize);
 212   //
 213   StdMeshers_TNode aTNode;
 214   gp_XYZ aCoords;
 215   gp_Pnt aP3D;
 216   //
 217   // 2. Fill the repers on base face (Z=0)
 218   i=0; j=0;
 219   // vertices
 220   for (k=0; k<aNbSIDs; ++k) {
 221     aSID=aSIDs[k];
 222     const TopoDS_Shape& aS = myBlock.Shape(aSID);
 223     SMDS_NodeIteratorPtr ite = pMesh->GetSubMeshContaining(aS)->GetSubMeshDS()->GetNodes();
 224     while(ite->more()) {
 225       const SMDS_MeshNode* aNode = ite->next();
 226       if(myTool->IsMedium(aNode))
 227         continue;
 228       aNodeID=aNode->GetID();
 229       //
 230       aTNode.SetNode(aNode);
 231       aTNode.SetShapeSupportID(aSID);
 232       aTNode.SetBaseNodeID(aNodeID);
 233       //
 234       if ( SMESH_Block::IsEdgeID (aSID)) {
 235         const SMDS_EdgePosition* epos =
 236           static_cast<const SMDS_EdgePosition*>(aNode->GetPosition().get());
 237         myBlock.ComputeParameters( epos->GetUParameter(), aS, aCoords );
 238       }
 239       else {
 240         aX=aNode->X();
 241         aY=aNode->Y();
 242         aZ=aNode->Z();
 243         aP3D.SetCoord(aX, aY, aZ);
 244         myBlock.ComputeParameters(aP3D, aS, aCoords);
 245       }
 246       iErr = myBlock.ErrorStatus();
 247       if (iErr) {
 248         MESSAGE("StdMeshers_Penta_3D::MakeNodes()," <<
 249                 "SMESHBlock: ComputeParameters operation failed");
 250         myErrorStatus=myBlock.GetError();
 251         return;
 252       }
 253       aTNode.SetNormCoord(aCoords);
 254       ij=i*myJSize+j;
 255       myTNodes[ij]=aTNode;
 256       ++j;
 257     }
 258   }
 259
 260   // 3.1 Fill maps of wall nodes
 261   SMESH_Block::TShapeID wallFaceID[ NB_WALL_FACES ] = {
 262     SMESH_Block::ID_Fx0z, SMESH_Block::ID_Fx1z,
 263     SMESH_Block::ID_F0yz, SMESH_Block::ID_F1yz
 264     };
 265   SMESH_Block::TShapeID baseEdgeID[ NB_WALL_FACES ] = {
 266     SMESH_Block::ID_Ex00, SMESH_Block::ID_Ex10,
 267     SMESH_Block::ID_E0y0, SMESH_Block::ID_E1y0
 268     };
 269   for ( i = 0; i < NB_WALL_FACES ; ++i ) {
 270     int fIndex = SMESH_Block::ShapeIndex( wallFaceID[ i ]);
 271     bool ok = LoadIJNodes (myWallNodesMaps[ fIndex ],
 272                            TopoDS::Face( myBlock.Shape( wallFaceID[ i ] )),
 273                            TopoDS::Edge( myBlock.Shape( baseEdgeID[ i ] )),
 274                            pMesh->GetMeshDS());
 275     if ( !ok ) {
 276       myErrorStatus->myName = COMPERR_BAD_INPUT_MESH;
 277       myErrorStatus->myComment = SMESH_Comment() <<
 278         "Can't find regular quadrangle mesh on a side face #" <<
 279         pMesh->GetMeshDS()->ShapeToIndex( myBlock.Shape( wallFaceID[ i ]));
 280       return;
 281     }
 282   }
 283
 284   // 3.2 find node columns for vertical edges and edge IDs
 285   vector<const SMDS_MeshNode*> * verticEdgeNodes[ NB_WALL_FACES ];
 286   SMESH_Block::TShapeID          verticEdgeID   [ NB_WALL_FACES ];
 287   for ( i = 0; i < NB_WALL_FACES ; ++i ) { // 4 first base nodes are nodes on vertices
 288     // edge ID
 289     SMESH_Block::TShapeID eID, vID = aSIDs[ i ];
 290     ShapeSupportID(false, vID, eID);
 291     verticEdgeID[ i ] = eID;
 292     // column nodes
 293     StdMeshers_TNode& aTNode = myTNodes[ i ];
 294     verticEdgeNodes[ i ] = 0;
 295     for ( j = 0; j < NB_WALL_FACES ; ++j ) { // loop on 4 wall faces
 296       int fIndex = SMESH_Block::ShapeIndex( wallFaceID[ j ]);
 297       StdMeshers_IJNodeMap & ijNodes= myWallNodesMaps[ fIndex ];
 298       if ( ijNodes.begin()->second[0] == aTNode.Node() )
 299         verticEdgeNodes[ i ] = & ijNodes.begin()->second;
 300       else if ( ijNodes.rbegin()->second[0] == aTNode.Node() )
 301         verticEdgeNodes[ i ] = & ijNodes.rbegin()->second;
 302       if ( verticEdgeNodes[ i ] )
 303         break;
 304     }
 305   }
 306
 307   // 3.3 set XYZ of vertices, and initialize of the rest
 308   SMESHDS_Mesh* aMesh = GetMesh()->GetMeshDS();
 309   for ( int id = SMESH_Block::ID_V000; id < SMESH_Block::ID_Shell; ++id ) {
 310     if ( SMESH_Block::IsVertexID( id )) {
 311       TopoDS_Shape V = myBlock.Shape( id );
 312       SMESHDS_SubMesh* sm = aMesh->MeshElements( V );
 313       const SMDS_MeshNode* n = sm->GetNodes()->next();
 314       myShapeXYZ[ id ].SetCoord( n->X(), n->Y(), n->Z() );
 315     }
 316     else
 317       myShapeXYZ[ id ].SetCoord( 0., 0., 0. );
 318   }
 319
 320
 321   // 4. Fill the rest repers
 322   bool bIsUpperLayer;
 323   int iBNID;
 324   SMESH_Block::TShapeID aSSID, aBNSSID;
 325   StdMeshers_TNode aTN;
 326   //
 327
 328   // create top face and find UV for it's corners
 329   const TopoDS_Face& TopFace = TopoDS::Face(myBlock.Shape(SMESH_Block::ID_Fxy1));
 330   SMESHDS_Mesh* meshDS = pMesh->GetMeshDS();
 331   int topfaceID = meshDS->ShapeToIndex(TopFace);
 332   const TopoDS_Vertex& v001 = TopoDS::Vertex(myBlock.Shape(SMESH_Block::ID_V001));
 333   SMDS_NodeIteratorPtr itn = pMesh->GetSubMeshContaining(v001)->GetSubMeshDS()->GetNodes();
 334   const SMDS_MeshNode* N = itn->next();
 335   gp_XY UV001 = myTool->GetNodeUV(TopFace,N);
 336   const TopoDS_Vertex& v101 = TopoDS::Vertex(myBlock.Shape(SMESH_Block::ID_V101));
 337   itn = pMesh->GetSubMeshContaining(v101)->GetSubMeshDS()->GetNodes();
 338   N = itn->next();
 339   gp_XY UV101 = myTool->GetNodeUV(TopFace,N);
 340   const TopoDS_Vertex& v011 = TopoDS::Vertex(myBlock.Shape(SMESH_Block::ID_V011));
 341   itn = pMesh->GetSubMeshContaining(v011)->GetSubMeshDS()->GetNodes();
 342   N = itn->next();
 343   gp_XY UV011 = myTool->GetNodeUV(TopFace,N);
 344   const TopoDS_Vertex& v111 = TopoDS::Vertex(myBlock.Shape(SMESH_Block::ID_V111));
 345   itn = pMesh->GetSubMeshContaining(v111)->GetSubMeshDS()->GetNodes();
 346   N = itn->next();
 347   gp_XY UV111 = myTool->GetNodeUV(TopFace,N);
 348
 349   for (j=0; j<myJSize; ++j) { // loop on all nodes of the base face (ID_Fxy0)
 350     // base node info
 351     const StdMeshers_TNode& aBN = myTNodes[j];
 352     aBNSSID = (SMESH_Block::TShapeID)aBN.ShapeSupportID();
 353     iBNID = aBN.BaseNodeID();
 354     const gp_XYZ& aBNXYZ = aBN.NormCoord();
 355     bool createNode = ( aBNSSID == SMESH_Block::ID_Fxy0 ); // if base node is inside a bottom face
 356     //
 357     // set XYZ on horizontal edges and get node columns of faces:
 358     // 2 columns for each face, between which a base node is located
 359     vector<const SMDS_MeshNode*>* nColumns[8];
 360     double ratio[ NB_WALL_FACES ]; // base node position between columns [0.-1.]
 361     if ( createNode ) {
 362       for ( k = 0; k < NB_WALL_FACES ; ++k ) {
 363         ratio[ k ] = SetHorizEdgeXYZ (aBNXYZ, wallFaceID[ k ],
 364                                       nColumns[k*2], nColumns[k*2+1]);
 365       }
 366     }
 367     //
 368     // XYZ on the bottom and top faces
 369     const SMDS_MeshNode* n = aBN.Node();
 370     myShapeXYZ[ SMESH_Block::ID_Fxy0 ].SetCoord( n->X(), n->Y(), n->Z() );
 371     myShapeXYZ[ SMESH_Block::ID_Fxy1 ].SetCoord( 0., 0., 0. );
 372     //
 373     // first create or find a top node, then the rest ones in a column
 374     for (i=myISize-1; i>0; --i) // vertical loop, from top to bottom
 375     {
 376       bIsUpperLayer = (i==(myISize-1));
 377       gp_XY UV_Ex01, UV_Ex11, UV_E0y1, UV_E1y1;
 378       if ( createNode ) // a base node is inside a top face
 379       {
 380         // set XYZ on vertical edges and faces
 381         for ( k = 0; k < NB_WALL_FACES ; ++k ) {
 382           // XYZ on a vertical edge
 383           const SMDS_MeshNode* n = (*verticEdgeNodes[ k ]) [ i ];
 384           myShapeXYZ[ verticEdgeID[ k ] ].SetCoord( n->X(), n->Y(), n->Z() );
 385           // XYZ on a face (part 1 from one column)
 386           n = (*nColumns[k*2]) [ i ];
 387           gp_XYZ xyz( n->X(), n->Y(), n->Z() );
 388           myShapeXYZ[ wallFaceID[ k ]] = ( 1. - ratio[ k ]) * xyz;
 389           gp_XY tmp1;
 390           if( bIsUpperLayer ) {
 391             tmp1 = myTool->GetNodeUV(TopFace,n);
 392             tmp1 = ( 1. - ratio[ k ]) * tmp1;
 393           }
 394           // XYZ on a face (part 2 from other column)
 395           n = (*nColumns[k*2+1]) [ i ];
 396           xyz.SetCoord( n->X(), n->Y(), n->Z() );
 397           myShapeXYZ[ wallFaceID[ k ]] += ratio[ k ] * xyz;
 398           if( bIsUpperLayer ) {
 399             gp_XY tmp2 = myTool->GetNodeUV(TopFace,n);
 400             tmp1 +=  ratio[ k ] * tmp2;
 401             if( k==0 )
 402               UV_Ex01 = tmp1;
 403             else if( k==1 )
 404               UV_Ex11 = tmp1;
 405             else if( k==2 )
 406               UV_E0y1 = tmp1;
 407             else
 408               UV_E1y1 = tmp1;
 409           }
 410         }
 411       }
 412       // fill current node info
 413       //   -index in aTNodes
 414       ij=i*myJSize+j;
 415       //   -normalized coordinates
 416       aX=aBNXYZ.X();
 417       aY=aBNXYZ.Y();
 418       //aZ=aZL[i];
 419       aZ=(double)i/(double)(myISize-1);
 420       aCoords.SetCoord(aX, aY, aZ);
 421       //
 422       //   suporting shape ID
 423       ShapeSupportID(bIsUpperLayer, aBNSSID, aSSID);
 424       if (!myErrorStatus->IsOK()) {
 425         MESSAGE("StdMeshers_Penta_3D::MakeNodes() ");
 426         return;
 427       }
 428       //
 429       aTN.SetShapeSupportID(aSSID);
 430       aTN.SetNormCoord(aCoords);
 431       aTN.SetBaseNodeID(iBNID);
 432       //
 433       if (aSSID!=SMESH_Block::ID_NONE){
 434         // try to find the node
 435         const TopoDS_Shape& aS=myBlock.Shape((int)aSSID);
 436         FindNodeOnShape(aS, aCoords, i, aTN);
 437       }
 438       else{
 439         // create node and get it id
 440         CreateNode (bIsUpperLayer, aCoords, aTN);
 441         //
 442         if ( bIsUpperLayer ) {
 443           const SMDS_MeshNode* n = aTN.Node();
 444           myShapeXYZ[ SMESH_Block::ID_Fxy1 ].SetCoord( n->X(), n->Y(), n->Z() );
 445           // set node on top face:
 446           // find UV parameter for this node
 447           //              UV_Ex11
 448           //   UV011+-----+----------+UV111
 449           //        |                |
 450           //        |                |
 451           // UV_E0y1+     +node      +UV_E1y1
 452           //        |                |
 453           //        |                |
 454           //        |                |
 455           //   UV001+-----+----------+UV101
 456           //              UV_Ex01
 457           gp_Pnt2d aP;
 458           double u = aCoords.X(), v = aCoords.Y();
 459           double u1 = ( 1. - u ), v1 = ( 1. - v );
 460           aP.ChangeCoord()  = UV_Ex01 * v1;
 461           aP.ChangeCoord() += UV_Ex11 * v;
 462           aP.ChangeCoord() += UV_E0y1 * u1;
 463           aP.ChangeCoord() += UV_E1y1 * u;
 464           aP.ChangeCoord() -= UV001 * u1 * v1;
 465           aP.ChangeCoord() -= UV101 * u  * v1;
 466           aP.ChangeCoord() -= UV011 * u1 * v;
 467           aP.ChangeCoord() -= UV111 * u  * v;
 468           meshDS->SetNodeOnFace((SMDS_MeshNode*)n, topfaceID, aP.X(), aP.Y());
 469         }
 470       }
 471       if (!myErrorStatus->IsOK()) {
 472         MESSAGE("StdMeshers_Penta_3D::MakeNodes() ");
 473         return;
 474       }
 475       //
 476       myTNodes[ij]=aTN;
 477     }
 478   }
 479   //DEB
 480   /*
 481   {
 482     int iSSID, iBNID, aID;
 483     //
 484     for (i=0; i<myISize; ++i) {
 485       printf(" Layer# %d\n", i);
 486       for (j=0; j<myJSize; ++j) {
 487         ij=i*myJSize+j;
 488         const StdMeshers_TNode& aTN=myTNodes[ij];
 489         //const StdMeshers_TNode& aTN=aTNodes[ij];
 490         const gp_XYZ& aXYZ=aTN.NormCoord();
 491         iSSID=aTN.ShapeSupportID();
 492         iBNID=aTN.BaseNodeID();
 493         //
 494         const SMDS_MeshNode* aNode=aTN.Node();
 495         aID=aNode->GetID();
 496         aX=aNode->X();
 497         aY=aNode->Y();
 498         aZ=aNode->Z();
 499         printf("*** j:%d BNID#%d iSSID:%d ID:%d { %lf %lf %lf },  { %lf %lf %lf }\n",
 500                j,  iBNID, iSSID, aID, aXYZ.X(),  aXYZ.Y(), aXYZ.Z(), aX, aY, aZ);
 501       }
 502     }
 503   }
 504   */
 505   //DEB t
 506 }
 507
 508
 509 //=======================================================================
 510 //function : FindNodeOnShape
 511 //purpose  :
 512 //=======================================================================
 513
 514 void StdMeshers_Penta_3D::FindNodeOnShape(const TopoDS_Shape& aS,
 515                                           const gp_XYZ&       aParams,
 516                                           const int           z,
 517                                           StdMeshers_TNode&   aTN)
 518 {
 519   double aX, aY, aZ, aD, aTol2, minD;
 520   gp_Pnt aP1, aP2;
 521   //
 522   SMESH_Mesh* pMesh = GetMesh();
 523   aTol2 = myTol3D*myTol3D;
 524   minD = 1.e100;
 525   SMDS_MeshNode* pNode = NULL;
 526   //
 527   if ( aS.ShapeType() == TopAbs_FACE ||
 528        aS.ShapeType() == TopAbs_EDGE ) {
 529     // find a face ID to which aTN belongs to
 530     int faceID;
 531     if ( aS.ShapeType() == TopAbs_FACE )
 532       faceID = myBlock.ShapeID( aS );
 533     else { // edge maybe vertical or top horizontal
 534       gp_XYZ aCoord = aParams;
 535       if ( aCoord.Z() == 1. )
 536         aCoord.SetZ( 0.5 ); // move from top down
 537       else
 538         aCoord.SetX( 0.5 ); // move along X
 539       faceID = SMESH_Block::GetShapeIDByParams( aCoord );
 540     }
 541     ASSERT( SMESH_Block::IsFaceID( faceID ));
 542     int fIndex = SMESH_Block::ShapeIndex( faceID );
 543     StdMeshers_IJNodeMap & ijNodes = myWallNodesMaps[ fIndex ];
 544     // look for a base node in ijNodes
 545     const SMDS_MeshNode* baseNode = pMesh->GetMeshDS()->FindNode( aTN.BaseNodeID() );
 546     StdMeshers_IJNodeMap::const_iterator par_nVec = ijNodes.begin();
 547     for ( ; par_nVec != ijNodes.end(); par_nVec++ )
 548       if ( par_nVec->second[ 0 ] == baseNode ) {
 549         pNode = (SMDS_MeshNode*)par_nVec->second.at( z );
 550         aTN.SetNode(pNode);
 551         return;
 552       }
 553   }
 554   //
 555   myBlock.Point(aParams, aS, aP1);
 556   //
 557   SMDS_NodeIteratorPtr ite=
 558     pMesh->GetSubMeshContaining(aS)->GetSubMeshDS()->GetNodes();
 559   while(ite->more()) {
 560     const SMDS_MeshNode* aNode = ite->next();
 561     if(myTool->IsMedium(aNode))
 562       continue;
 563     aX=aNode->X();
 564     aY=aNode->Y();
 565     aZ=aNode->Z();
 566     aP2.SetCoord(aX, aY, aZ);
 567     aD=(double)aP1.SquareDistance(aP2);
 568     //printf("** D=%lf ", aD, aTol2);
 569     if (aD < minD) {
 570       pNode=(SMDS_MeshNode*)aNode;
 571       aTN.SetNode(pNode);
 572       minD = aD;
 573       //printf(" Ok\n");
 574       if (aD<aTol2)
 575         return;
 576     }
 577   }
 578   //
 579   //printf(" KO\n");
 580   //aTN.SetNode(pNode);
 581   //MESSAGE("StdMeshers_Penta_3D::FindNodeOnShape(), can not find the node");
 582   //myErrorStatus=11; // can not find the node;
 583 }
 584
 585
 586 //=======================================================================
 587 //function : SetHorizEdgeXYZ
 588 //purpose  :
 589 //=======================================================================
 590
 591 double StdMeshers_Penta_3D::SetHorizEdgeXYZ(const gp_XYZ&                  aBaseNodeParams,
 592                                             const int                      aFaceID,
 593                                             vector<const SMDS_MeshNode*>*& aCol1,
 594                                             vector<const SMDS_MeshNode*>*& aCol2)
 595 {
 596   // find base and top edges of the face
 597   enum { BASE = 0, TOP };
 598   vector< int > edgeVec; // 0-base, 1-top
 599   SMESH_Block::GetFaceEdgesIDs( aFaceID, edgeVec );
 600   //
 601   int coord = SMESH_Block::GetCoordIndOnEdge( edgeVec[ BASE ] );
 602   bool isForward = myBlock.IsForwadEdge( edgeVec[ BASE ] );
 603
 604   double param = aBaseNodeParams.Coord( coord );
 605   if ( !isForward)
 606     param = 1. - param;
 607   //
 608   // look for columns around param
 609   StdMeshers_IJNodeMap & ijNodes =
 610     myWallNodesMaps[ SMESH_Block::ShapeIndex( aFaceID )];
 611   StdMeshers_IJNodeMap::iterator par_nVec_1 = ijNodes.begin();
 612   while ( par_nVec_1->first < param )
 613     par_nVec_1++;
 614   StdMeshers_IJNodeMap::iterator par_nVec_2 = par_nVec_1;
 615   //
 616   double r = 0;
 617   if ( par_nVec_1 != ijNodes.begin() ) {
 618     par_nVec_1--;
 619     r = ( param - par_nVec_1->first ) / ( par_nVec_2->first - par_nVec_1->first );
 620   }
 621   aCol1 = & par_nVec_1->second;
 622   aCol2 = & par_nVec_2->second;
 623
 624   // top edge
 625   if (1) {
 626     // this variant is better for cases with curved edges and
 627     // different nodes distribution on top and base edges
 628     const SMDS_MeshNode* n1 = aCol1->back();
 629     const SMDS_MeshNode* n2 = aCol2->back();
 630     gp_XYZ xyz1( n1->X(), n1->Y(), n1->Z() );
 631     gp_XYZ xyz2( n2->X(), n2->Y(), n2->Z() );
 632     myShapeXYZ[ edgeVec[ 1 ] ] = ( 1. - r ) * xyz1 + r * xyz2;
 633   }
 634   else {
 635     // this variant is better for other cases
 636 //   SMESH_MesherHelper helper( *GetMesh() );
 637 //   const TopoDS_Edge & edge = TopoDS::Edge( myBlock.Shape( edgeVec[ TOP ]));
 638 //   double u1 = helper.GetNodeU( edge, n1 );
 639 //   double u2 = helper.GetNodeU( edge, n2 );
 640 //   double u = ( 1. - r ) * u1 + r * u2;
 641 //   gp_XYZ topNodeParams;
 642 //   myBlock.Block().EdgeParameters( edgeVec[ TOP ], u, topNodeParams );
 643 //   myBlock.Block().EdgePoint( edgeVec[ TOP ],
 644 //                              topNodeParams,
 645 //                              myShapeXYZ[ edgeVec[ TOP ]]);
 646   }
 647
 648   // base edge
 649   myBlock.Block().EdgePoint( edgeVec[ BASE ],
 650                              aBaseNodeParams,
 651                              myShapeXYZ[ edgeVec[ BASE ]]);
 652   return r;
 653 }
 654
 655
 656 //=======================================================================
 657 //function : MakeVolumeMesh
 658 //purpose  :
 659 //=======================================================================
 660 void StdMeshers_Penta_3D::MakeVolumeMesh()
 661 {
 662   int i, j, ij, ik, i1, i2, aSSID;
 663   //
 664   SMESH_Mesh*   pMesh = GetMesh();
 665   SMESHDS_Mesh* meshDS = pMesh->GetMeshDS();
 666   //
 667   int shapeID = meshDS->ShapeToIndex( myShape );
 668   //
 669   // 1. Set Node In Volume
 670   ik = myISize-1;
 671   for (i=1; i<ik; ++i){
 672     for (j=0; j<myJSize; ++j){
 673       ij=i*myJSize+j;
 674       const StdMeshers_TNode& aTN = myTNodes[ij];
 675       aSSID=aTN.ShapeSupportID();
 676       if (aSSID==SMESH_Block::ID_NONE) {
 677         SMDS_MeshNode* aNode = (SMDS_MeshNode*)aTN.Node();
 678         meshDS->SetNodeInVolume(aNode, shapeID);
 679       }
 680     }
 681   }
 682   //
 683   // 2. Make pentahedrons
 684   int aID0, k , aJ[3];
 685   vector<const SMDS_MeshNode*> aN;
 686   //
 687   SMDS_ElemIteratorPtr itf, aItNodes;
 688   //
 689   const TopoDS_Face& aFxy0=
 690     TopoDS::Face(myBlock.Shape(SMESH_Block::ID_Fxy0));
 691   SMESH_subMesh *aSubMesh0 = pMesh->GetSubMeshContaining(aFxy0);
 692   SMESHDS_SubMesh *aSM0 = aSubMesh0->GetSubMeshDS();
 693   //
 694   itf = aSM0->GetElements();
 695   while(itf->more()) {
 696     const SMDS_MeshElement* pE0 = itf->next();
 697     //
 698     int nbFaceNodes = pE0->NbNodes();
 699     if(myCreateQuadratic)
 700       nbFaceNodes = nbFaceNodes/2;
 701     if ( aN.size() < nbFaceNodes * 2 )
 702       aN.resize( nbFaceNodes * 2 );
 703     //
 704     k=0;
 705     aItNodes=pE0->nodesIterator();
 706     while (aItNodes->more()) {
 707       //const SMDS_MeshElement* pNode = aItNodes->next();
 708       const SMDS_MeshNode* pNode =
 709         static_cast<const SMDS_MeshNode*> (aItNodes->next());
 710       if(myTool->IsMedium(pNode))
 711         continue;
 712       aID0 = pNode->GetID();
 713       aJ[k] = GetIndexOnLayer(aID0);
 714       if (!myErrorStatus->IsOK()) {
 715         MESSAGE("StdMeshers_Penta_3D::MakeVolumeMesh");
 716         return;
 717       }
 718       //
 719       ++k;
 720     }
 721     //
 722     bool forward = true;
 723     for (i=0; i<ik; ++i) {
 724       i1=i;
 725       i2=i+1;
 726       for(j=0; j<nbFaceNodes; ++j) {
 727         ij = i1*myJSize+aJ[j];
 728         const StdMeshers_TNode& aTN1 = myTNodes[ij];
 729         const SMDS_MeshNode* aN1 = aTN1.Node();
 730         aN[j]=aN1;
 731         //
 732         ij=i2*myJSize+aJ[j];
 733         const StdMeshers_TNode& aTN2 = myTNodes[ij];
 734         const SMDS_MeshNode* aN2 = aTN2.Node();
 735         aN[j+nbFaceNodes] = aN2;
 736       }
 737       // check if volume orientation will be ok
 738       if ( i == 0 ) {
 739         SMDS_VolumeTool vTool;
 740         switch ( nbFaceNodes ) {
 741         case 3: {
 742           SMDS_VolumeOfNodes tmpVol (aN[0], aN[1], aN[2],
 743                                      aN[3], aN[4], aN[5]);
 744           vTool.Set( &tmpVol );
 745           break;
 746         }
 747         case 4: {
 748           SMDS_VolumeOfNodes tmpVol(aN[0], aN[1], aN[2], aN[3],
 749                                     aN[4], aN[5], aN[6], aN[7]);
 750           vTool.Set( &tmpVol );
 751           break;
 752         }
 753         default:
 754           continue;
 755         }
 756         forward = vTool.IsForward();
 757       }
 758       // add volume
 759       SMDS_MeshVolume* aV = 0;
 760       switch ( nbFaceNodes ) {
 761       case 3:
 762         if ( forward ) {
 763           //aV = meshDS->AddVolume(aN[0], aN[1], aN[2],
 764           //                       aN[3], aN[4], aN[5]);
 765           aV = myTool->AddVolume(aN[0], aN[1], aN[2], aN[3], aN[4], aN[5]);
 766         }
 767         else {
 768           //aV = meshDS->AddVolume(aN[0], aN[2], aN[1],
 769           //                       aN[3], aN[5], aN[4]);
 770           aV = myTool->AddVolume(aN[0], aN[2], aN[1], aN[3], aN[5], aN[4]);
 771         }
 772         break;
 773       case 4:
 774         if ( forward ) {
 775           //aV = meshDS->AddVolume(aN[0], aN[1], aN[2], aN[3],
 776           //                       aN[4], aN[5], aN[6], aN[7]);
 777           aV = myTool->AddVolume(aN[0], aN[1], aN[2], aN[3],
 778                                  aN[4], aN[5], aN[6], aN[7]);
 779         }
 780         else {
 781           //aV = meshDS->AddVolume(aN[0], aN[3], aN[2], aN[1],
 782           //                       aN[4], aN[7], aN[6], aN[5]);
 783           aV = myTool->AddVolume(aN[0], aN[3], aN[2], aN[1],
 784                                  aN[4], aN[7], aN[6], aN[5]);
 785         }
 786         break;
 787       default:
 788         continue;
 789       }
 790       meshDS->SetMeshElementOnShape(aV, shapeID);
 791     }
 792   }
 793 }
 794
 795 //=======================================================================
 796 //function : MakeMeshOnFxy1
 797 //purpose  :
 798 //=======================================================================
 799 void StdMeshers_Penta_3D::MakeMeshOnFxy1()
 800 {
 801   int aID0, aJ, aLevel, ij, aNbNodes, k;
 802   //
 803   SMDS_NodeIteratorPtr itn;
 804   SMDS_ElemIteratorPtr itf, aItNodes;
 805   SMDSAbs_ElementType aElementType;
 806   //
 807   const TopoDS_Face& aFxy0=
 808     TopoDS::Face(myBlock.Shape(SMESH_Block::ID_Fxy0));
 809   const TopoDS_Face& aFxy1=
 810     TopoDS::Face(myBlock.Shape(SMESH_Block::ID_Fxy1));
 811   SMESH_MesherHelper faceHelper( *GetMesh() );
 812   faceHelper.IsQuadraticSubMesh(aFxy1);
 813   //
 814   SMESH_Mesh* pMesh = GetMesh();
 815   SMESHDS_Mesh * meshDS = pMesh->GetMeshDS();
 816   //
 817   SMESH_subMesh *aSubMesh1 = pMesh->GetSubMeshContaining(aFxy1);
 818   SMESH_subMesh *aSubMesh0 = pMesh->GetSubMeshContaining(aFxy0);
 819   SMESHDS_SubMesh *aSM0 = aSubMesh0->GetSubMeshDS();
 820   //
 821   // set nodes on aFxy1
 822   aLevel = myISize-1;
 823   itn = aSM0->GetNodes();
 824   aNbNodes = aSM0->NbNodes();
 825   //printf("** aNbNodes=%d\n", aNbNodes);
 826
 827   //
 828   // set elements on aFxy1
 829   vector<const SMDS_MeshNode*> aNodes1;
 830   //
 831   itf = aSM0->GetElements();
 832   while(itf->more()) {
 833     const SMDS_MeshElement* pE0 = itf->next();
 834     aElementType = pE0->GetType();
 835     if (!aElementType==SMDSAbs_Face) {
 836       continue;
 837     }
 838     aNbNodes = pE0->NbNodes();
 839     if(myCreateQuadratic)
 840       aNbNodes = aNbNodes/2;
 841     if ( aNodes1.size() < aNbNodes )
 842       aNodes1.resize( aNbNodes );
 843     //
 844     k = aNbNodes-1; // reverse a face
 845     aItNodes = pE0->nodesIterator();
 846     while (aItNodes->more()) {
 847       const SMDS_MeshNode* pNode =
 848         static_cast<const SMDS_MeshNode*> (aItNodes->next());
 849       if(myTool->IsMedium(pNode))
 850         continue;
 851       aID0 = pNode->GetID();
 852       aJ = GetIndexOnLayer(aID0);
 853       if (!myErrorStatus->IsOK()) {
 854         MESSAGE("StdMeshers_Penta_3D::MakeMeshOnFxy1() ");
 855         return;
 856       }
 857       //
 858       ij = aLevel*myJSize + aJ;
 859       const StdMeshers_TNode& aTN1 = myTNodes[ij];
 860       const SMDS_MeshNode* aN1 = aTN1.Node();
 861       aNodes1[k] = aN1;
 862       --k;
 863     }
 864     SMDS_MeshFace * face = 0;
 865     switch ( aNbNodes ) {
 866     case 3:
 867       face = faceHelper.AddFace(aNodes1[0], aNodes1[1], aNodes1[2]);
 868       break;
 869     case 4:
 870       face = faceHelper.AddFace(aNodes1[0], aNodes1[1], aNodes1[2], aNodes1[3]);
 871       break;
 872     default:
 873       continue;
 874     }
 875     meshDS->SetMeshElementOnShape(face, aFxy1);
 876   }
 877
 878   // update compute state of top face submesh
 879   aSubMesh1->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
 880
 881   // assure that mesh on the top face will be cleaned when it is cleaned
 882   // on the bottom face
 883   SMESH_subMesh* volSM = pMesh->GetSubMesh( myTool->GetSubShape() );
 884   volSM->SetEventListener( new SMESH_subMeshEventListener(true),
 885                            SMESH_subMeshEventListenerData::MakeData( aSubMesh1 ),
 886                            aSubMesh0 ); // translate CLEAN event of aSubMesh0 to aSubMesh1
 887 }
 888
 889 //=======================================================================
 890 //function : ClearMeshOnFxy1
 891 //purpose  :
 892 //=======================================================================
 893 void StdMeshers_Penta_3D::ClearMeshOnFxy1()
 894 {
 895   SMESH_subMesh* aSubMesh;
 896   SMESH_Mesh* pMesh=GetMesh();
 897   //
 898   const TopoDS_Shape& aFxy1=myBlock.Shape(SMESH_Block::ID_Fxy1);
 899   aSubMesh = pMesh->GetSubMeshContaining(aFxy1);
 900   if (aSubMesh)
 901     aSubMesh->ComputeStateEngine( SMESH_subMesh::CLEAN );
 902 }
 903
 904 //=======================================================================
 905 //function : GetIndexOnLayer
 906 //purpose  :
 907 //=======================================================================
 908 int StdMeshers_Penta_3D::GetIndexOnLayer(const int aID)
 909 {
 910   int j=-1;
 911   StdMeshers_IteratorOfDataMapOfIntegerInteger aMapIt;
 912   //
 913   aMapIt=myConnectingMap.find(aID);
 914   if (aMapIt==myConnectingMap.end()) {
 915     myErrorStatus->myName    = 200;
 916     myErrorStatus->myComment = "Internal error of StdMeshers_Penta_3D";
 917     return j;
 918   }
 919   j=(*aMapIt).second;
 920   return j;
 921 }
 922
 923 //=======================================================================
 924 //function : MakeConnectingMap
 925 //purpose  :
 926 //=======================================================================
 927 void StdMeshers_Penta_3D::MakeConnectingMap()
 928 {
 929   int j, aBNID;
 930   //
 931   for (j=0; j<myJSize; ++j) {
 932     const StdMeshers_TNode& aBN=myTNodes[j];
 933     aBNID=aBN.BaseNodeID();
 934     myConnectingMap[aBNID]=j;
 935   }
 936 }
 937
 938 //=======================================================================
 939 //function : CreateNode
 940 //purpose  :
 941 //=======================================================================
 942 void StdMeshers_Penta_3D::CreateNode(const bool bIsUpperLayer,
 943                                      const gp_XYZ& aParams,
 944                                      StdMeshers_TNode& aTN)
 945 {
 946   double aX, aY, aZ;
 947   //
 948   gp_Pnt aP;
 949   //
 950   SMDS_MeshNode* pNode=NULL;
 951   aTN.SetNode(pNode);
 952   //
 953 //   if (bIsUpperLayer) {
 954 //     // point on face Fxy1
 955 //     const TopoDS_Shape& aS=myBlock.Shape(SMESH_Block::ID_Fxy1);
 956 //     myBlock.Point(aParams, aS, aP);
 957 //   }
 958 //   else {
 959 //     // point inside solid
 960 //     myBlock.Point(aParams, aP);
 961 //   }
 962   if (bIsUpperLayer) {
 963     double u = aParams.X(), v = aParams.Y();
 964     double u1 = ( 1. - u ), v1 = ( 1. - v );
 965     aP.ChangeCoord()  = myShapeXYZ[ SMESH_Block::ID_Ex01 ] * v1;
 966     aP.ChangeCoord() += myShapeXYZ[ SMESH_Block::ID_Ex11 ] * v;
 967     aP.ChangeCoord() += myShapeXYZ[ SMESH_Block::ID_E0y1 ] * u1;
 968     aP.ChangeCoord() += myShapeXYZ[ SMESH_Block::ID_E1y1 ] * u;
 969
 970     aP.ChangeCoord() -= myShapeXYZ[ SMESH_Block::ID_V001 ] * u1 * v1;
 971     aP.ChangeCoord() -= myShapeXYZ[ SMESH_Block::ID_V101 ] * u  * v1;
 972     aP.ChangeCoord() -= myShapeXYZ[ SMESH_Block::ID_V011 ] * u1 * v;
 973     aP.ChangeCoord() -= myShapeXYZ[ SMESH_Block::ID_V111 ] * u  * v;
 974   }
 975   else {
 976     SMESH_Block::ShellPoint( aParams, myShapeXYZ, aP.ChangeCoord() );
 977   }
 978   //
 979 //   iErr=myBlock.ErrorStatus();
 980 //   if (iErr) {
 981 //     myErrorStatus=12; // can not find the node point;
 982 //     return;
 983 //   }
 984   //
 985   aX=aP.X(); aY=aP.Y(); aZ=aP.Z();
 986   //
 987   SMESH_Mesh* pMesh = GetMesh();
 988   SMESHDS_Mesh* pMeshDS = pMesh->GetMeshDS();
 989   //
 990   pNode = pMeshDS->AddNode(aX, aY, aZ);
 991
 992   aTN.SetNode(pNode);
 993 }
 994
 995 //=======================================================================
 996 //function : ShapeSupportID
 997 //purpose  :
 998 //=======================================================================
 999 void StdMeshers_Penta_3D::ShapeSupportID(const bool bIsUpperLayer,
1000                                          const SMESH_Block::TShapeID aBNSSID,
1001                                          SMESH_Block::TShapeID& aSSID)
1002 {
1003   switch (aBNSSID) {
1004     case SMESH_Block::ID_V000:
1005       aSSID=(bIsUpperLayer) ?  SMESH_Block::ID_V001 : SMESH_Block::ID_E00z;
1006       break;
1007     case SMESH_Block::ID_V100:
1008       aSSID=(bIsUpperLayer) ?  SMESH_Block::ID_V101 : SMESH_Block::ID_E10z;
1009       break;
1010     case SMESH_Block::ID_V110:
1011       aSSID=(bIsUpperLayer) ?  SMESH_Block::ID_V111 : SMESH_Block::ID_E11z;
1012       break;
1013     case SMESH_Block::ID_V010:
1014       aSSID=(bIsUpperLayer) ?  SMESH_Block::ID_V011 : SMESH_Block::ID_E01z;
1015       break;
1016     case SMESH_Block::ID_Ex00:
1017       aSSID=(bIsUpperLayer) ?  SMESH_Block::ID_Ex01 : SMESH_Block::ID_Fx0z;
1018       break;
1019     case SMESH_Block::ID_Ex10:
1020       aSSID=(bIsUpperLayer) ?  SMESH_Block::ID_Ex11 : SMESH_Block::ID_Fx1z;
1021       break;
1022     case SMESH_Block::ID_E0y0:
1023       aSSID=(bIsUpperLayer) ?  SMESH_Block::ID_E0y1 : SMESH_Block::ID_F0yz;
1024       break;
1025     case SMESH_Block::ID_E1y0:
1026       aSSID=(bIsUpperLayer) ?  SMESH_Block::ID_E1y1 : SMESH_Block::ID_F1yz;
1027       break;
1028     case SMESH_Block::ID_Fxy0:
1029       aSSID=SMESH_Block::ID_NONE;//(bIsUpperLayer) ?  Shape_ID_Fxy1 : Shape_ID_NONE;
1030       break;
1031     default:
1032       aSSID=SMESH_Block::ID_NONE;
1033       myErrorStatus->myName=10; // Can not find supporting shape ID
1034       myErrorStatus->myComment = "Internal error of StdMeshers_Penta_3D";
1035       break;
1036   }
1037   return;
1038 }
1039 //=======================================================================
1040 //function : MakeBlock
1041 //purpose  :
1042 //=======================================================================
1043 void StdMeshers_Penta_3D::MakeBlock()
1044 {
1045   bool bFound;
1046   int i, j, iNbEV, iNbE, iErr, iCnt, iNbNodes, iNbF;
1047   //
1048   TopoDS_Vertex aV000, aV001;
1049   TopoDS_Shape aFTr;
1050   TopTools_IndexedDataMapOfShapeListOfShape aMVES;
1051   TopTools_IndexedMapOfShape aME ,aMEV, aM;
1052   TopTools_ListIteratorOfListOfShape aIt;
1053   //
1054   TopExp::MapShapes(myShape, TopAbs_FACE, aM);
1055   //
1056   // 0. Find triangulated face aFTr
1057   SMDSAbs_ElementType aElementType;
1058   SMESH_Mesh* pMesh=GetMesh();
1059   //
1060   iCnt = 0;
1061   iNbF = aM.Extent();
1062   for (i=1; i<=iNbF; ++i) {
1063     const TopoDS_Shape& aF = aM(i);
1064     SMESH_subMesh *aSubMesh = pMesh->GetSubMeshContaining(aF);
1065     ASSERT(aSubMesh);
1066     SMESHDS_SubMesh *aSM = aSubMesh->GetSubMeshDS();
1067     SMDS_ElemIteratorPtr itf = aSM->GetElements();
1068     while(itf->more()) {
1069       const SMDS_MeshElement * pElement = itf->next();
1070       aElementType = pElement->GetType();
1071       if (aElementType==SMDSAbs_Face) {
1072         iNbNodes = pElement->NbNodes();
1073         if ( iNbNodes==3 || (myCreateQuadratic && iNbNodes==6) ) {
1074           aFTr = aF;
1075           ++iCnt;
1076           if (iCnt>1) {
1077             // \begin{E.A.}
1078             // The current algorithm fails if there is more that one
1079             // face wich contains triangles ...
1080             // In that case, replace return by break to try another
1081             // method (coded in "if (iCnt != 1) { ... }")
1082             //
1083             // MESSAGE("StdMeshers_Penta_3D::MakeBlock() ");
1084             // myErrorStatus=5; // more than one face has triangulation
1085             // return;
1086             break;
1087             // \end{E.A.}
1088           }
1089           break; // next face
1090         }
1091       }
1092     }
1093   }
1094   //
1095   // \begin{E.A.}
1096   // The current algorithm fails if "iCnt != 1", the case "iCnt == 0"
1097   // was not reached 'cause it was not called from Hexa_3D ... Now it
1098   // can occurs and in my opinion, it is the most common case.
1099   //
1100   if (iCnt != 1) {
1101     // The suggested algorithm is the following :
1102     //
1103     // o Check that nb_of_faces == 6 and nb_of_edges == 12
1104     //   then the shape is tologically equivalent to a box
1105     // o In a box, there are three set of four // edges ...
1106     //   In the cascade notation, it seems to be the edges
1107     //   numbered :
1108     //     - 1, 3, 5, 7
1109     //     - 2, 4, 6, 8
1110     //     - 9, 10, 11, 12
1111     // o For each one of this set, check if the four edges
1112     //   have the same number of element.
1113     // o If so, check if the "corresponding" // faces contains
1114     //   only quads. It's the faces numbered:
1115     //     - 1, 2, 3, 4
1116     //     - 1, 2, 5, 6
1117     //     - 3, 4, 5, 6
1118     // o If so, check if the opposite edges of each // faces
1119     //   have the same number of elements. It is the edges
1120     //   numbered :
1121     //     - 2 and 4, 6 and 8, 9 and 10, 11 and 12
1122     //     - 1 and 3, 5 and 7, 9 and 11, 10 and 12
1123     //     - 1 and 5, 3 and 7, 4 and 8, 2 and 6
1124     // o If so, check if the two other faces have the same
1125     //   number of elements. It is the faces numbered:
1126     //     - 5, 6
1127     //     - 3, 4
1128     //     - 1, 2
1129     //   This test should be improved to test if the nodes
1130     //   of the two faces are really "en face".
1131     // o If so, one of the two faces is a candidate to an extrusion,
1132     //   It is the faces numbered :
1133     //     - 5
1134     //     - 3
1135     //     - 1
1136     // o Finally, if there is only one candidate, let do the
1137     //   extrusion job for the corresponding face
1138     //
1139     int isOK = 0;
1140     //
1141     int iNbF = aM.Extent();
1142     if (iNbF == 6) {
1143       //
1144       int nb_f1 = pMesh->GetSubMeshContaining(aM(1))->GetSubMeshDS()->NbElements();
1145       int nb_f2 = pMesh->GetSubMeshContaining(aM(2))->GetSubMeshDS()->NbElements();
1146       int nb_f3 = pMesh->GetSubMeshContaining(aM(3))->GetSubMeshDS()->NbElements();
1147       int nb_f4 = pMesh->GetSubMeshContaining(aM(4))->GetSubMeshDS()->NbElements();
1148       int nb_f5 = pMesh->GetSubMeshContaining(aM(5))->GetSubMeshDS()->NbElements();
1149       int nb_f6 = pMesh->GetSubMeshContaining(aM(6))->GetSubMeshDS()->NbElements();
1150       //
1151       int has_only_quad_f1 = 1;
1152       int has_only_quad_f2 = 1;
1153       int has_only_quad_f3 = 1;
1154       int has_only_quad_f4 = 1;
1155       int has_only_quad_f5 = 1;
1156       int has_only_quad_f6 = 1;
1157       //
1158       for (i=1; i<=iNbF; ++i) {
1159         int ok = 1;
1160         const TopoDS_Shape& aF = aM(i);
1161         SMESH_subMesh *aSubMesh = pMesh->GetSubMeshContaining(aF);
1162         SMESHDS_SubMesh *aSM = aSubMesh->GetSubMeshDS();
1163         SMDS_ElemIteratorPtr itf = aSM->GetElements();
1164         while(itf->more()) {
1165           const SMDS_MeshElement * pElement = itf->next();
1166           aElementType = pElement->GetType();
1167           if (aElementType==SMDSAbs_Face) {
1168             iNbNodes = pElement->NbNodes();
1169             if ( iNbNodes!=4 ) {
1170               ok = 0;
1171               break ;
1172             }
1173           }
1174         }
1175         if (i==1) has_only_quad_f1 = ok ;
1176         if (i==2) has_only_quad_f2 = ok ;
1177         if (i==3) has_only_quad_f3 = ok ;
1178         if (i==4) has_only_quad_f4 = ok ;
1179         if (i==5) has_only_quad_f5 = ok ;
1180         if (i==6) has_only_quad_f6 = ok ;
1181       }
1182       //
1183       TopTools_IndexedMapOfShape aE;
1184       TopExp::MapShapes(myShape, TopAbs_EDGE, aE);
1185       int iNbE = aE.Extent();
1186       if (iNbE == 12) {
1187         //
1188         int nb_e01 = pMesh->GetSubMeshContaining(aE(1))->GetSubMeshDS()->NbElements();
1189         int nb_e02 = pMesh->GetSubMeshContaining(aE(2))->GetSubMeshDS()->NbElements();
1190         int nb_e03 = pMesh->GetSubMeshContaining(aE(3))->GetSubMeshDS()->NbElements();
1191         int nb_e04 = pMesh->GetSubMeshContaining(aE(4))->GetSubMeshDS()->NbElements();
1192         int nb_e05 = pMesh->GetSubMeshContaining(aE(5))->GetSubMeshDS()->NbElements();
1193         int nb_e06 = pMesh->GetSubMeshContaining(aE(6))->GetSubMeshDS()->NbElements();
1194         int nb_e07 = pMesh->GetSubMeshContaining(aE(7))->GetSubMeshDS()->NbElements();
1195         int nb_e08 = pMesh->GetSubMeshContaining(aE(8))->GetSubMeshDS()->NbElements();
1196         int nb_e09 = pMesh->GetSubMeshContaining(aE(9))->GetSubMeshDS()->NbElements();
1197         int nb_e10 = pMesh->GetSubMeshContaining(aE(10))->GetSubMeshDS()->NbElements();
1198         int nb_e11 = pMesh->GetSubMeshContaining(aE(11))->GetSubMeshDS()->NbElements();
1199         int nb_e12 = pMesh->GetSubMeshContaining(aE(12))->GetSubMeshDS()->NbElements();
1200         //
1201         int nb_ok = 0 ;
1202         //
1203         if ( (nb_e01==nb_e03) && (nb_e03==nb_e05) && (nb_e05==nb_e07) ) {
1204           if ( has_only_quad_f1 && has_only_quad_f2 && has_only_quad_f3 && has_only_quad_f4 ) {
1205             if ( (nb_e09==nb_e10) && (nb_e08==nb_e06) && (nb_e11==nb_e12) && (nb_e04==nb_e02) ) {
1206               if (nb_f5==nb_f6) {
1207                 nb_ok += 1;
1208                 aFTr = aM(5);
1209               }
1210             }
1211           }
1212         }
1213         if ( (nb_e02==nb_e04) && (nb_e04==nb_e06) && (nb_e06==nb_e08) ) {
1214           if ( has_only_quad_f1 && has_only_quad_f2 && has_only_quad_f5 && has_only_quad_f6 ) {
1215             if ( (nb_e01==nb_e03) && (nb_e10==nb_e12) && (nb_e05==nb_e07) && (nb_e09==nb_e11) ) {
1216               if (nb_f3==nb_f4) {
1217                 nb_ok += 1;
1218                 aFTr = aM(3);
1219               }
1220             }
1221           }
1222         }
1223         if ( (nb_e09==nb_e10) && (nb_e10==nb_e11) && (nb_e11==nb_e12) ) {
1224           if ( has_only_quad_f3 && has_only_quad_f4 && has_only_quad_f5 && has_only_quad_f6 ) {
1225             if ( (nb_e01==nb_e05) && (nb_e02==nb_e06) && (nb_e03==nb_e07) && (nb_e04==nb_e08) ) {
1226               if (nb_f1==nb_f2) {
1227                 nb_ok += 1;
1228                 aFTr = aM(1);
1229               }
1230             }
1231           }
1232         }
1233         //
1234         if ( nb_ok == 1 ) {
1235           isOK = 1;
1236         }
1237         //
1238       }
1239     }
1240     if (!isOK) {
1241       myErrorStatus->myName=5; // more than one face has triangulation
1242       myErrorStatus->myComment="Incorrect input mesh";
1243       return;
1244     }
1245   }
1246   // \end{E.A.}
1247   //
1248   // 1. Vetrices V00, V001;
1249   //
1250   TopExp::MapShapes(aFTr, TopAbs_EDGE, aME);
1251   TopExp::MapShapesAndAncestors(myShape, TopAbs_VERTEX, TopAbs_EDGE, aMVES);
1252   //
1253   // 1.1 Base vertex V000
1254   iNbE = aME.Extent();
1255   if (iNbE!= NB_WALL_FACES ){
1256     MESSAGE("StdMeshers_Penta_3D::MakeBlock() ");
1257     myErrorStatus->myName=7; // too few edges are in base face aFTr
1258     myErrorStatus->myComment=SMESH_Comment("Not a quadrilateral face #")
1259       <<pMesh->GetMeshDS()->ShapeToIndex( aFTr )<<": "<<iNbE<<" edges" ;
1260     return;
1261   }
1262   const TopoDS_Edge& aE1=TopoDS::Edge(aME(1));
1263   aV000=TopExp::FirstVertex(aE1);
1264   //
1265   const TopTools_ListOfShape& aLE=aMVES.FindFromKey(aV000);
1266   aIt.Initialize(aLE);
1267   for (; aIt.More(); aIt.Next()) {
1268     const TopoDS_Shape& aEx=aIt.Value();
1269     aMEV.Add(aEx);
1270   }
1271   iNbEV=aMEV.Extent();
1272   if (iNbEV!=3){
1273     MESSAGE("StdMeshers_Penta_3D::MakeBlock() ");
1274     myErrorStatus->myName=7; // too few edges meet in base vertex
1275     myErrorStatus->myComment=SMESH_Comment("3 edges must share vertex #")
1276       <<pMesh->GetMeshDS()->ShapeToIndex( aV000 )<<" but there are "<<iNbEV<<" edges";
1277     return;
1278   }
1279   //
1280   // 1.2 Vertex V001
1281   bFound=false;
1282   for (j=1; j<=iNbEV; ++j) {
1283     const TopoDS_Shape& aEx=aMEV(j);
1284     if (!aME.Contains(aEx)) {
1285       TopoDS_Vertex aV[2];
1286       //
1287       const TopoDS_Edge& aE=TopoDS::Edge(aEx);
1288       TopExp::Vertices(aE, aV[0], aV[1]);
1289       for (i=0; i<2; ++i) {
1290         if (!aV[i].IsSame(aV000)) {
1291           aV001=aV[i];
1292           bFound=!bFound;
1293           break;
1294         }
1295       }
1296     }
1297   }
1298   //
1299   if (!bFound) {
1300     MESSAGE("StdMeshers_Penta_3D::MakeBlock() ");
1301     myErrorStatus->myName=8; // can not find reper V001
1302     myErrorStatus->myComment=SMESH_Comment("Can't find opposite vertex for vertex #")
1303       <<pMesh->GetMeshDS()->ShapeToIndex( aV000 );
1304     return;
1305   }
1306   //DEB
1307   //gp_Pnt aP000, aP001;
1308   //
1309   //aP000=BRep_Tool::Pnt(TopoDS::Vertex(aV000));
1310   //printf("*** aP000 { %lf, %lf, %lf }\n", aP000.X(), aP000.Y(), aP000.Z());
1311   //aP001=BRep_Tool::Pnt(TopoDS::Vertex(aV001));
1312   //printf("*** aP001 { %lf, %lf, %lf }\n", aP001.X(), aP001.Y(), aP001.Z());
1313   //DEB
1314   //
1315   aME.Clear();
1316   TopExp::MapShapes(myShape, TopAbs_SHELL, aME);
1317   iNbE=aME.Extent();
1318   if (iNbE!=1) {
1319     MESSAGE("StdMeshers_Penta_3D::MakeBlock() ");
1320     myErrorStatus->myName=9; // number of shells in source shape !=1
1321     myErrorStatus->myComment=SMESH_Comment("Unexpected nb of shells ")<<iNbE;
1322     return;
1323   }
1324   //
1325   // 2. Load Block
1326   const TopoDS_Shell& aShell=TopoDS::Shell(aME(1));
1327   myBlock.Load(aShell, aV000, aV001);
1328   iErr = myBlock.ErrorStatus();
1329   if (iErr) {
1330     MESSAGE("StdMeshers_Penta_3D::MakeBlock() ");
1331     myErrorStatus=myBlock.GetError(); // SMESHBlock: Load operation failed
1332     return;
1333   }
1334 }
1335 //=======================================================================
1336 //function : CheckData
1337 //purpose  :
1338 //=======================================================================
1339 void StdMeshers_Penta_3D::CheckData()
1340 {
1341   int i, iNb;
1342   int iNbEx[]={8, 12, 6};
1343   //
1344   TopAbs_ShapeEnum aST;
1345   TopAbs_ShapeEnum aSTEx[]={
1346     TopAbs_VERTEX, TopAbs_EDGE, TopAbs_FACE
1347   };
1348   TopTools_IndexedMapOfShape aM;
1349   //
1350   if (myShape.IsNull()){
1351     MESSAGE("StdMeshers_Penta_3D::CheckData() ");
1352     myErrorStatus->myName=2; // null shape
1353     myErrorStatus->myComment="Null shape";
1354     return;
1355   }
1356   //
1357   aST=myShape.ShapeType();
1358   if (!(aST==TopAbs_SOLID || aST==TopAbs_SHELL)) {
1359     MESSAGE("StdMeshers_Penta_3D::CheckData() ");
1360     myErrorStatus->myName=3; // not compatible type of shape
1361     myErrorStatus->myComment=SMESH_Comment("Wrong shape type (TopAbs_ShapeEnum) ")<<aST;
1362     return;
1363   }
1364   //
1365   for (i=0; i<3; ++i) {
1366     aM.Clear();
1367     TopExp::MapShapes(myShape, aSTEx[i], aM);
1368     iNb=aM.Extent();
1369     if (iNb!=iNbEx[i]){
1370       MESSAGE("StdMeshers_Penta_3D::CheckData() ");
1371       myErrorStatus->myName=4; // number of subshape is not compatible
1372       myErrorStatus->myComment="Wrong number of subshapes of a block";
1373       return;
1374     }
1375   }
1376 }
1377
1378 //=======================================================================
1379 //function : LoadIJNodes
1380 //purpose  : Load nodes bound to theFace into column (vectors) and rows
1381 //           of theIJNodes.
1382 //           The value of theIJNodes map is a vector of ordered nodes so
1383 //           that the 0-the one lies on theBaseEdge.
1384 //           The key of theIJNodes map is a normalized parameter of each
1385 //           0-the node on theBaseEdge.
1386 //=======================================================================
1387
1388 bool StdMeshers_Penta_3D::LoadIJNodes(StdMeshers_IJNodeMap & theIJNodes,
1389                                       const TopoDS_Face&     theFace,
1390                                       const TopoDS_Edge&     theBaseEdge,
1391                                       SMESHDS_Mesh*          theMesh)
1392 {
1393   // get vertices of theBaseEdge
1394   TopoDS_Vertex vfb, vlb, vft; // first and last, bottom and top vertices
1395   TopoDS_Edge eFrw = TopoDS::Edge( theBaseEdge.Oriented( TopAbs_FORWARD ));
1396   TopExp::Vertices( eFrw, vfb, vlb );
1397
1398   // find the other edges of theFace and orientation of e1
1399   TopoDS_Edge e1, e2, eTop;
1400   bool rev1, CumOri = false;
1401   TopExp_Explorer exp( theFace, TopAbs_EDGE );
1402   int nbEdges = 0;
1403   for ( ; exp.More(); exp.Next() ) {
1404     if ( ++nbEdges > NB_WALL_FACES ) {
1405       return false; // more than 4 edges in theFace
1406     }
1407     TopoDS_Edge e = TopoDS::Edge( exp.Current() );
1408     if ( theBaseEdge.IsSame( e ))
1409       continue;
1410     TopoDS_Vertex vCommon;
1411     if ( !TopExp::CommonVertex( theBaseEdge, e, vCommon ))
1412       eTop = e;
1413     else if ( vCommon.IsSame( vfb )) {
1414       e1 = e;
1415       vft = TopExp::LastVertex( e1, CumOri );
1416       rev1 = vfb.IsSame( vft );
1417       if ( rev1 )
1418         vft = TopExp::FirstVertex( e1, CumOri );
1419     }
1420     else
1421       e2 = e;
1422   }
1423   if ( nbEdges < NB_WALL_FACES ) {
1424     return false; // less than 4 edges in theFace
1425   }
1426
1427   // submeshes corresponding to shapes
1428   SMESHDS_SubMesh* smFace = theMesh->MeshElements( theFace );
1429   SMESHDS_SubMesh* smb = theMesh->MeshElements( theBaseEdge );
1430   SMESHDS_SubMesh* smt = theMesh->MeshElements( eTop );
1431   SMESHDS_SubMesh* sm1 = theMesh->MeshElements( e1 );
1432   SMESHDS_SubMesh* sm2 = theMesh->MeshElements( e2 );
1433   SMESHDS_SubMesh* smVfb = theMesh->MeshElements( vfb );
1434   SMESHDS_SubMesh* smVlb = theMesh->MeshElements( vlb );
1435   SMESHDS_SubMesh* smVft = theMesh->MeshElements( vft );
1436   if (!smFace || !smb || !smt || !sm1 || !sm2 || !smVfb || !smVlb || !smVft ) {
1437     MESSAGE( "NULL submesh " <<smFace<<" "<<smb<<" "<<smt<<" "<<
1438             sm1<<" "<<sm2<<" "<<smVfb<<" "<<smVlb<<" "<<smVft);
1439     return false;
1440   }
1441   if ( smb->NbNodes() != smt->NbNodes() || sm1->NbNodes() != sm2->NbNodes() ) {
1442     MESSAGE(" Diff nb of nodes on opposite edges" );
1443     return false;
1444   }
1445   if (smVfb->NbNodes() != 1 || smVlb->NbNodes() != 1 || smVft->NbNodes() != 1) {
1446     MESSAGE("Empty submesh of vertex");
1447     return false;
1448   }
1449   if ( sm1->NbNodes() * smb->NbNodes() != smFace->NbNodes() ) {
1450     // check quadratic case
1451     if ( myCreateQuadratic ) {
1452       int n1 = sm1->NbNodes()/2;
1453       int n2 = smb->NbNodes()/2;
1454       int n3 = sm1->NbNodes() - n1;
1455       int n4 = smb->NbNodes() - n2;
1456       int nf = sm1->NbNodes()*smb->NbNodes() - n3*n4;
1457       if( nf != smFace->NbNodes() ) {
1458         MESSAGE( "Wrong nb face nodes: " <<
1459                 sm1->NbNodes()<<" "<<smb->NbNodes()<<" "<<smFace->NbNodes());
1460         return false;
1461       }
1462     }
1463     else {
1464       MESSAGE( "Wrong nb face nodes: " <<
1465               sm1->NbNodes()<<" "<<smb->NbNodes()<<" "<<smFace->NbNodes());
1466       return false;
1467     }
1468   }
1469   // IJ size
1470   int vsize = sm1->NbNodes() + 2;
1471   int hsize = smb->NbNodes() + 2;
1472   if(myCreateQuadratic) {
1473     vsize = vsize - sm1->NbNodes()/2 -1;
1474     hsize = hsize - smb->NbNodes()/2 -1;
1475   }
1476
1477   // load nodes from theBaseEdge
1478
1479   set<const SMDS_MeshNode*> loadedNodes;
1480   const SMDS_MeshNode* nullNode = 0;
1481
1482   vector<const SMDS_MeshNode*> & nVecf = theIJNodes[ 0.];
1483   nVecf.resize( vsize, nullNode );
1484   loadedNodes.insert( nVecf[ 0 ] = smVfb->GetNodes()->next() );
1485
1486   vector<const SMDS_MeshNode*> & nVecl = theIJNodes[ 1.];
1487   nVecl.resize( vsize, nullNode );
1488   loadedNodes.insert( nVecl[ 0 ] = smVlb->GetNodes()->next() );
1489
1490   double f, l;
1491   BRep_Tool::Range( eFrw, f, l );
1492   double range = l - f;
1493   SMDS_NodeIteratorPtr nIt = smb->GetNodes();
1494   const SMDS_MeshNode* node;
1495   while ( nIt->more() ) {
1496     node = nIt->next();
1497     if(myTool->IsMedium(node))
1498       continue;
1499     const SMDS_EdgePosition* pos =
1500       dynamic_cast<const SMDS_EdgePosition*>( node->GetPosition().get() );
1501     if ( !pos ) {
1502       return false;
1503     }
1504     double u = ( pos->GetUParameter() - f ) / range;
1505     vector<const SMDS_MeshNode*> & nVec = theIJNodes[ u ];
1506     nVec.resize( vsize, nullNode );
1507     loadedNodes.insert( nVec[ 0 ] = node );
1508   }
1509   if ( theIJNodes.size() != hsize ) {
1510     MESSAGE( "Wrong node positions on theBaseEdge" );
1511     return false;
1512   }
1513
1514   // load nodes from e1
1515
1516   map< double, const SMDS_MeshNode*> sortedNodes; // sort by param on edge
1517   nIt = sm1->GetNodes();
1518   while ( nIt->more() ) {
1519     node = nIt->next();
1520     if(myTool->IsMedium(node))
1521       continue;
1522     const SMDS_EdgePosition* pos =
1523       dynamic_cast<const SMDS_EdgePosition*>( node->GetPosition().get() );
1524     if ( !pos ) {
1525       return false;
1526     }
1527     sortedNodes.insert( make_pair( pos->GetUParameter(), node ));
1528   }
1529   loadedNodes.insert( nVecf[ vsize - 1 ] = smVft->GetNodes()->next() );
1530   map< double, const SMDS_MeshNode*>::iterator u_n = sortedNodes.begin();
1531   int row = rev1 ? vsize - 1 : 0;
1532   for ( ; u_n != sortedNodes.end(); u_n++ ) {
1533     if ( rev1 ) row--;
1534     else        row++;
1535     loadedNodes.insert( nVecf[ row ] = u_n->second );
1536   }
1537
1538   // try to load the rest nodes
1539
1540   // get all faces from theFace
1541   TIDSortedElemSet allFaces, foundFaces;
1542   SMDS_ElemIteratorPtr eIt = smFace->GetElements();
1543   while ( eIt->more() ) {
1544     const SMDS_MeshElement* e = eIt->next();
1545     if ( e->GetType() == SMDSAbs_Face )
1546       allFaces.insert( e );
1547   }
1548   // Starting from 2 neighbour nodes on theBaseEdge, look for a face
1549   // the nodes belong to, and between the nodes of the found face,
1550   // look for a not loaded node considering this node to be the next
1551   // in a column of the starting second node. Repeat, starting
1552   // from nodes next to the previous starting nodes in their columns,
1553   // and so on while a face can be found. Then go the the next pair
1554   // of nodes on theBaseEdge.
1555   StdMeshers_IJNodeMap::iterator par_nVec_1 = theIJNodes.begin();
1556   StdMeshers_IJNodeMap::iterator par_nVec_2 = par_nVec_1;
1557   // loop on columns
1558   int col = 0;
1559   for ( par_nVec_2++; par_nVec_2 != theIJNodes.end(); par_nVec_1++, par_nVec_2++ ) {
1560     col++;
1561     row = 0;
1562     const SMDS_MeshNode* n1 = par_nVec_1->second[ row ];
1563     const SMDS_MeshNode* n2 = par_nVec_2->second[ row ];
1564     const SMDS_MeshElement* face = 0;
1565     do {
1566       // look for a face by 2 nodes
1567       face = SMESH_MeshEditor::FindFaceInSet( n1, n2, allFaces, foundFaces );
1568       if ( face ) {
1569         int nbFaceNodes = face->NbNodes();
1570         if ( (!myCreateQuadratic && nbFaceNodes>4) ||
1571              (myCreateQuadratic && nbFaceNodes>8) ) {
1572           MESSAGE(" Too many nodes in a face: " << nbFaceNodes );
1573           return false;
1574         }
1575         // look for a not loaded node of the <face>
1576         bool found = false;
1577         const SMDS_MeshNode* n3 = 0; // a node defferent from n1 and n2
1578         eIt = face->nodesIterator() ;
1579         while ( !found && eIt->more() ) {
1580           node = static_cast<const SMDS_MeshNode*>( eIt->next() );
1581           if(myTool->IsMedium(node))
1582             continue;
1583           found = loadedNodes.insert( node ).second;
1584           if ( !found && node != n1 && node != n2 )
1585             n3 = node;
1586         }
1587         if ( found ) {
1588           if ( ++row > vsize - 1 ) {
1589             MESSAGE( "Too many nodes in column "<< col <<": "<< row+1);
1590             return false;
1591           }
1592           par_nVec_2->second[ row ] = node;
1593           foundFaces.insert( face );
1594           n2 = node;
1595           if ( nbFaceNodes==4 || (myCreateQuadratic && nbFaceNodes==8) ) {
1596             n1 = par_nVec_1->second[ row ];
1597           }
1598         }
1599         else if ( (nbFaceNodes==3 || (myCreateQuadratic && nbFaceNodes==6) )  &&
1600                  n3 == par_nVec_1->second[ row ] ) {
1601           n1 = n3;
1602         }
1603         else {
1604           MESSAGE( "Not quad mesh, column "<< col );
1605           return false;
1606         }
1607       }
1608     }
1609     while ( face && n1 && n2 );
1610
1611     if ( row < vsize - 1 ) {
1612       MESSAGE( "Too few nodes in column "<< col <<": "<< row+1);
1613       MESSAGE( "Base node 1: "<< par_nVec_1->second[0]);
1614       MESSAGE( "Base node 2: "<< par_nVec_2->second[0]);
1615       MESSAGE( "Current node 1: "<< n1);
1616       MESSAGE( "Current node 2: "<< n2);
1617       MESSAGE( "first base node: "<< theIJNodes.begin()->second[0]);
1618       MESSAGE( "last base node: "<< theIJNodes.rbegin()->second[0]);
1619       return false;
1620     }
1621   } // loop on columns
1622
1623   return true;
1624 }
1625
1626 //////////////////////////////////////////////////////////////////////////
1627 //
1628 //   StdMeshers_SMESHBlock
1629 //
1630 //////////////////////////////////////////////////////////////////////////
1631
1632 //=======================================================================
1633 //function : StdMeshers_SMESHBlock
1634 //purpose  :
1635 //=======================================================================
1636 StdMeshers_SMESHBlock::StdMeshers_SMESHBlock()
1637 {
1638   myErrorStatus=1;
1639   myIsEdgeForward.resize( SMESH_Block::NbEdges(), -1 );
1640 }
1641
1642 //=======================================================================
1643 //function : IsForwadEdge
1644 //purpose  :
1645 //=======================================================================
1646
1647 bool StdMeshers_SMESHBlock::IsForwadEdge(const int theEdgeID)
1648 {
1649   int index = myTBlock.ShapeIndex( theEdgeID );
1650   if ( !myTBlock.IsEdgeID( theEdgeID ))
1651     return false;
1652
1653   if ( myIsEdgeForward[ index ] < 0 )
1654     myIsEdgeForward[ index ] =
1655       myTBlock.IsForwardEdge( TopoDS::Edge( Shape( theEdgeID )), myShapeIDMap );
1656
1657   return myIsEdgeForward[ index ];
1658 }
1659
1660 //=======================================================================
1661 //function : ErrorStatus
1662 //purpose  :
1663 //=======================================================================
1664 int StdMeshers_SMESHBlock::ErrorStatus() const
1665 {
1666   return myErrorStatus;
1667 }
1668
1669 //================================================================================
1670 /*!
1671  * \brief Return problem description
1672  */
1673 //================================================================================
1674
1675 SMESH_ComputeErrorPtr StdMeshers_SMESHBlock::GetError() const
1676 {
1677   SMESH_ComputeErrorPtr err = SMESH_ComputeError::New();
1678   string & text = err->myComment;
1679   switch ( myErrorStatus ) {
1680   case 2:
1681   case 3: text = "Internal error of StdMeshers_Penta_3D"; break;
1682   case 4: text = "Can't compute normalized parameters of a point inside a block"; break;
1683   case 5: text = "Can't compute coordinates by normalized parameters inside a block"; break;
1684   case 6: text = "Can't detect block subshapes. Not a block?"; break;
1685   }
1686   if (!text.empty())
1687     err->myName = myErrorStatus;
1688   return err;
1689 }
1690
1691 //=======================================================================
1692 //function : Load
1693 //purpose  :
1694 //=======================================================================
1695 void StdMeshers_SMESHBlock::Load(const TopoDS_Shell& theShell)
1696 {
1697   TopoDS_Vertex aV000, aV001;
1698   //
1699   Load(theShell, aV000, aV001);
1700 }
1701
1702 //=======================================================================
1703 //function : Load
1704 //purpose  :
1705 //=======================================================================
1706 void StdMeshers_SMESHBlock::Load(const TopoDS_Shell& theShell,
1707                                  const TopoDS_Vertex& theV000,
1708                                  const TopoDS_Vertex& theV001)
1709 {
1710   myErrorStatus=0;
1711   //
1712   myShell=theShell;
1713   //
1714   bool bOk;
1715   //
1716   myShapeIDMap.Clear();
1717   bOk = myTBlock.LoadBlockShapes(myShell, theV000, theV001, myShapeIDMap);
1718   if (!bOk) {
1719     myErrorStatus=6;
1720     return;
1721   }
1722 }
1723
1724 //=======================================================================
1725 //function : ComputeParameters
1726 //purpose  :
1727 //=======================================================================
1728 void StdMeshers_SMESHBlock::ComputeParameters(const gp_Pnt& thePnt,
1729                                               gp_XYZ& theXYZ)
1730 {
1731   ComputeParameters(thePnt, myShell, theXYZ);
1732 }
1733
1734 //=======================================================================
1735 //function : ComputeParameters
1736 //purpose  :
1737 //=======================================================================
1738 void StdMeshers_SMESHBlock::ComputeParameters(const gp_Pnt& thePnt,
1739                                               const TopoDS_Shape& theShape,
1740                                               gp_XYZ& theXYZ)
1741 {
1742   myErrorStatus=0;
1743   //
1744   int aID;
1745   bool bOk;
1746   //
1747   aID = ShapeID(theShape);
1748   if (myErrorStatus) {
1749     return;
1750   }
1751   bOk = myTBlock.ComputeParameters(thePnt, theXYZ, aID);
1752   if (!bOk) {
1753     myErrorStatus=4; // problems with computation Parameters
1754     return;
1755   }
1756 }
1757
1758 //=======================================================================
1759 //function : ComputeParameters
1760 //purpose  :
1761 //=======================================================================
1762
1763 void StdMeshers_SMESHBlock::ComputeParameters(const double& theU,
1764                                               const TopoDS_Shape& theShape,
1765                                               gp_XYZ& theXYZ)
1766 {
1767   myErrorStatus=0;
1768   //
1769   int aID;
1770   bool bOk=false;
1771   //
1772   aID = ShapeID(theShape);
1773   if (myErrorStatus) {
1774     return;
1775   }
1776   if ( SMESH_Block::IsEdgeID( aID ))
1777       bOk = myTBlock.EdgeParameters( aID, theU, theXYZ );
1778   if (!bOk) {
1779     myErrorStatus=4; // problems with computation Parameters
1780     return;
1781   }
1782 }
1783
1784 //=======================================================================
1785 //function : Point
1786 //purpose  :
1787 //=======================================================================
1788  void StdMeshers_SMESHBlock::Point(const gp_XYZ& theParams,
1789                                    gp_Pnt& aP3D)
1790 {
1791   TopoDS_Shape aS;
1792   //
1793   Point(theParams, aS, aP3D);
1794 }
1795
1796 //=======================================================================
1797 //function : Point
1798 //purpose  :
1799 //=======================================================================
1800  void StdMeshers_SMESHBlock::Point(const gp_XYZ& theParams,
1801                                    const TopoDS_Shape& theShape,
1802                                    gp_Pnt& aP3D)
1803 {
1804   myErrorStatus = 0;
1805   //
1806   int aID;
1807   bool bOk = false;
1808   gp_XYZ aXYZ(99.,99.,99.);
1809   aP3D.SetXYZ(aXYZ);
1810   //
1811   if (theShape.IsNull()) {
1812     bOk = myTBlock.ShellPoint(theParams, aXYZ);
1813   }
1814   //
1815   else {
1816     aID=ShapeID(theShape);
1817     if (myErrorStatus) {
1818       return;
1819     }
1820     //
1821     if (SMESH_Block::IsVertexID(aID)) {
1822       bOk = myTBlock.VertexPoint(aID, aXYZ);
1823     }
1824     else if (SMESH_Block::IsEdgeID(aID)) {
1825       bOk = myTBlock.EdgePoint(aID, theParams, aXYZ);
1826     }
1827     //
1828     else if (SMESH_Block::IsFaceID(aID)) {
1829       bOk = myTBlock.FacePoint(aID, theParams, aXYZ);
1830     }
1831   }
1832   if (!bOk) {
1833     myErrorStatus=5; // problems with point computation
1834     return;
1835   }
1836   aP3D.SetXYZ(aXYZ);
1837 }
1838
1839 //=======================================================================
1840 //function : ShapeID
1841 //purpose  :
1842 //=======================================================================
1843 int StdMeshers_SMESHBlock::ShapeID(const TopoDS_Shape& theShape)
1844 {
1845   myErrorStatus=0;
1846   //
1847   int aID=-1;
1848   TopoDS_Shape aSF, aSR;
1849   //
1850   aSF=theShape;
1851   aSF.Orientation(TopAbs_FORWARD);
1852   aSR=theShape;
1853   aSR.Orientation(TopAbs_REVERSED);
1854   //
1855   if (myShapeIDMap.Contains(aSF)) {
1856     aID=myShapeIDMap.FindIndex(aSF);
1857     return aID;
1858   }
1859   if (myShapeIDMap.Contains(aSR)) {
1860     aID=myShapeIDMap.FindIndex(aSR);
1861     return aID;
1862   }
1863   myErrorStatus=2; // unknown shape;
1864   return aID;
1865 }
1866
1867 //=======================================================================
1868 //function : Shape
1869 //purpose  :
1870 //=======================================================================
1871 const TopoDS_Shape& StdMeshers_SMESHBlock::Shape(const int theID)
1872 {
1873   myErrorStatus=0;
1874   //
1875   int aNb;
1876   //
1877   aNb=myShapeIDMap.Extent();
1878   if (theID<1 || theID>aNb) {
1879     myErrorStatus=3; // ID is out of range
1880     return myEmptyShape;
1881   }
1882   //
1883   const TopoDS_Shape& aS=myShapeIDMap.FindKey(theID);
1884   return aS;
1885 }
1886
1887