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