Back to index

salome-smesh  6.5.0
StdMeshers_Penta_3D.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 
00023 //  SMESH StdMeshers_Penta_3D implementaion of SMESH idl descriptions
00024 //  File   : StdMeshers_Penta_3D.cxx
00025 //  Module : SMESH
00026 //
00027 #include "StdMeshers_Penta_3D.hxx"
00028 
00029 #include "utilities.h"
00030 #include "Utils_ExceptHandlers.hxx"
00031 
00032 #include "SMDS_EdgePosition.hxx"
00033 #include "SMDS_MeshElement.hxx"
00034 #include "SMDS_VolumeOfNodes.hxx"
00035 #include "SMDS_VolumeTool.hxx"
00036 #include "SMESHDS_SubMesh.hxx"
00037 #include "SMESH_Mesh.hxx"
00038 #include "SMESH_MesherHelper.hxx"
00039 #include "SMESH_subMesh.hxx"
00040 #include "SMESH_subMeshEventListener.hxx"
00041 #include "SMESH_Comment.hxx"
00042 
00043 #include <BRep_Tool.hxx>
00044 #include <TopExp.hxx>
00045 #include <TopExp_Explorer.hxx>
00046 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
00047 #include <TopTools_IndexedMapOfShape.hxx>
00048 #include <TopTools_ListIteratorOfListOfShape.hxx>
00049 #include <TopTools_ListOfShape.hxx>
00050 #include <TopTools_SequenceOfShape.hxx>
00051 #include <TopTools_MapOfShape.hxx>
00052 #include <TopoDS.hxx>
00053 #include <TopoDS_Edge.hxx>
00054 #include <TopoDS_Shell.hxx>
00055 #include <TopoDS_Vertex.hxx>
00056 #include <gp_Pnt.hxx>
00057 
00058 #include <stdio.h>
00059 #include <algorithm>
00060 
00061 using namespace std;
00062 
00063 typedef map < int, int, less<int> >::iterator   \
00064   StdMeshers_IteratorOfDataMapOfIntegerInteger;
00065 
00066 enum { NB_WALL_FACES = 4 };
00067 
00068 //=======================================================================
00069 //function : StdMeshers_Penta_3D
00070 //purpose  : 
00071 //=======================================================================
00072 StdMeshers_Penta_3D::StdMeshers_Penta_3D()
00073   : myErrorStatus(SMESH_ComputeError::New())
00074 {
00075   myTol3D=0.1;
00076   myWallNodesMaps.resize( SMESH_Block::NbFaces() );
00077   myShapeXYZ.resize( SMESH_Block::NbSubShapes() );
00078   myTool = 0;
00079 }
00080 
00081 //=======================================================================
00082 //function : ~StdMeshers_Penta_3D
00083 //purpose  : 
00084 //=======================================================================
00085 
00086 StdMeshers_Penta_3D::~StdMeshers_Penta_3D()
00087 {
00088 }
00089 
00090 //=======================================================================
00091 //function : Compute
00092 //purpose  : 
00093 //=======================================================================
00094 bool StdMeshers_Penta_3D::Compute(SMESH_Mesh& aMesh, 
00095                                   const TopoDS_Shape& aShape)
00096 {
00097   MESSAGE("StdMeshers_Penta_3D::Compute()");
00098   //
00099   bool bOK=false;
00100   //
00101   myShape=aShape;
00102   SetMesh(aMesh);
00103   //
00104   CheckData();
00105   if (!myErrorStatus->IsOK()) {
00106     return bOK;
00107   }
00108 
00109   //
00110   MakeBlock();
00111   if (!myErrorStatus->IsOK()) {
00112     return bOK;
00113   }
00114   //
00115   ClearMeshOnFxy1();
00116   if (!myErrorStatus->IsOK()) {
00117     return bOK;
00118   }
00119 
00120   // now unnecessary faces removed, we can load medium nodes
00121   SMESH_MesherHelper helper(aMesh);
00122   myTool = &helper;
00123   myCreateQuadratic = myTool->IsQuadraticSubMesh(aShape);
00124 
00125   //
00126   MakeNodes();
00127   if (!myErrorStatus->IsOK()) {
00128     return bOK;
00129   }
00130   //
00131   MakeConnectingMap();
00132   //
00133   MakeMeshOnFxy1();
00134   if (!myErrorStatus->IsOK()) {
00135     return bOK;
00136   }
00137   //
00138   MakeVolumeMesh();
00139   //
00140   return !bOK;
00141 }
00142 
00143 //=======================================================================
00144 //function : MakeNodes
00145 //purpose  : 
00146 //=======================================================================
00147 void StdMeshers_Penta_3D::MakeNodes()
00148 {
00149   const int aNbSIDs=9;
00150   int i, j, k, ij, iNbN, aNodeID, aSize, iErr;
00151   double aX, aY, aZ;
00152   SMESH_Block::TShapeID aSID, aSIDs[aNbSIDs]={
00153     SMESH_Block::ID_V000, SMESH_Block::ID_V100, 
00154     SMESH_Block::ID_V110, SMESH_Block::ID_V010,
00155     SMESH_Block::ID_Ex00, SMESH_Block::ID_E1y0, 
00156     SMESH_Block::ID_Ex10, SMESH_Block::ID_E0y0,
00157     SMESH_Block::ID_Fxy0
00158   }; 
00159   //
00160   SMESH_Mesh* pMesh=GetMesh();
00161   //
00162   // 1. Define the sizes of mesh
00163   //
00164   // 1.1 Horizontal size
00165   myJSize=0;
00166   for (i=0; i<aNbSIDs; ++i) {
00167     const TopoDS_Shape& aS = myBlock.Shape(aSIDs[i]);
00168     SMESH_subMesh *aSubMesh = pMesh->GetSubMeshContaining(aS);
00169     ASSERT(aSubMesh);
00170     SMESHDS_SubMesh *aSM = aSubMesh->GetSubMeshDS();
00171     if(!myCreateQuadratic) {
00172       iNbN = aSM->NbNodes();
00173     }
00174     else {
00175       iNbN = 0;
00176       SMDS_NodeIteratorPtr itn = aSM->GetNodes();
00177       while(itn->more()) {
00178         const SMDS_MeshNode* aNode = itn->next();
00179         if(myTool->IsMedium(aNode))
00180           continue;
00181         iNbN++;
00182       }
00183     }
00184     myJSize += iNbN;
00185   }
00186   //printf("***  Horizontal: number of nodes summary=%d\n", myJSize);
00187   //
00188   // 1.2 Vertical size
00189   myISize=2;
00190   {
00191     const TopoDS_Shape& aS=myBlock.Shape(SMESH_Block::ID_E00z);
00192     SMESH_subMesh *aSubMesh = pMesh->GetSubMeshContaining(aS);
00193     ASSERT(aSubMesh);
00194     SMESHDS_SubMesh *aSM = aSubMesh->GetSubMeshDS();
00195     if(!myCreateQuadratic) {
00196       iNbN = aSM->NbNodes();
00197     }
00198     else {
00199       iNbN = 0;
00200       SMDS_NodeIteratorPtr itn = aSM->GetNodes();
00201       while(itn->more()) {
00202         const SMDS_MeshNode* aNode = itn->next();
00203         if(myTool->IsMedium(aNode))
00204           continue;
00205         iNbN++;
00206       }
00207     }
00208     myISize += iNbN;
00209   }
00210   //printf("***  Vertical: number of nodes on edges and vertices=%d\n", myISize);
00211   //
00212   aSize=myISize*myJSize;
00213   myTNodes.resize(aSize);
00214   //
00215   StdMeshers_TNode aTNode;
00216   gp_XYZ aCoords;
00217   gp_Pnt aP3D;
00218   //
00219   // 2. Fill the repers on base face (Z=0)
00220   i=0; j=0;
00221   // vertices
00222   for (k=0; k<aNbSIDs; ++k) {
00223     aSID=aSIDs[k];
00224     const TopoDS_Shape& aS = myBlock.Shape(aSID);
00225     SMDS_NodeIteratorPtr ite = pMesh->GetSubMeshContaining(aS)->GetSubMeshDS()->GetNodes();
00226     while(ite->more()) {
00227       const SMDS_MeshNode* aNode = ite->next();
00228       if(myTool->IsMedium(aNode))
00229         continue;
00230       aNodeID=aNode->GetID();
00231       //
00232       aTNode.SetNode(aNode);
00233       aTNode.SetShapeSupportID(aSID);
00234       aTNode.SetBaseNodeID(aNodeID);
00235       //
00236       if ( SMESH_Block::IsEdgeID (aSID)) {
00237         const SMDS_EdgePosition* epos =
00238           static_cast<const SMDS_EdgePosition*>(aNode->GetPosition());
00239         myBlock.ComputeParameters( epos->GetUParameter(), aS, aCoords );
00240       }
00241       else {
00242         aX=aNode->X();
00243         aY=aNode->Y();
00244         aZ=aNode->Z();
00245         aP3D.SetCoord(aX, aY, aZ);
00246         myBlock.ComputeParameters(aP3D, aS, aCoords);
00247       }
00248       iErr = myBlock.ErrorStatus();
00249       if (iErr) {
00250         MESSAGE("StdMeshers_Penta_3D::MakeNodes()," <<
00251                 "SMESHBlock: ComputeParameters operation failed");
00252         myErrorStatus=myBlock.GetError();
00253         return;
00254       }
00255       aTNode.SetNormCoord(aCoords);
00256       ij=i*myJSize+j;
00257       myTNodes[ij]=aTNode;
00258       ++j;
00259     }
00260   }
00261 
00262   // 3.1 Fill maps of wall nodes
00263   SMESH_Block::TShapeID wallFaceID[ NB_WALL_FACES ] = {
00264     SMESH_Block::ID_Fx0z, SMESH_Block::ID_Fx1z,
00265     SMESH_Block::ID_F0yz, SMESH_Block::ID_F1yz
00266   };
00267   SMESH_Block::TShapeID baseEdgeID[ NB_WALL_FACES ] = {
00268     SMESH_Block::ID_Ex00, SMESH_Block::ID_Ex10,
00269     SMESH_Block::ID_E0y0, SMESH_Block::ID_E1y0
00270   };
00271   for ( i = 0; i < NB_WALL_FACES ; ++i ) {
00272     int fIndex = SMESH_Block::ShapeIndex( wallFaceID[ i ]);
00273     bool ok = LoadIJNodes (myWallNodesMaps[ fIndex ],
00274                            TopoDS::Face( myBlock.Shape( wallFaceID[ i ] )),
00275                            TopoDS::Edge( myBlock.Shape( baseEdgeID[ i ] )),
00276                            pMesh->GetMeshDS());
00277     if ( !ok ) {
00278       myErrorStatus->myName = COMPERR_BAD_INPUT_MESH;
00279       myErrorStatus->myComment = SMESH_Comment() <<
00280         "Can't find regular quadrangle mesh on a side face #" <<
00281         pMesh->GetMeshDS()->ShapeToIndex( myBlock.Shape( wallFaceID[ i ]));
00282       return;
00283     }
00284   }
00285 
00286   // 3.2 find node columns for vertical edges and edge IDs
00287   vector<const SMDS_MeshNode*> * verticEdgeNodes[ NB_WALL_FACES ];
00288   SMESH_Block::TShapeID          verticEdgeID   [ NB_WALL_FACES ];
00289   for ( i = 0; i < NB_WALL_FACES ; ++i ) { // 4 first base nodes are nodes on vertices
00290     // edge ID
00291     SMESH_Block::TShapeID eID, vID = aSIDs[ i ];
00292     ShapeSupportID(false, vID, eID);
00293     verticEdgeID[ i ] = eID;
00294     // column nodes
00295     StdMeshers_TNode& aTNode = myTNodes[ i ];
00296     verticEdgeNodes[ i ] = 0;
00297     for ( j = 0; j < NB_WALL_FACES ; ++j ) { // loop on 4 wall faces
00298       int fIndex = SMESH_Block::ShapeIndex( wallFaceID[ j ]);
00299       StdMeshers_IJNodeMap & ijNodes= myWallNodesMaps[ fIndex ];
00300       if ( ijNodes.begin()->second[0] == aTNode.Node() )
00301         verticEdgeNodes[ i ] = & ijNodes.begin()->second;
00302       else if ( ijNodes.rbegin()->second[0] == aTNode.Node() )
00303         verticEdgeNodes[ i ] = & ijNodes.rbegin()->second;
00304       if ( verticEdgeNodes[ i ] )
00305         break;
00306     }
00307   }
00308 
00309   // 3.3 set XYZ of vertices, and initialize of the rest
00310   SMESHDS_Mesh* aMesh = GetMesh()->GetMeshDS();
00311   for ( int id = SMESH_Block::ID_V000; id < SMESH_Block::ID_Shell; ++id ) {
00312     if ( SMESH_Block::IsVertexID( id )) {
00313       TopoDS_Shape V = myBlock.Shape( id );
00314       SMESHDS_SubMesh* sm = aMesh->MeshElements( V );
00315       const SMDS_MeshNode* n = sm->GetNodes()->next();
00316       myShapeXYZ[ id ].SetCoord( n->X(), n->Y(), n->Z() );
00317     }
00318     else
00319       myShapeXYZ[ id ].SetCoord( 0., 0., 0. );
00320   }
00321 
00322 
00323   // 4. Fill the rest repers
00324   bool bIsUpperLayer;
00325   int iBNID;
00326   SMESH_Block::TShapeID aSSID, aBNSSID;
00327   StdMeshers_TNode aTN;
00328   //
00329 
00330   // create top face and find UV for it's corners
00331   const TopoDS_Face& TopFace = TopoDS::Face(myBlock.Shape(SMESH_Block::ID_Fxy1));
00332   SMESHDS_Mesh* meshDS = pMesh->GetMeshDS();
00333   int topfaceID = meshDS->ShapeToIndex(TopFace);
00334   const TopoDS_Vertex& v001 = TopoDS::Vertex(myBlock.Shape(SMESH_Block::ID_V001));
00335   SMDS_NodeIteratorPtr itn = pMesh->GetSubMeshContaining(v001)->GetSubMeshDS()->GetNodes();
00336   const SMDS_MeshNode* N = itn->next();
00337   gp_XY UV001 = myTool->GetNodeUV(TopFace,N);
00338   const TopoDS_Vertex& v101 = TopoDS::Vertex(myBlock.Shape(SMESH_Block::ID_V101));
00339   itn = pMesh->GetSubMeshContaining(v101)->GetSubMeshDS()->GetNodes();
00340   N = itn->next();
00341   gp_XY UV101 = myTool->GetNodeUV(TopFace,N);
00342   const TopoDS_Vertex& v011 = TopoDS::Vertex(myBlock.Shape(SMESH_Block::ID_V011));
00343   itn = pMesh->GetSubMeshContaining(v011)->GetSubMeshDS()->GetNodes();
00344   N = itn->next();
00345   gp_XY UV011 = myTool->GetNodeUV(TopFace,N);
00346   const TopoDS_Vertex& v111 = TopoDS::Vertex(myBlock.Shape(SMESH_Block::ID_V111));
00347   itn = pMesh->GetSubMeshContaining(v111)->GetSubMeshDS()->GetNodes();
00348   N = itn->next();
00349   gp_XY UV111 = myTool->GetNodeUV(TopFace,N);
00350 
00351   for (j=0; j<myJSize; ++j) { // loop on all nodes of the base face (ID_Fxy0)
00352     // base node info
00353     const StdMeshers_TNode& aBN = myTNodes[j];
00354     aBNSSID = (SMESH_Block::TShapeID)aBN.ShapeSupportID();
00355     iBNID = aBN.BaseNodeID();
00356     const gp_XYZ& aBNXYZ = aBN.NormCoord();
00357     bool createNode = ( aBNSSID == SMESH_Block::ID_Fxy0 ); // if base node is inside a bottom face
00358     //
00359     // set XYZ on horizontal edges and get node columns of faces:
00360     // 2 columns for each face, between which a base node is located
00361     vector<const SMDS_MeshNode*>* nColumns[8];
00362     double ratio[ NB_WALL_FACES ]; // base node position between columns [0.-1.]
00363     if ( createNode ) {
00364       for ( k = 0; k < NB_WALL_FACES ; ++k ) {
00365         ratio[ k ] = SetHorizEdgeXYZ (aBNXYZ, wallFaceID[ k ],
00366                                       nColumns[k*2], nColumns[k*2+1]);
00367       }
00368     }
00369     //
00370     // XYZ on the bottom and top faces
00371     const SMDS_MeshNode* n = aBN.Node();
00372     myShapeXYZ[ SMESH_Block::ID_Fxy0 ].SetCoord( n->X(), n->Y(), n->Z() );
00373     myShapeXYZ[ SMESH_Block::ID_Fxy1 ].SetCoord( 0., 0., 0. );
00374     //
00375     // first create or find a top node, then the rest ones in a column
00376     for (i=myISize-1; i>0; --i) // vertical loop, from top to bottom
00377     {
00378       bIsUpperLayer = (i==(myISize-1));
00379       gp_XY UV_Ex01, UV_Ex11, UV_E0y1, UV_E1y1;
00380       if ( createNode ) // a base node is inside a top face
00381       {
00382         // set XYZ on vertical edges and faces
00383         for ( k = 0; k < NB_WALL_FACES ; ++k ) {
00384           // XYZ on a vertical edge 
00385           const SMDS_MeshNode* n = (*verticEdgeNodes[ k ]) [ i ];
00386           myShapeXYZ[ verticEdgeID[ k ] ].SetCoord( n->X(), n->Y(), n->Z() );
00387           // XYZ on a face (part 1 from one column)
00388           n = (*nColumns[k*2]) [ i ];
00389           gp_XYZ xyz( n->X(), n->Y(), n->Z() );
00390           myShapeXYZ[ wallFaceID[ k ]] = ( 1. - ratio[ k ]) * xyz;
00391           gp_XY tmp1;
00392           if( bIsUpperLayer ) {
00393             tmp1 = myTool->GetNodeUV(TopFace,n);
00394             tmp1 = ( 1. - ratio[ k ]) * tmp1;
00395           }
00396           // XYZ on a face (part 2 from other column)
00397           n = (*nColumns[k*2+1]) [ i ];
00398           xyz.SetCoord( n->X(), n->Y(), n->Z() );
00399           myShapeXYZ[ wallFaceID[ k ]] += ratio[ k ] * xyz;
00400           if( bIsUpperLayer ) {
00401             gp_XY tmp2 = myTool->GetNodeUV(TopFace,n);
00402             tmp1 +=  ratio[ k ] * tmp2;
00403             if( k==0 )
00404               UV_Ex01 = tmp1;
00405             else if( k==1 )
00406               UV_Ex11 = tmp1;
00407             else if( k==2 )
00408               UV_E0y1 = tmp1;
00409             else
00410               UV_E1y1 = tmp1;
00411           }
00412         }
00413       }
00414       // fill current node info
00415       //   -index in aTNodes
00416       ij=i*myJSize+j; 
00417       //   -normalized coordinates  
00418       aX=aBNXYZ.X();  
00419       aY=aBNXYZ.Y();
00420       //aZ=aZL[i];
00421       aZ=(double)i/(double)(myISize-1);
00422       aCoords.SetCoord(aX, aY, aZ);
00423       //
00424       //   suporting shape ID
00425       ShapeSupportID(bIsUpperLayer, aBNSSID, aSSID);
00426       if (!myErrorStatus->IsOK()) {
00427         MESSAGE("StdMeshers_Penta_3D::MakeNodes() ");
00428         return;
00429       }
00430       //
00431       aTN.SetShapeSupportID(aSSID);
00432       aTN.SetNormCoord(aCoords);
00433       aTN.SetBaseNodeID(iBNID);
00434       //
00435       if (aSSID!=SMESH_Block::ID_NONE){
00436         // try to find the node
00437         const TopoDS_Shape& aS=myBlock.Shape((int)aSSID);
00438         FindNodeOnShape(aS, aCoords, i, aTN);
00439       }
00440       else{
00441         // create node and get its id
00442         CreateNode (bIsUpperLayer, aCoords, aTN);
00443         //
00444         if ( bIsUpperLayer ) {
00445           const SMDS_MeshNode* n = aTN.Node();
00446           myShapeXYZ[ SMESH_Block::ID_Fxy1 ].SetCoord( n->X(), n->Y(), n->Z() );
00447           // set node on top face:
00448           // find UV parameter for this node
00449           //              UV_Ex11
00450           //   UV011+-----+----------+UV111
00451           //        |                |
00452           //        |                |
00453           // UV_E0y1+     +node      +UV_E1y1
00454           //        |                |
00455           //        |                |
00456           //        |                |
00457           //   UV001+-----+----------+UV101
00458           //              UV_Ex01
00459           gp_Pnt2d aP;
00460           double u = aCoords.X(), v = aCoords.Y();
00461           double u1 = ( 1. - u ), v1 = ( 1. - v );
00462           aP.ChangeCoord()  = UV_Ex01 * v1;
00463           aP.ChangeCoord() += UV_Ex11 * v;
00464           aP.ChangeCoord() += UV_E0y1 * u1;
00465           aP.ChangeCoord() += UV_E1y1 * u;
00466           aP.ChangeCoord() -= UV001 * u1 * v1;
00467           aP.ChangeCoord() -= UV101 * u  * v1;
00468           aP.ChangeCoord() -= UV011 * u1 * v;
00469           aP.ChangeCoord() -= UV111 * u  * v;
00470           meshDS->SetNodeOnFace((SMDS_MeshNode*)n, topfaceID, aP.X(), aP.Y());
00471         }
00472       }
00473       if (!myErrorStatus->IsOK()) {
00474         MESSAGE("StdMeshers_Penta_3D::MakeNodes() ");
00475         return;
00476       }
00477       //
00478       myTNodes[ij]=aTN;
00479     }
00480   }
00481 }
00482 
00483 
00484 //=======================================================================
00485 //function : FindNodeOnShape
00486 //purpose  : 
00487 //=======================================================================
00488 
00489 void StdMeshers_Penta_3D::FindNodeOnShape(const TopoDS_Shape& aS,
00490                                           const gp_XYZ&       aParams,
00491                                           const int           z,
00492                                           StdMeshers_TNode&   aTN)
00493 {
00494   double aX, aY, aZ, aD, aTol2, minD;
00495   gp_Pnt aP1, aP2;
00496   //
00497   SMESH_Mesh* pMesh = GetMesh();
00498   aTol2 = myTol3D*myTol3D;
00499   minD = 1.e100;
00500   SMDS_MeshNode* pNode = NULL;
00501   //
00502   if ( aS.ShapeType() == TopAbs_FACE ||
00503        aS.ShapeType() == TopAbs_EDGE ) {
00504     // find a face ID to which aTN belongs to
00505     int faceID;
00506     if ( aS.ShapeType() == TopAbs_FACE )
00507       faceID = myBlock.ShapeID( aS );
00508     else { // edge maybe vertical or top horizontal
00509       gp_XYZ aCoord = aParams;
00510       if ( aCoord.Z() == 1. )
00511         aCoord.SetZ( 0.5 ); // move from top down
00512       else
00513         aCoord.SetX( 0.5 ); // move along X
00514       faceID = SMESH_Block::GetShapeIDByParams( aCoord );
00515     }
00516     ASSERT( SMESH_Block::IsFaceID( faceID ));
00517     int fIndex = SMESH_Block::ShapeIndex( faceID );
00518     StdMeshers_IJNodeMap & ijNodes = myWallNodesMaps[ fIndex ];
00519     // look for a base node in ijNodes
00520     const SMDS_MeshNode* baseNode = pMesh->GetMeshDS()->FindNode( aTN.BaseNodeID() );
00521     StdMeshers_IJNodeMap::const_iterator par_nVec = ijNodes.begin();
00522     for ( ; par_nVec != ijNodes.end(); par_nVec++ )
00523       if ( par_nVec->second[ 0 ] == baseNode ) {
00524         pNode = (SMDS_MeshNode*)par_nVec->second.at( z );
00525         aTN.SetNode(pNode);
00526         return;
00527       }
00528   }
00529   //
00530   myBlock.Point(aParams, aS, aP1);
00531   //
00532   SMDS_NodeIteratorPtr ite=
00533     pMesh->GetSubMeshContaining(aS)->GetSubMeshDS()->GetNodes();
00534   while(ite->more()) {
00535     const SMDS_MeshNode* aNode = ite->next();
00536     if(myTool->IsMedium(aNode))
00537       continue;
00538     aX=aNode->X();
00539     aY=aNode->Y();
00540     aZ=aNode->Z();
00541     aP2.SetCoord(aX, aY, aZ);
00542     aD=(double)aP1.SquareDistance(aP2);
00543     //printf("** D=%lf ", aD, aTol2);
00544     if (aD < minD) {
00545       pNode=(SMDS_MeshNode*)aNode;
00546       aTN.SetNode(pNode);
00547       minD = aD;
00548       //printf(" Ok\n");
00549       if (aD<aTol2)
00550         return; 
00551     }
00552   }
00553   //
00554   //printf(" KO\n");
00555   //aTN.SetNode(pNode);
00556   //MESSAGE("StdMeshers_Penta_3D::FindNodeOnShape(), can not find the node");
00557   //myErrorStatus=11; // can not find the node;
00558 }
00559 
00560 
00561 //=======================================================================
00562 //function : SetHorizEdgeXYZ
00563 //purpose  : 
00564 //=======================================================================
00565 
00566 double StdMeshers_Penta_3D::SetHorizEdgeXYZ(const gp_XYZ&                  aBaseNodeParams,
00567                                             const int                      aFaceID,
00568                                             vector<const SMDS_MeshNode*>*& aCol1,
00569                                             vector<const SMDS_MeshNode*>*& aCol2)
00570 {
00571   // find base and top edges of the face
00572   enum { BASE = 0, TOP };
00573   vector< int > edgeVec; // 0-base, 1-top
00574   SMESH_Block::GetFaceEdgesIDs( aFaceID, edgeVec );
00575   //
00576   int coord = SMESH_Block::GetCoordIndOnEdge( edgeVec[ BASE ] );
00577   bool isForward = myBlock.IsForwadEdge( edgeVec[ BASE ] );
00578 
00579   double param = aBaseNodeParams.Coord( coord );
00580   if ( !isForward)
00581     param = 1. - param;
00582   //
00583   // look for columns around param
00584   StdMeshers_IJNodeMap & ijNodes =
00585     myWallNodesMaps[ SMESH_Block::ShapeIndex( aFaceID )];
00586   StdMeshers_IJNodeMap::iterator par_nVec_1 = ijNodes.begin();
00587   while ( par_nVec_1->first < param )
00588     par_nVec_1++;
00589   StdMeshers_IJNodeMap::iterator par_nVec_2 = par_nVec_1;
00590   //
00591   double r = 0;
00592   if ( par_nVec_1 != ijNodes.begin() ) {
00593     par_nVec_1--;
00594     r = ( param - par_nVec_1->first ) / ( par_nVec_2->first - par_nVec_1->first );
00595   }
00596   aCol1 = & par_nVec_1->second;
00597   aCol2 = & par_nVec_2->second;
00598 
00599   // top edge
00600   if (1) {
00601     // this variant is better for cases with curved edges and
00602     // different nodes distribution on top and base edges
00603     const SMDS_MeshNode* n1 = aCol1->back();
00604     const SMDS_MeshNode* n2 = aCol2->back();
00605     gp_XYZ xyz1( n1->X(), n1->Y(), n1->Z() );
00606     gp_XYZ xyz2( n2->X(), n2->Y(), n2->Z() );
00607     myShapeXYZ[ edgeVec[ 1 ] ] = ( 1. - r ) * xyz1 + r * xyz2;
00608   }
00609   else {
00610     // this variant is better for other cases
00611     //   SMESH_MesherHelper helper( *GetMesh() );
00612     //   const TopoDS_Edge & edge = TopoDS::Edge( myBlock.Shape( edgeVec[ TOP ]));
00613     //   double u1 = helper.GetNodeU( edge, n1 );
00614     //   double u2 = helper.GetNodeU( edge, n2 );
00615     //   double u = ( 1. - r ) * u1 + r * u2;
00616     //   gp_XYZ topNodeParams;
00617     //   myBlock.Block().EdgeParameters( edgeVec[ TOP ], u, topNodeParams );
00618     //   myBlock.Block().EdgePoint( edgeVec[ TOP ],
00619     //                              topNodeParams,
00620     //                              myShapeXYZ[ edgeVec[ TOP ]]);
00621   }
00622 
00623   // base edge
00624   myBlock.Block().EdgePoint( edgeVec[ BASE ],
00625                              aBaseNodeParams,
00626                              myShapeXYZ[ edgeVec[ BASE ]]);
00627   return r;
00628 }
00629 
00630 
00631 //=======================================================================
00632 //function : MakeVolumeMesh
00633 //purpose  : 
00634 //=======================================================================
00635 void StdMeshers_Penta_3D::MakeVolumeMesh()
00636 {
00637   int i, j, ij, ik, i1, i2, aSSID; 
00638   //
00639   SMESH_Mesh*   pMesh = GetMesh();
00640   SMESHDS_Mesh* meshDS = pMesh->GetMeshDS();
00641   //
00642   int shapeID = meshDS->ShapeToIndex( myShape );
00643   //
00644   // 1. Set Node In Volume
00645   ik = myISize-1;
00646   for (i=1; i<ik; ++i){
00647     for (j=0; j<myJSize; ++j){
00648       ij=i*myJSize+j;
00649       const StdMeshers_TNode& aTN = myTNodes[ij];
00650       aSSID=aTN.ShapeSupportID();
00651       if (aSSID==SMESH_Block::ID_NONE) {
00652         SMDS_MeshNode* aNode = (SMDS_MeshNode*)aTN.Node();
00653         meshDS->SetNodeInVolume(aNode, shapeID);
00654       }
00655     }
00656   }
00657   //
00658   // 2. Make pentahedrons
00659   int aID0, k , aJ[3];
00660   vector<const SMDS_MeshNode*> aN;
00661   //
00662   SMDS_ElemIteratorPtr itf, aItNodes;
00663   //
00664   const TopoDS_Face& aFxy0=
00665     TopoDS::Face(myBlock.Shape(SMESH_Block::ID_Fxy0));
00666   SMESH_subMesh *aSubMesh0 = pMesh->GetSubMeshContaining(aFxy0);
00667   SMESHDS_SubMesh *aSM0 = aSubMesh0->GetSubMeshDS();
00668   //
00669   itf = aSM0->GetElements();
00670   while(itf->more()) {
00671     const SMDS_MeshElement* pE0 = itf->next();
00672     //
00673     int nbFaceNodes = pE0->NbNodes();
00674     if(myCreateQuadratic)
00675       nbFaceNodes = nbFaceNodes/2;
00676     if ( aN.size() < nbFaceNodes * 2 )
00677       aN.resize( nbFaceNodes * 2 );
00678     //
00679     for ( k=0; k<nbFaceNodes; ++k ) {
00680       const SMDS_MeshNode* pNode = pE0->GetNode(k);
00681 //       if(myTool->IsMedium(pNode))
00682 //         continue;
00683       aID0 = pNode->GetID();
00684       aJ[k] = GetIndexOnLayer(aID0);
00685       if (!myErrorStatus->IsOK()) {
00686         MESSAGE("StdMeshers_Penta_3D::MakeVolumeMesh");
00687         return;
00688       }
00689     }
00690     //
00691     bool forward = true;
00692     for (i=0; i<ik; ++i) {
00693       i1=i;
00694       i2=i+1;
00695       for(j=0; j<nbFaceNodes; ++j) {
00696         ij = i1*myJSize+aJ[j];
00697         const StdMeshers_TNode& aTN1 = myTNodes[ij];
00698         const SMDS_MeshNode* aN1 = aTN1.Node();
00699         aN[j]=aN1;
00700         //
00701         ij=i2*myJSize+aJ[j];
00702         const StdMeshers_TNode& aTN2 = myTNodes[ij];
00703         const SMDS_MeshNode* aN2 = aTN2.Node();
00704         aN[j+nbFaceNodes] = aN2;
00705       }
00706       // check if volume orientation will be ok
00707       if ( i == 0 ) {
00708         SMDS_VolumeTool vTool;
00709         switch ( nbFaceNodes ) {
00710         case 3: {
00711           SMDS_VolumeOfNodes tmpVol (aN[0], aN[1], aN[2],
00712                                      aN[3], aN[4], aN[5]);
00713           vTool.Set( &tmpVol );
00714           break;
00715         }
00716         case 4: {
00717           SMDS_VolumeOfNodes tmpVol(aN[0], aN[1], aN[2], aN[3],
00718                                     aN[4], aN[5], aN[6], aN[7]);
00719           vTool.Set( &tmpVol );
00720           break;
00721         }
00722         default:
00723           continue;
00724         }
00725         forward = vTool.IsForward();
00726       }
00727       // add volume
00728       SMDS_MeshVolume* aV = 0;
00729       switch ( nbFaceNodes ) {
00730       case 3:
00731         if ( forward ) {
00732           //aV = meshDS->AddVolume(aN[0], aN[1], aN[2],
00733           //                       aN[3], aN[4], aN[5]);
00734           aV = myTool->AddVolume(aN[0], aN[1], aN[2], aN[3], aN[4], aN[5]);
00735         }
00736         else {
00737           //aV = meshDS->AddVolume(aN[0], aN[2], aN[1],
00738           //                       aN[3], aN[5], aN[4]);
00739           aV = myTool->AddVolume(aN[0], aN[2], aN[1], aN[3], aN[5], aN[4]);
00740         }
00741         break;
00742       case 4:
00743         if ( forward ) {
00744           //aV = meshDS->AddVolume(aN[0], aN[1], aN[2], aN[3],
00745           //                       aN[4], aN[5], aN[6], aN[7]);
00746           aV = myTool->AddVolume(aN[0], aN[1], aN[2], aN[3],
00747                                  aN[4], aN[5], aN[6], aN[7]);
00748         }
00749         else {
00750           //aV = meshDS->AddVolume(aN[0], aN[3], aN[2], aN[1],
00751           //                       aN[4], aN[7], aN[6], aN[5]);
00752           aV = myTool->AddVolume(aN[0], aN[3], aN[2], aN[1],
00753                                  aN[4], aN[7], aN[6], aN[5]);
00754         }
00755         break;
00756       default:
00757         continue;
00758       }
00759       meshDS->SetMeshElementOnShape(aV, shapeID);
00760     }
00761   }
00762 }
00763 
00764 //=======================================================================
00765 //function : MakeMeshOnFxy1
00766 //purpose  : 
00767 //=======================================================================
00768 void StdMeshers_Penta_3D::MakeMeshOnFxy1()
00769 {
00770   int aID0, aJ, aLevel, ij, aNbNodes, k;
00771   //
00772   SMDS_NodeIteratorPtr itn;
00773   SMDS_ElemIteratorPtr itf, aItNodes;
00774   SMDSAbs_ElementType aElementType;
00775   //
00776   const TopoDS_Face& aFxy0=
00777     TopoDS::Face(myBlock.Shape(SMESH_Block::ID_Fxy0));
00778   const TopoDS_Face& aFxy1=
00779     TopoDS::Face(myBlock.Shape(SMESH_Block::ID_Fxy1));
00780   //
00781   SMESH_Mesh* pMesh = GetMesh();
00782   SMESHDS_Mesh * meshDS = pMesh->GetMeshDS();
00783   //
00784   SMESH_subMesh *aSubMesh1 = pMesh->GetSubMeshContaining(aFxy1);
00785   SMESH_subMesh *aSubMesh0 = pMesh->GetSubMeshContaining(aFxy0);
00786   SMESHDS_SubMesh *aSM0 = aSubMesh0->GetSubMeshDS();
00787   //
00788   // set nodes on aFxy1
00789   aLevel = myISize-1;
00790   itn = aSM0->GetNodes();
00791   aNbNodes = aSM0->NbNodes();
00792   //printf("** aNbNodes=%d\n", aNbNodes);
00793   myTool->SetSubShape( aFxy1 ); // to set medium nodes to aFxy1
00794   //
00795   // set elements on aFxy1
00796   vector<const SMDS_MeshNode*> aNodes1;
00797   //
00798   itf = aSM0->GetElements();
00799   while(itf->more()) {
00800     const SMDS_MeshElement* pE0 = itf->next();
00801     aElementType = pE0->GetType();
00802     if (!aElementType==SMDSAbs_Face) {
00803       continue;
00804     }
00805     aNbNodes = pE0->NbNodes();
00806     if(myCreateQuadratic)
00807       aNbNodes = aNbNodes/2;
00808     if ( aNodes1.size() < aNbNodes )
00809       aNodes1.resize( aNbNodes );
00810     //
00811     k = aNbNodes-1; // reverse a face
00812     aItNodes = pE0->nodesIterator();
00813     while (aItNodes->more()) {
00814       const SMDS_MeshNode* pNode =
00815         static_cast<const SMDS_MeshNode*> (aItNodes->next());
00816       if(myTool->IsMedium(pNode))
00817         continue;
00818       aID0 = pNode->GetID();
00819       aJ = GetIndexOnLayer(aID0);
00820       if (!myErrorStatus->IsOK()) {
00821         MESSAGE("StdMeshers_Penta_3D::MakeMeshOnFxy1() ");
00822         return;
00823       }
00824       //
00825       ij = aLevel*myJSize + aJ;
00826       const StdMeshers_TNode& aTN1 = myTNodes[ij];
00827       const SMDS_MeshNode* aN1 = aTN1.Node();
00828       aNodes1[k] = aN1;
00829       --k;
00830     }
00831     SMDS_MeshFace * face = 0;
00832     switch ( aNbNodes ) {
00833     case 3:
00834       face = myTool->AddFace(aNodes1[0], aNodes1[1], aNodes1[2]);
00835       break;
00836     case 4:
00837       face = myTool->AddFace(aNodes1[0], aNodes1[1], aNodes1[2], aNodes1[3]);
00838       break;
00839     default:
00840       continue;
00841     }
00842     meshDS->SetMeshElementOnShape(face, aFxy1);
00843   }
00844   myTool->SetSubShape( myShape );
00845 
00846   // update compute state of top face submesh
00847   aSubMesh1->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE );
00848 
00849   // assure that mesh on the top face will be cleaned when it is cleaned
00850   // on the bottom face
00851   SMESH_subMesh* volSM = pMesh->GetSubMesh( myTool->GetSubShape() );
00852   volSM->SetEventListener( new SMESH_subMeshEventListener(true, // deletable by SMESH_subMesh
00853                                                           "StdMeshers_Penta_3D"),
00854                            SMESH_subMeshEventListenerData::MakeData( aSubMesh1 ),
00855                            aSubMesh0 ); // translate CLEAN event of aSubMesh0 to aSubMesh1
00856 }
00857 
00858 //=======================================================================
00859 //function : ClearMeshOnFxy1
00860 //purpose  : 
00861 //=======================================================================
00862 void StdMeshers_Penta_3D::ClearMeshOnFxy1()
00863 {
00864   SMESH_subMesh* aSubMesh;
00865   SMESH_Mesh* pMesh=GetMesh();
00866   //
00867   const TopoDS_Shape& aFxy1=myBlock.Shape(SMESH_Block::ID_Fxy1);
00868   aSubMesh = pMesh->GetSubMeshContaining(aFxy1);
00869   if (aSubMesh)
00870     aSubMesh->ComputeStateEngine( SMESH_subMesh::CLEAN );
00871 }
00872 
00873 //=======================================================================
00874 //function : GetIndexOnLayer
00875 //purpose  : 
00876 //=======================================================================
00877 int StdMeshers_Penta_3D::GetIndexOnLayer(const int aID)
00878 {
00879   int j=-1;
00880   StdMeshers_IteratorOfDataMapOfIntegerInteger aMapIt;
00881   //
00882   aMapIt=myConnectingMap.find(aID);
00883   if (aMapIt==myConnectingMap.end()) {
00884     myErrorStatus->myName    = 200;
00885     myErrorStatus->myComment = "Internal error of StdMeshers_Penta_3D";
00886     return j;
00887   }
00888   j=(*aMapIt).second;
00889   return j;
00890 }
00891 
00892 //=======================================================================
00893 //function : MakeConnectingMap
00894 //purpose  : 
00895 //=======================================================================
00896 void StdMeshers_Penta_3D::MakeConnectingMap()
00897 {
00898   int j, aBNID;
00899   //
00900   for (j=0; j<myJSize; ++j) {
00901     const StdMeshers_TNode& aBN=myTNodes[j];
00902     aBNID=aBN.BaseNodeID();
00903     myConnectingMap[aBNID]=j;
00904   }
00905 }
00906 
00907 //=======================================================================
00908 //function : CreateNode
00909 //purpose  : 
00910 //=======================================================================
00911 void StdMeshers_Penta_3D::CreateNode(const bool bIsUpperLayer,
00912                                      const gp_XYZ& aParams,
00913                                      StdMeshers_TNode& aTN)
00914 {
00915   double aX, aY, aZ;
00916   //
00917   gp_Pnt aP;
00918   //
00919   SMDS_MeshNode* pNode=NULL; 
00920   aTN.SetNode(pNode);  
00921   //
00922   //   if (bIsUpperLayer) {
00923   //     // point on face Fxy1
00924   //     const TopoDS_Shape& aS=myBlock.Shape(SMESH_Block::ID_Fxy1);
00925   //     myBlock.Point(aParams, aS, aP);
00926   //   }
00927   //   else {
00928   //     // point inside solid
00929   //     myBlock.Point(aParams, aP);
00930   //   }
00931   if (bIsUpperLayer) {
00932     double u = aParams.X(), v = aParams.Y();
00933     double u1 = ( 1. - u ), v1 = ( 1. - v );
00934     aP.ChangeCoord()  = myShapeXYZ[ SMESH_Block::ID_Ex01 ] * v1;
00935     aP.ChangeCoord() += myShapeXYZ[ SMESH_Block::ID_Ex11 ] * v;
00936     aP.ChangeCoord() += myShapeXYZ[ SMESH_Block::ID_E0y1 ] * u1;
00937     aP.ChangeCoord() += myShapeXYZ[ SMESH_Block::ID_E1y1 ] * u;
00938 
00939     aP.ChangeCoord() -= myShapeXYZ[ SMESH_Block::ID_V001 ] * u1 * v1;
00940     aP.ChangeCoord() -= myShapeXYZ[ SMESH_Block::ID_V101 ] * u  * v1;
00941     aP.ChangeCoord() -= myShapeXYZ[ SMESH_Block::ID_V011 ] * u1 * v;
00942     aP.ChangeCoord() -= myShapeXYZ[ SMESH_Block::ID_V111 ] * u  * v;
00943   }
00944   else {
00945     SMESH_Block::ShellPoint( aParams, myShapeXYZ, aP.ChangeCoord() );
00946   }
00947   //
00948   //   iErr=myBlock.ErrorStatus();
00949   //   if (iErr) {
00950   //     myErrorStatus=12; // can not find the node point;
00951   //     return;
00952   //   }
00953   //
00954   aX=aP.X(); aY=aP.Y(); aZ=aP.Z(); 
00955   //
00956   SMESH_Mesh* pMesh = GetMesh();
00957   SMESHDS_Mesh* pMeshDS = pMesh->GetMeshDS();
00958   //
00959   pNode = pMeshDS->AddNode(aX, aY, aZ);
00960 
00961   aTN.SetNode(pNode);
00962 }
00963 
00964 //=======================================================================
00965 //function : ShapeSupportID
00966 //purpose  : 
00967 //=======================================================================
00968 void StdMeshers_Penta_3D::ShapeSupportID(const bool bIsUpperLayer,
00969                                          const SMESH_Block::TShapeID aBNSSID,
00970                                          SMESH_Block::TShapeID& aSSID)
00971 {
00972   switch (aBNSSID) {
00973   case SMESH_Block::ID_V000:
00974     aSSID=(bIsUpperLayer) ?  SMESH_Block::ID_V001 : SMESH_Block::ID_E00z;
00975     break;
00976   case SMESH_Block::ID_V100:
00977     aSSID=(bIsUpperLayer) ?  SMESH_Block::ID_V101 : SMESH_Block::ID_E10z;
00978     break; 
00979   case SMESH_Block::ID_V110:
00980     aSSID=(bIsUpperLayer) ?  SMESH_Block::ID_V111 : SMESH_Block::ID_E11z;
00981     break;
00982   case SMESH_Block::ID_V010:
00983     aSSID=(bIsUpperLayer) ?  SMESH_Block::ID_V011 : SMESH_Block::ID_E01z;
00984     break;
00985   case SMESH_Block::ID_Ex00:
00986     aSSID=(bIsUpperLayer) ?  SMESH_Block::ID_Ex01 : SMESH_Block::ID_Fx0z;
00987     break;
00988   case SMESH_Block::ID_Ex10:
00989     aSSID=(bIsUpperLayer) ?  SMESH_Block::ID_Ex11 : SMESH_Block::ID_Fx1z;
00990     break; 
00991   case SMESH_Block::ID_E0y0:
00992     aSSID=(bIsUpperLayer) ?  SMESH_Block::ID_E0y1 : SMESH_Block::ID_F0yz;
00993     break; 
00994   case SMESH_Block::ID_E1y0:
00995     aSSID=(bIsUpperLayer) ?  SMESH_Block::ID_E1y1 : SMESH_Block::ID_F1yz;
00996     break; 
00997   case SMESH_Block::ID_Fxy0:
00998     aSSID=SMESH_Block::ID_NONE;//(bIsUpperLayer) ?  Shape_ID_Fxy1 : Shape_ID_NONE;
00999     break;   
01000   default:
01001     aSSID=SMESH_Block::ID_NONE;
01002     myErrorStatus->myName=10; // Can not find supporting shape ID
01003     myErrorStatus->myComment = "Internal error of StdMeshers_Penta_3D";
01004     break;
01005   }
01006   return;
01007 }
01008 //=======================================================================
01009 //function : MakeBlock
01010 //purpose  : 
01011 //=======================================================================
01012 void StdMeshers_Penta_3D::MakeBlock()
01013 {
01014   bool bFound;
01015   int i, j, iNbEV, iNbE, iErr, iCnt, iNbNodes, iNbF;
01016   //
01017   TopoDS_Vertex aV000, aV001;
01018   TopoDS_Shape aFTr;
01019   TopTools_IndexedDataMapOfShapeListOfShape aMVES;
01020   TopTools_IndexedMapOfShape aME ,aMEV, aM;
01021   TopTools_ListIteratorOfListOfShape aIt;
01022   //
01023   TopExp::MapShapes(myShape, TopAbs_FACE, aM);
01024   //
01025   // 0. Find triangulated face aFTr
01026   SMDSAbs_ElementType aElementType;
01027   SMESH_Mesh* pMesh=GetMesh();
01028   //
01029   iCnt = 0;
01030   iNbF = aM.Extent();
01031   for (i=1; i<=iNbF; ++i) {
01032     const TopoDS_Shape& aF = aM(i);
01033     SMESH_subMesh *aSubMesh = pMesh->GetSubMeshContaining(aF);
01034     ASSERT(aSubMesh);
01035     SMESHDS_SubMesh *aSM = aSubMesh->GetSubMeshDS();
01036     SMDS_ElemIteratorPtr itf = aSM->GetElements();
01037     while(itf->more()) {
01038       const SMDS_MeshElement * pElement = itf->next();
01039       aElementType = pElement->GetType();
01040       if (aElementType==SMDSAbs_Face) {
01041         iNbNodes = pElement->NbNodes();
01042         if ( iNbNodes==3 || (pElement->IsQuadratic() && iNbNodes==6) ) {
01043           aFTr = aF;
01044           ++iCnt;
01045           if (iCnt>1) {
01046             // \begin{E.A.}
01047             // The current algorithm fails if there is more that one
01048             // face wich contains triangles ...
01049             // In that case, replace return by break to try another
01050             // method (coded in "if (iCnt != 1) { ... }")
01051             //
01052             // MESSAGE("StdMeshers_Penta_3D::MakeBlock() ");
01053             // myErrorStatus=5; // more than one face has triangulation
01054             // return;
01055             break;
01056             // \end{E.A.}
01057           }
01058           break; // next face
01059         }
01060       }
01061     }
01062   }
01063   //
01064   // \begin{E.A.}
01065   // The current algorithm fails if "iCnt != 1", the case "iCnt == 0"
01066   // was not reached 'cause it was not called from Hexa_3D ... Now it
01067   // can occurs and in my opinion, it is the most common case.
01068   //
01069   if (iCnt != 1) {
01070     // The suggested algorithm is the following :
01071     //
01072     // o Check that nb_of_faces == 6 and nb_of_edges == 12
01073     //   then the shape is tologically equivalent to a box
01074     // o In a box, there are three set of four // edges ...
01075     //   In the cascade notation, it seems to be the edges
01076     //   numbered : 
01077     //     - 1, 3, 5, 7
01078     //     - 2, 4, 6, 8
01079     //     - 9, 10, 11, 12
01080     // o For each one of this set, check if the four edges
01081     //   have the same number of element.
01082     // o If so, check if the "corresponding" // faces contains
01083     //   only quads. It's the faces numbered:
01084     //     - 1, 2, 3, 4
01085     //     - 1, 2, 5, 6
01086     //     - 3, 4, 5, 6
01087     // o If so, check if the opposite edges of each // faces
01088     //   have the same number of elements. It is the edges
01089     //   numbered :
01090     //     - 2 and 4, 6 and 8, 9 and 10, 11 and 12
01091     //     - 1 and 3, 5 and 7, 9 and 11, 10 and 12
01092     //     - 1 and 5, 3 and 7, 4 and 8, 2 and 6
01093     // o If so, check if the two other faces have the same
01094     //   number of elements. It is the faces numbered:
01095     //     - 5, 6
01096     //     - 3, 4
01097     //     - 1, 2
01098     //   This test should be improved to test if the nodes
01099     //   of the two faces are really "en face".
01100     // o If so, one of the two faces is a candidate to an extrusion,
01101     //   It is the faces numbered :
01102     //     - 5
01103     //     - 3
01104     //     - 1
01105     // o Finally, if there is only one candidate, let do the
01106     //   extrusion job for the corresponding face
01107     //
01108     int isOK = 0;
01109     //
01110     int iNbF = aM.Extent();
01111     if (iNbF == 6) {
01112       //
01113       int nb_f1 = pMesh->GetSubMeshContaining(aM(1))->GetSubMeshDS()->NbElements();
01114       int nb_f2 = pMesh->GetSubMeshContaining(aM(2))->GetSubMeshDS()->NbElements();
01115       int nb_f3 = pMesh->GetSubMeshContaining(aM(3))->GetSubMeshDS()->NbElements();
01116       int nb_f4 = pMesh->GetSubMeshContaining(aM(4))->GetSubMeshDS()->NbElements();
01117       int nb_f5 = pMesh->GetSubMeshContaining(aM(5))->GetSubMeshDS()->NbElements();
01118       int nb_f6 = pMesh->GetSubMeshContaining(aM(6))->GetSubMeshDS()->NbElements();
01119       //
01120       int has_only_quad_f1 = 1;
01121       int has_only_quad_f2 = 1;
01122       int has_only_quad_f3 = 1;
01123       int has_only_quad_f4 = 1;
01124       int has_only_quad_f5 = 1;
01125       int has_only_quad_f6 = 1;
01126       //
01127       for (i=1; i<=iNbF; ++i) {
01128         int ok = 1;
01129         const TopoDS_Shape& aF = aM(i);
01130         SMESH_subMesh *aSubMesh = pMesh->GetSubMeshContaining(aF);
01131         SMESHDS_SubMesh *aSM = aSubMesh->GetSubMeshDS();
01132         SMDS_ElemIteratorPtr itf = aSM->GetElements();
01133         while(itf->more()) {
01134           const SMDS_MeshElement * pElement = itf->next();
01135           aElementType = pElement->GetType();
01136           if (aElementType==SMDSAbs_Face) {
01137             iNbNodes = pElement->NbNodes();
01138             if ( iNbNodes!=4 ) {
01139               ok = 0;
01140               break ;
01141             }
01142           }
01143         }
01144         if (i==1) has_only_quad_f1 = ok ;
01145         if (i==2) has_only_quad_f2 = ok ;
01146         if (i==3) has_only_quad_f3 = ok ;
01147         if (i==4) has_only_quad_f4 = ok ;
01148         if (i==5) has_only_quad_f5 = ok ;
01149         if (i==6) has_only_quad_f6 = ok ;
01150       }
01151       //
01152       TopTools_IndexedMapOfShape aE;
01153       TopExp::MapShapes(myShape, TopAbs_EDGE, aE);
01154       int iNbE = aE.Extent();
01155       if (iNbE == 12) {
01156         //
01157         int nb_e01 = pMesh->GetSubMeshContaining(aE(1))->GetSubMeshDS()->NbElements();
01158         int nb_e02 = pMesh->GetSubMeshContaining(aE(2))->GetSubMeshDS()->NbElements();
01159         int nb_e03 = pMesh->GetSubMeshContaining(aE(3))->GetSubMeshDS()->NbElements();
01160         int nb_e04 = pMesh->GetSubMeshContaining(aE(4))->GetSubMeshDS()->NbElements();
01161         int nb_e05 = pMesh->GetSubMeshContaining(aE(5))->GetSubMeshDS()->NbElements();
01162         int nb_e06 = pMesh->GetSubMeshContaining(aE(6))->GetSubMeshDS()->NbElements();
01163         int nb_e07 = pMesh->GetSubMeshContaining(aE(7))->GetSubMeshDS()->NbElements();
01164         int nb_e08 = pMesh->GetSubMeshContaining(aE(8))->GetSubMeshDS()->NbElements();
01165         int nb_e09 = pMesh->GetSubMeshContaining(aE(9))->GetSubMeshDS()->NbElements();
01166         int nb_e10 = pMesh->GetSubMeshContaining(aE(10))->GetSubMeshDS()->NbElements();
01167         int nb_e11 = pMesh->GetSubMeshContaining(aE(11))->GetSubMeshDS()->NbElements();
01168         int nb_e12 = pMesh->GetSubMeshContaining(aE(12))->GetSubMeshDS()->NbElements();
01169         //
01170         int nb_ok = 0 ;
01171         //
01172         if ( (nb_e01==nb_e03) && (nb_e03==nb_e05) && (nb_e05==nb_e07) ) {
01173           if ( has_only_quad_f1 && has_only_quad_f2 && has_only_quad_f3 && has_only_quad_f4 ) {
01174             if ( (nb_e09==nb_e10) && (nb_e08==nb_e06) && (nb_e11==nb_e12) && (nb_e04==nb_e02) ) {
01175               if (nb_f5==nb_f6) {
01176                 nb_ok += 1;
01177                 aFTr = aM(5);
01178               }
01179             }
01180           }
01181         }
01182         if ( (nb_e02==nb_e04) && (nb_e04==nb_e06) && (nb_e06==nb_e08) ) {
01183           if ( has_only_quad_f1 && has_only_quad_f2 && has_only_quad_f5 && has_only_quad_f6 ) {
01184             if ( (nb_e01==nb_e03) && (nb_e10==nb_e12) && (nb_e05==nb_e07) && (nb_e09==nb_e11) ) {
01185               if (nb_f3==nb_f4) {
01186                 nb_ok += 1;
01187                 aFTr = aM(3);
01188               }
01189             }
01190           }
01191         }
01192         if ( (nb_e09==nb_e10) && (nb_e10==nb_e11) && (nb_e11==nb_e12) ) {
01193           if ( has_only_quad_f3 && has_only_quad_f4 && has_only_quad_f5 && has_only_quad_f6 ) {
01194             if ( (nb_e01==nb_e05) && (nb_e02==nb_e06) && (nb_e03==nb_e07) && (nb_e04==nb_e08) ) {
01195               if (nb_f1==nb_f2) {
01196                 nb_ok += 1;
01197                 aFTr = aM(1);
01198               }
01199             }
01200           }
01201         }
01202         //
01203         if ( nb_ok == 1 ) {
01204           isOK = 1;
01205         }
01206         //
01207       }
01208     }
01209     if (!isOK) {
01210       myErrorStatus->myName=5; // more than one face has triangulation
01211       myErrorStatus->myComment="Incorrect input mesh";
01212       return;
01213     }
01214   }
01215   // \end{E.A.}
01216   // 
01217   // 1. Vetrices V00, V001;
01218   //
01219   TopExp::MapShapes(aFTr, TopAbs_EDGE, aME);
01220   TopExp::MapShapesAndAncestors(myShape, TopAbs_VERTEX, TopAbs_EDGE, aMVES);
01221   //
01222   // 1.1 Base vertex V000
01223   iNbE = aME.Extent();
01224   if (iNbE!= NB_WALL_FACES ){
01225     MESSAGE("StdMeshers_Penta_3D::MakeBlock() ");
01226     myErrorStatus->myName=7; // too few edges are in base face aFTr
01227     myErrorStatus->myComment=SMESH_Comment("Not a quadrilateral face #")
01228       <<pMesh->GetMeshDS()->ShapeToIndex( aFTr )<<": "<<iNbE<<" edges" ;
01229     return;
01230   }
01231   const TopoDS_Edge& aE1=TopoDS::Edge(aME(1));
01232   aV000=TopExp::FirstVertex(aE1);
01233   //
01234   const TopTools_ListOfShape& aLE=aMVES.FindFromKey(aV000);
01235   aIt.Initialize(aLE);
01236   for (; aIt.More(); aIt.Next()) {
01237     const TopoDS_Shape& aEx=aIt.Value();
01238     aMEV.Add(aEx);
01239   }
01240   iNbEV=aMEV.Extent();
01241   if (iNbEV!=3){
01242     MESSAGE("StdMeshers_Penta_3D::MakeBlock() ");
01243     myErrorStatus->myName=7; // too few edges meet in base vertex 
01244     myErrorStatus->myComment=SMESH_Comment("3 edges must share vertex #")
01245       <<pMesh->GetMeshDS()->ShapeToIndex( aV000 )<<" but there are "<<iNbEV<<" edges";
01246     return;
01247   }
01248   //
01249   // 1.2 Vertex V001
01250   bFound=false;
01251   for (j=1; j<=iNbEV; ++j) {
01252     const TopoDS_Shape& aEx=aMEV(j);
01253     if (!aME.Contains(aEx)) {
01254       TopoDS_Vertex aV[2];
01255       //
01256       const TopoDS_Edge& aE=TopoDS::Edge(aEx);
01257       TopExp::Vertices(aE, aV[0], aV[1]);
01258       for (i=0; i<2; ++i) {
01259         if (!aV[i].IsSame(aV000)) {
01260           aV001=aV[i];
01261           bFound=!bFound;
01262           break;
01263         }
01264       }
01265     }
01266   }
01267   //
01268   if (!bFound) {
01269     MESSAGE("StdMeshers_Penta_3D::MakeBlock() ");
01270     myErrorStatus->myName=8; // can not find reper V001
01271     myErrorStatus->myComment=SMESH_Comment("Can't find opposite vertex for vertex #")
01272       <<pMesh->GetMeshDS()->ShapeToIndex( aV000 );
01273     return;
01274   }
01275   //DEB
01276   //gp_Pnt aP000, aP001;
01277   //
01278   //aP000=BRep_Tool::Pnt(TopoDS::Vertex(aV000));
01279   //printf("*** aP000 { %lf, %lf, %lf }\n", aP000.X(), aP000.Y(), aP000.Z());
01280   //aP001=BRep_Tool::Pnt(TopoDS::Vertex(aV001));
01281   //printf("*** aP001 { %lf, %lf, %lf }\n", aP001.X(), aP001.Y(), aP001.Z());
01282   //DEB
01283   //
01284   aME.Clear();
01285   TopExp::MapShapes(myShape, TopAbs_SHELL, aME);
01286   iNbE=aME.Extent();
01287   if (iNbE!=1) {
01288     MESSAGE("StdMeshers_Penta_3D::MakeBlock() ");
01289     myErrorStatus->myName=9; // number of shells in source shape !=1
01290     myErrorStatus->myComment=SMESH_Comment("Unexpected nb of shells ")<<iNbE;
01291     return;
01292   }
01293   //
01294   // 2. Load Block
01295   const TopoDS_Shell& aShell=TopoDS::Shell(aME(1));
01296   myBlock.Load(aShell, aV000, aV001);
01297   iErr = myBlock.ErrorStatus();
01298   if (iErr) {
01299     MESSAGE("StdMeshers_Penta_3D::MakeBlock() ");
01300     myErrorStatus=myBlock.GetError(); // SMESHBlock: Load operation failed
01301     return;
01302   }
01303 }
01304 //=======================================================================
01305 //function : CheckData
01306 //purpose  : 
01307 //=======================================================================
01308 void StdMeshers_Penta_3D::CheckData()
01309 {
01310   int i, iNb;
01311   int iNbEx[]={8, 12, 6};
01312   //
01313   TopAbs_ShapeEnum aST;
01314   TopAbs_ShapeEnum aSTEx[]={
01315     TopAbs_VERTEX, TopAbs_EDGE, TopAbs_FACE
01316   }; 
01317   TopTools_IndexedMapOfShape aM;
01318   //
01319   if (myShape.IsNull()){
01320     MESSAGE("StdMeshers_Penta_3D::CheckData() ");
01321     myErrorStatus->myName=2; // null shape
01322     myErrorStatus->myComment="Null shape";
01323     return;
01324   }
01325   //
01326   aST=myShape.ShapeType();
01327   if (!(aST==TopAbs_SOLID || aST==TopAbs_SHELL)) {
01328     MESSAGE("StdMeshers_Penta_3D::CheckData() ");
01329     myErrorStatus->myName=3; // not compatible type of shape
01330     myErrorStatus->myComment=SMESH_Comment("Wrong shape type (TopAbs_ShapeEnum) ")<<aST;
01331     return;
01332   }
01333   //
01334   for (i=0; i<3; ++i) {
01335     aM.Clear();
01336     TopExp::MapShapes(myShape, aSTEx[i], aM);
01337     iNb=aM.Extent();
01338     if (iNb!=iNbEx[i]){
01339       MESSAGE("StdMeshers_Penta_3D::CheckData() ");
01340       myErrorStatus->myName=4; // number of sub-shape is not compatible
01341       myErrorStatus->myComment="Wrong number of sub-shapes of a block";
01342       return;
01343     }
01344   }
01345 }
01346 
01347 //=======================================================================
01348 //function : LoadIJNodes
01349 //purpose  : Load nodes bound to theFace into column (vectors) and rows
01350 //           of theIJNodes.
01351 //           The value of theIJNodes map is a vector of ordered nodes so
01352 //           that the 0-the one lies on theBaseEdge.
01353 //           The key of theIJNodes map is a normalized parameter of each
01354 //           0-the node on theBaseEdge.
01355 //=======================================================================
01356 
01357 bool StdMeshers_Penta_3D::LoadIJNodes(StdMeshers_IJNodeMap & theIJNodes,
01358                                       const TopoDS_Face&     theFace,
01359                                       const TopoDS_Edge&     theBaseEdge,
01360                                       SMESHDS_Mesh*          theMesh)
01361 {
01362   // get vertices of theBaseEdge
01363   TopoDS_Vertex vfb, vlb, vft; // first and last, bottom and top vertices
01364   TopoDS_Edge eFrw = TopoDS::Edge( theBaseEdge.Oriented( TopAbs_FORWARD ));
01365   TopExp::Vertices( eFrw, vfb, vlb );
01366 
01367   // find the other edges of theFace and orientation of e1
01368   TopoDS_Edge e1, e2, eTop;
01369   bool rev1, CumOri = false;
01370   TopExp_Explorer exp( theFace, TopAbs_EDGE );
01371   int nbEdges = 0;
01372   for ( ; exp.More(); exp.Next() ) {
01373     if ( ++nbEdges > NB_WALL_FACES ) {
01374       return false; // more than 4 edges in theFace
01375     }
01376     TopoDS_Edge e = TopoDS::Edge( exp.Current() );
01377     if ( theBaseEdge.IsSame( e ))
01378       continue;
01379     TopoDS_Vertex vCommon;
01380     if ( !TopExp::CommonVertex( theBaseEdge, e, vCommon ))
01381       eTop = e;
01382     else if ( vCommon.IsSame( vfb )) {
01383       e1 = e;
01384       vft = TopExp::LastVertex( e1, CumOri );
01385       rev1 = vfb.IsSame( vft );
01386       if ( rev1 )
01387         vft = TopExp::FirstVertex( e1, CumOri );
01388     }
01389     else
01390       e2 = e;
01391   }
01392   if ( nbEdges < NB_WALL_FACES ) {
01393     return false; // less than 4 edges in theFace
01394   }
01395 
01396   // submeshes corresponding to shapes
01397   SMESHDS_SubMesh* smFace = theMesh->MeshElements( theFace );
01398   SMESHDS_SubMesh* smb = theMesh->MeshElements( theBaseEdge );
01399   SMESHDS_SubMesh* smt = theMesh->MeshElements( eTop );
01400   SMESHDS_SubMesh* sm1 = theMesh->MeshElements( e1 );
01401   SMESHDS_SubMesh* sm2 = theMesh->MeshElements( e2 );
01402   SMESHDS_SubMesh* smVfb = theMesh->MeshElements( vfb );
01403   SMESHDS_SubMesh* smVlb = theMesh->MeshElements( vlb );
01404   SMESHDS_SubMesh* smVft = theMesh->MeshElements( vft );
01405   if (!smFace || !smb || !smt || !sm1 || !sm2 || !smVfb || !smVlb || !smVft ) {
01406     MESSAGE( "NULL submesh " <<smFace<<" "<<smb<<" "<<smt<<" "<<
01407              sm1<<" "<<sm2<<" "<<smVfb<<" "<<smVlb<<" "<<smVft);
01408     return false;
01409   }
01410   if ( smb->NbNodes() != smt->NbNodes() || sm1->NbNodes() != sm2->NbNodes() ) {
01411     MESSAGE(" Diff nb of nodes on opposite edges" );
01412     return false;
01413   }
01414   if (smVfb->NbNodes() != 1 || smVlb->NbNodes() != 1 || smVft->NbNodes() != 1) {
01415     MESSAGE("Empty submesh of vertex");
01416     return false;
01417   }
01418   if ( sm1->NbNodes() * smb->NbNodes() != smFace->NbNodes() ) {
01419     // check quadratic case
01420     if ( myCreateQuadratic ) {
01421       int n1 = sm1->NbNodes()/2;
01422       int n2 = smb->NbNodes()/2;
01423       int n3 = sm1->NbNodes() - n1;
01424       int n4 = smb->NbNodes() - n2;
01425       int nf = sm1->NbNodes()*smb->NbNodes() - n3*n4;
01426       if( nf != smFace->NbNodes() ) {
01427         MESSAGE( "Wrong nb face nodes: " <<
01428                  sm1->NbNodes()<<" "<<smb->NbNodes()<<" "<<smFace->NbNodes());
01429         return false;
01430       }
01431     }
01432     else {
01433       MESSAGE( "Wrong nb face nodes: " <<
01434                sm1->NbNodes()<<" "<<smb->NbNodes()<<" "<<smFace->NbNodes());
01435       return false;
01436     }
01437   }
01438   // IJ size
01439   int vsize = sm1->NbNodes() + 2;
01440   int hsize = smb->NbNodes() + 2;
01441   if(myCreateQuadratic) {
01442     vsize = vsize - sm1->NbNodes()/2 -1;
01443     hsize = hsize - smb->NbNodes()/2 -1;
01444   }
01445 
01446   // load nodes from theBaseEdge
01447 
01448   set<const SMDS_MeshNode*> loadedNodes;
01449   const SMDS_MeshNode* nullNode = 0;
01450 
01451   vector<const SMDS_MeshNode*> & nVecf = theIJNodes[ 0.];
01452   nVecf.resize( vsize, nullNode );
01453   loadedNodes.insert( nVecf[ 0 ] = smVfb->GetNodes()->next() );
01454 
01455   vector<const SMDS_MeshNode*> & nVecl = theIJNodes[ 1.];
01456   nVecl.resize( vsize, nullNode );
01457   loadedNodes.insert( nVecl[ 0 ] = smVlb->GetNodes()->next() );
01458 
01459   double f, l;
01460   BRep_Tool::Range( eFrw, f, l );
01461   double range = l - f;
01462   SMDS_NodeIteratorPtr nIt = smb->GetNodes();
01463   const SMDS_MeshNode* node;
01464   while ( nIt->more() ) {
01465     node = nIt->next();
01466     if(myTool->IsMedium(node))
01467       continue;
01468     const SMDS_EdgePosition* pos =
01469       dynamic_cast<const SMDS_EdgePosition*>( node->GetPosition() );
01470     if ( !pos ) {
01471       return false;
01472     }
01473     double u = ( pos->GetUParameter() - f ) / range;
01474     vector<const SMDS_MeshNode*> & nVec = theIJNodes[ u ];
01475     nVec.resize( vsize, nullNode );
01476     loadedNodes.insert( nVec[ 0 ] = node );
01477   }
01478   if ( theIJNodes.size() != hsize ) {
01479     MESSAGE( "Wrong node positions on theBaseEdge" );
01480     return false;
01481   }
01482 
01483   // load nodes from e1
01484 
01485   map< double, const SMDS_MeshNode*> sortedNodes; // sort by param on edge
01486   nIt = sm1->GetNodes();
01487   while ( nIt->more() ) {
01488     node = nIt->next();
01489     if(myTool->IsMedium(node))
01490       continue;
01491     const SMDS_EdgePosition* pos =
01492       dynamic_cast<const SMDS_EdgePosition*>( node->GetPosition() );
01493     if ( !pos ) {
01494       return false;
01495     }
01496     sortedNodes.insert( make_pair( pos->GetUParameter(), node ));
01497   }
01498   loadedNodes.insert( nVecf[ vsize - 1 ] = smVft->GetNodes()->next() );
01499   map< double, const SMDS_MeshNode*>::iterator u_n = sortedNodes.begin();
01500   int row = rev1 ? vsize - 1 : 0;
01501   for ( ; u_n != sortedNodes.end(); u_n++ ) {
01502     if ( rev1 ) row--;
01503     else        row++;
01504     loadedNodes.insert( nVecf[ row ] = u_n->second );
01505   }
01506 
01507   // try to load the rest nodes
01508 
01509   // get all faces from theFace
01510   TIDSortedElemSet allFaces, foundFaces;
01511   SMDS_ElemIteratorPtr eIt = smFace->GetElements();
01512   while ( eIt->more() ) {
01513     const SMDS_MeshElement* e = eIt->next();
01514     if ( e->GetType() == SMDSAbs_Face )
01515       allFaces.insert( e );
01516   }
01517   // Starting from 2 neighbour nodes on theBaseEdge, look for a face
01518   // the nodes belong to, and between the nodes of the found face,
01519   // look for a not loaded node considering this node to be the next
01520   // in a column of the starting second node. Repeat, starting
01521   // from nodes next to the previous starting nodes in their columns,
01522   // and so on while a face can be found. Then go the the next pair
01523   // of nodes on theBaseEdge.
01524   StdMeshers_IJNodeMap::iterator par_nVec_1 = theIJNodes.begin();
01525   StdMeshers_IJNodeMap::iterator par_nVec_2 = par_nVec_1;
01526   // loop on columns
01527   int col = 0;
01528   for ( par_nVec_2++; par_nVec_2 != theIJNodes.end(); par_nVec_1++, par_nVec_2++ ) {
01529     col++;
01530     row = 0;
01531     const SMDS_MeshNode* n1 = par_nVec_1->second[ row ];
01532     const SMDS_MeshNode* n2 = par_nVec_2->second[ row ];
01533     const SMDS_MeshElement* face = 0;
01534     do {
01535       // look for a face by 2 nodes
01536       face = SMESH_MeshEditor::FindFaceInSet( n1, n2, allFaces, foundFaces );
01537       if ( face ) {
01538         int nbFaceNodes = face->NbNodes();
01539         if ( (!myCreateQuadratic && nbFaceNodes>4) ||
01540              (myCreateQuadratic && nbFaceNodes>8) ) {
01541           MESSAGE(" Too many nodes in a face: " << nbFaceNodes );
01542           return false;
01543         }
01544         // look for a not loaded node of the <face>
01545         bool found = false;
01546         const SMDS_MeshNode* n3 = 0; // a node defferent from n1 and n2
01547         eIt = face->nodesIterator() ;
01548         while ( !found && eIt->more() ) {
01549           node = static_cast<const SMDS_MeshNode*>( eIt->next() );
01550           if(myTool->IsMedium(node))
01551             continue;
01552           found = loadedNodes.insert( node ).second;
01553           if ( !found && node != n1 && node != n2 )
01554             n3 = node;
01555         }
01556         if ( found ) {
01557           if ( ++row > vsize - 1 ) {
01558             MESSAGE( "Too many nodes in column "<< col <<": "<< row+1);
01559             return false;
01560           }
01561           par_nVec_2->second[ row ] = node;
01562           foundFaces.insert( face );
01563           n2 = node;
01564           if ( nbFaceNodes==4 || (myCreateQuadratic && nbFaceNodes==8) ) {
01565             n1 = par_nVec_1->second[ row ];
01566           }
01567         }
01568         else if ( (nbFaceNodes==3 || (myCreateQuadratic && nbFaceNodes==6) )  &&
01569                   n3 == par_nVec_1->second[ row ] ) {
01570           n1 = n3;
01571         }
01572         else {
01573           MESSAGE( "Not quad mesh, column "<< col );
01574           return false;
01575         }
01576       }
01577     }
01578     while ( face && n1 && n2 );
01579 
01580     if ( row < vsize - 1 ) {
01581       MESSAGE( "Too few nodes in column "<< col <<": "<< row+1);
01582       MESSAGE( "Base node 1: "<< par_nVec_1->second[0]);
01583       MESSAGE( "Base node 2: "<< par_nVec_2->second[0]);
01584       MESSAGE( "Current node 1: "<< n1);
01585       MESSAGE( "Current node 2: "<< n2);
01586       MESSAGE( "first base node: "<< theIJNodes.begin()->second[0]);
01587       MESSAGE( "last base node: "<< theIJNodes.rbegin()->second[0]);
01588       return false;
01589     }
01590   } // loop on columns
01591 
01592   return true;
01593 }
01594 
01596 //
01597 //   StdMeshers_SMESHBlock
01598 //
01600 
01601 //=======================================================================
01602 //function : StdMeshers_SMESHBlock
01603 //purpose  : 
01604 //=======================================================================
01605 StdMeshers_SMESHBlock::StdMeshers_SMESHBlock()
01606 {
01607   myErrorStatus=1;
01608   myIsEdgeForward.resize( SMESH_Block::NbEdges(), -1 );
01609 }
01610 
01611 //=======================================================================
01612 //function : IsForwadEdge
01613 //purpose  : 
01614 //=======================================================================
01615 
01616 bool StdMeshers_SMESHBlock::IsForwadEdge(const int theEdgeID)
01617 {
01618   int index = myTBlock.ShapeIndex( theEdgeID );
01619   if ( !myTBlock.IsEdgeID( theEdgeID ))
01620     return false;
01621 
01622   if ( myIsEdgeForward[ index ] < 0 )
01623     myIsEdgeForward[ index ] =
01624       myTBlock.IsForwardEdge( TopoDS::Edge( Shape( theEdgeID )), myShapeIDMap );
01625 
01626   return myIsEdgeForward[ index ];
01627 }
01628 
01629 //=======================================================================
01630 //function : ErrorStatus
01631 //purpose  : 
01632 //=======================================================================
01633 int StdMeshers_SMESHBlock::ErrorStatus() const
01634 {
01635   return myErrorStatus;
01636 }
01637 
01638 //================================================================================
01642 //================================================================================
01643 
01644 SMESH_ComputeErrorPtr StdMeshers_SMESHBlock::GetError() const
01645 {
01646   SMESH_ComputeErrorPtr err = SMESH_ComputeError::New();
01647   string & text = err->myComment;
01648   switch ( myErrorStatus ) {
01649   case 2:
01650   case 3: text = "Internal error of StdMeshers_Penta_3D"; break; 
01651   case 4: text = "Can't compute normalized parameters of a point inside a block"; break;
01652   case 5: text = "Can't compute coordinates by normalized parameters inside a block"; break;
01653   case 6: text = "Can't detect block sub-shapes. Not a block?"; break;
01654   }
01655   if (!text.empty())
01656     err->myName = myErrorStatus;
01657   return err;
01658 }
01659 
01660 //=======================================================================
01661 //function : Load
01662 //purpose  : 
01663 //=======================================================================
01664 void StdMeshers_SMESHBlock::Load(const TopoDS_Shell& theShell)
01665 {
01666   TopoDS_Vertex aV000, aV001;
01667   //
01668   Load(theShell, aV000, aV001);
01669 }
01670 
01671 //=======================================================================
01672 //function : Load
01673 //purpose  : 
01674 //=======================================================================
01675 void StdMeshers_SMESHBlock::Load(const TopoDS_Shell& theShell,
01676                                  const TopoDS_Vertex& theV000,
01677                                  const TopoDS_Vertex& theV001)
01678 {
01679   myErrorStatus=0;
01680   //
01681   myShell=theShell;
01682   //
01683   bool bOk;
01684   //
01685   myShapeIDMap.Clear();  
01686   bOk = myTBlock.LoadBlockShapes(myShell, theV000, theV001, myShapeIDMap);
01687   if (!bOk) {
01688     myErrorStatus=6;
01689     return;
01690   }
01691 }
01692 
01693 //=======================================================================
01694 //function : ComputeParameters
01695 //purpose  : 
01696 //=======================================================================
01697 void StdMeshers_SMESHBlock::ComputeParameters(const gp_Pnt& thePnt, 
01698                                               gp_XYZ& theXYZ)
01699 {
01700   ComputeParameters(thePnt, myShell, theXYZ);
01701 }
01702 
01703 //=======================================================================
01704 //function : ComputeParameters
01705 //purpose  : 
01706 //=======================================================================
01707 void StdMeshers_SMESHBlock::ComputeParameters(const gp_Pnt& thePnt,
01708                                               const TopoDS_Shape& theShape,
01709                                               gp_XYZ& theXYZ)
01710 {
01711   myErrorStatus=0;
01712   //
01713   int aID;
01714   bool bOk;
01715   //
01716   aID = ShapeID(theShape);
01717   if (myErrorStatus) {
01718     return;
01719   }
01720   bOk = myTBlock.ComputeParameters(thePnt, theXYZ, aID);
01721   if (!bOk) {
01722     myErrorStatus=4; // problems with computation Parameters 
01723     return;
01724   }
01725 }
01726 
01727 //=======================================================================
01728 //function : ComputeParameters
01729 //purpose  : 
01730 //=======================================================================
01731 
01732 void StdMeshers_SMESHBlock::ComputeParameters(const double& theU,
01733                                               const TopoDS_Shape& theShape,
01734                                               gp_XYZ& theXYZ)
01735 {
01736   myErrorStatus=0;
01737   //
01738   int aID;
01739   bool bOk=false;
01740   //
01741   aID = ShapeID(theShape);
01742   if (myErrorStatus) {
01743     return;
01744   }
01745   if ( SMESH_Block::IsEdgeID( aID ))
01746     bOk = myTBlock.EdgeParameters( aID, theU, theXYZ );
01747   if (!bOk) {
01748     myErrorStatus=4; // problems with computation Parameters 
01749     return;
01750   }
01751 }
01752 
01753 //=======================================================================
01754 //function : Point
01755 //purpose  : 
01756 //=======================================================================
01757 void StdMeshers_SMESHBlock::Point(const gp_XYZ& theParams, gp_Pnt& aP3D)
01758 {
01759   TopoDS_Shape aS;
01760   //
01761   Point(theParams, aS, aP3D);
01762 }
01763 
01764 //=======================================================================
01765 //function : Point
01766 //purpose  : 
01767 //=======================================================================
01768 void StdMeshers_SMESHBlock::Point(const gp_XYZ& theParams,
01769                                   const TopoDS_Shape& theShape,
01770                                   gp_Pnt& aP3D)
01771 {
01772   myErrorStatus = 0;
01773   //
01774   int aID;
01775   bool bOk = false;
01776   gp_XYZ aXYZ(99.,99.,99.);
01777   aP3D.SetXYZ(aXYZ);
01778   //
01779   if (theShape.IsNull()) {
01780     bOk = myTBlock.ShellPoint(theParams, aXYZ);
01781   }
01782   //
01783   else {
01784     aID=ShapeID(theShape);
01785     if (myErrorStatus) {
01786       return;
01787     }
01788     //
01789     if (SMESH_Block::IsVertexID(aID)) {
01790       bOk = myTBlock.VertexPoint(aID, aXYZ);
01791     }
01792     else if (SMESH_Block::IsEdgeID(aID)) {
01793       bOk = myTBlock.EdgePoint(aID, theParams, aXYZ);
01794     }
01795     //
01796     else if (SMESH_Block::IsFaceID(aID)) {
01797       bOk = myTBlock.FacePoint(aID, theParams, aXYZ);
01798     }
01799   }
01800   if (!bOk) {
01801     myErrorStatus=5; // problems with point computation 
01802     return;
01803   }
01804   aP3D.SetXYZ(aXYZ);
01805 }
01806 
01807 //=======================================================================
01808 //function : ShapeID
01809 //purpose  : 
01810 //=======================================================================
01811 int StdMeshers_SMESHBlock::ShapeID(const TopoDS_Shape& theShape)
01812 {
01813   myErrorStatus=0;
01814   //
01815   int aID=-1;
01816   TopoDS_Shape aSF, aSR;
01817   //
01818   aSF=theShape;
01819   aSF.Orientation(TopAbs_FORWARD);
01820   aSR=theShape;
01821   aSR.Orientation(TopAbs_REVERSED);
01822   //
01823   if (myShapeIDMap.Contains(aSF)) {
01824     aID=myShapeIDMap.FindIndex(aSF);
01825     return aID;
01826   }
01827   if (myShapeIDMap.Contains(aSR)) {
01828     aID=myShapeIDMap.FindIndex(aSR);
01829     return aID;
01830   }
01831   myErrorStatus=2; // unknown shape;
01832   return aID;
01833 }
01834 
01835 //=======================================================================
01836 //function : Shape
01837 //purpose  : 
01838 //=======================================================================
01839 const TopoDS_Shape& StdMeshers_SMESHBlock::Shape(const int theID)
01840 {
01841   myErrorStatus=0;
01842   //
01843   int aNb;
01844   //
01845   aNb=myShapeIDMap.Extent();
01846   if (theID<1 || theID>aNb) {
01847     myErrorStatus=3; // ID is out of range
01848     return myEmptyShape;
01849   }
01850   //
01851   const TopoDS_Shape& aS=myShapeIDMap.FindKey(theID);
01852   return aS;
01853 }
01854 
01855 
01856 //=======================================================================
01857 //function : Evaluate
01858 //purpose  : 
01859 //=======================================================================
01860 bool StdMeshers_Penta_3D::Evaluate(SMESH_Mesh& aMesh, 
01861                                    const TopoDS_Shape& aShape,
01862                                    MapShapeNbElems& aResMap)
01863 {
01864   MESSAGE("StdMeshers_Penta_3D::Evaluate()");
01865 
01866   // find face contains only triangles
01867   vector < SMESH_subMesh * >meshFaces;
01868   TopTools_SequenceOfShape aFaces;
01869   int NumBase = 0, i = 0;
01870   for (TopExp_Explorer exp(aShape, TopAbs_FACE); exp.More(); exp.Next()) {
01871     i++;
01872     aFaces.Append(exp.Current());
01873     SMESH_subMesh *aSubMesh = aMesh.GetSubMesh(exp.Current());
01874     meshFaces.push_back(aSubMesh);
01875     MapShapeNbElemsItr anIt = aResMap.find(meshFaces[i]);
01876     if( anIt == aResMap.end() ) {
01877       NumBase = 0;
01878       break;
01879     }
01880     std::vector<int> aVec = (*anIt).second;
01881     int nbtri = Max(aVec[SMDSEntity_Triangle],aVec[SMDSEntity_Quad_Triangle]);
01882     int nbqua = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
01883     if( nbtri>0 && nbqua==0 ) {
01884       NumBase = i;
01885     }
01886   }
01887 
01888   if(NumBase==0) {
01889     std::vector<int> aResVec(SMDSEntity_Last);
01890     for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
01891     SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
01892     aResMap.insert(std::make_pair(sm,aResVec));
01893     myErrorStatus->myName    = COMPERR_ALGO_FAILED;
01894     myErrorStatus->myComment = "Submesh can not be evaluated";
01895     return false;
01896   }
01897 
01898   // find number of 1d elems for base face
01899   int nb1d = 0;
01900   TopTools_MapOfShape Edges1;
01901   for (TopExp_Explorer exp(aFaces.Value(NumBase), TopAbs_EDGE); exp.More(); exp.Next()) {
01902     Edges1.Add(exp.Current());
01903     SMESH_subMesh *sm = aMesh.GetSubMesh(exp.Current());
01904     if( sm ) {
01905       MapShapeNbElemsItr anIt = aResMap.find(sm);
01906       if( anIt == aResMap.end() ) continue;
01907       std::vector<int> aVec = (*anIt).second;
01908       nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
01909     }
01910   }
01911   // find face opposite to base face
01912   int OppNum = 0;
01913   for(i=1; i<=6; i++) {
01914     if(i==NumBase) continue;
01915     bool IsOpposite = true;
01916     for(TopExp_Explorer exp(aFaces.Value(i), TopAbs_EDGE); exp.More(); exp.Next()) {
01917       if( Edges1.Contains(exp.Current()) ) {
01918         IsOpposite = false;
01919         break;
01920       }
01921     }
01922     if(IsOpposite) {
01923       OppNum = i;
01924       break;
01925     }
01926   }
01927   // find number of 2d elems on side faces
01928   int nb2d = 0;
01929   for(i=1; i<=6; i++) {
01930     if( i==OppNum || i==NumBase ) continue;
01931     MapShapeNbElemsItr anIt = aResMap.find( meshFaces[i-1] );
01932     if( anIt == aResMap.end() ) continue;
01933     std::vector<int> aVec = (*anIt).second;
01934     nb2d += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
01935   }
01936 
01937   MapShapeNbElemsItr anIt = aResMap.find( meshFaces[NumBase-1] );
01938   std::vector<int> aVec = (*anIt).second;
01939   int nb2d_face0 = Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
01940   int nb0d_face0 = aVec[SMDSEntity_Node];
01941 
01942   anIt = aResMap.find( meshFaces[OppNum-1] );
01943   for(i=SMDSEntity_Node; i<SMDSEntity_Last; i++)
01944     (*anIt).second[i] = aVec[i];
01945 
01946   SMESH_MesherHelper aTool (aMesh);
01947   bool _quadraticMesh = aTool.IsQuadraticSubMesh(aShape);
01948 
01949   std::vector<int> aResVec(SMDSEntity_Last);
01950   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
01951   if(_quadraticMesh) {
01952     aResVec[SMDSEntity_Quad_Penta] = nb2d_face0 * ( nb2d/nb1d );
01953     aResVec[SMDSEntity_Node] = nb0d_face0 * ( 2*nb2d/nb1d - 1 );
01954   }
01955   else {
01956     aResVec[SMDSEntity_Node] = nb0d_face0 * ( nb2d/nb1d - 1 );
01957     aResVec[SMDSEntity_Penta] = nb2d_face0 * ( nb2d/nb1d );
01958   }
01959   SMESH_subMesh * sm = aMesh.GetSubMesh(aShape);
01960   aResMap.insert(std::make_pair(sm,aResVec));
01961 
01962   return true;
01963 }
01964