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