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_MeshEditor.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   SMESH_MesherHelper helper(aMesh);
 107   myTool = &helper;
 108   myCreateQuadratic = myTool->IsQuadraticSubMesh(aShape);
 109
 110   //
 111   MakeBlock();
 112   if (!myErrorStatus->IsOK()) {
 113     return bOK;
 114   }
 115   //
 116   ClearMeshOnFxy1();
 117   if (!myErrorStatus->IsOK()) {
 118     return bOK;
 119   }
 120   //
 121   MakeNodes();
 122   if (!myErrorStatus->IsOK()) {
 123     return bOK;
 124   }
 125   //
 126   MakeConnectingMap();
 127   //
 128   MakeMeshOnFxy1();
 129   if (!myErrorStatus->IsOK()) {
 130     return bOK;
 131   }
 132   //
 133   MakeVolumeMesh();
 134   //
 135   return !bOK;
 136 }
 137
 138 //=======================================================================
 139 //function : MakeNodes
 140 //purpose  :
 141 //=======================================================================
 142 void StdMeshers_Penta_3D::MakeNodes()
 143 {
 144   const int aNbSIDs=9;
 145   int i, j, k, ij, iNbN, aNodeID, aSize, iErr;
 146   double aX, aY, aZ;
 147   SMESH_Block::TShapeID aSID, aSIDs[aNbSIDs]={
 148     SMESH_Block::ID_V000, SMESH_Block::ID_V100,
 149     SMESH_Block::ID_V110, SMESH_Block::ID_V010,
 150     SMESH_Block::ID_Ex00, SMESH_Block::ID_E1y0,
 151     SMESH_Block::ID_Ex10, SMESH_Block::ID_E0y0,
 152     SMESH_Block::ID_Fxy0
 153   };
 154   //
 155   SMESH_Mesh* pMesh=GetMesh();
 156   //
 157   // 1. Define the sizes of mesh
 158   //
 159   // 1.1 Horizontal size
 160   myJSize=0;
 161   for (i=0; i<aNbSIDs; ++i) {
 162     const TopoDS_Shape& aS = myBlock.Shape(aSIDs[i]);
 163     SMESH_subMesh *aSubMesh = pMesh->GetSubMeshContaining(aS);
 164     ASSERT(aSubMesh);
 165     SMESHDS_SubMesh *aSM = aSubMesh->GetSubMeshDS();
 166     if(!myCreateQuadratic) {
 167       iNbN = aSM->NbNodes();
 168     }
 169     else {
 170       iNbN = 0;
 171       SMDS_NodeIteratorPtr itn = aSM->GetNodes();
 172       while(itn->more()) {
 173         const SMDS_MeshNode* aNode = itn->next();
 174         if(myTool->IsMedium(aNode))
 175           continue;
 176         iNbN++;
 177       }
 178     }
 179     myJSize += iNbN;
 180   }
 181   //printf("***  Horizontal: number of nodes summary=%d\n", myJSize);
 182   //
 183   // 1.2 Vertical size
 184   myISize=2;
 185   {
 186     const TopoDS_Shape& aS=myBlock.Shape(SMESH_Block::ID_E00z);
 187     SMESH_subMesh *aSubMesh = pMesh->GetSubMeshContaining(aS);
 188     ASSERT(aSubMesh);
 189     SMESHDS_SubMesh *aSM = aSubMesh->GetSubMeshDS();
 190     if(!myCreateQuadratic) {
 191       iNbN = aSM->NbNodes();
 192     }
 193     else {
 194       iNbN = 0;
 195       SMDS_NodeIteratorPtr itn = aSM->GetNodes();
 196       while(itn->more()) {
 197         const SMDS_MeshNode* aNode = itn->next();
 198         if(myTool->IsMedium(aNode))
 199           continue;
 200         iNbN++;
 201       }
 202     }
 203     myISize += iNbN;
 204   }
 205   //printf("***  Vertical: number of nodes on edges and vertices=%d\n", myISize);
 206   //
 207   aSize=myISize*myJSize;
 208   myTNodes.resize(aSize);
 209   //
 210   StdMeshers_TNode aTNode;
 211   gp_XYZ aCoords;
 212   gp_Pnt aP3D;
 213   //
 214   // 2. Fill the repers on base face (Z=0)
 215   i=0; j=0;
 216   // vertices
 217   for (k=0; k<aNbSIDs; ++k) {
 218     aSID=aSIDs[k];
 219     const TopoDS_Shape& aS = myBlock.Shape(aSID);
 220     SMDS_NodeIteratorPtr ite = pMesh->GetSubMeshContaining(aS)->GetSubMeshDS()->GetNodes();
 221     while(ite->more()) {
 222       const SMDS_MeshNode* aNode = ite->next();
 223       if(myTool->IsMedium(aNode))
 224         continue;
 225       aNodeID=aNode->GetID();
 226       //
 227       aTNode.SetNode(aNode);
 228       aTNode.SetShapeSupportID(aSID);
 229       aTNode.SetBaseNodeID(aNodeID);
 230       //
 231       if ( SMESH_Block::IsEdgeID (aSID)) {
 232         const SMDS_EdgePosition* epos =
 233           static_cast<const SMDS_EdgePosition*>(aNode->GetPosition().get());
 234         myBlock.ComputeParameters( epos->GetUParameter(), aS, aCoords );
 235       }
 236       else {
 237         aX=aNode->X();
 238         aY=aNode->Y();
 239         aZ=aNode->Z();
 240         aP3D.SetCoord(aX, aY, aZ);
 241         myBlock.ComputeParameters(aP3D, aS, aCoords);
 242       }
 243       iErr = myBlock.ErrorStatus();
 244       if (iErr) {
 245         MESSAGE("StdMeshers_Penta_3D::MakeNodes()," <<
 246                 "SMESHBlock: ComputeParameters operation failed");
 247         myErrorStatus=myBlock.GetError();
 248         return;
 249       }
 250       aTNode.SetNormCoord(aCoords);
 251       ij=i*myJSize+j;
 252       myTNodes[ij]=aTNode;
 253       ++j;
 254     }
 255   }
 256
 257   // 3.1 Fill maps of wall nodes
 258   SMESH_Block::TShapeID wallFaceID[ NB_WALL_FACES ] = {
 259     SMESH_Block::ID_Fx0z, SMESH_Block::ID_Fx1z,
 260     SMESH_Block::ID_F0yz, SMESH_Block::ID_F1yz
 261     };
 262   SMESH_Block::TShapeID baseEdgeID[ NB_WALL_FACES ] = {
 263     SMESH_Block::ID_Ex00, SMESH_Block::ID_Ex10,
 264     SMESH_Block::ID_E0y0, SMESH_Block::ID_E1y0
 265     };
 266   for ( i = 0; i < NB_WALL_FACES ; ++i ) {
 267     int fIndex = SMESH_Block::ShapeIndex( wallFaceID[ i ]);
 268     bool ok = LoadIJNodes (myWallNodesMaps[ fIndex ],
 269                            TopoDS::Face( myBlock.Shape( wallFaceID[ i ] )),
 270                            TopoDS::Edge( myBlock.Shape( baseEdgeID[ i ] )),
 271                            pMesh->GetMeshDS());
 272     if ( !ok ) {
 273       myErrorStatus->myName = COMPERR_BAD_INPUT_MESH;
 274       myErrorStatus->myComment = SMESH_Comment() <<
 275         "Can't find regular quadrangle mesh on a side face #" <<
 276         pMesh->GetMeshDS()->ShapeToIndex( myBlock.Shape( wallFaceID[ i ]));
 277       return;
 278     }
 279   }
 280
 281   // 3.2 find node columns for vertical edges and edge IDs
 282   vector<const SMDS_MeshNode*> * verticEdgeNodes[ NB_WALL_FACES ];
 283   SMESH_Block::TShapeID          verticEdgeID   [ NB_WALL_FACES ];
 284   for ( i = 0; i < NB_WALL_FACES ; ++i ) { // 4 first base nodes are nodes on vertices
 285     // edge ID
 286     SMESH_Block::TShapeID eID, vID = aSIDs[ i ];
 287     ShapeSupportID(false, vID, eID);
 288     verticEdgeID[ i ] = eID;
 289     // column nodes
 290     StdMeshers_TNode& aTNode = myTNodes[ i ];
 291     verticEdgeNodes[ i ] = 0;
 292     for ( j = 0; j < NB_WALL_FACES ; ++j ) { // loop on 4 wall faces
 293       int fIndex = SMESH_Block::ShapeIndex( wallFaceID[ j ]);
 294       StdMeshers_IJNodeMap & ijNodes= myWallNodesMaps[ fIndex ];
 295       if ( ijNodes.begin()->second[0] == aTNode.Node() )
 296         verticEdgeNodes[ i ] = & ijNodes.begin()->second;
 297       else if ( ijNodes.rbegin()->second[0] == aTNode.Node() )
 298         verticEdgeNodes[ i ] = & ijNodes.rbegin()->second;
 299       if ( verticEdgeNodes[ i ] )
 300         break;
 301     }
 302   }
 303
 304   // 3.3 set XYZ of vertices, and initialize of the rest
 305   SMESHDS_Mesh* aMesh = GetMesh()->GetMeshDS();
 306   for ( int id = SMESH_Block::ID_V000; id < SMESH_Block::ID_Shell; ++id ) {
 307     if ( SMESH_Block::IsVertexID( id )) {
 308       TopoDS_Shape V = myBlock.Shape( id );
 309       SMESHDS_SubMesh* sm = aMesh->MeshElements( V );
 310       const SMDS_MeshNode* n = sm->GetNodes()->next();
 311       myShapeXYZ[ id ].SetCoord( n->X(), n->Y(), n->Z() );
 312     }
 313     else
 314       myShapeXYZ[ id ].SetCoord( 0., 0., 0. );
 315   }
 316
 317
 318   // 4. Fill the rest repers
 319   bool bIsUpperLayer;
 320   int iBNID;
 321   SMESH_Block::TShapeID aSSID, aBNSSID;
 322   StdMeshers_TNode aTN;
 323   //
 324
 325   // create top face and find UV for it's corners
 326   const TopoDS_Face& TopFace = TopoDS::Face(myBlock.Shape(SMESH_Block::ID_Fxy1));
 327   SMESHDS_Mesh* meshDS = pMesh->GetMeshDS();
 328   int topfaceID = meshDS->ShapeToIndex(TopFace);
 329   const TopoDS_Vertex& v001 = TopoDS::Vertex(myBlock.Shape(SMESH_Block::ID_V001));
 330   SMDS_NodeIteratorPtr itn = pMesh->GetSubMeshContaining(v001)->GetSubMeshDS()->GetNodes();
 331   const SMDS_MeshNode* N = itn->next();
 332   gp_XY UV001 = myTool->GetNodeUV(TopFace,N);
 333   const TopoDS_Vertex& v101 = TopoDS::Vertex(myBlock.Shape(SMESH_Block::ID_V101));
 334   itn = pMesh->GetSubMeshContaining(v101)->GetSubMeshDS()->GetNodes();
 335   N = itn->next();
 336   gp_XY UV101 = myTool->GetNodeUV(TopFace,N);
 337   const TopoDS_Vertex& v011 = TopoDS::Vertex(myBlock.Shape(SMESH_Block::ID_V011));
 338   itn = pMesh->GetSubMeshContaining(v011)->GetSubMeshDS()->GetNodes();
 339   N = itn->next();
 340   gp_XY UV011 = myTool->GetNodeUV(TopFace,N);
 341   const TopoDS_Vertex& v111 = TopoDS::Vertex(myBlock.Shape(SMESH_Block::ID_V111));
 342   itn = pMesh->GetSubMeshContaining(v111)->GetSubMeshDS()->GetNodes();
 343   N = itn->next();
 344   gp_XY UV111 = myTool->GetNodeUV(TopFace,N);
 345
 346   for (j=0; j<myJSize; ++j) { // loop on all nodes of the base face (ID_Fxy0)
 347     // base node info
 348     const StdMeshers_TNode& aBN = myTNodes[j];
 349     aBNSSID = (SMESH_Block::TShapeID)aBN.ShapeSupportID();
 350     iBNID = aBN.BaseNodeID();
 351     const gp_XYZ& aBNXYZ = aBN.NormCoord();
 352     bool createNode = ( aBNSSID == SMESH_Block::ID_Fxy0 ); // if base node is inside a bottom face
 353     //
 354     // set XYZ on horizontal edges and get node columns of faces:
 355     // 2 columns for each face, between which a base node is located
 356     vector<const SMDS_MeshNode*>* nColumns[8];
 357     double ratio[ NB_WALL_FACES ]; // base node position between columns [0.-1.]
 358     if ( createNode ) {
 359       for ( k = 0; k < NB_WALL_FACES ; ++k ) {
 360         ratio[ k ] = SetHorizEdgeXYZ (aBNXYZ, wallFaceID[ k ],
 361                                       nColumns[k*2], nColumns[k*2+1]);
 362       }
 363     }
 364     //
 365     // XYZ on the bottom and top faces
 366     const SMDS_MeshNode* n = aBN.Node();
 367     myShapeXYZ[ SMESH_Block::ID_Fxy0 ].SetCoord( n->X(), n->Y(), n->Z() );
 368     myShapeXYZ[ SMESH_Block::ID_Fxy1 ].SetCoord( 0., 0., 0. );
 369     //
 370     // first create or find a top node, then the rest ones in a column
 371     for (i=myISize-1; i>0; --i) // vertical loop, from top to bottom
 372     {
 373       bIsUpperLayer = (i==(myISize-1));
 374       gp_XY UV_Ex01, UV_Ex11, UV_E0y1, UV_E1y1;
 375       if ( createNode ) // a base node is inside a top face
 376       {
 377         // set XYZ on vertical edges and faces
 378         for ( k = 0; k < NB_WALL_FACES ; ++k ) {
 379           // XYZ on a vertical edge
 380           const SMDS_MeshNode* n = (*verticEdgeNodes[ k ]) [ i ];
 381           myShapeXYZ[ verticEdgeID[ k ] ].SetCoord( n->X(), n->Y(), n->Z() );
 382           // XYZ on a face (part 1 from one column)
 383           n = (*nColumns[k*2]) [ i ];
 384           gp_XYZ xyz( n->X(), n->Y(), n->Z() );
 385           myShapeXYZ[ wallFaceID[ k ]] = ( 1. - ratio[ k ]) * xyz;
 386           gp_XY tmp1;
 387           if( bIsUpperLayer ) {
 388             tmp1 = myTool->GetNodeUV(TopFace,n);
 389             tmp1 = ( 1. - ratio[ k ]) * tmp1;
 390           }
 391           // XYZ on a face (part 2 from other column)
 392           n = (*nColumns[k*2+1]) [ i ];
 393           xyz.SetCoord( n->X(), n->Y(), n->Z() );
 394           myShapeXYZ[ wallFaceID[ k ]] += ratio[ k ] * xyz;
 395           if( bIsUpperLayer ) {
 396             gp_XY tmp2 = myTool->GetNodeUV(TopFace,n);
 397             tmp1 +=  ratio[ k ] * tmp2;
 398             if( k==0 )
 399               UV_Ex01 = tmp1;
 400             else if( k==1 )
 401               UV_Ex11 = tmp1;
 402             else if( k==2 )
 403               UV_E0y1 = tmp1;
 404             else
 405               UV_E1y1 = tmp1;
 406           }
 407         }
 408       }
 409       // fill current node info
 410       //   -index in aTNodes
 411       ij=i*myJSize+j;
 412       //   -normalized coordinates
 413       aX=aBNXYZ.X();
 414       aY=aBNXYZ.Y();
 415       //aZ=aZL[i];
 416       aZ=(double)i/(double)(myISize-1);
 417       aCoords.SetCoord(aX, aY, aZ);
 418       //
 419       //   suporting shape ID
 420       ShapeSupportID(bIsUpperLayer, aBNSSID, aSSID);
 421       if (!myErrorStatus->IsOK()) {
 422         MESSAGE("StdMeshers_Penta_3D::MakeNodes() ");
 423         return;
 424       }
 425       //
 426       aTN.SetShapeSupportID(aSSID);
 427       aTN.SetNormCoord(aCoords);
 428       aTN.SetBaseNodeID(iBNID);
 429       //
 430       if (aSSID!=SMESH_Block::ID_NONE){
 431         // try to find the node
 432         const TopoDS_Shape& aS=myBlock.Shape((int)aSSID);
 433         FindNodeOnShape(aS, aCoords, i, aTN);
 434       }
 435       else{
 436         // create node and get it id
 437         CreateNode (bIsUpperLayer, aCoords, aTN);
 438         //
 439         if ( bIsUpperLayer ) {
 440           const SMDS_MeshNode* n = aTN.Node();
 441           myShapeXYZ[ SMESH_Block::ID_Fxy1 ].SetCoord( n->X(), n->Y(), n->Z() );
 442           // set node on top face:
 443           // find UV parameter for this node
 444           //              UV_Ex11
 445           //   UV011+-----+----------+UV111
 446           //        |                |
 447           //        |                |
 448           // UV_E0y1+     +node      +UV_E1y1
 449           //        |                |
 450           //        |                |
 451           //        |                |
 452           //   UV001+-----+----------+UV101
 453           //              UV_Ex01
 454           gp_Pnt2d aP;
 455           double u = aCoords.X(), v = aCoords.Y();
 456           double u1 = ( 1. - u ), v1 = ( 1. - v );
 457           aP.ChangeCoord()  = UV_Ex01 * v1;
 458           aP.ChangeCoord() += UV_Ex11 * v;
 459           aP.ChangeCoord() += UV_E0y1 * u1;
 460           aP.ChangeCoord() += UV_E1y1 * u;
 461           aP.ChangeCoord() -= UV001 * u1 * v1;
 462           aP.ChangeCoord() -= UV101 * u  * v1;
 463           aP.ChangeCoord() -= UV011 * u1 * v;
 464           aP.ChangeCoord() -= UV111 * u  * v;
 465           meshDS->SetNodeOnFace((SMDS_MeshNode*)n, topfaceID, aP.X(), aP.Y());
 466         }
 467       }
 468       if (!myErrorStatus->IsOK()) {
 469         MESSAGE("StdMeshers_Penta_3D::MakeNodes() ");
 470         return;
 471       }
 472       //
 473       myTNodes[ij]=aTN;
 474     }
 475   }
 476   //DEB
 477   /*
 478   {
 479     int iSSID, iBNID, aID;
 480     //
 481     for (i=0; i<myISize; ++i) {
 482       printf(" Layer# %d\n", i);
 483       for (j=0; j<myJSize; ++j) {
 484         ij=i*myJSize+j;
 485         const StdMeshers_TNode& aTN=myTNodes[ij];
 486         //const StdMeshers_TNode& aTN=aTNodes[ij];
 487         const gp_XYZ& aXYZ=aTN.NormCoord();
 488         iSSID=aTN.ShapeSupportID();
 489         iBNID=aTN.BaseNodeID();
 490         //
 491         const SMDS_MeshNode* aNode=aTN.Node();
 492         aID=aNode->GetID();
 493         aX=aNode->X();
 494         aY=aNode->Y();
 495         aZ=aNode->Z();
 496         printf("*** j:%d BNID#%d iSSID:%d ID:%d { %lf %lf %lf },  { %lf %lf %lf }\n",
 497                j,  iBNID, iSSID, aID, aXYZ.X(),  aXYZ.Y(), aXYZ.Z(), aX, aY, aZ);
 498       }
 499     }
 500   }
 501   */
 502   //DEB t
 503 }
 504
 505
 506 //=======================================================================
 507 //function : FindNodeOnShape
 508 //purpose  :
 509 //=======================================================================
 510
 511 void StdMeshers_Penta_3D::FindNodeOnShape(const TopoDS_Shape& aS,
 512                                           const gp_XYZ&       aParams,
 513                                           const int           z,
 514                                           StdMeshers_TNode&   aTN)
 515 {
 516   double aX, aY, aZ, aD, aTol2, minD;
 517   gp_Pnt aP1, aP2;
 518   //
 519   SMESH_Mesh* pMesh = GetMesh();
 520   aTol2 = myTol3D*myTol3D;
 521   minD = 1.e100;
 522   SMDS_MeshNode* pNode = NULL;
 523   //
 524   if ( aS.ShapeType() == TopAbs_FACE ||
 525        aS.ShapeType() == TopAbs_EDGE ) {
 526     // find a face ID to which aTN belongs to
 527     int faceID;
 528     if ( aS.ShapeType() == TopAbs_FACE )
 529       faceID = myBlock.ShapeID( aS );
 530     else { // edge maybe vertical or top horizontal
 531       gp_XYZ aCoord = aParams;
 532       if ( aCoord.Z() == 1. )
 533         aCoord.SetZ( 0.5 ); // move from top down
 534       else
 535         aCoord.SetX( 0.5 ); // move along X
 536       faceID = SMESH_Block::GetShapeIDByParams( aCoord );
 537     }
 538     ASSERT( SMESH_Block::IsFaceID( faceID ));
 539     int fIndex = SMESH_Block::ShapeIndex( faceID );
 540     StdMeshers_IJNodeMap & ijNodes = myWallNodesMaps[ fIndex ];
 541     // look for a base node in ijNodes
 542     const SMDS_MeshNode* baseNode = pMesh->GetMeshDS()->FindNode( aTN.BaseNodeID() );
 543     StdMeshers_IJNodeMap::const_iterator par_nVec = ijNodes.begin();
 544     for ( ; par_nVec != ijNodes.end(); par_nVec++ )
 545       if ( par_nVec->second[ 0 ] == baseNode ) {
 546         pNode = (SMDS_MeshNode*)par_nVec->second.at( z );
 547         aTN.SetNode(pNode);
 548         return;
 549       }
 550   }
 551   //
 552   myBlock.Point(aParams, aS, aP1);
 553   //
 554   SMDS_NodeIteratorPtr ite=
 555     pMesh->GetSubMeshContaining(aS)->GetSubMeshDS()->GetNodes();
 556   while(ite->more()) {
 557     const SMDS_MeshNode* aNode = ite->next();
 558     if(myTool->IsMedium(aNode))
 559       continue;
 560     aX=aNode->X();
 561     aY=aNode->Y();
 562     aZ=aNode->Z();
 563     aP2.SetCoord(aX, aY, aZ);
 564     aD=(double)aP1.SquareDistance(aP2);
 565     //printf("** D=%lf ", aD, aTol2);
 566     if (aD < minD) {
 567       pNode=(SMDS_MeshNode*)aNode;
 568       aTN.SetNode(pNode);
 569       minD = aD;
 570       //printf(" Ok\n");
 571       if (aD<aTol2)
 572         return;
 573     }
 574   }
 575   //
 576   //printf(" KO\n");
 577   //aTN.SetNode(pNode);
 578   //MESSAGE("StdMeshers_Penta_3D::FindNodeOnShape(), can not find the node");
 579   //myErrorStatus=11; // can not find the node;
 580 }
 581
 582
 583 //=======================================================================
 584 //function : SetHorizEdgeXYZ
 585 //purpose  :
 586 //=======================================================================
 587
 588 double StdMeshers_Penta_3D::SetHorizEdgeXYZ(const gp_XYZ&                  aBaseNodeParams,
 589                                             const int                      aFaceID,
 590                                             vector<const SMDS_MeshNode*>*& aCol1,
 591                                             vector<const SMDS_MeshNode*>*& aCol2)
 592 {
 593   // find base and top edges of the face
 594   enum { BASE = 0, TOP };
 595   vector< int > edgeVec; // 0-base, 1-top
 596   SMESH_Block::GetFaceEdgesIDs( aFaceID, edgeVec );
 597   //
 598   int coord = SMESH_Block::GetCoordIndOnEdge( edgeVec[ BASE ] );
 599   bool isForward = myBlock.IsForwadEdge( edgeVec[ BASE ] );
 600
 601   double param = aBaseNodeParams.Coord( coord );
 602   if ( !isForward)
 603     param = 1. - param;
 604   //
 605   // look for columns around param
 606   StdMeshers_IJNodeMap & ijNodes =
 607     myWallNodesMaps[ SMESH_Block::ShapeIndex( aFaceID )];
 608   StdMeshers_IJNodeMap::iterator par_nVec_1 = ijNodes.begin();
 609   while ( par_nVec_1->first < param )
 610     par_nVec_1++;
 611   StdMeshers_IJNodeMap::iterator par_nVec_2 = par_nVec_1;
 612   //
 613   double r = 0;
 614   if ( par_nVec_1 != ijNodes.begin() ) {
 615     par_nVec_1--;
 616     r = ( param - par_nVec_1->first ) / ( par_nVec_2->first - par_nVec_1->first );
 617   }
 618   aCol1 = & par_nVec_1->second;
 619   aCol2 = & par_nVec_2->second;
 620
 621   // top edge
 622   if (1) {
 623     // this variant is better for cases with curved edges and
 624     // different nodes distribution on top and base edges
 625     const SMDS_MeshNode* n1 = aCol1->back();
 626     const SMDS_MeshNode* n2 = aCol2->back();
 627     gp_XYZ xyz1( n1->X(), n1->Y(), n1->Z() );
 628     gp_XYZ xyz2( n2->X(), n2->Y(), n2->Z() );
 629     myShapeXYZ[ edgeVec[ 1 ] ] = ( 1. - r ) * xyz1 + r * xyz2;
 630   }
 631   else {
 632     // this variant is better for other cases
 633 //   SMESH_MesherHelper helper( *GetMesh() );
 634 //   const TopoDS_Edge & edge = TopoDS::Edge( myBlock.Shape( edgeVec[ TOP ]));
 635 //   double u1 = helper.GetNodeU( edge, n1 );
 636 //   double u2 = helper.GetNodeU( edge, n2 );
 637 //   double u = ( 1. - r ) * u1 + r * u2;
 638 //   gp_XYZ topNodeParams;
 639 //   myBlock.Block().EdgeParameters( edgeVec[ TOP ], u, topNodeParams );
 640 //   myBlock.Block().EdgePoint( edgeVec[ TOP ],
 641 //                              topNodeParams,
 642 //                              myShapeXYZ[ edgeVec[ TOP ]]);
 643   }
 644
 645   // base edge
 646   myBlock.Block().EdgePoint( edgeVec[ BASE ],
 647                              aBaseNodeParams,
 648                              myShapeXYZ[ edgeVec[ BASE ]]);
 649   return r;
 650 }
 651
 652
 653 //=======================================================================
 654 //function : MakeVolumeMesh
 655 //purpose  :
 656 //=======================================================================
 657 void StdMeshers_Penta_3D::MakeVolumeMesh()
 658 {
 659   int i, j, ij, ik, i1, i2, aSSID;
 660   //
 661   SMESH_Mesh*   pMesh = GetMesh();
 662   SMESHDS_Mesh* meshDS = pMesh->GetMeshDS();
 663   //
 664   int shapeID = meshDS->ShapeToIndex( myShape );
 665   //
 666   // 1. Set Node In Volume
 667   ik = myISize-1;
 668   for (i=1; i<ik; ++i){
 669     for (j=0; j<myJSize; ++j){
 670       ij=i*myJSize+j;
 671       const StdMeshers_TNode& aTN = myTNodes[ij];
 672       aSSID=aTN.ShapeSupportID();
 673       if (aSSID==SMESH_Block::ID_NONE) {
 674         SMDS_MeshNode* aNode = (SMDS_MeshNode*)aTN.Node();
 675         meshDS->SetNodeInVolume(aNode, shapeID);
 676       }
 677     }
 678   }
 679   //
 680   // 2. Make pentahedrons
 681   int aID0, k , aJ[3];
 682   vector<const SMDS_MeshNode*> aN;
 683   //
 684   SMDS_ElemIteratorPtr itf, aItNodes;
 685   //
 686   const TopoDS_Face& aFxy0=
 687     TopoDS::Face(myBlock.Shape(SMESH_Block::ID_Fxy0));
 688   SMESH_subMesh *aSubMesh0 = pMesh->GetSubMeshContaining(aFxy0);
 689   SMESHDS_SubMesh *aSM0 = aSubMesh0->GetSubMeshDS();
 690   //
 691   itf = aSM0->GetElements();
 692   while(itf->more()) {
 693     const SMDS_MeshElement* pE0 = itf->next();
 694     //
 695     int nbFaceNodes = pE0->NbNodes();
 696     if(myCreateQuadratic)
 697       nbFaceNodes = nbFaceNodes/2;
 698     if ( aN.size() < nbFaceNodes * 2 )
 699       aN.resize( nbFaceNodes * 2 );
 700     //
 701     k=0;
 702     aItNodes=pE0->nodesIterator();
 703     while (aItNodes->more()) {
 704       //const SMDS_MeshElement* pNode = aItNodes->next();
 705       const SMDS_MeshNode* pNode =
 706         static_cast<const SMDS_MeshNode*> (aItNodes->next());
 707       if(myTool->IsMedium(pNode))
 708         continue;
 709       aID0 = pNode->GetID();
 710       aJ[k] = GetIndexOnLayer(aID0);
 711       if (!myErrorStatus->IsOK()) {
 712         MESSAGE("StdMeshers_Penta_3D::MakeVolumeMesh");
 713         return;
 714       }
 715       //
 716       ++k;
 717     }
 718     //
 719     bool forward = true;
 720     for (i=0; i<ik; ++i) {
 721       i1=i;
 722       i2=i+1;
 723       for(j=0; j<nbFaceNodes; ++j) {
 724         ij = i1*myJSize+aJ[j];
 725         const StdMeshers_TNode& aTN1 = myTNodes[ij];
 726         const SMDS_MeshNode* aN1 = aTN1.Node();
 727         aN[j]=aN1;
 728         //
 729         ij=i2*myJSize+aJ[j];
 730         const StdMeshers_TNode& aTN2 = myTNodes[ij];
 731         const SMDS_MeshNode* aN2 = aTN2.Node();
 732         aN[j+nbFaceNodes] = aN2;
 733       }
 734       // check if volume orientation will be ok
 735       if ( i == 0 ) {
 736         SMDS_VolumeTool vTool;
 737         switch ( nbFaceNodes ) {
 738         case 3: {
 739           SMDS_VolumeOfNodes tmpVol (aN[0], aN[1], aN[2],
 740                                      aN[3], aN[4], aN[5]);
 741           vTool.Set( &tmpVol );
 742           break;
 743         }
 744         case 4: {
 745           SMDS_VolumeOfNodes tmpVol(aN[0], aN[1], aN[2], aN[3],
 746                                     aN[4], aN[5], aN[6], aN[7]);
 747           vTool.Set( &tmpVol );
 748           break;
 749         }
 750         default:
 751           continue;
 752         }
 753         forward = vTool.IsForward();
 754       }
 755       // add volume
 756       SMDS_MeshVolume* aV = 0;
 757       switch ( nbFaceNodes ) {
 758       case 3:
 759         if ( forward ) {
 760           //aV = meshDS->AddVolume(aN[0], aN[1], aN[2],
 761           //                       aN[3], aN[4], aN[5]);
 762           aV = myTool->AddVolume(aN[0], aN[1], aN[2], aN[3], aN[4], aN[5]);
 763         }
 764         else {
 765           //aV = meshDS->AddVolume(aN[0], aN[2], aN[1],
 766           //                       aN[3], aN[5], aN[4]);
 767           aV = myTool->AddVolume(aN[0], aN[2], aN[1], aN[3], aN[5], aN[4]);
 768         }
 769         break;
 770       case 4:
 771         if ( forward ) {
 772           //aV = meshDS->AddVolume(aN[0], aN[1], aN[2], aN[3],
 773           //                       aN[4], aN[5], aN[6], aN[7]);
 774           aV = myTool->AddVolume(aN[0], aN[1], aN[2], aN[3],
 775                                  aN[4], aN[5], aN[6], aN[7]);
 776         }
 777         else {
 778           //aV = meshDS->AddVolume(aN[0], aN[3], aN[2], aN[1],
 779           //                       aN[4], aN[7], aN[6], aN[5]);
 780           aV = myTool->AddVolume(aN[0], aN[3], aN[2], aN[1],
 781                                  aN[4], aN[7], aN[6], aN[5]);
 782         }
 783         break;
 784       default:
 785         continue;
 786       }
 787       meshDS->SetMeshElementOnShape(aV, shapeID);
 788     }
 789   }
 790 }
 791
 792 //=======================================================================
 793 //function : MakeMeshOnFxy1
 794 //purpose  :
 795 //=======================================================================
 796 void StdMeshers_Penta_3D::MakeMeshOnFxy1()
 797 {
 798   int aID0, aJ, aLevel, ij, aNbNodes, k;
 799   //
 800   SMDS_NodeIteratorPtr itn;
 801   SMDS_ElemIteratorPtr itf, aItNodes;
 802   SMDSAbs_ElementType aElementType;
 803   //
 804   const TopoDS_Face& aFxy0=
 805     TopoDS::Face(myBlock.Shape(SMESH_Block::ID_Fxy0));
 806   const TopoDS_Face& aFxy1=
 807     TopoDS::Face(myBlock.Shape(SMESH_Block::ID_Fxy1));
 808   SMESH_MesherHelper faceHelper( *GetMesh() );
 809   faceHelper.IsQuadraticSubMesh(aFxy1);
 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 = faceHelper.AddFace(aNodes1[0], aNodes1[1], aNodes1[2]);
 865       break;
 866     case 4:
 867       face = faceHelper.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 || (myCreateQuadratic && 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