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