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