Back to index

salome-smesh  6.5.0
StdMeshers_HexaFromSkin_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 // File      : StdMeshers_HexaFromSkin_3D.cxx
00021 // Created   : Wed Jan 27 12:28:07 2010
00022 // Author    : Edward AGAPOV (eap)
00023 
00024 #include "StdMeshers_HexaFromSkin_3D.hxx"
00025 
00026 #include "SMDS_VolumeOfNodes.hxx"
00027 #include "SMDS_VolumeTool.hxx"
00028 #include "SMESH_Block.hxx"
00029 #include "SMESH_MesherHelper.hxx"
00030 
00031 #include <gp_Ax2.hxx>
00032 
00033 //#include "utilities.h"
00034 #include <limits>
00035 
00036 // Define error message and _MYDEBUG_ if needed
00037 #ifdef _DEBUG_
00038 #define BAD_MESH_ERR \
00039   error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh.\n" \
00040                       __FILE__ ":" )<<__LINE__)
00041 //#define _MYDEBUG_
00042 #else
00043 #define BAD_MESH_ERR \
00044   error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh"))
00045 #endif
00046 
00047 
00048 // Debug output
00049 #ifdef _MYDEBUG_
00050 #define _DUMP_(msg) cout << msg << endl
00051 #else
00052 #define _DUMP_(msg)
00053 #endif
00054 
00055 
00056 namespace
00057 {
00058   enum EBoxSides 
00059     {
00060       B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, NB_BLOCK_SIDES
00061     };
00062 #ifdef _MYDEBUG_
00063   const char* SBoxSides[] = 
00064     {
00065       "BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED"
00066     };
00067 #endif
00068   enum EQuadEdge 
00069     {
00070       Q_BOTTOM = 0, Q_RIGHT, Q_TOP, Q_LEFT, NB_QUAD_SIDES
00071     };
00072 
00073 
00074   //================================================================================
00078   //================================================================================
00079 
00080   bool getEdgeEnds(EQuadEdge edge, bool& xMax1, bool& yMax1, bool& xMax2, bool& yMax2 )
00081   {
00082     xMax1=0, yMax1=0, xMax2=1, yMax2=1;
00083     switch( edge )
00084     {
00085     case Q_BOTTOM: yMax2 = 0; break;
00086     case Q_RIGHT:  xMax1 = 1; break;
00087     case Q_TOP:    yMax1 = 1; break;
00088     case Q_LEFT:   xMax2 = 0; break;
00089     default:
00090       return false;
00091     }
00092     return true;
00093   }
00094 
00095   //================================================================================
00101   //================================================================================
00102 
00103   bool isCornerNode( const SMDS_MeshNode* n )
00104   {
00105     int nbF = n ? n->NbInverseElements( SMDSAbs_Face ) : 1;
00106     if ( nbF % 2 )
00107       return true;
00108 
00109     set<const SMDS_MeshNode*> nodesInInverseFaces;
00110     SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face );
00111     while ( fIt->more() )
00112     {
00113       const SMDS_MeshElement* face = fIt->next();
00114       nodesInInverseFaces.insert( face->begin_nodes(), face->end_nodes() );
00115     }
00116 
00117     return nodesInInverseFaces.size() != ( 6 + (nbF/2-1)*3 );
00118   }
00119 
00120   //================================================================================
00124   //================================================================================
00125 
00126   bool isQuadrangle(const SMDS_MeshElement* e)
00127   {
00128     return ( e && e->NbCornerNodes() == 4 );
00129   }
00130 
00131   //================================================================================
00135   //================================================================================
00136 
00137   const SMDS_MeshNode* oppositeNode(const SMDS_MeshElement* quad, int iNode)
00138   {
00139     return quad->GetNode( (iNode+2) % 4 );
00140   }
00141 
00142   //================================================================================
00146   struct _Indexer
00147   {
00148     int _xSize, _ySize;
00149     _Indexer( int xSize=0, int ySize=0 ): _xSize(xSize), _ySize(ySize) {}
00150     int size() const { return _xSize * _ySize; }
00151     int operator()(int x, int y) const { return y * _xSize + x; }
00152   };
00153   //================================================================================
00157   class _OrientedIndexer : public _Indexer
00158   {
00159   public:
00160     enum OriFlags 
00161       {
00162         REV_X = 1, REV_Y = 2, SWAP_XY = 4, MAX_ORI = REV_X|REV_Y|SWAP_XY
00163       };
00164     _OrientedIndexer( const _Indexer& indexer, const int oriFlags ):
00165       _Indexer( indexer._xSize, indexer._ySize ),
00166       _xSize (indexer._xSize), _ySize(indexer._ySize),
00167       _xRevFun((oriFlags & REV_X) ? & reverse : & lazy),
00168       _yRevFun((oriFlags & REV_Y) ? & reverse : & lazy),
00169       _swapFun((oriFlags & SWAP_XY ) ? & swap : & lazy)
00170     {
00171       (*_swapFun)( _xSize, _ySize );
00172     }
00174     int operator()(int x, int y) const
00175     {
00176       (*_xRevFun)( x, const_cast<int&>( _xSize ));
00177       (*_yRevFun)( y, const_cast<int&>( _ySize ));
00178       (*_swapFun)( x, y );
00179       return _Indexer::operator()(x,y);
00180     }
00182     int corner(bool xMax, bool yMax) const
00183     {
00184       int x = xMax, y = yMax, size = 2;
00185       (*_xRevFun)( x, size );
00186       (*_yRevFun)( y, size );
00187       (*_swapFun)( x, y );
00188       return _Indexer::operator()(x ? _Indexer::_xSize-1 : 0 , y ? _Indexer::_ySize-1 : 0);
00189     }
00190     int xSize() const { return _xSize; }
00191     int ySize() const { return _ySize; }
00192   private:
00193     _Indexer _indexer;
00194     int _xSize, _ySize;
00195 
00196     typedef void (*TFun)(int& x, int& y);
00197     TFun _xRevFun, _yRevFun, _swapFun;
00198     
00199     static void lazy   (int&, int&) {}
00200     static void reverse(int& x, int& size) { x = size - x - 1; }
00201     static void swap   (int& x, int& y) { std::swap(x,y); }
00202   };
00203   //================================================================================
00207   struct _BlockSide
00208   {
00209     vector<const SMDS_MeshNode*> _grid;
00210     _Indexer                     _index;
00211     int                          _nbBlocksExpected;
00212     int                          _nbBlocksFound;
00213 
00214 #ifdef _DEBUG_ // want to get SIGSEGV in case of invalid index
00215 #define _grid_access_(pobj, i) pobj->_grid[ ((i) < pobj->_grid.size()) ? i : int(1e100)]
00216 #else
00217 #define _grid_access_(pobj, i) pobj->_grid[ i ]
00218 #endif
00219 
00220     const SMDS_MeshNode* getNode(int x, int y) const { return _grid_access_(this, _index( x,y ));}
00222     void setNode(int x, int y, const SMDS_MeshNode* n) { _grid_access_(this, _index( x,y )) = n; }
00224     SMESH_OrientedLink getEdge(EQuadEdge edge) const
00225     {
00226       bool x1, y1, x2, y2; getEdgeEnds( edge, x1, y1, x2, y2 );
00227       return SMESH_OrientedLink( getCornerNode ( x1, y1 ), getCornerNode( x2, y2 ));
00228     }
00230     const SMDS_MeshNode* getCornerNode(bool isXMax, bool isYMax) const
00231     {
00232       return getNode( isXMax ? _index._xSize-1 : 0 , isYMax ? _index._ySize-1 : 0 );
00233     }
00234     const SMDS_MeshElement* getCornerFace(const SMDS_MeshNode* cornerNode) const;
00236     bool isBound() const { return _nbBlocksExpected <= _nbBlocksFound; }
00238     gp_XYZ getXYZ(int x, int y) const { return SMESH_TNodeXYZ( getNode( x, y )); }
00240     gp_XYZ getGC() const
00241     {
00242       gp_XYZ xyz =
00243         getXYZ( 0, 0 ) +
00244         getXYZ( _index._xSize-1, 0 ) +
00245         getXYZ( 0, _index._ySize-1 ) +
00246         getXYZ( _index._xSize-1, _index._ySize-1 ) +
00247         getXYZ( _index._xSize/2, _index._ySize/2 );
00248       return xyz / 5;
00249     }
00251     int getNbFaces() const { return (_index._xSize-1) * (_index._ySize-1); }
00252   };
00253   //================================================================================
00257   struct _OrientedBlockSide
00258   {
00259     _BlockSide*       _side;
00260     _OrientedIndexer  _index;
00261 
00262     _OrientedBlockSide( _BlockSide* side=0, const int oriFlags=0 ):
00263       _side(side), _index(side ? side->_index : _Indexer(), oriFlags ) {}
00265     gp_XYZ xyz(int x, int y) const
00266     {
00267       return SMESH_TNodeXYZ( _grid_access_(_side, _index( x, y )) );
00268     }
00270     const SMDS_MeshNode* node(int x, int y) const
00271     {
00272       int i = _index( x, y );
00273       return ( i < 0 || i >= _side->_grid.size()) ? 0 : _side->_grid[i];
00274     }
00276     SMESH_OrientedLink edge(EQuadEdge edge) const
00277     {
00278       bool x1, y1, x2, y2; getEdgeEnds( edge, x1, y1, x2, y2 );
00279       return SMESH_OrientedLink( cornerNode ( x1, y1 ), cornerNode( x2, y2 ));
00280     }
00282     const SMDS_MeshNode* cornerNode(bool isXMax, bool isYMax) const
00283     {
00284       return _grid_access_(_side, _index.corner( isXMax, isYMax ));
00285     }
00287     int getHoriSize() const { return _index.xSize(); }
00288     int getVertSize() const  { return _index.ySize(); }
00290     operator bool() const { return _side; }
00292     const _BlockSide* operator->() const { return _side; }
00293     _BlockSide* operator->() { return _side; }
00294   };
00295   //================================================================================
00299   struct _Block
00300   {
00301     _OrientedBlockSide        _side[6]; // 6 sides of a sub-block
00302     set<const SMDS_MeshNode*> _corners;
00303 
00304     const _OrientedBlockSide& getSide(int i) const { return _side[i]; }
00305     bool setSide( int i, const _OrientedBlockSide& s)
00306     {
00307       if (( _side[i] = s ))
00308       {
00309         _corners.insert( s.cornerNode(0,0));
00310         _corners.insert( s.cornerNode(1,0));
00311         _corners.insert( s.cornerNode(0,1));
00312         _corners.insert( s.cornerNode(1,1));
00313       }
00314       return s;
00315     }
00316     void clear() { for (int i=0;i<6;++i) _side[i]=0; _corners.clear(); }
00317     bool hasSide( const _OrientedBlockSide& s) const
00318     {
00319       if ( s ) for (int i=0;i<6;++i) if ( _side[i] && _side[i]._side == s._side ) return true;
00320       return false;
00321     }
00322     int nbSides() const { int n=0; for (int i=0;i<6;++i) if ( _side[i] ) ++n; return n; }
00323     bool isValid() const;
00324   };
00325   //================================================================================
00329   class _Skin
00330   {
00331   public:
00332 
00333     int findBlocks(SMESH_Mesh& mesh);
00335     const _Block& getBlock(int i) const { return _blocks[i]; }
00337     const SMESH_Comment& error() const { return _error; }
00338 
00339   private:
00340     bool fillSide( _BlockSide&             side,
00341                    const SMDS_MeshElement* cornerQuad,
00342                    const SMDS_MeshNode*    cornerNode);
00343     bool fillRowsUntilCorner(const SMDS_MeshElement* quad,
00344                              const SMDS_MeshNode*    n1,
00345                              const SMDS_MeshNode*    n2,
00346                              vector<const SMDS_MeshNode*>& verRow1,
00347                              vector<const SMDS_MeshNode*>& verRow2,
00348                              bool alongN1N2 );
00349     _OrientedBlockSide findBlockSide( EBoxSides           startBlockSide,
00350                                       EQuadEdge           sharedSideEdge1,
00351                                       EQuadEdge           sharedSideEdge2,
00352                                       bool                withGeometricAnalysis,
00353                                       set< _BlockSide* >& sidesAround);
00355     void setSideBoundToBlock( _BlockSide& side )
00356     {
00357       if ( side._nbBlocksFound++, side.isBound() )
00358         for ( int e = 0; e < int(NB_QUAD_SIDES); ++e )
00359           _edge2sides[ side.getEdge( (EQuadEdge) e ) ].erase( &side );
00360     }
00362     int error(const SMESH_Comment& reason) { _error = reason; return 0; }
00363 
00364     SMESH_Comment      _error;
00365 
00366     list< _BlockSide > _allSides;
00367     vector< _Block >   _blocks;
00368 
00369     //map< const SMDS_MeshNode*, set< _BlockSide* > > _corner2sides;
00370     map< SMESH_OrientedLink, set< _BlockSide* > > _edge2sides;
00371   };
00372 
00373   //================================================================================
00377   //================================================================================
00378 
00379   int _Skin::findBlocks(SMESH_Mesh& mesh)
00380   {
00381     SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
00382 
00383     // Find a node at any block corner
00384 
00385     SMDS_NodeIteratorPtr nIt = meshDS->nodesIterator(/*idInceasingOrder=*/true);
00386     if ( !nIt->more() ) return error("Empty mesh");
00387 
00388     const SMDS_MeshNode* nCorner = 0;
00389     while ( nIt->more() )
00390     {
00391       nCorner = nIt->next();
00392       if ( isCornerNode( nCorner ))
00393         break;
00394       else
00395         nCorner = 0;
00396     }
00397     if ( !nCorner )
00398       return BAD_MESH_ERR;
00399 
00400     // --------------------------------------------------------------------
00401     // Find all block sides starting from mesh faces sharing the corner node
00402     // --------------------------------------------------------------------
00403 
00404     int nbFacesOnSides = 0;
00405     TIDSortedElemSet cornerFaces; // corner faces of found _BlockSide's
00406     list< const SMDS_MeshNode* > corners( 1, nCorner );
00407     list< const SMDS_MeshNode* >::iterator corner = corners.begin();
00408     while ( corner != corners.end() )
00409     {
00410       SMDS_ElemIteratorPtr faceIt = (*corner)->GetInverseElementIterator( SMDSAbs_Face );
00411       while ( faceIt->more() )
00412       {
00413         const SMDS_MeshElement* face = faceIt->next();
00414         if ( !cornerFaces.insert( face ).second )
00415           continue; // already loaded block side
00416 
00417         if ( !isQuadrangle( face ))
00418           return error("Non-quadrangle elements in the input mesh");
00419 
00420         if ( _allSides.empty() || !_allSides.back()._grid.empty() )
00421           _allSides.push_back( _BlockSide() );
00422 
00423         _BlockSide& side = _allSides.back();
00424         if ( !fillSide( side, face, *corner ) )
00425         {
00426           if ( !_error.empty() )
00427             return false;
00428         }
00429         else
00430         {
00431           for ( int isXMax = 0; isXMax < 2; ++isXMax )
00432             for ( int isYMax = 0; isYMax < 2; ++isYMax )
00433             {
00434               const SMDS_MeshNode* nCorner = side.getCornerNode(isXMax,isYMax );
00435               corners.push_back( nCorner );
00436               cornerFaces.insert( side.getCornerFace( nCorner ));
00437             }
00438           for ( int e = 0; e < int(NB_QUAD_SIDES); ++e )
00439             _edge2sides[ side.getEdge( (EQuadEdge) e ) ].insert( &side );
00440 
00441           nbFacesOnSides += side.getNbFaces();
00442         }
00443       }
00444       ++corner;
00445 
00446       // find block sides of other domains if any
00447       if ( corner == corners.end() && nbFacesOnSides < mesh.NbQuadrangles() )
00448       {
00449         while ( nIt->more() )
00450         {
00451           nCorner = nIt->next();
00452           if ( isCornerNode( nCorner ))
00453             corner = corners.insert( corner, nCorner );
00454         }
00455         nbFacesOnSides = mesh.NbQuadrangles();
00456       }
00457     }
00458     
00459     if ( _allSides.empty() )
00460       return BAD_MESH_ERR;
00461     if ( _allSides.back()._grid.empty() )
00462       _allSides.pop_back();
00463     _DUMP_("Nb detected sides "<< _allSides.size());
00464 
00465     // ---------------------------
00466     // Organize sides into blocks
00467     // ---------------------------
00468 
00469     // analyse sharing of sides by blocks and sort sides by nb of adjacent sides
00470     int nbBlockSides = 0; // total nb of block sides taking into account their sharing
00471     multimap<int, _BlockSide* > sortedSides;
00472     {
00473       list < _BlockSide >::iterator sideIt = _allSides.begin();
00474       for ( ; sideIt != _allSides.end(); ++sideIt )
00475       {
00476         _BlockSide& side = *sideIt;
00477         bool isSharedSide = true;
00478         int nbAdjacent = 0;
00479         for ( int e = 0; e < int(NB_QUAD_SIDES) && isSharedSide; ++e )
00480         {
00481           int nbAdj = _edge2sides[ side.getEdge( (EQuadEdge) e ) ].size();
00482           nbAdjacent += nbAdj;
00483           isSharedSide = ( nbAdj > 2 );
00484         }
00485         side._nbBlocksFound = 0;
00486         side._nbBlocksExpected = isSharedSide ? 2 : 1;
00487         nbBlockSides += side._nbBlocksExpected;
00488         sortedSides.insert( make_pair( nbAdjacent, & side ));
00489       }
00490     }
00491 
00492     // find sides of each block
00493     int nbBlocks = 0;
00494     while ( nbBlockSides >= 6 )
00495     {
00496       // get any side not bound to all blocks it belongs to
00497       multimap<int, _BlockSide*>::iterator i_side = sortedSides.begin();
00498       while ( i_side != sortedSides.end() && i_side->second->isBound())
00499         ++i_side;
00500 
00501       // start searching for block sides from the got side
00502       bool ok = true;
00503       if ( _blocks.empty() || _blocks.back()._side[B_FRONT] )
00504         _blocks.resize( _blocks.size() + 1 );
00505 
00506       _Block& block = _blocks.back();
00507       block.setSide( B_FRONT, i_side->second );
00508       setSideBoundToBlock( *i_side->second );
00509       nbBlockSides--;
00510 
00511       // edges of adjacent sides of B_FRONT corresponding to front's edges
00512       EQuadEdge edgeOfFront[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
00513       EQuadEdge edgeOfAdj  [4] = { Q_BOTTOM, Q_LEFT, Q_BOTTOM, Q_LEFT };
00514       // first find all sides detectable w/o advanced analysis,
00515       // then repeat the search, which then may pass without advanced analysis
00516       set< _BlockSide* > sidesAround;
00517       for ( int advAnalys = 0; advAnalys < 2; ++advAnalys )
00518       {
00519         // try to find 4 sides adjacent to a FRONT side
00520         for ( int i = 0; (ok || !advAnalys) && i < NB_QUAD_SIDES; ++i )
00521           if ( !block._side[i] )
00522             ok = block.setSide( i, findBlockSide( B_FRONT, edgeOfFront[i], edgeOfAdj[i],
00523                                                   advAnalys, sidesAround));
00524         // try to find a BACK side by a TOP one
00525         if ( ok || !advAnalys)
00526           if ( !block._side[B_BACK] && block._side[B_TOP] )
00527             ok = block.setSide( B_BACK, findBlockSide( B_TOP, Q_TOP, Q_TOP,
00528                                                        advAnalys, sidesAround ));
00529         if ( !advAnalys ) ok = true;
00530       }
00531       ok = block.isValid();
00532       if ( ok )
00533       {
00534         // check if just found block is same as one of previously found blocks
00535         bool isSame = false;
00536         for ( int i = 1; i < _blocks.size() && !isSame; ++i )
00537           isSame = ( block._corners == _blocks[i-1]._corners );
00538         ok = !isSame;
00539       }
00540 
00541       // count the found sides
00542       _DUMP_(endl << "** Block " << _blocks.size() << " valid: " << block.isValid());
00543       for (int i = 0; i < NB_BLOCK_SIDES; ++i )
00544       {
00545         _DUMP_("\tSide "<< SBoxSides[i] <<" "<< block._side[ i ]._side);
00546         if ( block._side[ i ] )
00547         {
00548           if ( ok && i != B_FRONT)
00549           {
00550             setSideBoundToBlock( *block._side[ i ]._side );
00551             nbBlockSides--;
00552           }
00553           _DUMP_("\t corners "<<
00554                  block._side[ i ].cornerNode(0,0)->GetID() << ", " <<
00555                  block._side[ i ].cornerNode(1,0)->GetID() << ", " <<
00556                  block._side[ i ].cornerNode(1,1)->GetID() << ", " <<
00557                  block._side[ i ].cornerNode(0,1)->GetID() << ", "<<endl);
00558         }
00559         else
00560         {
00561           _DUMP_("\t not found"<<endl);
00562         }
00563       }
00564       if ( !ok )
00565         block.clear();
00566       else
00567         nbBlocks++;
00568     }
00569     _DUMP_("Nb found blocks "<< nbBlocks <<endl);
00570 
00571     if ( nbBlocks == 0 && _error.empty() )
00572       return BAD_MESH_ERR;
00573 
00574     return nbBlocks;
00575   }
00576 
00577   //================================================================================
00581   //================================================================================
00582 
00583   bool _Skin::fillSide( _BlockSide&             side,
00584                         const SMDS_MeshElement* cornerQuad,
00585                         const SMDS_MeshNode*    nCorner)
00586   {
00587     // Find out size of block side mesured in nodes and by the way find two rows
00588     // of nodes in two directions.
00589 
00590     int x, y, nbX, nbY;
00591     const SMDS_MeshElement* firstQuad = cornerQuad;
00592     {
00593       // get a node on block edge
00594       int iCorner = firstQuad->GetNodeIndex( nCorner );
00595       const SMDS_MeshNode* nOnEdge = firstQuad->GetNode( (iCorner+1) % 4);
00596 
00597       // find out size of block side
00598       vector<const SMDS_MeshNode*> horRow1, horRow2, verRow1, verRow2;
00599       if ( !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, horRow1, horRow2, true ) ||
00600            !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, verRow1, verRow2, false ))
00601         return false;
00602       nbX = horRow1.size(), nbY = verRow1.size();
00603 
00604       // store found nodes
00605       side._index._xSize = horRow1.size();
00606       side._index._ySize = verRow1.size();
00607       side._grid.resize( side._index.size(), NULL );
00608 
00609       for ( x = 0; x < horRow1.size(); ++x )
00610       {
00611         side.setNode( x, 0, horRow1[x] );
00612         side.setNode( x, 1, horRow2[x] );
00613       }
00614       for ( y = 0; y < verRow1.size(); ++y )
00615       {
00616         side.setNode( 0, y, verRow1[y] );
00617         side.setNode( 1, y, verRow2[y] );
00618       }
00619     }
00620     // Find the rest nodes
00621 
00622     y = 1; // y of the row to fill
00623     TIDSortedElemSet emptySet, avoidSet;
00624     while ( ++y < nbY )
00625     {
00626       // get next firstQuad in the next row of quadrangles
00627       //
00628       //          n2up
00629       //     o---o               <- y row
00630       //     |   |
00631       //     o---o  o  o  o  o   <- found nodes
00632       //n1down    n2down       
00633       //
00634       int i1down, i2down, i2up;
00635       const SMDS_MeshNode* n1down = side.getNode( 0, y-1 );
00636       const SMDS_MeshNode* n2down = side.getNode( 1, y-1 );
00637       avoidSet.clear(); avoidSet.insert( firstQuad );
00638       firstQuad = SMESH_MeshEditor::FindFaceInSet( n1down, n2down, emptySet, avoidSet,
00639                                                    &i1down, &i2down);
00640       if ( !isQuadrangle( firstQuad ))
00641         return BAD_MESH_ERR;
00642 
00643       const SMDS_MeshNode* n2up = oppositeNode( firstQuad, i1down );
00644       avoidSet.clear(); avoidSet.insert( firstQuad );
00645 
00646       // find the rest nodes in the y-th row by faces in the row
00647 
00648       x = 1; 
00649       while ( ++x < nbX )
00650       {
00651         const SMDS_MeshElement* quad = SMESH_MeshEditor::FindFaceInSet( n2up, n2down, emptySet,
00652                                                                         avoidSet, &i2up, &i2down);
00653         if ( !isQuadrangle( quad ))
00654           return BAD_MESH_ERR;
00655 
00656         n2up   = oppositeNode( quad, i2down );
00657         n2down = oppositeNode( quad, i2up );
00658         avoidSet.clear(); avoidSet.insert( quad );
00659 
00660         side.setNode( x, y, n2up );
00661       }
00662     }
00663 
00664     // check side validity
00665     bool ok =
00666       side.getCornerFace( side.getCornerNode( 0, 0 )) &&
00667       side.getCornerFace( side.getCornerNode( 1, 0 )) &&
00668       side.getCornerFace( side.getCornerNode( 0, 1 )) &&
00669       side.getCornerFace( side.getCornerNode( 1, 1 ));
00670 
00671     return ok;
00672   }
00673 
00674   //================================================================================
00679   //================================================================================
00680 
00681   bool isClosedChainOfSides( _BlockSide*                                        startSide,
00682                              map< const SMDS_MeshNode*, list< _BlockSide* > > & corner2Sides )
00683   {
00684     // get start and end nodes
00685     const SMDS_MeshNode *n1 = 0, *n2 = 0, *n;
00686     for ( int y = 0; y < 2; ++y )
00687       for ( int x = 0; x < 2; ++x )
00688       {
00689         n = startSide->getCornerNode(x,y);
00690         if ( !corner2Sides.count( n )) continue;
00691         if ( n1 )
00692           n2 = n;
00693         else
00694           n1 = n;
00695       }
00696     if ( !n2 ) return false;
00697 
00698     map< const SMDS_MeshNode*, list< _BlockSide* > >::iterator
00699       c2sides = corner2Sides.find( n1 );
00700     if ( c2sides == corner2Sides.end() ) return false;
00701 
00702     int nbChainLinks = 1;
00703     n = n1;
00704     _BlockSide* prevSide = startSide;
00705     while ( n != n2 )
00706     {
00707       // get the next side sharing n
00708       list< _BlockSide* > & sides = c2sides->second;
00709       _BlockSide* nextSide = ( sides.back() == prevSide ? sides.front() : sides.back() );
00710       if ( nextSide == prevSide ) return false;
00711 
00712       // find the next corner of the nextSide being in corner2Sides
00713       n1 = n;
00714       n = 0;
00715       for ( int y = 0; y < 2 && !n; ++y )
00716         for ( int x = 0; x < 2; ++x )
00717         {
00718           n = nextSide->getCornerNode(x,y);
00719           c2sides = corner2Sides.find( n );
00720           if ( n == n1 || c2sides == corner2Sides.end() )
00721             n = 0;
00722           else
00723             break;
00724         }
00725       if ( !n ) return false;
00726 
00727       prevSide = nextSide;
00728       nbChainLinks++;
00729     }
00730 
00731     return ( n == n2 && nbChainLinks == NB_QUAD_SIDES );
00732   }
00733 
00734   //================================================================================
00738   //================================================================================
00739 
00740   _OrientedBlockSide _Skin::findBlockSide( EBoxSides           startBlockSide,
00741                                            EQuadEdge           sharedSideEdge1,
00742                                            EQuadEdge           sharedSideEdge2,
00743                                            bool                withGeometricAnalysis,
00744                                            set< _BlockSide* >& sidesAround)
00745   {
00746     _Block& block = _blocks.back();
00747     _OrientedBlockSide& side1 = block._side[ startBlockSide ];
00748 
00749     // get corner nodes of the given block edge
00750     SMESH_OrientedLink edge = side1.edge( sharedSideEdge1 );
00751     const SMDS_MeshNode* n1 = edge.node1();
00752     const SMDS_MeshNode* n2 = edge.node2();
00753     if ( edge._reversed ) swap( n1, n2 );
00754 
00755     // find all sides sharing both nodes n1 and n2
00756     set< _BlockSide* > sidesOnEdge = _edge2sides[ edge ]; // copy a set
00757 
00758     // exclude loaded sides of block from sidesOnEdge
00759     for (int i = 0; i < NB_BLOCK_SIDES; ++i )
00760       if ( block._side[ i ] )
00761         sidesOnEdge.erase( block._side[ i ]._side );
00762 
00763     int nbSidesOnEdge = sidesOnEdge.size();
00764     _DUMP_("nbSidesOnEdge "<< nbSidesOnEdge << " " << n1->GetID() << "-" << n2->GetID() );
00765     if ( nbSidesOnEdge == 0 )
00766       return 0;
00767 
00768     _BlockSide* foundSide = 0;
00769     if ( nbSidesOnEdge == 1 )
00770     {
00771       foundSide = *sidesOnEdge.begin();
00772     }
00773     else
00774     {
00775       set< _BlockSide* >::iterator sideIt = sidesOnEdge.begin();
00776       int nbLoadedSides = block.nbSides();
00777       if ( nbLoadedSides > 1 )
00778       {
00779         // Find the side having more than 2 corners common with already loaded sides
00780         for (; !foundSide && sideIt != sidesOnEdge.end(); ++sideIt )
00781         {
00782           _BlockSide* sideI = *sideIt;
00783           int nbCommonCorners =
00784             block._corners.count( sideI->getCornerNode(0,0)) +
00785             block._corners.count( sideI->getCornerNode(1,0)) +
00786             block._corners.count( sideI->getCornerNode(0,1)) +
00787             block._corners.count( sideI->getCornerNode(1,1));
00788           if ( nbCommonCorners > 2 )
00789             foundSide = sideI;
00790         }
00791       }
00792 
00793       if ( !foundSide )
00794       {
00795         if ( !withGeometricAnalysis )
00796         {
00797           sidesAround.insert( sidesOnEdge.begin(), sidesOnEdge.end() );
00798           return 0;
00799         }
00800         if ( nbLoadedSides == 1 )
00801         {
00802           // Issue 0021529. There are at least 2 sides by each edge and
00803           // position of block gravity center is undefined.
00804           // Find a side starting from which we can walk around the startBlockSide
00805 
00806           // fill in corner2Sides
00807           map< const SMDS_MeshNode*, list< _BlockSide* > > corner2Sides;
00808           for ( sideIt = sidesAround.begin(); sideIt != sidesAround.end(); ++sideIt )
00809           {
00810             _BlockSide* sideI = *sideIt;
00811             corner2Sides[ sideI->getCornerNode(0,0) ].push_back( sideI );
00812             corner2Sides[ sideI->getCornerNode(1,0) ].push_back( sideI );
00813             corner2Sides[ sideI->getCornerNode(0,1) ].push_back( sideI );
00814             corner2Sides[ sideI->getCornerNode(1,1) ].push_back( sideI );
00815           }
00816           // remove corners of startBlockSide from corner2Sides
00817           set<const SMDS_MeshNode*>::iterator nIt = block._corners.begin();
00818           for ( ; nIt != block._corners.end(); ++nIt )
00819             corner2Sides.erase( *nIt );
00820 
00821           // select a side
00822           for ( sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
00823           {
00824             if ( isClosedChainOfSides( *sideIt, corner2Sides ))
00825             {
00826               foundSide = *sideIt;
00827               break;
00828             }
00829           }
00830           if ( !foundSide )
00831             return 0;
00832         }
00833         else
00834         {
00835           // Select one of found sides most close to startBlockSide
00836 
00837           gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()),  p2 (n2->X(),n2->Y(),n2->Z());
00838           gp_Vec p1p2( p1, p2 );
00839 
00840           const SMDS_MeshElement* face1 = side1->getCornerFace( n1 );
00841           gp_XYZ p1Op = SMESH_TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(n1)));
00842           gp_Vec side1Dir( p1, p1Op );
00843           gp_Ax2 pln( p1, p1p2, side1Dir ); // plane with normal p1p2 and X dir side1Dir
00844           _DUMP_("  Select adjacent for "<< side1._side << " - side dir ("
00845                  << side1Dir.X() << ", " << side1Dir.Y() << ", " << side1Dir.Z() << ")" );
00846 
00847           map < double , _BlockSide* > angleOfSide;
00848           for (sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
00849           {
00850             _BlockSide* sideI = *sideIt;
00851             const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 );
00852             gp_XYZ p1Op = SMESH_TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(n1)));
00853             gp_Vec sideIDir( p1, p1Op );
00854             // compute angle of (sideIDir projection to pln) and (X dir of pln)
00855             gp_Vec2d sideIDirProj( sideIDir * pln.XDirection(), sideIDir * pln.YDirection());
00856             double angle = sideIDirProj.Angle( gp::DX2d() );
00857             if ( angle < 0 ) angle += 2. * M_PI; // angle [0-2*PI]
00858             angleOfSide.insert( make_pair( angle, sideI ));
00859             _DUMP_("  "<< sideI << " - side dir ("
00860                    << sideIDir.X() << ", " << sideIDir.Y() << ", " << sideIDir.Z() << ")"
00861                    << " angle " << angle);
00862           }
00863 
00864           gp_XYZ gc(0,0,0); // gravity center of already loaded block sides
00865           for (int i = 0; i < NB_BLOCK_SIDES; ++i )
00866             if ( block._side[ i ] )
00867               gc += block._side[ i ]._side->getGC();
00868           gc /= nbLoadedSides;
00869 
00870           gp_Vec gcDir( p1, gc );
00871           gp_Vec2d gcDirProj( gcDir * pln.XDirection(), gcDir * pln.YDirection());
00872           double gcAngle = gcDirProj.Angle( gp::DX2d() );
00873           foundSide = gcAngle < 0 ? angleOfSide.rbegin()->second : angleOfSide.begin()->second;
00874         }
00875       }
00876       _DUMP_("  selected "<< foundSide );
00877     }
00878 
00879     // Orient the found side correctly
00880 
00881     // corners of found side corresponding to nodes n1 and n2
00882     bool xMax1, yMax1, xMax2, yMax2;
00883     if ( !getEdgeEnds( sharedSideEdge2, xMax1, yMax1, xMax2, yMax2 ))
00884       return error(SMESH_Comment("Internal error at ")<<__FILE__<<":"<<__LINE__),
00885         _OrientedBlockSide(0);
00886 
00887     for ( int ori = 0; ori < _OrientedIndexer::MAX_ORI+1; ++ori )
00888     {
00889       _OrientedBlockSide orientedSide( foundSide, ori );
00890       const SMDS_MeshNode* n12 = orientedSide.cornerNode( xMax1, yMax1);
00891       const SMDS_MeshNode* n22 = orientedSide.cornerNode( xMax2, yMax2);
00892       if ( n1 == n12 && n2 == n22 )
00893         return orientedSide;
00894     }
00895     error(SMESH_Comment("Failed to orient a block side found by edge ")<<sharedSideEdge1
00896           << " of side " << startBlockSide
00897           << " of block " << _blocks.size());
00898     return 0;
00899   }
00900 
00901   //================================================================================
00907   //================================================================================
00908 
00909   bool _Skin::fillRowsUntilCorner(const SMDS_MeshElement*       quad,
00910                                   const SMDS_MeshNode*          n1,
00911                                   const SMDS_MeshNode*          n2,
00912                                   vector<const SMDS_MeshNode*>& row1,
00913                                   vector<const SMDS_MeshNode*>& row2,
00914                                   const bool                    alongN1N2 )
00915   {
00916     const SMDS_MeshNode* corner1 = n1;
00917 
00918     // Store nodes of quad in the rows and find new n1 and n2 to get
00919     // the next face so that new n2 is on block edge
00920     int i1 = quad->GetNodeIndex( n1 );
00921     int i2 = quad->GetNodeIndex( n2 );
00922     row1.clear(); row2.clear();
00923     row1.push_back( n1 );
00924     if ( alongN1N2 )
00925     {
00926       row1.push_back( n2 );
00927       row2.push_back( oppositeNode( quad, i2 ));
00928       row2.push_back( n1 = oppositeNode( quad, i1 ));
00929     }
00930     else
00931     {
00932       row2.push_back( n2 );
00933       row1.push_back( n2 = oppositeNode( quad, i2 ));
00934       row2.push_back( n1 = oppositeNode( quad, i1 ));
00935     }
00936 
00937     if ( isCornerNode( row1[1] ))
00938       return true;
00939 
00940     // Find the rest nodes
00941     TIDSortedElemSet emptySet, avoidSet;
00942     while ( !isCornerNode( n2 ) )
00943     {
00944       avoidSet.clear(); avoidSet.insert( quad );
00945       quad = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2 );
00946       if ( !isQuadrangle( quad ))
00947         return BAD_MESH_ERR;
00948 
00949       row1.push_back( n2 = oppositeNode( quad, i1 ));
00950       row2.push_back( n1 = oppositeNode( quad, i2 ));
00951     }
00952     return n1 != corner1;
00953   }
00954 
00955   //================================================================================
00959   //================================================================================
00960 
00961   const SMDS_MeshElement* _BlockSide::getCornerFace(const SMDS_MeshNode* cornerNode) const
00962   {
00963     int x, y, isXMax, isYMax, found = 0;
00964     for ( isXMax = 0; isXMax < 2; ++isXMax )
00965     {
00966       for ( isYMax = 0; isYMax < 2; ++isYMax )
00967       {
00968         x = isXMax ? _index._xSize-1 : 0;
00969         y = isYMax ? _index._ySize-1 : 0;
00970         found = ( getNode(x,y) == cornerNode );
00971         if ( found ) break;
00972       }
00973       if ( found ) break;
00974     }
00975     if ( !found ) return 0;
00976     int dx = isXMax ? -1 : +1;
00977     int dy = isYMax ? -1 : +1;
00978     const SMDS_MeshNode* n1 = getNode(x,y);
00979     const SMDS_MeshNode* n2 = getNode(x+dx,y);
00980     const SMDS_MeshNode* n3 = getNode(x,y+dy);
00981     const SMDS_MeshNode* n4 = getNode(x+dx,y+dy);
00982     return SMDS_Mesh::FindFace(n1, n2, n3, n4 );
00983   }
00984 
00985   //================================================================================
00989   //================================================================================
00990 
00991   bool _Block::isValid() const
00992   {
00993     bool ok = ( nbSides() == 6 );
00994 
00995     // check only corners depending on side selection
00996     EBoxSides adjacent[4] = { B_BOTTOM, B_RIGHT, B_TOP, B_LEFT };
00997     EQuadEdge edgeAdj [4] = { Q_TOP,    Q_RIGHT, Q_TOP, Q_RIGHT };
00998     EQuadEdge edgeBack[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
00999 
01000     for ( int i=0; ok && i < NB_QUAD_SIDES; ++i )
01001     { 
01002       SMESH_OrientedLink eBack = _side[ B_BACK      ].edge( edgeBack[i] );
01003       SMESH_OrientedLink eAdja = _side[ adjacent[i] ].edge( edgeAdj[i] );
01004       ok = ( eBack == eAdja );
01005     }
01006     return ok;
01007   }
01008 
01009 } // namespace
01010 
01011 //=======================================================================
01012 //function : StdMeshers_HexaFromSkin_3D
01013 //purpose  : 
01014 //=======================================================================
01015 
01016 StdMeshers_HexaFromSkin_3D::StdMeshers_HexaFromSkin_3D(int hypId, int studyId, SMESH_Gen* gen)
01017   :SMESH_3D_Algo(hypId, studyId, gen)
01018 {
01019   MESSAGE("StdMeshers_HexaFromSkin_3D::StdMeshers_HexaFromSkin_3D");
01020   _name = "HexaFromSkin_3D";
01021 }
01022 
01023 StdMeshers_HexaFromSkin_3D::~StdMeshers_HexaFromSkin_3D()
01024 {
01025   MESSAGE("StdMeshers_HexaFromSkin_3D::~StdMeshers_HexaFromSkin_3D");
01026 }
01027 
01028 //================================================================================
01032 //================================================================================
01033 
01034 bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
01035 {
01036   _Skin skin;
01037   int nbBlocks = skin.findBlocks(aMesh);
01038   if ( nbBlocks == 0 )
01039     return error( skin.error());
01040 
01041   vector< vector< const SMDS_MeshNode* > > columns;
01042   int x, xSize, y, ySize, z, zSize;
01043   _Indexer colIndex;
01044 
01045   for ( int i = 0; i < nbBlocks; ++i )
01046   {
01047     const _Block& block = skin.getBlock( i );
01048 
01049     // ------------------------------------------
01050     // Fill columns of nodes with existing nodes
01051     // ------------------------------------------
01052 
01053     xSize = block.getSide(B_BOTTOM).getHoriSize();
01054     ySize = block.getSide(B_BOTTOM).getVertSize();
01055     zSize = block.getSide(B_FRONT ).getVertSize();
01056     int X = xSize - 1, Y = ySize - 1, Z = zSize - 1;
01057     colIndex = _Indexer( xSize, ySize );
01058     columns.resize( colIndex.size() );
01059 
01060     // fill node columns by front and back box sides
01061     for ( x = 0; x < xSize; ++x ) {
01062       vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
01063       vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
01064       column0.resize( zSize );
01065       column1.resize( zSize );
01066       for ( z = 0; z < zSize; ++z ) {
01067         column0[ z ] = block.getSide(B_FRONT).node( x, z );
01068         column1[ z ] = block.getSide(B_BACK) .node( x, z );
01069       }
01070     }
01071     // fill node columns by left and right box sides
01072     for ( y = 1; y < ySize-1; ++y ) {
01073       vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
01074       vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
01075       column0.resize( zSize );
01076       column1.resize( zSize );
01077       for ( z = 0; z < zSize; ++z ) {
01078         column0[ z ] = block.getSide(B_LEFT) .node( y, z );
01079         column1[ z ] = block.getSide(B_RIGHT).node( y, z );
01080       }
01081     }
01082     // get nodes from top and bottom box sides
01083     for ( x = 1; x < xSize-1; ++x ) {
01084       for ( y = 1; y < ySize-1; ++y ) {
01085         vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
01086         column.resize( zSize );
01087         column.front() = block.getSide(B_BOTTOM).node( x, y );
01088         column.back()  = block.getSide(B_TOP)   .node( x, y );
01089       }
01090     }
01091 
01092     // ----------------------------
01093     // Add internal nodes of a box
01094     // ----------------------------
01095     // projection points of internal nodes on box sub-shapes by which
01096     // coordinates of internal nodes are computed
01097     vector<gp_XYZ> pointOnShape( SMESH_Block::ID_Shell );
01098 
01099     // projections on vertices are constant
01100     pointOnShape[ SMESH_Block::ID_V000 ] = block.getSide(B_BOTTOM).xyz( 0, 0 );
01101     pointOnShape[ SMESH_Block::ID_V100 ] = block.getSide(B_BOTTOM).xyz( X, 0 );
01102     pointOnShape[ SMESH_Block::ID_V010 ] = block.getSide(B_BOTTOM).xyz( 0, Y );
01103     pointOnShape[ SMESH_Block::ID_V110 ] = block.getSide(B_BOTTOM).xyz( X, Y );
01104     pointOnShape[ SMESH_Block::ID_V001 ] = block.getSide(B_TOP).xyz( 0, 0 );
01105     pointOnShape[ SMESH_Block::ID_V101 ] = block.getSide(B_TOP).xyz( X, 0 );
01106     pointOnShape[ SMESH_Block::ID_V011 ] = block.getSide(B_TOP).xyz( 0, Y );
01107     pointOnShape[ SMESH_Block::ID_V111 ] = block.getSide(B_TOP).xyz( X, Y );
01108 
01109     for ( x = 1; x < xSize-1; ++x )
01110     {
01111       gp_XYZ params; // normalized parameters of internal node within a unit box
01112       params.SetCoord( 1, x / double(X) );
01113       for ( y = 1; y < ySize-1; ++y )
01114       {
01115         params.SetCoord( 2, y / double(Y) );
01116         // column to fill during z loop
01117         vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
01118         // projections on horizontal edges
01119         pointOnShape[ SMESH_Block::ID_Ex00 ] = block.getSide(B_BOTTOM).xyz( x, 0 );
01120         pointOnShape[ SMESH_Block::ID_Ex10 ] = block.getSide(B_BOTTOM).xyz( x, Y );
01121         pointOnShape[ SMESH_Block::ID_E0y0 ] = block.getSide(B_BOTTOM).xyz( 0, y );
01122         pointOnShape[ SMESH_Block::ID_E1y0 ] = block.getSide(B_BOTTOM).xyz( X, y );
01123         pointOnShape[ SMESH_Block::ID_Ex01 ] = block.getSide(B_TOP).xyz( x, 0 );
01124         pointOnShape[ SMESH_Block::ID_Ex11 ] = block.getSide(B_TOP).xyz( x, Y );
01125         pointOnShape[ SMESH_Block::ID_E0y1 ] = block.getSide(B_TOP).xyz( 0, y );
01126         pointOnShape[ SMESH_Block::ID_E1y1 ] = block.getSide(B_TOP).xyz( X, y );
01127         // projections on horizontal sides
01128         pointOnShape[ SMESH_Block::ID_Fxy0 ] = block.getSide(B_BOTTOM).xyz( x, y );
01129         pointOnShape[ SMESH_Block::ID_Fxy1 ] = block.getSide(B_TOP)   .xyz( x, y );
01130         for ( z = 1; z < zSize-1; ++z ) // z loop
01131         {
01132           params.SetCoord( 3, z / double(Z) );
01133           // projections on vertical edges
01134           pointOnShape[ SMESH_Block::ID_E00z ] = block.getSide(B_FRONT).xyz( 0, z );    
01135           pointOnShape[ SMESH_Block::ID_E10z ] = block.getSide(B_FRONT).xyz( X, z );    
01136           pointOnShape[ SMESH_Block::ID_E01z ] = block.getSide(B_BACK).xyz( 0, z );    
01137           pointOnShape[ SMESH_Block::ID_E11z ] = block.getSide(B_BACK).xyz( X, z );
01138           // projections on vertical sides
01139           pointOnShape[ SMESH_Block::ID_Fx0z ] = block.getSide(B_FRONT).xyz( x, z );    
01140           pointOnShape[ SMESH_Block::ID_Fx1z ] = block.getSide(B_BACK) .xyz( x, z );    
01141           pointOnShape[ SMESH_Block::ID_F0yz ] = block.getSide(B_LEFT) .xyz( y, z );    
01142           pointOnShape[ SMESH_Block::ID_F1yz ] = block.getSide(B_RIGHT).xyz( y, z );
01143 
01144           // compute internal node coordinates
01145           gp_XYZ coords;
01146           SMESH_Block::ShellPoint( params, pointOnShape, coords );
01147           column[ z ] = aHelper->AddNode( coords.X(), coords.Y(), coords.Z() );
01148 
01149 #ifdef DEB_GRID
01150           // debug
01151           //cout << "----------------------------------------------------------------------"<<endl;
01152           //for ( int id = SMESH_Block::ID_V000; id < SMESH_Block::ID_Shell; ++id)
01153           //{
01154           //  gp_XYZ p = pointOnShape[ id ];
01155           //  SMESH_Block::DumpShapeID( id,cout)<<" ( "<<p.X()<<", "<<p.Y()<<", "<<p.Z()<<" )"<<endl;
01156           //}
01157           //cout << "Params: ( "<< params.X()<<", "<<params.Y()<<", "<<params.Z()<<" )"<<endl;
01158           //cout << "coords: ( "<< coords.X()<<", "<<coords.Y()<<", "<<coords.Z()<<" )"<<endl;
01159 #endif
01160         }
01161       }
01162     }
01163     // ----------------
01164     // Add hexahedrons
01165     // ----------------
01166 
01167     // find out orientation by a least distorted hexahedron (issue 0020855);
01168     // the last is defined by evaluating sum of face normals of 8 corner hexahedrons
01169     double badness = numeric_limits<double>::max();
01170     bool isForw = true;
01171     for ( int xMax = 0; xMax < 2; ++xMax )
01172       for ( int yMax = 0; yMax < 2; ++yMax )
01173         for ( int zMax = 0; zMax < 2; ++zMax )
01174         {
01175           x = xMax ? xSize-1 : 1;
01176           y = yMax ? ySize-1 : 1;
01177           z = zMax ? zSize-1 : 1;
01178           vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x-1, y-1 )];
01179           vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x  , y-1 )];
01180           vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x-1, y   )];
01181           vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x  , y )];
01182           
01183           const SMDS_MeshNode* n000 = col00[z-1];
01184           const SMDS_MeshNode* n100 = col10[z-1];
01185           const SMDS_MeshNode* n010 = col01[z-1];
01186           const SMDS_MeshNode* n110 = col11[z-1];
01187           const SMDS_MeshNode* n001 = col00[z];
01188           const SMDS_MeshNode* n101 = col10[z];
01189           const SMDS_MeshNode* n011 = col01[z];
01190           const SMDS_MeshNode* n111 = col11[z];
01191           SMDS_VolumeOfNodes probeVolume (n000,n010,n110,n100,
01192                                           n001,n011,n111,n101);
01193           SMDS_VolumeTool volTool( &probeVolume );
01194           double Nx=0.,Ny=0.,Nz=0.;
01195           for ( int iFace = 0; iFace < volTool.NbFaces(); ++iFace )
01196           {
01197             double nx,ny,nz;
01198             volTool.GetFaceNormal( iFace, nx,ny,nz );
01199             Nx += nx;
01200             Ny += ny;
01201             Nz += nz;
01202           }
01203           double quality = Nx*Nx + Ny*Ny + Nz*Nz;
01204           if ( quality < badness )
01205           {
01206             badness = quality;
01207             isForw = volTool.IsForward();
01208           }
01209         }
01210 
01211     // add elements
01212     for ( x = 0; x < xSize-1; ++x ) {
01213       for ( y = 0; y < ySize-1; ++y ) {
01214         vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
01215         vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
01216         vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
01217         vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
01218         // bottom face normal of a hexa mush point outside the volume
01219         if ( isForw )
01220           for ( z = 0; z < zSize-1; ++z )
01221             aHelper->AddVolume(col00[z],   col01[z],   col11[z],   col10[z],
01222                                col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
01223         else
01224           for ( z = 0; z < zSize-1; ++z )
01225             aHelper->AddVolume(col00[z],   col10[z],   col11[z],   col01[z],
01226                                col00[z+1], col10[z+1], col11[z+1], col01[z+1]);
01227       }
01228     }
01229   } // loop on blocks
01230 
01231   return true;
01232 }
01233 
01234 //================================================================================
01238 //================================================================================
01239 
01240 bool StdMeshers_HexaFromSkin_3D::Evaluate(SMESH_Mesh &         aMesh,
01241                                           const TopoDS_Shape & aShape,
01242                                           MapShapeNbElems&     aResMap)
01243 {
01244   _Skin skin;
01245   int nbBlocks = skin.findBlocks(aMesh);
01246   if ( nbBlocks == 0 )
01247     return error( skin.error());
01248 
01249   bool secondOrder = aMesh.NbFaces( ORDER_QUADRATIC );
01250 
01251   int entity = secondOrder ? SMDSEntity_Quad_Hexa : SMDSEntity_Hexa;
01252   vector<int>& nbByType = aResMap[ aMesh.GetSubMesh( aShape )];
01253   if ( entity >= nbByType.size() )
01254     nbByType.resize( SMDSEntity_Last, 0 );
01255 
01256   for ( int i = 0; i < nbBlocks; ++i )
01257   {
01258     const _Block& block = skin.getBlock( i );
01259 
01260     int nbX = block.getSide(B_BOTTOM).getHoriSize();
01261     int nbY = block.getSide(B_BOTTOM).getVertSize();
01262     int nbZ = block.getSide(B_FRONT ).getVertSize();
01263 
01264     int nbHexa  = (nbX-1) * (nbY-1) * (nbZ-1);
01265     int nbNodes = (nbX-2) * (nbY-2) * (nbZ-2);
01266     if ( secondOrder )
01267       nbNodes +=
01268         (nbX-2) * (nbY-2) * (nbZ-1) +
01269         (nbX-2) * (nbY-1) * (nbZ-2) +
01270         (nbX-1) * (nbY-2) * (nbZ-2);
01271 
01272 
01273     nbByType[ entity ] += nbHexa;
01274     nbByType[ SMDSEntity_Node ] += nbNodes;
01275   }
01276 
01277   return true;
01278 }
01279 
01280 //================================================================================
01284 //================================================================================
01285 
01286 bool StdMeshers_HexaFromSkin_3D::CheckHypothesis(SMESH_Mesh&, const TopoDS_Shape&,
01287                                                  Hypothesis_Status& aStatus)
01288 {
01289   aStatus = SMESH_Hypothesis::HYP_OK;
01290   return true;
01291 }
01292 
01293 //================================================================================
01298 //================================================================================
01299 
01300 bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh&, const TopoDS_Shape&)
01301 {
01302   return error("Algorithm can't work with geometrical shapes");
01303 }