src/StdMeshers/StdMeshers_Penta_3D.cxx

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