Back to index

salome-smesh  6.5.0
StdMeshers_CompositeHexa_3D.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 //  SMESH SMESH : implementaion of SMESH idl descriptions
00021 // File      : StdMeshers_CompositeHexa_3D.cxx
00022 // Module    : SMESH
00023 // Created   : Tue Nov 25 11:04:59 2008
00024 // Author    : Edward AGAPOV (eap)
00025 
00026 #include "StdMeshers_CompositeHexa_3D.hxx"
00027 
00028 #include "SMDS_Mesh.hxx"
00029 #include "SMDS_MeshNode.hxx"
00030 #include "SMDS_SetIterator.hxx"
00031 #include "SMESH_Block.hxx"
00032 #include "SMESH_Comment.hxx"
00033 #include "SMESH_ComputeError.hxx"
00034 #include "SMESH_Mesh.hxx"
00035 #include "SMESH_MesherHelper.hxx"
00036 #include "SMESH_subMesh.hxx"
00037 
00038 #include <BRepAdaptor_Surface.hxx>
00039 #include <BRep_Tool.hxx>
00040 #include <Standard_ErrorHandler.hxx>
00041 #include <Standard_Failure.hxx>
00042 #include <TopExp_Explorer.hxx>
00043 #include <TopTools_MapIteratorOfMapOfShape.hxx>
00044 #include <TopTools_MapOfShape.hxx>
00045 #include <TopTools_SequenceOfShape.hxx>
00046 #include <TopoDS.hxx>
00047 #include <TopoDS_Edge.hxx>
00048 #include <TopoDS_Face.hxx>
00049 #include <TopoDS_Vertex.hxx>
00050 #include <gp_Pnt.hxx>
00051 #include <gp_Pnt2d.hxx>
00052 #include <gp_Vec.hxx>
00053 #include <gp_XYZ.hxx>
00054 
00055 #include <list>
00056 #include <set>
00057 #include <vector>
00058 
00059 
00060 #ifdef _DEBUG_
00061 
00062 // #define DEB_FACES
00063 // #define DEB_GRID
00064 #define DUMP_VERT(msg,V) \
00065 // { TopoDS_Vertex v = V; gp_Pnt p = BRep_Tool::Pnt(v);\
00066 //   cout << msg << "( "<< p.X()<<", "<<p.Y()<<", "<<p.Z()<<" )"<<endl;}
00067 
00068 #else
00069 
00070 #define DUMP_VERT(msg,v)
00071 
00072 #endif
00073 
00074 //================================================================================
00075 // text for message about an internal error
00076 #define ERR_LI(txt) SMESH_Comment(txt) << ":" << __LINE__
00077 
00078 // order corresponds to right order of edges in CASCADE face
00079 enum EQuadSides{ Q_BOTTOM=0, Q_RIGHT, Q_TOP, Q_LEFT,   Q_CHILD, Q_PARENT };
00080 
00081 enum EBoxSides{ B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_UNDEFINED };
00082 
00083 //================================================================================
00087 struct _Indexer
00088 {
00089   int _xSize, _ySize;
00090   _Indexer( int xSize, int ySize ): _xSize(xSize), _ySize(ySize) {}
00091   int size() const { return _xSize * _ySize; }
00092   int operator()(const int x, const int y) const { return y * _xSize + x; }
00093 };
00094 
00095 //================================================================================
00099 class _FaceSide
00100 {
00101 public:
00102   _FaceSide(const _FaceSide& other);
00103   _FaceSide(const TopoDS_Edge& edge=TopoDS_Edge());
00104   _FaceSide(const list<TopoDS_Edge>& edges);
00105   _FaceSide* GetSide(const int i);
00106   const _FaceSide* GetSide(const int i) const;
00107   int size() { return myChildren.size(); }
00108   int NbVertices() const;
00109   TopoDS_Vertex FirstVertex() const;
00110   TopoDS_Vertex LastVertex() const;
00111   TopoDS_Vertex Vertex(int i) const;
00112   bool Contain( const _FaceSide& side, int* which=0 ) const;
00113   bool Contain( const TopoDS_Vertex& vertex ) const;
00114   void AppendSide( const _FaceSide& side );
00115   void SetBottomSide( int i );
00116   int GetNbSegments(SMESH_Mesh& mesh) const;
00117   bool StoreNodes(SMESH_Mesh& mesh, vector<const SMDS_MeshNode*>& myGrid, bool reverse );
00118   void SetID(EQuadSides id) { myID = id; }
00119   static inline const TopoDS_TShape* ptr(const TopoDS_Shape& theShape)
00120   { return theShape.TShape().operator->(); }
00121   void Dump() const;
00122 
00123 private:
00124 
00125 
00126   TopoDS_Edge       myEdge;
00127   list< _FaceSide > myChildren;
00128   int               myNbChildren;
00129 
00130   //set<const TopoDS_TShape*> myVertices;
00131   TopTools_MapOfShape myVertices;
00132 
00133   EQuadSides        myID; // debug
00134 };
00135 //================================================================================
00140 class _QuadFaceGrid
00141 {
00142   typedef list< _QuadFaceGrid > TChildren;
00143 public:
00144   _QuadFaceGrid();
00145 
00146 public: //** Methods to find and orient faces of 6 sides of the box **//
00147   
00149   bool Init(const TopoDS_Face& f);
00150 
00152   bool AddContinuousFace( const _QuadFaceGrid& f );
00153 
00155   bool SetBottomSide(const _FaceSide& side, int* sideIndex=0);
00156 
00158   _QuadFaceGrid* FindAdjacentForSide(int i, vector<_QuadFaceGrid>& faces) const; // (0<i<4)
00159 
00161   void ReverseEdges(/*int e1, int e2*/);
00162 
00163   bool IsComplex() const { return !myChildren.empty(); }
00164 
00165   typedef SMDS_SetIterator< const _QuadFaceGrid&, TChildren::const_iterator > TChildIterator;
00166 
00167   TChildIterator GetChildren() const
00168   { return TChildIterator( myChildren.begin(), myChildren.end()); }
00169 
00170 public: //** Loading and access to mesh **//
00171 
00173   bool LoadGrid( SMESH_Mesh& mesh );
00174 
00176   int GetNbHoriSegments(SMESH_Mesh& mesh, bool withBrothers=false) const;
00177 
00179   int GetNbVertSegments(SMESH_Mesh& mesh, bool withBrothers=false) const;
00180 
00182   const SMDS_MeshNode* GetNode(int iHori, int iVert) const;
00183 
00185   gp_XYZ GetXYZ(int iHori, int iVert) const;
00186 
00187 public: //** Access to member fields **//
00188 
00190   const _FaceSide& GetSide(int i) const;
00191 
00193   TopoDS_Face GetFace() const { return myFace; }
00194 
00196   bool GetNormal( const TopoDS_Vertex& v, gp_Vec& n ) const;
00197 
00198   SMESH_ComputeErrorPtr GetError() const { return myError; }
00199 
00200   void SetID(EBoxSides id) { myID = id; }
00201 
00202   void DumpGrid() const;
00203 
00204   void DumpVertices() const;
00205 
00206 private:
00207 
00208   bool error(const std::string& text, int code = COMPERR_ALGO_FAILED)
00209   { myError = SMESH_ComputeError::New( code, text ); return false; }
00210 
00211   bool error(const SMESH_ComputeErrorPtr& err)
00212   { myError = err; return ( !myError || myError->IsOK() ); }
00213 
00214   bool loadCompositeGrid(SMESH_Mesh& mesh);
00215 
00216   bool fillGrid(SMESH_Mesh&                    theMesh,
00217                 vector<const SMDS_MeshNode*> & theGrid,
00218                 const _Indexer&                theIndexer,
00219                 int                            theX,
00220                 int                            theY);
00221 
00222   bool locateChildren();
00223 
00224   void setBrothers( set< _QuadFaceGrid* >& notLocatedBrothers );
00225 
00226   TopoDS_Face myFace;
00227   _FaceSide   mySides;
00228   bool        myReverse;
00229 
00230   TChildren   myChildren;
00231 
00232   _QuadFaceGrid* myLeftBottomChild;
00233   _QuadFaceGrid* myRightBrother;
00234   _QuadFaceGrid* myUpBrother;
00235 
00236   _Indexer    myIndexer;
00237   vector<const SMDS_MeshNode*>  myGrid;
00238 
00239   SMESH_ComputeErrorPtr         myError;
00240 
00241   EBoxSides   myID; // debug
00242 };
00243 
00244 //================================================================================
00248 //================================================================================
00249 
00250 StdMeshers_CompositeHexa_3D::StdMeshers_CompositeHexa_3D(int hypId, int studyId, SMESH_Gen* gen)
00251   :SMESH_3D_Algo(hypId, studyId, gen)
00252 {
00253   _name = "CompositeHexa_3D";
00254   _shapeType = (1 << TopAbs_SHELL) | (1 << TopAbs_SOLID);       // 1 bit /shape type
00255 }
00256 
00257 //================================================================================
00261 //================================================================================
00262 
00263 bool StdMeshers_CompositeHexa_3D::CheckHypothesis(SMESH_Mesh&         aMesh,
00264                                                   const TopoDS_Shape& aShape,
00265                                                   Hypothesis_Status&  aStatus)
00266 {
00267   aStatus = HYP_OK;
00268   return true;
00269 }
00270 
00271 //================================================================================
00278 //================================================================================
00279 
00280 bool StdMeshers_CompositeHexa_3D::Compute(SMESH_Mesh&         theMesh,
00281                                           const TopoDS_Shape& theShape)
00282 {
00283   SMESH_MesherHelper helper( theMesh );
00284   _quadraticMesh = helper.IsQuadraticSubMesh( theShape );
00285   helper.SetElementsOnShape( true );
00286 
00287   // -------------------------
00288   // Try to find 6 side faces
00289   // -------------------------
00290   vector< _QuadFaceGrid > boxFaces; boxFaces.reserve( 6 );
00291   TopExp_Explorer exp;
00292   int iFace, nbFaces = 0;
00293   for ( exp.Init(theShape, TopAbs_FACE); exp.More(); exp.Next(), ++nbFaces )
00294   {
00295     _QuadFaceGrid f;
00296     if ( !f.Init( TopoDS::Face( exp.Current() )))
00297       return error (COMPERR_BAD_SHAPE);
00298     bool isContinuous = false;
00299     for ( int i=0; i < boxFaces.size() && !isContinuous; ++i )
00300       isContinuous = boxFaces[ i ].AddContinuousFace( f );
00301     if ( !isContinuous )
00302       boxFaces.push_back( f );
00303   }
00304   // Check what we have
00305   if ( boxFaces.size() != 6 && nbFaces != 6)
00306     return error
00307       (COMPERR_BAD_SHAPE,
00308        SMESH_Comment("Can't find 6 sides of a box. Number of found sides - ")<<boxFaces.size());
00309 
00310   if ( boxFaces.size() != 6 && nbFaces == 6 ) { // strange ordinary box with continuous faces
00311     boxFaces.resize( 6 );
00312     iFace = 0;
00313     for ( exp.Init(theShape, TopAbs_FACE); exp.More(); exp.Next(), ++iFace )
00314       boxFaces[ iFace ].Init( TopoDS::Face( exp.Current() ) );
00315   }
00316   // ----------------------------------------
00317   // Find out position of faces within a box
00318   // ----------------------------------------
00319 
00320   _QuadFaceGrid *fBottom, *fTop, *fFront, *fBack, *fLeft, *fRight;
00321   // start from a bottom face
00322   fBottom = &boxFaces[0];
00323   // find vertical faces
00324   fFront = fBottom->FindAdjacentForSide( Q_BOTTOM, boxFaces );
00325   fLeft  = fBottom->FindAdjacentForSide( Q_RIGHT, boxFaces );
00326   fBack  = fBottom->FindAdjacentForSide( Q_TOP, boxFaces );
00327   fRight = fBottom->FindAdjacentForSide( Q_LEFT, boxFaces );
00328   // check the found
00329   if ( !fFront || !fBack || !fLeft || !fRight )
00330     return error(COMPERR_BAD_SHAPE);
00331   // top face
00332   fTop = 0;
00333   for ( int i=1; i < boxFaces.size() && !fTop; ++i ) {
00334     fTop = & boxFaces[ i ];
00335     if ( fTop==fFront || fTop==fLeft || fTop==fBack || fTop==fRight )
00336       fTop = 0;
00337   }
00338   // set bottom of the top side
00339   if ( !fTop->SetBottomSide( fFront->GetSide( Q_TOP ) )) {
00340     if ( !fFront->IsComplex() )
00341       return error( ERR_LI("Error in StdMeshers_CompositeHexa_3D::Compute()"));
00342     else {
00343       _QuadFaceGrid::TChildIterator chIt = fFront->GetChildren();
00344       while ( chIt.more() ) {
00345         const _QuadFaceGrid& frontChild = chIt.next();
00346         if ( fTop->SetBottomSide( frontChild.GetSide( Q_TOP )))
00347           break;
00348       }
00349     }
00350   }
00351   if ( !fTop )
00352     return error(COMPERR_BAD_SHAPE);
00353 
00354   fBottom->SetID( B_BOTTOM );
00355   fBack  ->SetID( B_BACK );
00356   fLeft  ->SetID( B_LEFT );
00357   fFront ->SetID( B_FRONT );
00358   fRight ->SetID( B_RIGHT );
00359   fTop   ->SetID( B_TOP );
00360 
00361   // orient bottom egde of faces along axes of the unit box
00362   fBottom->ReverseEdges();
00363   fBack  ->ReverseEdges();
00364   fLeft  ->ReverseEdges();
00365 
00366   // ------------------------------------------
00367   // Fill columns of nodes with existing nodes
00368   // ------------------------------------------
00369 
00370   // let faces load their grids
00371   if ( !fBottom->LoadGrid( theMesh )) return error( fBottom->GetError() );
00372   if ( !fBack  ->LoadGrid( theMesh )) return error( fBack  ->GetError() );
00373   if ( !fLeft  ->LoadGrid( theMesh )) return error( fLeft  ->GetError() );
00374   if ( !fFront ->LoadGrid( theMesh )) return error( fFront ->GetError() );
00375   if ( !fRight ->LoadGrid( theMesh )) return error( fRight ->GetError() );
00376   if ( !fTop   ->LoadGrid( theMesh )) return error( fTop   ->GetError() );
00377 
00378   int x, xSize = fBottom->GetNbHoriSegments(theMesh) + 1, X = xSize - 1;
00379   int y, ySize = fBottom->GetNbVertSegments(theMesh) + 1, Y = ySize - 1;
00380   int z, zSize = fFront ->GetNbVertSegments(theMesh) + 1, Z = zSize - 1;
00381   _Indexer colIndex( xSize, ySize );
00382   vector< vector< const SMDS_MeshNode* > > columns( colIndex.size() );
00383 
00384   // fill node columns by front and back box sides
00385   for ( x = 0; x < xSize; ++x ) {
00386     vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
00387     vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
00388     column0.resize( zSize );
00389     column1.resize( zSize );
00390     for ( z = 0; z < zSize; ++z ) {
00391       column0[ z ] = fFront->GetNode( x, z );
00392       column1[ z ] = fBack ->GetNode( x, z );
00393     }
00394   }
00395   // fill node columns by left and right box sides
00396   for ( y = 1; y < ySize-1; ++y ) {
00397     vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
00398     vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
00399     column0.resize( zSize );
00400     column1.resize( zSize );
00401     for ( z = 0; z < zSize; ++z ) {
00402       column0[ z ] = fLeft ->GetNode( y, z );
00403       column1[ z ] = fRight->GetNode( y, z );
00404     }
00405   }
00406   // get nodes from top and bottom box sides
00407   for ( x = 1; x < xSize-1; ++x ) {
00408     for ( y = 1; y < ySize-1; ++y ) {
00409       vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
00410       column.resize( zSize );
00411       column.front() = fBottom->GetNode( x, y );
00412       column.back()  = fTop   ->GetNode( x, y );
00413     }
00414   }
00415 
00416   // ----------------------------
00417   // Add internal nodes of a box
00418   // ----------------------------
00419   // projection points of internal nodes on box sub-shapes by which
00420   // coordinates of internal nodes are computed
00421   vector<gp_XYZ> pointsOnShapes( SMESH_Block::ID_Shell );
00422 
00423   // projections on vertices are constant
00424   pointsOnShapes[ SMESH_Block::ID_V000 ] = fBottom->GetXYZ( 0, 0 );
00425   pointsOnShapes[ SMESH_Block::ID_V100 ] = fBottom->GetXYZ( X, 0 );
00426   pointsOnShapes[ SMESH_Block::ID_V010 ] = fBottom->GetXYZ( 0, Y );
00427   pointsOnShapes[ SMESH_Block::ID_V110 ] = fBottom->GetXYZ( X, Y );
00428   pointsOnShapes[ SMESH_Block::ID_V001 ] = fTop->GetXYZ( 0, 0 );
00429   pointsOnShapes[ SMESH_Block::ID_V101 ] = fTop->GetXYZ( X, 0 );
00430   pointsOnShapes[ SMESH_Block::ID_V011 ] = fTop->GetXYZ( 0, Y );
00431   pointsOnShapes[ SMESH_Block::ID_V111 ] = fTop->GetXYZ( X, Y );
00432 
00433   for ( x = 1; x < xSize-1; ++x )
00434   {
00435     gp_XYZ params; // normalized parameters of internal node within a unit box
00436     params.SetCoord( 1, x / double(X) );
00437     for ( y = 1; y < ySize-1; ++y )
00438     {
00439       params.SetCoord( 2, y / double(Y) );
00440       // column to fill during z loop
00441       vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
00442       // points projections on horizontal edges
00443       pointsOnShapes[ SMESH_Block::ID_Ex00 ] = fBottom->GetXYZ( x, 0 );
00444       pointsOnShapes[ SMESH_Block::ID_Ex10 ] = fBottom->GetXYZ( x, Y );
00445       pointsOnShapes[ SMESH_Block::ID_E0y0 ] = fBottom->GetXYZ( 0, y );
00446       pointsOnShapes[ SMESH_Block::ID_E1y0 ] = fBottom->GetXYZ( X, y );
00447       pointsOnShapes[ SMESH_Block::ID_Ex01 ] = fTop->GetXYZ( x, 0 );
00448       pointsOnShapes[ SMESH_Block::ID_Ex11 ] = fTop->GetXYZ( x, Y );
00449       pointsOnShapes[ SMESH_Block::ID_E0y1 ] = fTop->GetXYZ( 0, y );
00450       pointsOnShapes[ SMESH_Block::ID_E1y1 ] = fTop->GetXYZ( X, y );
00451       // points projections on horizontal faces
00452       pointsOnShapes[ SMESH_Block::ID_Fxy0 ] = fBottom->GetXYZ( x, y );
00453       pointsOnShapes[ SMESH_Block::ID_Fxy1 ] = fTop   ->GetXYZ( x, y );
00454       for ( z = 1; z < zSize-1; ++z ) // z loop
00455       {
00456         params.SetCoord( 3, z / double(Z) );
00457         // point projections on vertical edges
00458         pointsOnShapes[ SMESH_Block::ID_E00z ] = fFront->GetXYZ( 0, z );    
00459         pointsOnShapes[ SMESH_Block::ID_E10z ] = fFront->GetXYZ( X, z );    
00460         pointsOnShapes[ SMESH_Block::ID_E01z ] = fBack->GetXYZ( 0, z );    
00461         pointsOnShapes[ SMESH_Block::ID_E11z ] = fBack->GetXYZ( X, z );
00462         // point projections on vertical faces
00463         pointsOnShapes[ SMESH_Block::ID_Fx0z ] = fFront->GetXYZ( x, z );    
00464         pointsOnShapes[ SMESH_Block::ID_Fx1z ] = fBack ->GetXYZ( x, z );    
00465         pointsOnShapes[ SMESH_Block::ID_F0yz ] = fLeft ->GetXYZ( y, z );    
00466         pointsOnShapes[ SMESH_Block::ID_F1yz ] = fRight->GetXYZ( y, z );
00467 
00468         // compute internal node coordinates
00469         gp_XYZ coords;
00470         SMESH_Block::ShellPoint( params, pointsOnShapes, coords );
00471         column[ z ] = helper.AddNode( coords.X(), coords.Y(), coords.Z() );
00472 
00473 #ifdef DEB_GRID
00474         // debug
00475         //cout << "----------------------------------------------------------------------"<<endl;
00476         //for ( int id = SMESH_Block::ID_V000; id < SMESH_Block::ID_Shell; ++id)
00477         //{
00478         //  gp_XYZ p = pointsOnShapes[ id ];
00479         //  SMESH_Block::DumpShapeID( id,cout)<<" ( "<<p.X()<<", "<<p.Y()<<", "<<p.Z()<<" )"<<endl;
00480         //}
00481         //cout << "Params: ( "<< params.X()<<", "<<params.Y()<<", "<<params.Z()<<" )"<<endl;
00482         //cout << "coords: ( "<< coords.X()<<", "<<coords.Y()<<", "<<coords.Z()<<" )"<<endl;
00483 #endif
00484       }
00485     }
00486   }
00487   // faces no more needed, free memory
00488   boxFaces.clear();
00489 
00490   // ----------------
00491   // Add hexahedrons
00492   // ----------------
00493   for ( x = 0; x < xSize-1; ++x ) {
00494     for ( y = 0; y < ySize-1; ++y ) {
00495       vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
00496       vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
00497       vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
00498       vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
00499       for ( z = 0; z < zSize-1; ++z )
00500       {
00501         // bottom face normal of a hexa mush point outside the volume
00502         helper.AddVolume(col00[z],   col01[z],   col11[z],   col10[z],
00503                          col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
00504       }
00505     }
00506   }
00507   return true;
00508 }
00509 
00510 
00511 //=======================================================================
00512 //function : GetNb2d
00513 //purpose  : auxilary for Evaluate
00514 //=======================================================================
00515 int GetNb2d(_QuadFaceGrid* QFG, SMESH_Mesh& theMesh,
00516             MapShapeNbElems& aResMap)
00517 {
00518   int nb2d = 0;
00519   _QuadFaceGrid::TChildIterator aCI = QFG->GetChildren();
00520   while( aCI.more() ) {
00521     const _QuadFaceGrid& currChild = aCI.next();
00522     SMESH_subMesh *sm = theMesh.GetSubMesh(currChild.GetFace());
00523     if( sm ) {
00524       MapShapeNbElemsItr anIt = aResMap.find(sm);
00525       if( anIt == aResMap.end() ) continue;
00526       std::vector<int> aVec = (*anIt).second;
00527       nb2d += Max(aVec[SMDSEntity_Quadrangle],aVec[SMDSEntity_Quad_Quadrangle]);
00528     }
00529   }
00530   return nb2d;
00531 }
00532 
00533 
00534 //================================================================================
00538 //================================================================================
00539 
00540 bool StdMeshers_CompositeHexa_3D::Evaluate(SMESH_Mesh& theMesh,
00541                                            const TopoDS_Shape& theShape,
00542                                            MapShapeNbElems& aResMap)
00543 {
00544   SMESH_MesherHelper aTool(theMesh);
00545   bool _quadraticMesh = aTool.IsQuadraticSubMesh(theShape);
00546 
00547 
00548   // -------------------------
00549   // Try to find 6 side faces
00550   // -------------------------
00551   vector< _QuadFaceGrid > boxFaces; boxFaces.reserve( 6 );
00552   TopExp_Explorer exp;
00553   int iFace, nbFaces = 0;
00554   for ( exp.Init(theShape, TopAbs_FACE); exp.More(); exp.Next(), ++nbFaces )
00555   {
00556     _QuadFaceGrid f;
00557     if ( !f.Init( TopoDS::Face( exp.Current() )))
00558       //return error (COMPERR_BAD_SHAPE);
00559       return false;
00560     bool isContinuous = false;
00561     for ( int i=0; i < boxFaces.size() && !isContinuous; ++i )
00562       isContinuous = boxFaces[ i ].AddContinuousFace( f );
00563     if ( !isContinuous )
00564       boxFaces.push_back( f );
00565   }
00566   // Check what we have
00567   if ( boxFaces.size() != 6 && nbFaces != 6)
00568     //return error
00569     //  (COMPERR_BAD_SHAPE,
00570     //   SMESH_Comment("Can't find 6 sides of a box. Number of found sides - ")<<boxFaces.size());
00571     return false;
00572 
00573   if ( boxFaces.size() != 6 && nbFaces == 6 ) { // strange ordinary box with continuous faces
00574     boxFaces.resize( 6 );
00575     iFace = 0;
00576     for ( exp.Init(theShape, TopAbs_FACE); exp.More(); exp.Next(), ++iFace )
00577       boxFaces[ iFace ].Init( TopoDS::Face( exp.Current() ) );
00578   }
00579 
00580   // ----------------------------------------
00581   // Find out position of faces within a box
00582   // ----------------------------------------
00583 
00584   _QuadFaceGrid *fBottom, *fTop, *fFront, *fBack, *fLeft, *fRight;
00585   // start from a bottom face
00586   fBottom = &boxFaces[0];
00587   // find vertical faces
00588   fFront = fBottom->FindAdjacentForSide( Q_BOTTOM, boxFaces );
00589   fLeft  = fBottom->FindAdjacentForSide( Q_RIGHT, boxFaces );
00590   fBack  = fBottom->FindAdjacentForSide( Q_TOP, boxFaces );
00591   fRight = fBottom->FindAdjacentForSide( Q_LEFT, boxFaces );
00592   // check the found
00593   if ( !fFront || !fBack || !fLeft || !fRight )
00594     //return error(COMPERR_BAD_SHAPE);
00595     return false;
00596   // top face
00597   fTop = 0;
00598   int i = 1;
00599   for(; i < boxFaces.size() && !fTop; ++i ) {
00600     fTop = & boxFaces[ i ];
00601     if ( fTop==fFront || fTop==fLeft || fTop==fBack || fTop==fRight )
00602       fTop = 0;
00603   }
00604   // set bottom of the top side
00605   if ( !fTop->SetBottomSide( fFront->GetSide( Q_TOP ) )) {
00606     if ( !fFront->IsComplex() )
00607       //return error( ERR_LI("Error in StdMeshers_CompositeHexa_3D::Compute()"));
00608       return false;
00609     else {
00610       _QuadFaceGrid::TChildIterator chIt = fFront->GetChildren();
00611       while ( chIt.more() ) {
00612         const _QuadFaceGrid& frontChild = chIt.next();
00613         if ( fTop->SetBottomSide( frontChild.GetSide( Q_TOP )))
00614           break;
00615       }
00616     }
00617   }
00618   if ( !fTop )
00619     //return error(COMPERR_BAD_SHAPE);
00620     return false;
00621 
00622 
00623   TopTools_SequenceOfShape BottomFaces;
00624   _QuadFaceGrid::TChildIterator aCI = fBottom->GetChildren();
00625   while( aCI.more() ) {
00626     const _QuadFaceGrid& currChild = aCI.next();
00627     BottomFaces.Append(currChild.GetFace());
00628   }
00629   // find boundary edges and internal nodes for bottom face
00630   TopTools_SequenceOfShape BndEdges;
00631   int nb0d_in = 0;
00632   //TopTools_MapOfShape BndEdges;
00633   for(i=1; i<=BottomFaces.Length(); i++) {
00634     for (TopExp_Explorer exp(BottomFaces.Value(i), TopAbs_EDGE); exp.More(); exp.Next()) {
00635       int nb0 = 0;
00636       SMESH_subMesh *sm = theMesh.GetSubMesh(exp.Current());
00637       if( sm ) {
00638         MapShapeNbElemsItr anIt = aResMap.find(sm);
00639         if( anIt == aResMap.end() ) continue;
00640         std::vector<int> aVec = (*anIt).second;
00641         nb0 = aVec[SMDSEntity_Node];
00642       }
00643       int j = 1;
00644       for(; j<=BndEdges.Length(); j++) {
00645         if( BndEdges.Value(j) == exp.Current() ) {
00646           // internal edge => remove it
00647           BndEdges.Remove(j);
00648           nb0d_in += nb0;
00649           break;
00650         }
00651       }
00652       if( j > BndEdges.Length() ) {
00653         BndEdges.Append(exp.Current());
00654       }
00655       //if( BndEdges.Contains(exp.Current()) ) {
00656       //BndEdges.Remove( exp.Current() );
00657       //}
00658       //else {
00659       //BndEdges.Add( exp.Current() );
00660       //}
00661     }
00662   }
00663 
00664   // find number of 1d elems for bottom face
00665   int nb1d = 0;
00666   for(i=1; i<=BndEdges.Length(); i++) {
00667     SMESH_subMesh *sm = theMesh.GetSubMesh(BndEdges.Value(i));
00668     if( sm ) {
00669       MapShapeNbElemsItr anIt = aResMap.find(sm);
00670       if( anIt == aResMap.end() ) continue;
00671       std::vector<int> aVec = (*anIt).second;
00672       nb1d += Max(aVec[SMDSEntity_Edge],aVec[SMDSEntity_Quad_Edge]);
00673     }
00674   }
00675 
00676   // find number of 2d elems on side faces
00677   int nb2d = 0;
00678   nb2d += GetNb2d(fFront, theMesh, aResMap);
00679   nb2d += GetNb2d(fRight, theMesh, aResMap);
00680   nb2d += GetNb2d(fBack, theMesh, aResMap);
00681   nb2d += GetNb2d(fLeft, theMesh, aResMap);
00682 
00683   // find number of 2d elems and nodes on bottom faces
00684   int nb0d=0, nb2d_3=0, nb2d_4=0;
00685   for(i=1; i<=BottomFaces.Length(); i++) {
00686     SMESH_subMesh *sm = theMesh.GetSubMesh(BottomFaces.Value(i));
00687     if( sm ) {
00688       MapShapeNbElemsItr anIt = aResMap.find(sm);
00689       if( anIt == aResMap.end() ) continue;
00690       std::vector<int> aVec = (*anIt).second;
00691       nb0d += aVec[SMDSEntity_Node];
00692       nb2d_3 += Max(aVec[SMDSEntity_Triangle],   aVec[SMDSEntity_Quad_Triangle]);
00693       nb2d_4 += Max(aVec[SMDSEntity_Quadrangle], aVec[SMDSEntity_Quad_Quadrangle]);
00694     }
00695   }
00696   nb0d += nb0d_in;
00697 
00698   std::vector<int> aResVec(SMDSEntity_Last);
00699   for(int i=SMDSEntity_Node; i<SMDSEntity_Last; i++) aResVec[i] = 0;
00700   if(_quadraticMesh) {
00701     aResVec[SMDSEntity_Quad_Penta] = nb2d_3 * ( nb2d/nb1d );
00702     aResVec[SMDSEntity_Quad_Hexa]  = nb2d_4 * ( nb2d/nb1d );
00703     aResVec[SMDSEntity_Node] = nb0d * ( 2*nb2d/nb1d - 1 );
00704   }
00705   else {
00706     aResVec[SMDSEntity_Node]  = nb0d * ( nb2d/nb1d - 1 );
00707     aResVec[SMDSEntity_Penta] = nb2d_3 * ( nb2d/nb1d );
00708     aResVec[SMDSEntity_Hexa]  = nb2d_4 * ( nb2d/nb1d );
00709   }
00710   SMESH_subMesh * sm = theMesh.GetSubMesh(theShape);
00711   aResMap.insert(std::make_pair(sm,aResVec));
00712 
00713   return true;
00714 }
00715 
00716 
00717 //================================================================================
00721 //================================================================================
00722 
00723 _QuadFaceGrid::_QuadFaceGrid():
00724   myReverse(false), myRightBrother(0), myUpBrother(0), myIndexer(0,0), myID(B_UNDEFINED)
00725 {
00726 }
00727 
00728 //================================================================================
00732 //================================================================================
00733 
00734 bool _QuadFaceGrid::Init(const TopoDS_Face& f)
00735 {
00736   myFace         = f;
00737   mySides        = _FaceSide();
00738   myReverse      = false;
00739   myLeftBottomChild = myRightBrother = myUpBrother = 0;
00740   myChildren.clear();
00741   myGrid.clear();
00742   //if ( myFace.Orientation() != TopAbs_FORWARD )
00743     //myFace.Reverse();
00744 
00745   TopoDS_Vertex V;
00746   list< TopoDS_Edge > edges;
00747   list< int > nbEdgesInWire;
00748   int nbWire = SMESH_Block::GetOrderedEdges (myFace, V, edges, nbEdgesInWire);
00749   if ( nbWire != 1 )
00750     return false;
00751 
00752   list< TopoDS_Edge >::iterator edgeIt = edges.begin();
00753   if ( nbEdgesInWire.front() == 4 ) // exactly 4 edges
00754   {
00755     for ( ; edgeIt != edges.end(); ++edgeIt )
00756       mySides.AppendSide( _FaceSide( *edgeIt ));
00757   }
00758   else if ( nbEdgesInWire.front() > 4 ) { // more than 4 edges - try to unite some
00759     list< TopoDS_Edge > sideEdges;
00760     while ( !edges.empty()) {
00761       sideEdges.clear();
00762       sideEdges.splice( sideEdges.end(), edges, edges.begin());// edges.front()->sideEdges.back()
00763       while ( !edges.empty() ) {
00764         if ( SMESH_Algo::IsContinuous( sideEdges.back(), edges.front() )) {
00765           sideEdges.splice( sideEdges.end(), edges, edges.begin());
00766         }
00767         else if ( SMESH_Algo::IsContinuous( sideEdges.front(), edges.back() )) {
00768           sideEdges.splice( sideEdges.begin(), edges, --edges.end());
00769         }
00770         else {
00771           break;
00772         }
00773       }
00774       mySides.AppendSide( _FaceSide( sideEdges ));
00775     }
00776   }
00777   if (mySides.size() != 4)
00778     return false;
00779 
00780 #ifdef _DEBUG_
00781   mySides.GetSide( Q_BOTTOM )->SetID( Q_BOTTOM );
00782   mySides.GetSide( Q_RIGHT  )->SetID( Q_RIGHT );
00783   mySides.GetSide( Q_TOP    )->SetID( Q_TOP );
00784   mySides.GetSide( Q_LEFT   )->SetID( Q_LEFT );
00785 #endif
00786 
00787   return true;
00788 }
00789 
00790 //================================================================================
00794 //================================================================================
00795 
00796 bool _QuadFaceGrid::AddContinuousFace( const _QuadFaceGrid& other )
00797 {
00798   for ( int i = 0; i < 4; ++i ) {
00799     const _FaceSide& otherSide = other.GetSide( i );
00800     int iMyCommon;
00801     if ( mySides.Contain( otherSide, &iMyCommon ) ) {
00802       // check if normals of two faces are collinear at all vertices of a otherSide
00803       const double angleTol = M_PI / 180. / 2.;
00804       int iV, nbV = otherSide.NbVertices(), nbCollinear = 0;
00805       for ( iV = 0; iV < nbV; ++iV )
00806       {
00807         TopoDS_Vertex v = otherSide.Vertex( iV );
00808         gp_Vec n1, n2;
00809         if ( !GetNormal( v, n1 ) || !other.GetNormal( v, n2 ))
00810           continue;
00811         if ( n1 * n2 < 0 )
00812           n1.Reverse();
00813         if ( n1.Angle(n2) < angleTol )
00814           nbCollinear++;
00815         else
00816           break;
00817       }
00818       if ( nbCollinear > 1 ) { // this face becomes composite if not yet is
00819         DUMP_VERT("Cont 1", mySides.GetSide(iMyCommon)->FirstVertex());
00820         DUMP_VERT("Cont 2", mySides.GetSide(iMyCommon)->LastVertex());
00821         DUMP_VERT("Cont 3", otherSide.FirstVertex());
00822         DUMP_VERT("Cont 4", otherSide.LastVertex());
00823         if ( myChildren.empty() ) {
00824           myChildren.push_back( *this );
00825           myFace.Nullify();
00826         }
00827         myChildren.push_back( other );
00828         int otherBottomIndex = ( 4 + i - iMyCommon + 2 ) % 4;
00829         myChildren.back().SetBottomSide( other.GetSide( otherBottomIndex ));
00830         // collect vertices in mySides
00831         mySides.AppendSide( other.GetSide(0) );
00832         mySides.AppendSide( other.GetSide(1) );
00833         mySides.AppendSide( other.GetSide(2) );
00834         mySides.AppendSide( other.GetSide(3) );
00835         return true;
00836       }
00837     }
00838   }
00839   return false;
00840 }
00841 
00842 //================================================================================
00846 //================================================================================
00847 
00848 bool _QuadFaceGrid::SetBottomSide(const _FaceSide& bottom, int* sideIndex)
00849 {
00850   myLeftBottomChild = myRightBrother = myUpBrother = 0;
00851 
00852   int myBottomIndex;
00853   if ( myChildren.empty() )
00854   {
00855     if ( mySides.Contain( bottom, &myBottomIndex )) {
00856       mySides.SetBottomSide( myBottomIndex );
00857       if ( sideIndex )
00858         *sideIndex = myBottomIndex;
00859       return true;
00860     }
00861   }
00862   else
00863   {
00864     TChildren::iterator childFace = myChildren.begin(), childEnd = myChildren.end();
00865     for ( ; childFace != childEnd; ++childFace )
00866     {
00867       if ( childFace->SetBottomSide( bottom, &myBottomIndex ))
00868       {
00869         TChildren::iterator orientedCild = childFace;
00870         for ( childFace = myChildren.begin(); childFace != childEnd; ++childFace ) {
00871           if ( childFace != orientedCild )
00872             childFace->SetBottomSide( childFace->GetSide( myBottomIndex ));
00873         }
00874         if ( sideIndex )
00875           *sideIndex = myBottomIndex;
00876         return true;
00877       }
00878     }
00879   }
00880   return false;
00881 }
00882 
00883 //================================================================================
00887 //================================================================================
00888 
00889 _QuadFaceGrid* _QuadFaceGrid::FindAdjacentForSide(int i, vector<_QuadFaceGrid>& faces) const
00890 {
00891   for ( int iF = 0; iF < faces.size(); ++iF ) {
00892     _QuadFaceGrid* f  = &faces[ iF ];
00893     if ( f != this && f->SetBottomSide( GetSide( i )))
00894       return f;
00895   }
00896   return (_QuadFaceGrid*) 0;
00897 }
00898 
00899 //================================================================================
00903 //================================================================================
00904 
00905 const _FaceSide& _QuadFaceGrid::GetSide(int i) const
00906 {
00907   if ( myChildren.empty() )
00908     return *mySides.GetSide(i);
00909 
00910   _QuadFaceGrid* me = const_cast<_QuadFaceGrid*>(this);
00911   if ( !me->locateChildren() || !myLeftBottomChild )
00912     return *mySides.GetSide(i);
00913 
00914   const _QuadFaceGrid* child = myLeftBottomChild;
00915   switch ( i ){
00916   case Q_BOTTOM:
00917   case Q_LEFT:
00918     break;
00919   case Q_RIGHT:
00920     while ( child->myRightBrother )
00921       child = child->myRightBrother;
00922     break;
00923   case Q_TOP:
00924     while ( child->myUpBrother )
00925       child = child->myUpBrother;
00926     break;
00927   default: ;
00928   }
00929   return child->GetSide( i );
00930 }
00931 
00932 //================================================================================
00936 //================================================================================
00937 
00938 void _QuadFaceGrid::ReverseEdges(/*int e1, int e2*/)
00939 {
00940   myReverse = !myReverse;
00941 
00942 // #ifdef DEB_FACES
00943 //   if ( !myFace.IsNull() )
00944 //     TopAbs::Print(myFace.Orientation(), cout);
00945 // #endif
00946 
00947   if ( myChildren.empty() )
00948   {
00949 //     mySides.GetSide( e1 )->Reverse();
00950 //     mySides.GetSide( e2 )->Reverse();
00951     DumpVertices();
00952   }
00953   else
00954   {
00955     DumpVertices();
00956     TChildren::iterator child = myChildren.begin(), childEnd = myChildren.end();
00957     for ( ; child != childEnd; ++child )
00958       child->ReverseEdges( /*e1, e2*/ );
00959   }
00960 }
00961 
00962 //================================================================================
00966 //================================================================================
00967 
00968 bool _QuadFaceGrid::LoadGrid( SMESH_Mesh& mesh )
00969 {
00970   if ( !myChildren.empty() )
00971   {
00972     // Let child faces load their grids
00973     TChildren::iterator child = myChildren.begin(), childEnd = myChildren.end();
00974     for ( ; child != childEnd; ++child ) {
00975       child->SetID( myID );
00976       if ( !child->LoadGrid( mesh ) )
00977         return error( child->GetError() );
00978     }
00979     // Fill myGrid with nodes of patches
00980     return loadCompositeGrid( mesh );
00981   }
00982 
00983   // ---------------------------------------
00984   // Fill myGrid with nodes bound to myFace
00985   // ---------------------------------------
00986 
00987   if ( !myGrid.empty() )
00988     return true;
00989 
00990   SMESHDS_SubMesh* faceSubMesh = mesh.GetSubMesh( myFace )->GetSubMeshDS();
00991   // check that all faces are quadrangular
00992   SMDS_ElemIteratorPtr fIt = faceSubMesh->GetElements();
00993   while ( fIt->more() )
00994     if ( fIt->next()->NbNodes() % 4 > 0 )
00995       return error("Non-quadrangular mesh faces are not allowed on sides of a composite block");
00996   
00997   myIndexer._xSize = 1 + mySides.GetSide( Q_BOTTOM )->GetNbSegments( mesh );
00998   myIndexer._ySize = 1 + mySides.GetSide( Q_LEFT   )->GetNbSegments( mesh );
00999 
01000   myGrid.resize( myIndexer.size() );
01001 
01002   // strore nodes bound to the bottom edge
01003   mySides.GetSide( Q_BOTTOM )->StoreNodes( mesh, myGrid, myReverse );
01004 
01005   // store the rest nodes row by row
01006 
01007   const SMDS_MeshNode* dummy = mesh.GetMeshDS()->AddNode(0,0,0);
01008   const SMDS_MeshElement* firstQuad = dummy; // most left face above the last row of found nodes
01009   
01010   int nbFoundNodes = myIndexer._xSize;
01011   while ( nbFoundNodes != myGrid.size() )
01012   {
01013     // first and last nodes of the last filled row of nodes
01014     const SMDS_MeshNode* n1down = myGrid[ nbFoundNodes - myIndexer._xSize ];
01015     const SMDS_MeshNode* n2down = myGrid[ nbFoundNodes - myIndexer._xSize + 1];
01016     const SMDS_MeshNode* n1downLast = myGrid[ nbFoundNodes-1 ];
01017 
01018     // find the first face above the row by the first two left nodes
01019     //
01020     // n1up     n2up
01021     //     o---o
01022     //     |   |
01023     //     o---o  o  o  o  o
01024     //n1down    n2down
01025     //
01026     TIDSortedElemSet emptySet, avoidSet;
01027     avoidSet.insert( firstQuad );
01028     firstQuad = SMESH_MeshEditor::FindFaceInSet( n1down, n2down, emptySet, avoidSet);
01029     while ( firstQuad && !faceSubMesh->Contains( firstQuad )) {
01030       avoidSet.insert( firstQuad );
01031       firstQuad = SMESH_MeshEditor::FindFaceInSet( n1down, n2down, emptySet, avoidSet);
01032     }
01033     if ( !firstQuad || !faceSubMesh->Contains( firstQuad ))
01034       return error(ERR_LI("Error in _QuadFaceGrid::LoadGrid()"));
01035 
01036     // find the node of quad bound to the left geom edge
01037     int i2down = firstQuad->GetNodeIndex( n2down );
01038     const SMDS_MeshNode* n1up = firstQuad->GetNode(( i2down+2 ) % 4 );
01039     myGrid[ nbFoundNodes++ ] = n1up;
01040     // the 4-the node of the first quad
01041     int i1down = firstQuad->GetNodeIndex( n1down );
01042     const SMDS_MeshNode* n2up = firstQuad->GetNode(( i1down+2 ) % 4 );
01043     myGrid[ nbFoundNodes++ ] = n2up;
01044 
01045     n1down = n2down;
01046     n1up   = n2up;
01047     const SMDS_MeshElement* quad = firstQuad;
01048 
01049     // find the rest nodes by remaining faces above the row
01050     //
01051     //             n1up
01052     //     o---o--o
01053     //     |   |  | ->
01054     //     o---o--o  o  o  o
01055     //                      n1downLast
01056     //
01057     while ( n1down != n1downLast )
01058     {
01059       // next face
01060       avoidSet.clear(); avoidSet.insert( quad );
01061       quad = SMESH_MeshEditor::FindFaceInSet( n1down, n1up, emptySet, avoidSet );
01062       if ( !quad || quad->NbNodes() % 4 > 0)
01063         return error(ERR_LI("Error in _QuadFaceGrid::LoadGrid()"));
01064 
01065       // next node
01066       if ( quad->GetNode( i1down ) != n1down ) // check already found index
01067         i1down = quad->GetNodeIndex( n1down );
01068       n2up = quad->GetNode(( i1down+2 ) % 4 );
01069       myGrid[ nbFoundNodes++ ] = n2up;
01070 
01071       n1down = myGrid[ nbFoundNodes - myIndexer._xSize - 1 ];
01072       n1up   = n2up;
01073     }
01074   }
01075   mesh.GetMeshDS()->RemoveNode(dummy);
01076   DumpGrid(); // debug
01077 
01078   return true;
01079 }
01080 
01081 //================================================================================
01085 //================================================================================
01086 
01087 bool _QuadFaceGrid::locateChildren()
01088 {
01089   if ( myLeftBottomChild )
01090     return true;
01091 
01092   TChildren::iterator child = myChildren.begin(), childEnd = myChildren.end();
01093 
01094   // find a child sharing it's first bottom vertex with no other brother
01095   myLeftBottomChild = 0;
01096   for ( ; !myLeftBottomChild && child != childEnd; ++child )
01097   {
01098     TopoDS_Vertex leftVertex = child->GetSide( Q_BOTTOM ).FirstVertex();
01099     bool sharedVertex = false;
01100     TChildren::iterator otherChild = myChildren.begin();
01101     for ( ; otherChild != childEnd && !sharedVertex; ++otherChild )
01102       if ( otherChild != child )
01103         sharedVertex = otherChild->mySides.Contain( leftVertex );
01104     if ( !sharedVertex ) {
01105       myLeftBottomChild = & (*child);
01106       DUMP_VERT("0 left bottom Vertex: ",leftVertex );
01107     }
01108   }
01109   if (!myLeftBottomChild)
01110     return error(ERR_LI("Error in locateChildren()"));
01111 
01112   set< _QuadFaceGrid* > notLocatedChilren;
01113   for (child = myChildren.begin() ; child != childEnd; ++child )
01114     notLocatedChilren.insert( & (*child));
01115 
01116   // connect myLeftBottomChild to it's right and upper brothers
01117   notLocatedChilren.erase( myLeftBottomChild );
01118   myLeftBottomChild->setBrothers( notLocatedChilren );
01119   if ( !notLocatedChilren.empty() )
01120     return error(ERR_LI("Error in locateChildren()"));
01121 
01122   return true;
01123 }
01124 
01125 //================================================================================
01129 //================================================================================
01130 
01131 bool _QuadFaceGrid::loadCompositeGrid(SMESH_Mesh& mesh)
01132 {
01133   // Find out mutual location of children: find their right and up brothers
01134   if ( !locateChildren() )
01135     return false;
01136 
01137   // Load nodes according to mutual location of children
01138 
01139   // grid size
01140   myIndexer._xSize = 1 + myLeftBottomChild->GetNbHoriSegments(mesh, /*withBrothers=*/true);
01141   myIndexer._ySize = 1 + myLeftBottomChild->GetNbVertSegments(mesh, /*withBrothers=*/true);
01142 
01143   myGrid.resize( myIndexer.size() );
01144 
01145   int fromX = myReverse ? myIndexer._xSize : 0;
01146   if (!myLeftBottomChild->fillGrid( mesh, myGrid, myIndexer, fromX, 0 ))
01147     return error( myLeftBottomChild->GetError() );
01148 
01149   DumpGrid();
01150 
01151   return true;
01152 }
01153 
01154 //================================================================================
01158 //================================================================================
01159 
01160 void _QuadFaceGrid::setBrothers( set< _QuadFaceGrid* >& notLocatedBrothers )
01161 {
01162   if ( !notLocatedBrothers.empty() )
01163   {
01164     // find right brother
01165     TopoDS_Vertex rightVertex = GetSide( Q_BOTTOM ).LastVertex();
01166     DUMP_VERT("1 right bottom Vertex: ",rightVertex );
01167     set< _QuadFaceGrid* >::iterator brIt, brEnd = notLocatedBrothers.end();
01168     for ( brIt = notLocatedBrothers.begin(); !myRightBrother && brIt != brEnd; ++brIt )
01169     {
01170       _QuadFaceGrid* brother = *brIt;
01171       TopoDS_Vertex brotherLeftVertex = brother->GetSide( Q_BOTTOM ).FirstVertex();
01172       DUMP_VERT( "brother left bottom: ", brotherLeftVertex );
01173       if ( rightVertex.IsSame( brotherLeftVertex )) {
01174         myRightBrother = brother;
01175         notLocatedBrothers.erase( myRightBrother );
01176       }
01177     }
01178     // find upper brother
01179     TopoDS_Vertex upVertex = GetSide( Q_LEFT ).FirstVertex();
01180     DUMP_VERT("1 left up Vertex: ",upVertex);
01181     brIt = notLocatedBrothers.begin(), brEnd = notLocatedBrothers.end();
01182     for ( ; !myUpBrother && brIt != brEnd; ++brIt )
01183     {
01184       _QuadFaceGrid* brother = *brIt;
01185       TopoDS_Vertex brotherLeftVertex = brother->GetSide( Q_BOTTOM ).FirstVertex();
01186       DUMP_VERT("brother left bottom: ", brotherLeftVertex);
01187       if ( upVertex.IsSame( brotherLeftVertex )) {
01188         myUpBrother = brother;
01189         notLocatedBrothers.erase( myUpBrother );
01190       }
01191     }
01192     // recursive call
01193     if ( myRightBrother )
01194       myRightBrother->setBrothers( notLocatedBrothers );
01195     if ( myUpBrother )
01196       myUpBrother->setBrothers( notLocatedBrothers );
01197   }
01198 }
01199 
01200 //================================================================================
01204 //================================================================================
01205 
01206 bool _QuadFaceGrid::fillGrid(SMESH_Mesh&                    theMesh,
01207                              vector<const SMDS_MeshNode*> & theGrid,
01208                              const _Indexer&                theIndexer,
01209                              int                            theX,
01210                              int                            theY)
01211 {
01212   if ( myGrid.empty() && !LoadGrid( theMesh ))
01213     return false;
01214 
01215   // store my own grid in the global grid
01216 
01217   int fromX = myReverse ? theX - myIndexer._xSize: theX;
01218 
01219   for ( int i = 0, x = fromX; i < myIndexer._xSize; ++i, ++x )
01220     for ( int j = 0, y = theY; j < myIndexer._ySize; ++j, ++y )
01221       theGrid[ theIndexer( x, y )] = myGrid[ myIndexer( i, j )];
01222 
01223   // store grids of my right and upper brothers
01224 
01225   if ( myRightBrother )
01226   {
01227     if ( myReverse )
01228       fromX += 1;
01229     else
01230       fromX += myIndexer._xSize - 1;
01231     if ( !myRightBrother->fillGrid( theMesh, theGrid, theIndexer, fromX, theY ))
01232       return error( myRightBrother->GetError() );
01233   }
01234   if ( myUpBrother )
01235   {
01236     if ( !myUpBrother->fillGrid( theMesh, theGrid, theIndexer,
01237                                  theX, theY + myIndexer._ySize - 1))
01238       return error( myUpBrother->GetError() );
01239   }
01240   return true;
01241 }
01242 
01243 //================================================================================
01247 //================================================================================
01248 
01249 int _QuadFaceGrid::GetNbHoriSegments(SMESH_Mesh& mesh, bool withBrothers) const
01250 {
01251   int nbSegs = 0;
01252   if ( myLeftBottomChild )
01253   {
01254     nbSegs += myLeftBottomChild->GetNbHoriSegments( mesh, true );
01255   }
01256   else
01257   {
01258     nbSegs = mySides.GetSide( Q_BOTTOM )->GetNbSegments(mesh);
01259     if ( withBrothers && myRightBrother )
01260       nbSegs += myRightBrother->GetNbHoriSegments( mesh, withBrothers );
01261   }
01262   return nbSegs;
01263 }
01264 
01265 //================================================================================
01269 //================================================================================
01270 
01271 int _QuadFaceGrid::GetNbVertSegments(SMESH_Mesh& mesh, bool withBrothers) const
01272 {
01273   int nbSegs = 0;
01274   if ( myLeftBottomChild )
01275   {
01276     nbSegs += myLeftBottomChild->GetNbVertSegments( mesh, true );
01277   }
01278   else
01279   {
01280     nbSegs = mySides.GetSide( Q_LEFT )->GetNbSegments(mesh);
01281     if ( withBrothers && myUpBrother )
01282       nbSegs += myUpBrother->GetNbVertSegments( mesh, withBrothers );
01283   }
01284   return nbSegs;
01285 }
01286 
01287 //================================================================================
01291 //================================================================================
01292 
01293 const SMDS_MeshNode* _QuadFaceGrid::GetNode(int iHori, int iVert) const
01294 {
01295   return myGrid[ myIndexer( iHori, iVert )];
01296 }
01297 
01298 //================================================================================
01302 //================================================================================
01303 
01304 gp_XYZ _QuadFaceGrid::GetXYZ(int iHori, int iVert) const
01305 {
01306   const SMDS_MeshNode* n = myGrid[ myIndexer( iHori, iVert )];
01307   return gp_XYZ( n->X(), n->Y(), n->Z() );
01308 }
01309 
01310 //================================================================================
01314 //================================================================================
01315 
01316 bool _QuadFaceGrid::GetNormal( const TopoDS_Vertex& v, gp_Vec& n ) const
01317 {
01318   if ( myChildren.empty() )
01319   {
01320     if ( mySides.Contain( v )) {
01321       try {
01322         gp_Pnt2d uv = BRep_Tool::Parameters( v, myFace );
01323         BRepAdaptor_Surface surface( myFace );
01324         gp_Pnt p; gp_Vec d1u, d1v;
01325         surface.D1( uv.X(), uv.Y(), p, d1u, d1v );
01326         n = d1u.Crossed( d1v );
01327         return true;
01328       }
01329       catch (Standard_Failure) {
01330         return false;
01331       }
01332     }
01333   }
01334   else
01335   {
01336     TChildren::const_iterator child = myChildren.begin(), childEnd = myChildren.end();
01337     for ( ; child != childEnd; ++child )
01338       if ( child->GetNormal( v, n ))
01339         return true;
01340   }
01341   return false;
01342 }
01343 
01344 //================================================================================
01348 //================================================================================
01349 
01350 void _QuadFaceGrid::DumpGrid() const
01351 {
01352 #ifdef DEB_GRID
01353   const char* names[] = { "B_BOTTOM", "B_RIGHT", "B_TOP", "B_LEFT", "B_FRONT", "B_BACK" };
01354   cout << "****** Face " << names[ myID ] << endl;
01355 
01356   if ( myChildren.empty() || !myGrid.empty() )
01357   {
01358     cout << "x size: " << myIndexer._xSize << "; y size: " << myIndexer._ySize << endl;
01359     for ( int y = 0; y < myIndexer._ySize; ++y ) {
01360       cout << "-- row " << y << endl;
01361       for ( int x = 0; x < myIndexer._xSize; ++x ) {
01362         const SMDS_MeshNode* n = myGrid[ myIndexer( x, y ) ];
01363         cout << x << " ( " << n->X() << ", " << n->Y() << ", " << n->Z() << " )" << endl;
01364       }
01365     }
01366   }
01367   else
01368   {
01369     cout << "Nb children: " << myChildren.size() << endl;
01370     TChildren::const_iterator child = myChildren.begin(), childEnd = myChildren.end();
01371     for ( int i=0; child != childEnd; ++child, ++i ) {
01372       cout << "   *** SUBFACE " << i+1 << endl;
01373       ((_QuadFaceGrid&)(*child)).SetID( myID );
01374       child->DumpGrid();
01375     }
01376   }
01377 #endif
01378 }
01379 
01380 //================================================================================
01384 //================================================================================
01385 
01386 void _QuadFaceGrid::DumpVertices() const
01387 {
01388 #ifdef DEB_FACES
01389   cout << "****** Face ";
01390   const char* names[] = { "B_BOTTOM", "B_RIGHT", "B_TOP", "B_LEFT", "B_FRONT", "B_BACK" };
01391   if ( myID >= B_BOTTOM && myID < B_BACK )
01392     cout << names[ myID ] << endl;
01393   else
01394     cout << "UNDEFINED" << endl;
01395 
01396   if ( myChildren.empty() )
01397   {
01398     for ( int i = 0; i < 4; ++i )
01399     {
01400       cout << "  Side "; mySides.GetSide( i )->Dump();
01401     }
01402   }
01403   else
01404   {
01405     cout << "-- Nb children: " << myChildren.size() << endl;
01406     TChildren::const_iterator child = myChildren.begin(), childEnd = myChildren.end();
01407     for ( int i=0; child != childEnd; ++child, ++i ) {
01408       cout << "   *** SUBFACE " << i+1 << endl;
01409       ((_QuadFaceGrid&)(*child)).SetID( myID );
01410       child->DumpVertices();
01411     }
01412   }
01413 #endif
01414 }
01415 
01416 //=======================================================================
01417 //function : _FaceSide
01418 //purpose  : copy constructor
01419 //=======================================================================
01420 
01421 _FaceSide::_FaceSide(const _FaceSide& other)
01422 {
01423   myEdge = other.myEdge;
01424   myChildren = other.myChildren;
01425   myNbChildren = other.myNbChildren;
01426   myVertices.Assign( other.myVertices );
01427   myID = other.myID;
01428 }
01429 
01430 //================================================================================
01434 //================================================================================
01435 
01436 _FaceSide::_FaceSide(const TopoDS_Edge& edge):
01437   myEdge( edge ), myNbChildren(0)
01438 {
01439   if ( !edge.IsNull() )
01440     for ( TopExp_Explorer exp( edge, TopAbs_VERTEX ); exp.More(); exp.Next() )
01441       //myVertices.insert( ptr ( exp.Current() ));
01442       myVertices.Add( exp.Current() );
01443 }
01444 
01445 //================================================================================
01449 //================================================================================
01450 
01451 _FaceSide::_FaceSide(const list<TopoDS_Edge>& edges):
01452   myNbChildren(0)
01453 {
01454   list<TopoDS_Edge>::const_iterator edge = edges.begin(), eEnd = edges.end();
01455   for ( ; edge != eEnd; ++edge ) {
01456     myChildren.push_back( _FaceSide( *edge ));
01457     myNbChildren++;
01458 //     myVertices.insert( myChildren.back().myVertices.begin(),
01459 //                        myChildren.back().myVertices.end() );
01460     myVertices.Add( myChildren.back().FirstVertex() );
01461     myVertices.Add( myChildren.back().LastVertex() );
01462     myChildren.back().SetID( Q_CHILD ); // not to splice them
01463   }
01464 }
01465 
01466 //=======================================================================
01467 //function : GetSide
01468 //purpose  : 
01469 //=======================================================================
01470 
01471 _FaceSide* _FaceSide::GetSide(const int i)
01472 {
01473   if ( i >= myNbChildren )
01474     return 0;
01475 
01476   list< _FaceSide >::iterator side = myChildren.begin();
01477   if ( i )
01478     std::advance( side, i );
01479   return & (*side);
01480 }
01481 
01482 //=======================================================================
01483 //function : GetSide
01484 //purpose  : 
01485 //=======================================================================
01486 
01487 const _FaceSide* _FaceSide::GetSide(const int i) const
01488 {
01489   return const_cast< _FaceSide* >(this)->GetSide(i);
01490 }
01491 
01492 //=======================================================================
01493 //function : NbVertices
01494 //purpose  : return nb of vertices in the side
01495 //=======================================================================
01496 
01497 int _FaceSide::NbVertices() const
01498 {
01499   if ( myChildren.empty() )
01500     return myVertices.Extent();
01501 //     return myVertices.size();
01502 
01503   return myNbChildren + 1;
01504 }
01505 
01506 //=======================================================================
01507 //function : FirstVertex
01508 //purpose  : 
01509 //=======================================================================
01510 
01511 TopoDS_Vertex _FaceSide::FirstVertex() const
01512 {
01513   if ( myChildren.empty() )
01514     return TopExp::FirstVertex( myEdge, Standard_True );
01515 
01516   return myChildren.front().FirstVertex();
01517 }
01518 
01519 //=======================================================================
01520 //function : LastVertex
01521 //purpose  : 
01522 //=======================================================================
01523 
01524 TopoDS_Vertex _FaceSide::LastVertex() const
01525 {
01526   if ( myChildren.empty() )
01527     return TopExp::LastVertex( myEdge, Standard_True );
01528 
01529   return myChildren.back().LastVertex();
01530 }
01531 
01532 //=======================================================================
01533 //function : Vertex
01534 //purpose  : 
01535 //=======================================================================
01536 
01537 TopoDS_Vertex _FaceSide::Vertex(int i) const
01538 {
01539   if ( myChildren.empty() )
01540     return i ? LastVertex() : FirstVertex();
01541       
01542   if ( i >= myNbChildren )
01543     return myChildren.back().LastVertex();
01544   
01545   return GetSide(i)->FirstVertex();
01546 }
01547 
01548 //=======================================================================
01549 //function : Contain
01550 //purpose  : 
01551 //=======================================================================
01552 
01553 bool _FaceSide::Contain( const _FaceSide& side, int* which ) const
01554 {
01555   if ( !which || myChildren.empty() )
01556   {
01557     if ( which )
01558       *which = 0;
01559     int nbCommon = 0;
01560 //     set<const TopoDS_TShape*>::iterator v, vEnd = side.myVertices.end();
01561 //     for ( v = side.myVertices.begin(); v != vEnd; ++v )
01562 //       nbCommon += ( myVertices.find( *v ) != myVertices.end() );
01563     TopTools_MapIteratorOfMapOfShape vIt ( side.myVertices );
01564     for ( ; vIt.More(); vIt.Next() )
01565       nbCommon += ( myVertices.Contains( vIt.Key() ));
01566     return (nbCommon > 1);
01567   }
01568   list< _FaceSide >::const_iterator mySide = myChildren.begin(), sideEnd = myChildren.end();
01569   for ( int i = 0; mySide != sideEnd; ++mySide, ++i ) {
01570     if ( mySide->Contain( side )) {
01571       *which = i;
01572       return true;
01573     }
01574   }
01575   return false;
01576 }
01577 
01578 //=======================================================================
01579 //function : Contain
01580 //purpose  : 
01581 //=======================================================================
01582 
01583 bool _FaceSide::Contain( const TopoDS_Vertex& vertex ) const
01584 {
01585   return myVertices.Contains( vertex );
01586 //   return myVertices.find( ptr( vertex )) != myVertices.end();
01587 }
01588 
01589 //=======================================================================
01590 //function : AppendSide
01591 //purpose  : 
01592 //=======================================================================
01593 
01594 void _FaceSide::AppendSide( const _FaceSide& side )
01595 {
01596   if ( !myEdge.IsNull() )
01597   {
01598     myChildren.push_back( *this );
01599     myNbChildren = 1;
01600     myEdge.Nullify();
01601   }
01602   myChildren.push_back( side );
01603   myNbChildren++;
01604   //myVertices.insert( side.myVertices.begin(), side.myVertices.end() );
01605   TopTools_MapIteratorOfMapOfShape vIt ( side.myVertices );
01606   for ( ; vIt.More(); vIt.Next() )
01607     myVertices.Add( vIt.Key() );
01608 
01609   myID = Q_PARENT;
01610   myChildren.back().SetID( EQuadSides( myNbChildren-1 ));
01611 }
01612 
01613 //=======================================================================
01614 //function : SetBottomSide
01615 //purpose  : 
01616 //=======================================================================
01617 
01618 void _FaceSide::SetBottomSide( int i )
01619 {
01620   if ( i > 0 && myID == Q_PARENT ) {
01621     list< _FaceSide >::iterator sideEnd, side = myChildren.begin();
01622     std::advance( side, i );
01623     myChildren.splice( myChildren.begin(), myChildren, side, myChildren.end() );
01624 
01625     side = myChildren.begin(), sideEnd = myChildren.end();
01626     for ( int i = 0; side != sideEnd; ++side, ++i ) {
01627       side->SetID( EQuadSides(i) );
01628       side->SetBottomSide(i);
01629     }
01630   }
01631 }
01632 
01633 //=======================================================================
01634 //function : GetNbSegments
01635 //purpose  : 
01636 //=======================================================================
01637 
01638 int _FaceSide::GetNbSegments(SMESH_Mesh& mesh) const
01639 {
01640   int nb = 0;
01641   if ( myChildren.empty() )
01642   {
01643     nb = mesh.GetSubMesh(myEdge)->GetSubMeshDS()->NbElements();
01644   }
01645   else
01646   {
01647     list< _FaceSide >::const_iterator side = myChildren.begin(), sideEnd = myChildren.end();
01648     for ( ; side != sideEnd; ++side )
01649       nb += side->GetNbSegments(mesh);
01650   }
01651   return nb;
01652 }
01653 
01654 //=======================================================================
01655 //function : StoreNodes
01656 //purpose  : 
01657 //=======================================================================
01658 
01659 bool _FaceSide::StoreNodes(SMESH_Mesh&                   mesh,
01660                            vector<const SMDS_MeshNode*>& myGrid,
01661                            bool                          reverse )
01662 {
01663   list< TopoDS_Edge > edges;
01664   if ( myChildren.empty() )
01665   {
01666     edges.push_back( myEdge );
01667   }
01668   else
01669   {
01670     list< _FaceSide >::const_iterator side = myChildren.begin(), sideEnd = myChildren.end();
01671     for ( ; side != sideEnd; ++side )
01672       if ( reverse )
01673         edges.push_front( side->myEdge );
01674       else
01675         edges.push_back ( side->myEdge );
01676   }
01677   int nbNodes = 0;
01678   list< TopoDS_Edge >::iterator edge = edges.begin(), eEnd = edges.end();
01679   for ( ; edge != eEnd; ++edge )
01680   {
01681     map< double, const SMDS_MeshNode* > nodes;
01682     bool ok = SMESH_Algo::GetSortedNodesOnEdge( mesh.GetMeshDS(),
01683                                                 *edge,
01684                                                 /*ignoreMediumNodes=*/true,
01685                                                 nodes);
01686     if ( !ok ) return false;
01687 
01688     bool forward = ( edge->Orientation() == TopAbs_FORWARD );
01689     if ( reverse ) forward = !forward;
01690     if ( forward )
01691     {
01692       map< double, const SMDS_MeshNode* >::iterator u_node, nEnd = nodes.end();
01693       for ( u_node = nodes.begin(); u_node != nEnd; ++u_node )
01694         myGrid[ nbNodes++ ] = u_node->second;
01695     }
01696     else 
01697     {
01698       map< double, const SMDS_MeshNode* >::reverse_iterator u_node, nEnd = nodes.rend();
01699       for ( u_node = nodes.rbegin(); u_node != nEnd; ++u_node )
01700         myGrid[ nbNodes++ ] = u_node->second;
01701     }
01702     nbNodes--; // node on vertex present in two adjacent edges
01703   }
01704   return nbNodes > 0;
01705 }
01706 
01707 //=======================================================================
01708 //function : Dump
01709 //purpose  : dump end vertices
01710 //=======================================================================
01711 
01712 void _FaceSide::Dump() const
01713 {
01714   if ( myChildren.empty() )
01715   {
01716     const char* sideNames[] = { "Q_BOTTOM", "Q_RIGHT", "Q_TOP", "Q_LEFT", "Q_CHILD", "Q_PARENT" };
01717     if ( myID >= Q_BOTTOM && myID < Q_PARENT )
01718       cout << sideNames[ myID ] << endl;
01719     else
01720       cout << "<UNDEFINED ID>" << endl;
01721     TopoDS_Vertex f = FirstVertex();
01722     TopoDS_Vertex l = LastVertex();
01723     gp_Pnt pf = BRep_Tool::Pnt(f);
01724     gp_Pnt pl = BRep_Tool::Pnt(l);
01725     cout << "\t ( "<< ptr( f ) << " - " << ptr( l )<< " )"
01726          << "\t ( "<< pf.X()<<", "<<pf.Y()<<", "<<pf.Z()<<" ) - "
01727          << " ( "<< pl.X()<<", "<<pl.Y()<<", "<<pl.Z()<<" )"<<endl;
01728   }
01729   else
01730   {
01731     list< _FaceSide >::const_iterator side = myChildren.begin(), sideEnd = myChildren.end();
01732     for ( ; side != sideEnd; ++side ) {
01733       side->Dump();
01734       cout << "\t";
01735     }
01736   }
01737 }