src/StdMeshers/StdMeshers_Penta_3D.cxx

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