]> SALOME platform Git repositories - modules/smesh.git/blob - src/StdMeshers/StdMeshers_HexaFromSkin_3D.cxx
Salome HOME
0021825: Error in the example of "Projection Algorithms" in the user's guide
[modules/smesh.git] / src / StdMeshers / StdMeshers_HexaFromSkin_3D.cxx
1 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
3 // This library is free software; you can redistribute it and/or
4 // modify it under the terms of the GNU Lesser General Public
5 // License as published by the Free Software Foundation; either
6 // version 2.1 of the License.
7 //
8 // This library is distributed in the hope that it will be useful,
9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11 // Lesser General Public License for more details.
12 //
13 // You should have received a copy of the GNU Lesser General Public
14 // License along with this library; if not, write to the Free Software
15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
16 //
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
18 //
19
20 // File      : StdMeshers_HexaFromSkin_3D.cxx
21 // Created   : Wed Jan 27 12:28:07 2010
22 // Author    : Edward AGAPOV (eap)
23
24 #include "StdMeshers_HexaFromSkin_3D.hxx"
25
26 #include "SMDS_VolumeOfNodes.hxx"
27 #include "SMDS_VolumeTool.hxx"
28 #include "SMESH_Block.hxx"
29 #include "SMESH_MesherHelper.hxx"
30
31 #include <gp_Ax2.hxx>
32
33 //#include "utilities.h"
34 #include <limits>
35
36 // Define error message and _MYDEBUG_ if needed
37 #ifdef _DEBUG_
38 #define BAD_MESH_ERR \
39   error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh.\n" \
40                       __FILE__ ":" )<<__LINE__)
41 //#define _MYDEBUG_
42 #else
43 #define BAD_MESH_ERR \
44   error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh"))
45 #endif
46
47
48 // Debug output
49 #ifdef _MYDEBUG_
50 #define _DUMP_(msg) cout << msg << endl
51 #else
52 #define _DUMP_(msg)
53 #endif
54
55
56 namespace
57 {
58   enum EBoxSides //!< sides of the block
59     {
60       B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, NB_BLOCK_SIDES
61     };
62 #ifdef _MYDEBUG_
63   const char* SBoxSides[] = //!< names of block sides -- needed for DEBUG only
64     {
65       "BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED"
66     };
67 #endif
68   enum EQuadEdge //!< edges of quadrangle side
69     {
70       Q_BOTTOM = 0, Q_RIGHT, Q_TOP, Q_LEFT, NB_QUAD_SIDES
71     };
72
73
74   //================================================================================
75   /*!
76    * \brief return logical coordinates (i.e. min or max) of ends of edge
77    */
78   //================================================================================
79
80   bool getEdgeEnds(EQuadEdge edge, bool& xMax1, bool& yMax1, bool& xMax2, bool& yMax2 )
81   {
82     xMax1=0, yMax1=0, xMax2=1, yMax2=1;
83     switch( edge )
84     {
85     case Q_BOTTOM: yMax2 = 0; break;
86     case Q_RIGHT:  xMax1 = 1; break;
87     case Q_TOP:    yMax1 = 1; break;
88     case Q_LEFT:   xMax2 = 0; break;
89     default:
90       return false;
91     }
92     return true;
93   }
94
95   //================================================================================
96   /*!
97    * \brief return true if a node is at block corner
98    *
99    * This check is valid for simple cases only
100    */
101   //================================================================================
102
103   bool isCornerNode( const SMDS_MeshNode* n )
104   {
105     int nbF = n ? n->NbInverseElements( SMDSAbs_Face ) : 1;
106     if ( nbF % 2 )
107       return true;
108
109     set<const SMDS_MeshNode*> nodesInInverseFaces;
110     SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face );
111     while ( fIt->more() )
112     {
113       const SMDS_MeshElement* face = fIt->next();
114       nodesInInverseFaces.insert( face->begin_nodes(), face->end_nodes() );
115     }
116
117     return nodesInInverseFaces.size() != ( 6 + (nbF/2-1)*3 );
118   }
119
120   //================================================================================
121   /*!
122    * \brief check element type
123    */
124   //================================================================================
125
126   bool isQuadrangle(const SMDS_MeshElement* e)
127   {
128     return ( e && e->NbCornerNodes() == 4 );
129   }
130
131   //================================================================================
132   /*!
133    * \brief return opposite node of a quadrangle face
134    */
135   //================================================================================
136
137   const SMDS_MeshNode* oppositeNode(const SMDS_MeshElement* quad, int iNode)
138   {
139     return quad->GetNode( (iNode+2) % 4 );
140   }
141
142   //================================================================================
143   /*!
144    * \brief Convertor of a pair of integers to a sole index
145    */
146   struct _Indexer
147   {
148     int _xSize, _ySize;
149     _Indexer( int xSize=0, int ySize=0 ): _xSize(xSize), _ySize(ySize) {}
150     int size() const { return _xSize * _ySize; }
151     int operator()(int x, int y) const { return y * _xSize + x; }
152   };
153   //================================================================================
154   /*!
155    * \brief Oriented convertor of a pair of integers to a sole index 
156    */
157   class _OrientedIndexer : public _Indexer
158   {
159   public:
160     enum OriFlags //!< types of block side orientation
161       {
162         REV_X = 1, REV_Y = 2, SWAP_XY = 4, MAX_ORI = REV_X|REV_Y|SWAP_XY
163       };
164     _OrientedIndexer( const _Indexer& indexer, const int oriFlags ):
165       _Indexer( indexer._xSize, indexer._ySize ),
166       _xSize (indexer._xSize), _ySize(indexer._ySize),
167       _xRevFun((oriFlags & REV_X) ? & reverse : & lazy),
168       _yRevFun((oriFlags & REV_Y) ? & reverse : & lazy),
169       _swapFun((oriFlags & SWAP_XY ) ? & swap : & lazy)
170     {
171       (*_swapFun)( _xSize, _ySize );
172     }
173     //!< Return index by XY
174     int operator()(int x, int y) const
175     {
176       (*_xRevFun)( x, const_cast<int&>( _xSize ));
177       (*_yRevFun)( y, const_cast<int&>( _ySize ));
178       (*_swapFun)( x, y );
179       return _Indexer::operator()(x,y);
180     }
181     //!< Return index for a corner
182     int corner(bool xMax, bool yMax) const
183     {
184       int x = xMax, y = yMax, size = 2;
185       (*_xRevFun)( x, size );
186       (*_yRevFun)( y, size );
187       (*_swapFun)( x, y );
188       return _Indexer::operator()(x ? _Indexer::_xSize-1 : 0 , y ? _Indexer::_ySize-1 : 0);
189     }
190     int xSize() const { return _xSize; }
191     int ySize() const { return _ySize; }
192   private:
193     _Indexer _indexer;
194     int _xSize, _ySize;
195
196     typedef void (*TFun)(int& x, int& y);
197     TFun _xRevFun, _yRevFun, _swapFun;
198     
199     static void lazy   (int&, int&) {}
200     static void reverse(int& x, int& size) { x = size - x - 1; }
201     static void swap   (int& x, int& y) { std::swap(x,y); }
202   };
203   //================================================================================
204   /*!
205    * \brief Structure corresponding to the meshed side of block
206    */
207   struct _BlockSide
208   {
209     vector<const SMDS_MeshNode*> _grid;
210     _Indexer                     _index;
211     int                          _nbBlocksExpected;
212     int                          _nbBlocksFound;
213
214 #ifdef _DEBUG_ // want to get SIGSEGV in case of invalid index
215 #define _grid_access_(pobj, i) pobj->_grid[ ((i) < pobj->_grid.size()) ? i : int(1e100)]
216 #else
217 #define _grid_access_(pobj, i) pobj->_grid[ i ]
218 #endif
219     //!< Return node at XY
220     const SMDS_MeshNode* getNode(int x, int y) const { return _grid_access_(this, _index( x,y ));}
221     //!< Set node at XY
222     void setNode(int x, int y, const SMDS_MeshNode* n) { _grid_access_(this, _index( x,y )) = n; }
223     //!< Return an edge
224     SMESH_OrientedLink getEdge(EQuadEdge edge) const
225     {
226       bool x1, y1, x2, y2; getEdgeEnds( edge, x1, y1, x2, y2 );
227       return SMESH_OrientedLink( getCornerNode ( x1, y1 ), getCornerNode( x2, y2 ));
228     }
229     //!< Return a corner node
230     const SMDS_MeshNode* getCornerNode(bool isXMax, bool isYMax) const
231     {
232       return getNode( isXMax ? _index._xSize-1 : 0 , isYMax ? _index._ySize-1 : 0 );
233     }
234     const SMDS_MeshElement* getCornerFace(const SMDS_MeshNode* cornerNode) const;
235     //!< True if all blocks this side belongs to have been found
236     bool isBound() const { return _nbBlocksExpected <= _nbBlocksFound; }
237     //!< Return coordinates of node at XY
238     gp_XYZ getXYZ(int x, int y) const { return SMESH_TNodeXYZ( getNode( x, y )); }
239     //!< Return gravity center of the four corners and the middle node
240     gp_XYZ getGC() const
241     {
242       gp_XYZ xyz =
243         getXYZ( 0, 0 ) +
244         getXYZ( _index._xSize-1, 0 ) +
245         getXYZ( 0, _index._ySize-1 ) +
246         getXYZ( _index._xSize-1, _index._ySize-1 ) +
247         getXYZ( _index._xSize/2, _index._ySize/2 );
248       return xyz / 5;
249     }
250     //!< Return number of mesh faces
251     int getNbFaces() const { return (_index._xSize-1) * (_index._ySize-1); }
252   };
253   //================================================================================
254   /*!
255    * \brief _BlockSide with changed orientation
256    */
257   struct _OrientedBlockSide
258   {
259     _BlockSide*       _side;
260     _OrientedIndexer  _index;
261
262     _OrientedBlockSide( _BlockSide* side=0, const int oriFlags=0 ):
263       _side(side), _index(side ? side->_index : _Indexer(), oriFlags ) {}
264     //!< return coordinates by XY
265     gp_XYZ xyz(int x, int y) const
266     {
267       return SMESH_TNodeXYZ( _grid_access_(_side, _index( x, y )) );
268     }
269     //!< safely return a node by XY
270     const SMDS_MeshNode* node(int x, int y) const
271     {
272       int i = _index( x, y );
273       return ( i < 0 || i >= _side->_grid.size()) ? 0 : _side->_grid[i];
274     }
275     //!< Return an edge
276     SMESH_OrientedLink edge(EQuadEdge edge) const
277     {
278       bool x1, y1, x2, y2; getEdgeEnds( edge, x1, y1, x2, y2 );
279       return SMESH_OrientedLink( cornerNode ( x1, y1 ), cornerNode( x2, y2 ));
280     }
281     //!< Return a corner node
282     const SMDS_MeshNode* cornerNode(bool isXMax, bool isYMax) const
283     {
284       return _grid_access_(_side, _index.corner( isXMax, isYMax ));
285     }
286     //!< return its size in nodes
287     int getHoriSize() const { return _index.xSize(); }
288     int getVertSize() const  { return _index.ySize(); }
289     //!< True if _side has been initialized
290     operator bool() const { return _side; }
291     //! Direct access to _side
292     const _BlockSide* operator->() const { return _side; }
293     _BlockSide* operator->() { return _side; }
294   };
295   //================================================================================
296   /*!
297    * \brief Meshed skin of block
298    */
299   struct _Block
300   {
301     _OrientedBlockSide        _side[6]; // 6 sides of a sub-block
302     set<const SMDS_MeshNode*> _corners;
303
304     const _OrientedBlockSide& getSide(int i) const { return _side[i]; }
305     bool setSide( int i, const _OrientedBlockSide& s)
306     {
307       if (( _side[i] = s ))
308       {
309         _corners.insert( s.cornerNode(0,0));
310         _corners.insert( s.cornerNode(1,0));
311         _corners.insert( s.cornerNode(0,1));
312         _corners.insert( s.cornerNode(1,1));
313       }
314       return s;
315     }
316     void clear() { for (int i=0;i<6;++i) _side[i]=0; _corners.clear(); }
317     bool hasSide( const _OrientedBlockSide& s) const
318     {
319       if ( s ) for (int i=0;i<6;++i) if ( _side[i] && _side[i]._side == s._side ) return true;
320       return false;
321     }
322     int nbSides() const { int n=0; for (int i=0;i<6;++i) if ( _side[i] ) ++n; return n; }
323     bool isValid() const;
324   };
325   //================================================================================
326   /*!
327    * \brief Skin mesh possibly containing several meshed blocks
328    */
329   class _Skin
330   {
331   public:
332
333     int findBlocks(SMESH_Mesh& mesh);
334     //!< return i-th block
335     const _Block& getBlock(int i) const { return _blocks[i]; }
336     //!< return error description
337     const SMESH_Comment& error() const { return _error; }
338
339   private:
340     bool fillSide( _BlockSide&             side,
341                    const SMDS_MeshElement* cornerQuad,
342                    const SMDS_MeshNode*    cornerNode);
343     bool fillRowsUntilCorner(const SMDS_MeshElement* quad,
344                              const SMDS_MeshNode*    n1,
345                              const SMDS_MeshNode*    n2,
346                              vector<const SMDS_MeshNode*>& verRow1,
347                              vector<const SMDS_MeshNode*>& verRow2,
348                              bool alongN1N2 );
349     _OrientedBlockSide findBlockSide( EBoxSides           startBlockSide,
350                                       EQuadEdge           sharedSideEdge1,
351                                       EQuadEdge           sharedSideEdge2,
352                                       bool                withGeometricAnalysis,
353                                       set< _BlockSide* >& sidesAround);
354     //!< update own data and data of the side bound to block
355     void setSideBoundToBlock( _BlockSide& side )
356     {
357       if ( side._nbBlocksFound++, side.isBound() )
358         for ( int e = 0; e < int(NB_QUAD_SIDES); ++e )
359           _edge2sides[ side.getEdge( (EQuadEdge) e ) ].erase( &side );
360     }
361     //!< store reason of error
362     int error(const SMESH_Comment& reason) { _error = reason; return 0; }
363
364     SMESH_Comment      _error;
365
366     list< _BlockSide > _allSides;
367     vector< _Block >   _blocks;
368
369     //map< const SMDS_MeshNode*, set< _BlockSide* > > _corner2sides;
370     map< SMESH_OrientedLink, set< _BlockSide* > > _edge2sides;
371   };
372
373   //================================================================================
374   /*!
375    * \brief Find and return number of submeshes corresponding to blocks
376    */
377   //================================================================================
378
379   int _Skin::findBlocks(SMESH_Mesh& mesh)
380   {
381     SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
382
383     // Find a node at any block corner
384
385     SMDS_NodeIteratorPtr nIt = meshDS->nodesIterator(/*idInceasingOrder=*/true);
386     if ( !nIt->more() ) return error("Empty mesh");
387
388     const SMDS_MeshNode* nCorner = 0;
389     while ( nIt->more() )
390     {
391       nCorner = nIt->next();
392       if ( isCornerNode( nCorner ))
393         break;
394       else
395         nCorner = 0;
396     }
397     if ( !nCorner )
398       return BAD_MESH_ERR;
399
400     // --------------------------------------------------------------------
401     // Find all block sides starting from mesh faces sharing the corner node
402     // --------------------------------------------------------------------
403
404     int nbFacesOnSides = 0;
405     TIDSortedElemSet cornerFaces; // corner faces of found _BlockSide's
406     list< const SMDS_MeshNode* > corners( 1, nCorner );
407     list< const SMDS_MeshNode* >::iterator corner = corners.begin();
408     while ( corner != corners.end() )
409     {
410       SMDS_ElemIteratorPtr faceIt = (*corner)->GetInverseElementIterator( SMDSAbs_Face );
411       while ( faceIt->more() )
412       {
413         const SMDS_MeshElement* face = faceIt->next();
414         if ( !cornerFaces.insert( face ).second )
415           continue; // already loaded block side
416
417         if ( !isQuadrangle( face ))
418           return error("Non-quadrangle elements in the input mesh");
419
420         if ( _allSides.empty() || !_allSides.back()._grid.empty() )
421           _allSides.push_back( _BlockSide() );
422
423         _BlockSide& side = _allSides.back();
424         if ( !fillSide( side, face, *corner ) )
425         {
426           if ( !_error.empty() )
427             return false;
428         }
429         else
430         {
431           for ( int isXMax = 0; isXMax < 2; ++isXMax )
432             for ( int isYMax = 0; isYMax < 2; ++isYMax )
433             {
434               const SMDS_MeshNode* nCorner = side.getCornerNode(isXMax,isYMax );
435               corners.push_back( nCorner );
436               cornerFaces.insert( side.getCornerFace( nCorner ));
437             }
438           for ( int e = 0; e < int(NB_QUAD_SIDES); ++e )
439             _edge2sides[ side.getEdge( (EQuadEdge) e ) ].insert( &side );
440
441           nbFacesOnSides += side.getNbFaces();
442         }
443       }
444       ++corner;
445
446       // find block sides of other domains if any
447       if ( corner == corners.end() && nbFacesOnSides < mesh.NbQuadrangles() )
448       {
449         while ( nIt->more() )
450         {
451           nCorner = nIt->next();
452           if ( isCornerNode( nCorner ))
453             corner = corners.insert( corner, nCorner );
454         }
455         nbFacesOnSides = mesh.NbQuadrangles();
456       }
457     }
458     
459     if ( _allSides.empty() )
460       return BAD_MESH_ERR;
461     if ( _allSides.back()._grid.empty() )
462       _allSides.pop_back();
463     _DUMP_("Nb detected sides "<< _allSides.size());
464
465     // ---------------------------
466     // Organize sides into blocks
467     // ---------------------------
468
469     // analyse sharing of sides by blocks and sort sides by nb of adjacent sides
470     int nbBlockSides = 0; // total nb of block sides taking into account their sharing
471     multimap<int, _BlockSide* > sortedSides;
472     {
473       list < _BlockSide >::iterator sideIt = _allSides.begin();
474       for ( ; sideIt != _allSides.end(); ++sideIt )
475       {
476         _BlockSide& side = *sideIt;
477         bool isSharedSide = true;
478         int nbAdjacent = 0;
479         for ( int e = 0; e < int(NB_QUAD_SIDES) && isSharedSide; ++e )
480         {
481           int nbAdj = _edge2sides[ side.getEdge( (EQuadEdge) e ) ].size();
482           nbAdjacent += nbAdj;
483           isSharedSide = ( nbAdj > 2 );
484         }
485         side._nbBlocksFound = 0;
486         side._nbBlocksExpected = isSharedSide ? 2 : 1;
487         nbBlockSides += side._nbBlocksExpected;
488         sortedSides.insert( make_pair( nbAdjacent, & side ));
489       }
490     }
491
492     // find sides of each block
493     int nbBlocks = 0;
494     while ( nbBlockSides >= 6 )
495     {
496       // get any side not bound to all blocks it belongs to
497       multimap<int, _BlockSide*>::iterator i_side = sortedSides.begin();
498       while ( i_side != sortedSides.end() && i_side->second->isBound())
499         ++i_side;
500
501       // start searching for block sides from the got side
502       bool ok = true;
503       if ( _blocks.empty() || _blocks.back()._side[B_FRONT] )
504         _blocks.resize( _blocks.size() + 1 );
505
506       _Block& block = _blocks.back();
507       block.setSide( B_FRONT, i_side->second );
508       setSideBoundToBlock( *i_side->second );
509       nbBlockSides--;
510
511       // edges of adjacent sides of B_FRONT corresponding to front's edges
512       EQuadEdge edgeOfFront[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
513       EQuadEdge edgeOfAdj  [4] = { Q_BOTTOM, Q_LEFT, Q_BOTTOM, Q_LEFT };
514       // first find all sides detectable w/o advanced analysis,
515       // then repeat the search, which then may pass without advanced analysis
516       set< _BlockSide* > sidesAround;
517       for ( int advAnalys = 0; advAnalys < 2; ++advAnalys )
518       {
519         // try to find 4 sides adjacent to a FRONT side
520         for ( int i = 0; (ok || !advAnalys) && i < NB_QUAD_SIDES; ++i )
521           if ( !block._side[i] )
522             ok = block.setSide( i, findBlockSide( B_FRONT, edgeOfFront[i], edgeOfAdj[i],
523                                                   advAnalys, sidesAround));
524         // try to find a BACK side by a TOP one
525         if ( ok || !advAnalys)
526           if ( !block._side[B_BACK] && block._side[B_TOP] )
527             ok = block.setSide( B_BACK, findBlockSide( B_TOP, Q_TOP, Q_TOP,
528                                                        advAnalys, sidesAround ));
529         if ( !advAnalys ) ok = true;
530       }
531       ok = block.isValid();
532       if ( ok )
533       {
534         // check if just found block is same as one of previously found blocks
535         bool isSame = false;
536         for ( int i = 1; i < _blocks.size() && !isSame; ++i )
537           isSame = ( block._corners == _blocks[i-1]._corners );
538         ok = !isSame;
539       }
540
541       // count the found sides
542       _DUMP_(endl << "** Block " << _blocks.size() << " valid: " << block.isValid());
543       for (int i = 0; i < NB_BLOCK_SIDES; ++i )
544       {
545         _DUMP_("\tSide "<< SBoxSides[i] <<" "<< block._side[ i ]._side);
546         if ( block._side[ i ] )
547         {
548           if ( ok && i != B_FRONT)
549           {
550             setSideBoundToBlock( *block._side[ i ]._side );
551             nbBlockSides--;
552           }
553           _DUMP_("\t corners "<<
554                  block._side[ i ].cornerNode(0,0)->GetID() << ", " <<
555                  block._side[ i ].cornerNode(1,0)->GetID() << ", " <<
556                  block._side[ i ].cornerNode(1,1)->GetID() << ", " <<
557                  block._side[ i ].cornerNode(0,1)->GetID() << ", "<<endl);
558         }
559         else
560         {
561           _DUMP_("\t not found"<<endl);
562         }
563       }
564       if ( !ok )
565         block.clear();
566       else
567         nbBlocks++;
568     }
569     _DUMP_("Nb found blocks "<< nbBlocks <<endl);
570
571     if ( nbBlocks == 0 && _error.empty() )
572       return BAD_MESH_ERR;
573
574     return nbBlocks;
575   }
576
577   //================================================================================
578   /*!
579    * \brief Fill block side data starting from its corner quadrangle
580    */
581   //================================================================================
582
583   bool _Skin::fillSide( _BlockSide&             side,
584                         const SMDS_MeshElement* cornerQuad,
585                         const SMDS_MeshNode*    nCorner)
586   {
587     // Find out size of block side mesured in nodes and by the way find two rows
588     // of nodes in two directions.
589
590     int x, y, nbX, nbY;
591     const SMDS_MeshElement* firstQuad = cornerQuad;
592     {
593       // get a node on block edge
594       int iCorner = firstQuad->GetNodeIndex( nCorner );
595       const SMDS_MeshNode* nOnEdge = firstQuad->GetNode( (iCorner+1) % 4);
596
597       // find out size of block side
598       vector<const SMDS_MeshNode*> horRow1, horRow2, verRow1, verRow2;
599       if ( !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, horRow1, horRow2, true ) ||
600            !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, verRow1, verRow2, false ))
601         return false;
602       nbX = horRow1.size(), nbY = verRow1.size();
603
604       // store found nodes
605       side._index._xSize = horRow1.size();
606       side._index._ySize = verRow1.size();
607       side._grid.resize( side._index.size(), NULL );
608
609       for ( x = 0; x < horRow1.size(); ++x )
610       {
611         side.setNode( x, 0, horRow1[x] );
612         side.setNode( x, 1, horRow2[x] );
613       }
614       for ( y = 0; y < verRow1.size(); ++y )
615       {
616         side.setNode( 0, y, verRow1[y] );
617         side.setNode( 1, y, verRow2[y] );
618       }
619     }
620     // Find the rest nodes
621
622     y = 1; // y of the row to fill
623     TIDSortedElemSet emptySet, avoidSet;
624     while ( ++y < nbY )
625     {
626       // get next firstQuad in the next row of quadrangles
627       //
628       //          n2up
629       //     o---o               <- y row
630       //     |   |
631       //     o---o  o  o  o  o   <- found nodes
632       //n1down    n2down       
633       //
634       int i1down, i2down, i2up;
635       const SMDS_MeshNode* n1down = side.getNode( 0, y-1 );
636       const SMDS_MeshNode* n2down = side.getNode( 1, y-1 );
637       avoidSet.clear(); avoidSet.insert( firstQuad );
638       firstQuad = SMESH_MeshEditor::FindFaceInSet( n1down, n2down, emptySet, avoidSet,
639                                                    &i1down, &i2down);
640       if ( !isQuadrangle( firstQuad ))
641         return BAD_MESH_ERR;
642
643       const SMDS_MeshNode* n2up = oppositeNode( firstQuad, i1down );
644       avoidSet.clear(); avoidSet.insert( firstQuad );
645
646       // find the rest nodes in the y-th row by faces in the row
647
648       x = 1; 
649       while ( ++x < nbX )
650       {
651         const SMDS_MeshElement* quad = SMESH_MeshEditor::FindFaceInSet( n2up, n2down, emptySet,
652                                                                         avoidSet, &i2up, &i2down);
653         if ( !isQuadrangle( quad ))
654           return BAD_MESH_ERR;
655
656         n2up   = oppositeNode( quad, i2down );
657         n2down = oppositeNode( quad, i2up );
658         avoidSet.clear(); avoidSet.insert( quad );
659
660         side.setNode( x, y, n2up );
661       }
662     }
663
664     // check side validity
665     bool ok =
666       side.getCornerFace( side.getCornerNode( 0, 0 )) &&
667       side.getCornerFace( side.getCornerNode( 1, 0 )) &&
668       side.getCornerFace( side.getCornerNode( 0, 1 )) &&
669       side.getCornerFace( side.getCornerNode( 1, 1 ));
670
671     return ok;
672   }
673
674   //================================================================================
675   /*!
676    * \brief Return true if it's possible to make a loop over corner2Sides starting
677    * from the startSide
678    */
679   //================================================================================
680
681   bool isClosedChainOfSides( _BlockSide*                                        startSide,
682                              map< const SMDS_MeshNode*, list< _BlockSide* > > & corner2Sides )
683   {
684     // get start and end nodes
685     const SMDS_MeshNode *n1 = 0, *n2 = 0, *n;
686     for ( int y = 0; y < 2; ++y )
687       for ( int x = 0; x < 2; ++x )
688       {
689         n = startSide->getCornerNode(x,y);
690         if ( !corner2Sides.count( n )) continue;
691         if ( n1 )
692           n2 = n;
693         else
694           n1 = n;
695       }
696     if ( !n2 ) return false;
697
698     map< const SMDS_MeshNode*, list< _BlockSide* > >::iterator
699       c2sides = corner2Sides.find( n1 );
700     if ( c2sides == corner2Sides.end() ) return false;
701
702     int nbChainLinks = 1;
703     n = n1;
704     _BlockSide* prevSide = startSide;
705     while ( n != n2 )
706     {
707       // get the next side sharing n
708       list< _BlockSide* > & sides = c2sides->second;
709       _BlockSide* nextSide = ( sides.back() == prevSide ? sides.front() : sides.back() );
710       if ( nextSide == prevSide ) return false;
711
712       // find the next corner of the nextSide being in corner2Sides
713       n1 = n;
714       n = 0;
715       for ( int y = 0; y < 2 && !n; ++y )
716         for ( int x = 0; x < 2; ++x )
717         {
718           n = nextSide->getCornerNode(x,y);
719           c2sides = corner2Sides.find( n );
720           if ( n == n1 || c2sides == corner2Sides.end() )
721             n = 0;
722           else
723             break;
724         }
725       if ( !n ) return false;
726
727       prevSide = nextSide;
728       nbChainLinks++;
729     }
730
731     return ( n == n2 && nbChainLinks == NB_QUAD_SIDES );
732   }
733
734   //================================================================================
735   /*!
736    * \brief Try to find a block side adjacent to the given side by given edge
737    */
738   //================================================================================
739
740   _OrientedBlockSide _Skin::findBlockSide( EBoxSides           startBlockSide,
741                                            EQuadEdge           sharedSideEdge1,
742                                            EQuadEdge           sharedSideEdge2,
743                                            bool                withGeometricAnalysis,
744                                            set< _BlockSide* >& sidesAround)
745   {
746     _Block& block = _blocks.back();
747     _OrientedBlockSide& side1 = block._side[ startBlockSide ];
748
749     // get corner nodes of the given block edge
750     SMESH_OrientedLink edge = side1.edge( sharedSideEdge1 );
751     const SMDS_MeshNode* n1 = edge.node1();
752     const SMDS_MeshNode* n2 = edge.node2();
753     if ( edge._reversed ) swap( n1, n2 );
754
755     // find all sides sharing both nodes n1 and n2
756     set< _BlockSide* > sidesOnEdge = _edge2sides[ edge ]; // copy a set
757
758     // exclude loaded sides of block from sidesOnEdge
759     for (int i = 0; i < NB_BLOCK_SIDES; ++i )
760       if ( block._side[ i ] )
761         sidesOnEdge.erase( block._side[ i ]._side );
762
763     int nbSidesOnEdge = sidesOnEdge.size();
764     _DUMP_("nbSidesOnEdge "<< nbSidesOnEdge << " " << n1->GetID() << "-" << n2->GetID() );
765     if ( nbSidesOnEdge == 0 )
766       return 0;
767
768     _BlockSide* foundSide = 0;
769     if ( nbSidesOnEdge == 1 )
770     {
771       foundSide = *sidesOnEdge.begin();
772     }
773     else
774     {
775       set< _BlockSide* >::iterator sideIt = sidesOnEdge.begin();
776       int nbLoadedSides = block.nbSides();
777       if ( nbLoadedSides > 1 )
778       {
779         // Find the side having more than 2 corners common with already loaded sides
780         for (; !foundSide && sideIt != sidesOnEdge.end(); ++sideIt )
781         {
782           _BlockSide* sideI = *sideIt;
783           int nbCommonCorners =
784             block._corners.count( sideI->getCornerNode(0,0)) +
785             block._corners.count( sideI->getCornerNode(1,0)) +
786             block._corners.count( sideI->getCornerNode(0,1)) +
787             block._corners.count( sideI->getCornerNode(1,1));
788           if ( nbCommonCorners > 2 )
789             foundSide = sideI;
790         }
791       }
792
793       if ( !foundSide )
794       {
795         if ( !withGeometricAnalysis )
796         {
797           sidesAround.insert( sidesOnEdge.begin(), sidesOnEdge.end() );
798           return 0;
799         }
800         if ( nbLoadedSides == 1 )
801         {
802           // Issue 0021529. There are at least 2 sides by each edge and
803           // position of block gravity center is undefined.
804           // Find a side starting from which we can walk around the startBlockSide
805
806           // fill in corner2Sides
807           map< const SMDS_MeshNode*, list< _BlockSide* > > corner2Sides;
808           for ( sideIt = sidesAround.begin(); sideIt != sidesAround.end(); ++sideIt )
809           {
810             _BlockSide* sideI = *sideIt;
811             corner2Sides[ sideI->getCornerNode(0,0) ].push_back( sideI );
812             corner2Sides[ sideI->getCornerNode(1,0) ].push_back( sideI );
813             corner2Sides[ sideI->getCornerNode(0,1) ].push_back( sideI );
814             corner2Sides[ sideI->getCornerNode(1,1) ].push_back( sideI );
815           }
816           // remove corners of startBlockSide from corner2Sides
817           set<const SMDS_MeshNode*>::iterator nIt = block._corners.begin();
818           for ( ; nIt != block._corners.end(); ++nIt )
819             corner2Sides.erase( *nIt );
820
821           // select a side
822           for ( sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
823           {
824             if ( isClosedChainOfSides( *sideIt, corner2Sides ))
825             {
826               foundSide = *sideIt;
827               break;
828             }
829           }
830           if ( !foundSide )
831             return 0;
832         }
833         else
834         {
835           // Select one of found sides most close to startBlockSide
836
837           gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()),  p2 (n2->X(),n2->Y(),n2->Z());
838           gp_Vec p1p2( p1, p2 );
839
840           const SMDS_MeshElement* face1 = side1->getCornerFace( n1 );
841           gp_XYZ p1Op = SMESH_TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(n1)));
842           gp_Vec side1Dir( p1, p1Op );
843           gp_Ax2 pln( p1, p1p2, side1Dir ); // plane with normal p1p2 and X dir side1Dir
844           _DUMP_("  Select adjacent for "<< side1._side << " - side dir ("
845                  << side1Dir.X() << ", " << side1Dir.Y() << ", " << side1Dir.Z() << ")" );
846
847           map < double , _BlockSide* > angleOfSide;
848           for (sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
849           {
850             _BlockSide* sideI = *sideIt;
851             const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 );
852             gp_XYZ p1Op = SMESH_TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(n1)));
853             gp_Vec sideIDir( p1, p1Op );
854             // compute angle of (sideIDir projection to pln) and (X dir of pln)
855             gp_Vec2d sideIDirProj( sideIDir * pln.XDirection(), sideIDir * pln.YDirection());
856             double angle = sideIDirProj.Angle( gp::DX2d() );
857             if ( angle < 0 ) angle += 2. * M_PI; // angle [0-2*PI]
858             angleOfSide.insert( make_pair( angle, sideI ));
859             _DUMP_("  "<< sideI << " - side dir ("
860                    << sideIDir.X() << ", " << sideIDir.Y() << ", " << sideIDir.Z() << ")"
861                    << " angle " << angle);
862           }
863
864           gp_XYZ gc(0,0,0); // gravity center of already loaded block sides
865           for (int i = 0; i < NB_BLOCK_SIDES; ++i )
866             if ( block._side[ i ] )
867               gc += block._side[ i ]._side->getGC();
868           gc /= nbLoadedSides;
869
870           gp_Vec gcDir( p1, gc );
871           gp_Vec2d gcDirProj( gcDir * pln.XDirection(), gcDir * pln.YDirection());
872           double gcAngle = gcDirProj.Angle( gp::DX2d() );
873           foundSide = gcAngle < 0 ? angleOfSide.rbegin()->second : angleOfSide.begin()->second;
874         }
875       }
876       _DUMP_("  selected "<< foundSide );
877     }
878
879     // Orient the found side correctly
880
881     // corners of found side corresponding to nodes n1 and n2
882     bool xMax1, yMax1, xMax2, yMax2;
883     if ( !getEdgeEnds( sharedSideEdge2, xMax1, yMax1, xMax2, yMax2 ))
884       return error(SMESH_Comment("Internal error at ")<<__FILE__<<":"<<__LINE__),
885         _OrientedBlockSide(0);
886
887     for ( int ori = 0; ori < _OrientedIndexer::MAX_ORI+1; ++ori )
888     {
889       _OrientedBlockSide orientedSide( foundSide, ori );
890       const SMDS_MeshNode* n12 = orientedSide.cornerNode( xMax1, yMax1);
891       const SMDS_MeshNode* n22 = orientedSide.cornerNode( xMax2, yMax2);
892       if ( n1 == n12 && n2 == n22 )
893         return orientedSide;
894     }
895     error(SMESH_Comment("Failed to orient a block side found by edge ")<<sharedSideEdge1
896           << " of side " << startBlockSide
897           << " of block " << _blocks.size());
898     return 0;
899   }
900
901   //================================================================================
902   /*!
903    * \brief: Fill rows (which are actually columns,if !alongN1N2) of nodes starting
904    * from the given quadrangle until another block corner encounters.
905    *  n1 and n2 are at bottom of quad, n1 is at block corner.
906    */
907   //================================================================================
908
909   bool _Skin::fillRowsUntilCorner(const SMDS_MeshElement*       quad,
910                                   const SMDS_MeshNode*          n1,
911                                   const SMDS_MeshNode*          n2,
912                                   vector<const SMDS_MeshNode*>& row1,
913                                   vector<const SMDS_MeshNode*>& row2,
914                                   const bool                    alongN1N2 )
915   {
916     const SMDS_MeshNode* corner1 = n1;
917
918     // Store nodes of quad in the rows and find new n1 and n2 to get
919     // the next face so that new n2 is on block edge
920     int i1 = quad->GetNodeIndex( n1 );
921     int i2 = quad->GetNodeIndex( n2 );
922     row1.clear(); row2.clear();
923     row1.push_back( n1 );
924     if ( alongN1N2 )
925     {
926       row1.push_back( n2 );
927       row2.push_back( oppositeNode( quad, i2 ));
928       row2.push_back( n1 = oppositeNode( quad, i1 ));
929     }
930     else
931     {
932       row2.push_back( n2 );
933       row1.push_back( n2 = oppositeNode( quad, i2 ));
934       row2.push_back( n1 = oppositeNode( quad, i1 ));
935     }
936
937     if ( isCornerNode( row1[1] ))
938       return true;
939
940     // Find the rest nodes
941     TIDSortedElemSet emptySet, avoidSet;
942     while ( !isCornerNode( n2 ) )
943     {
944       avoidSet.clear(); avoidSet.insert( quad );
945       quad = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2 );
946       if ( !isQuadrangle( quad ))
947         return BAD_MESH_ERR;
948
949       row1.push_back( n2 = oppositeNode( quad, i1 ));
950       row2.push_back( n1 = oppositeNode( quad, i2 ));
951     }
952     return n1 != corner1;
953   }
954
955   //================================================================================
956   /*!
957    * \brief Return a corner face by a corner node
958    */
959   //================================================================================
960
961   const SMDS_MeshElement* _BlockSide::getCornerFace(const SMDS_MeshNode* cornerNode) const
962   {
963     int x, y, isXMax, isYMax, found = 0;
964     for ( isXMax = 0; isXMax < 2; ++isXMax )
965     {
966       for ( isYMax = 0; isYMax < 2; ++isYMax )
967       {
968         x = isXMax ? _index._xSize-1 : 0;
969         y = isYMax ? _index._ySize-1 : 0;
970         found = ( getNode(x,y) == cornerNode );
971         if ( found ) break;
972       }
973       if ( found ) break;
974     }
975     if ( !found ) return 0;
976     int dx = isXMax ? -1 : +1;
977     int dy = isYMax ? -1 : +1;
978     const SMDS_MeshNode* n1 = getNode(x,y);
979     const SMDS_MeshNode* n2 = getNode(x+dx,y);
980     const SMDS_MeshNode* n3 = getNode(x,y+dy);
981     const SMDS_MeshNode* n4 = getNode(x+dx,y+dy);
982     return SMDS_Mesh::FindFace(n1, n2, n3, n4 );
983   }
984
985   //================================================================================
986   /*!
987    * \brief Checks own validity
988    */
989   //================================================================================
990
991   bool _Block::isValid() const
992   {
993     bool ok = ( nbSides() == 6 );
994
995     // check only corners depending on side selection
996     EBoxSides adjacent[4] = { B_BOTTOM, B_RIGHT, B_TOP, B_LEFT };
997     EQuadEdge edgeAdj [4] = { Q_TOP,    Q_RIGHT, Q_TOP, Q_RIGHT };
998     EQuadEdge edgeBack[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
999
1000     for ( int i=0; ok && i < NB_QUAD_SIDES; ++i )
1001     { 
1002       SMESH_OrientedLink eBack = _side[ B_BACK      ].edge( edgeBack[i] );
1003       SMESH_OrientedLink eAdja = _side[ adjacent[i] ].edge( edgeAdj[i] );
1004       ok = ( eBack == eAdja );
1005     }
1006     return ok;
1007   }
1008
1009 } // namespace
1010
1011 //=======================================================================
1012 //function : StdMeshers_HexaFromSkin_3D
1013 //purpose  : 
1014 //=======================================================================
1015
1016 StdMeshers_HexaFromSkin_3D::StdMeshers_HexaFromSkin_3D(int hypId, int studyId, SMESH_Gen* gen)
1017   :SMESH_3D_Algo(hypId, studyId, gen)
1018 {
1019   MESSAGE("StdMeshers_HexaFromSkin_3D::StdMeshers_HexaFromSkin_3D");
1020   _name = "HexaFromSkin_3D";
1021 }
1022
1023 StdMeshers_HexaFromSkin_3D::~StdMeshers_HexaFromSkin_3D()
1024 {
1025   MESSAGE("StdMeshers_HexaFromSkin_3D::~StdMeshers_HexaFromSkin_3D");
1026 }
1027
1028 //================================================================================
1029 /*!
1030  * \brief Main method, which generates hexaheda
1031  */
1032 //================================================================================
1033
1034 bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
1035 {
1036   _Skin skin;
1037   int nbBlocks = skin.findBlocks(aMesh);
1038   if ( nbBlocks == 0 )
1039     return error( skin.error());
1040
1041   vector< vector< const SMDS_MeshNode* > > columns;
1042   int x, xSize, y, ySize, z, zSize;
1043   _Indexer colIndex;
1044
1045   for ( int i = 0; i < nbBlocks; ++i )
1046   {
1047     const _Block& block = skin.getBlock( i );
1048
1049     // ------------------------------------------
1050     // Fill columns of nodes with existing nodes
1051     // ------------------------------------------
1052
1053     xSize = block.getSide(B_BOTTOM).getHoriSize();
1054     ySize = block.getSide(B_BOTTOM).getVertSize();
1055     zSize = block.getSide(B_FRONT ).getVertSize();
1056     int X = xSize - 1, Y = ySize - 1, Z = zSize - 1;
1057     colIndex = _Indexer( xSize, ySize );
1058     columns.resize( colIndex.size() );
1059
1060     // fill node columns by front and back box sides
1061     for ( x = 0; x < xSize; ++x ) {
1062       vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
1063       vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
1064       column0.resize( zSize );
1065       column1.resize( zSize );
1066       for ( z = 0; z < zSize; ++z ) {
1067         column0[ z ] = block.getSide(B_FRONT).node( x, z );
1068         column1[ z ] = block.getSide(B_BACK) .node( x, z );
1069       }
1070     }
1071     // fill node columns by left and right box sides
1072     for ( y = 1; y < ySize-1; ++y ) {
1073       vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
1074       vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
1075       column0.resize( zSize );
1076       column1.resize( zSize );
1077       for ( z = 0; z < zSize; ++z ) {
1078         column0[ z ] = block.getSide(B_LEFT) .node( y, z );
1079         column1[ z ] = block.getSide(B_RIGHT).node( y, z );
1080       }
1081     }
1082     // get nodes from top and bottom box sides
1083     for ( x = 1; x < xSize-1; ++x ) {
1084       for ( y = 1; y < ySize-1; ++y ) {
1085         vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
1086         column.resize( zSize );
1087         column.front() = block.getSide(B_BOTTOM).node( x, y );
1088         column.back()  = block.getSide(B_TOP)   .node( x, y );
1089       }
1090     }
1091
1092     // ----------------------------
1093     // Add internal nodes of a box
1094     // ----------------------------
1095     // projection points of internal nodes on box sub-shapes by which
1096     // coordinates of internal nodes are computed
1097     vector<gp_XYZ> pointOnShape( SMESH_Block::ID_Shell );
1098
1099     // projections on vertices are constant
1100     pointOnShape[ SMESH_Block::ID_V000 ] = block.getSide(B_BOTTOM).xyz( 0, 0 );
1101     pointOnShape[ SMESH_Block::ID_V100 ] = block.getSide(B_BOTTOM).xyz( X, 0 );
1102     pointOnShape[ SMESH_Block::ID_V010 ] = block.getSide(B_BOTTOM).xyz( 0, Y );
1103     pointOnShape[ SMESH_Block::ID_V110 ] = block.getSide(B_BOTTOM).xyz( X, Y );
1104     pointOnShape[ SMESH_Block::ID_V001 ] = block.getSide(B_TOP).xyz( 0, 0 );
1105     pointOnShape[ SMESH_Block::ID_V101 ] = block.getSide(B_TOP).xyz( X, 0 );
1106     pointOnShape[ SMESH_Block::ID_V011 ] = block.getSide(B_TOP).xyz( 0, Y );
1107     pointOnShape[ SMESH_Block::ID_V111 ] = block.getSide(B_TOP).xyz( X, Y );
1108
1109     for ( x = 1; x < xSize-1; ++x )
1110     {
1111       gp_XYZ params; // normalized parameters of internal node within a unit box
1112       params.SetCoord( 1, x / double(X) );
1113       for ( y = 1; y < ySize-1; ++y )
1114       {
1115         params.SetCoord( 2, y / double(Y) );
1116         // column to fill during z loop
1117         vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
1118         // projections on horizontal edges
1119         pointOnShape[ SMESH_Block::ID_Ex00 ] = block.getSide(B_BOTTOM).xyz( x, 0 );
1120         pointOnShape[ SMESH_Block::ID_Ex10 ] = block.getSide(B_BOTTOM).xyz( x, Y );
1121         pointOnShape[ SMESH_Block::ID_E0y0 ] = block.getSide(B_BOTTOM).xyz( 0, y );
1122         pointOnShape[ SMESH_Block::ID_E1y0 ] = block.getSide(B_BOTTOM).xyz( X, y );
1123         pointOnShape[ SMESH_Block::ID_Ex01 ] = block.getSide(B_TOP).xyz( x, 0 );
1124         pointOnShape[ SMESH_Block::ID_Ex11 ] = block.getSide(B_TOP).xyz( x, Y );
1125         pointOnShape[ SMESH_Block::ID_E0y1 ] = block.getSide(B_TOP).xyz( 0, y );
1126         pointOnShape[ SMESH_Block::ID_E1y1 ] = block.getSide(B_TOP).xyz( X, y );
1127         // projections on horizontal sides
1128         pointOnShape[ SMESH_Block::ID_Fxy0 ] = block.getSide(B_BOTTOM).xyz( x, y );
1129         pointOnShape[ SMESH_Block::ID_Fxy1 ] = block.getSide(B_TOP)   .xyz( x, y );
1130         for ( z = 1; z < zSize-1; ++z ) // z loop
1131         {
1132           params.SetCoord( 3, z / double(Z) );
1133           // projections on vertical edges
1134           pointOnShape[ SMESH_Block::ID_E00z ] = block.getSide(B_FRONT).xyz( 0, z );    
1135           pointOnShape[ SMESH_Block::ID_E10z ] = block.getSide(B_FRONT).xyz( X, z );    
1136           pointOnShape[ SMESH_Block::ID_E01z ] = block.getSide(B_BACK).xyz( 0, z );    
1137           pointOnShape[ SMESH_Block::ID_E11z ] = block.getSide(B_BACK).xyz( X, z );
1138           // projections on vertical sides
1139           pointOnShape[ SMESH_Block::ID_Fx0z ] = block.getSide(B_FRONT).xyz( x, z );    
1140           pointOnShape[ SMESH_Block::ID_Fx1z ] = block.getSide(B_BACK) .xyz( x, z );    
1141           pointOnShape[ SMESH_Block::ID_F0yz ] = block.getSide(B_LEFT) .xyz( y, z );    
1142           pointOnShape[ SMESH_Block::ID_F1yz ] = block.getSide(B_RIGHT).xyz( y, z );
1143
1144           // compute internal node coordinates
1145           gp_XYZ coords;
1146           SMESH_Block::ShellPoint( params, pointOnShape, coords );
1147           column[ z ] = aHelper->AddNode( coords.X(), coords.Y(), coords.Z() );
1148
1149 #ifdef DEB_GRID
1150           // debug
1151           //cout << "----------------------------------------------------------------------"<<endl;
1152           //for ( int id = SMESH_Block::ID_V000; id < SMESH_Block::ID_Shell; ++id)
1153           //{
1154           //  gp_XYZ p = pointOnShape[ id ];
1155           //  SMESH_Block::DumpShapeID( id,cout)<<" ( "<<p.X()<<", "<<p.Y()<<", "<<p.Z()<<" )"<<endl;
1156           //}
1157           //cout << "Params: ( "<< params.X()<<", "<<params.Y()<<", "<<params.Z()<<" )"<<endl;
1158           //cout << "coords: ( "<< coords.X()<<", "<<coords.Y()<<", "<<coords.Z()<<" )"<<endl;
1159 #endif
1160         }
1161       }
1162     }
1163     // ----------------
1164     // Add hexahedrons
1165     // ----------------
1166
1167     // find out orientation by a least distorted hexahedron (issue 0020855);
1168     // the last is defined by evaluating sum of face normals of 8 corner hexahedrons
1169     double badness = numeric_limits<double>::max();
1170     bool isForw = true;
1171     for ( int xMax = 0; xMax < 2; ++xMax )
1172       for ( int yMax = 0; yMax < 2; ++yMax )
1173         for ( int zMax = 0; zMax < 2; ++zMax )
1174         {
1175           x = xMax ? xSize-1 : 1;
1176           y = yMax ? ySize-1 : 1;
1177           z = zMax ? zSize-1 : 1;
1178           vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x-1, y-1 )];
1179           vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x  , y-1 )];
1180           vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x-1, y   )];
1181           vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x  , y )];
1182           
1183           const SMDS_MeshNode* n000 = col00[z-1];
1184           const SMDS_MeshNode* n100 = col10[z-1];
1185           const SMDS_MeshNode* n010 = col01[z-1];
1186           const SMDS_MeshNode* n110 = col11[z-1];
1187           const SMDS_MeshNode* n001 = col00[z];
1188           const SMDS_MeshNode* n101 = col10[z];
1189           const SMDS_MeshNode* n011 = col01[z];
1190           const SMDS_MeshNode* n111 = col11[z];
1191           SMDS_VolumeOfNodes probeVolume (n000,n010,n110,n100,
1192                                           n001,n011,n111,n101);
1193           SMDS_VolumeTool volTool( &probeVolume );
1194           double Nx=0.,Ny=0.,Nz=0.;
1195           for ( int iFace = 0; iFace < volTool.NbFaces(); ++iFace )
1196           {
1197             double nx,ny,nz;
1198             volTool.GetFaceNormal( iFace, nx,ny,nz );
1199             Nx += nx;
1200             Ny += ny;
1201             Nz += nz;
1202           }
1203           double quality = Nx*Nx + Ny*Ny + Nz*Nz;
1204           if ( quality < badness )
1205           {
1206             badness = quality;
1207             isForw = volTool.IsForward();
1208           }
1209         }
1210
1211     // add elements
1212     for ( x = 0; x < xSize-1; ++x ) {
1213       for ( y = 0; y < ySize-1; ++y ) {
1214         vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
1215         vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
1216         vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
1217         vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
1218         // bottom face normal of a hexa mush point outside the volume
1219         if ( isForw )
1220           for ( z = 0; z < zSize-1; ++z )
1221             aHelper->AddVolume(col00[z],   col01[z],   col11[z],   col10[z],
1222                                col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
1223         else
1224           for ( z = 0; z < zSize-1; ++z )
1225             aHelper->AddVolume(col00[z],   col10[z],   col11[z],   col01[z],
1226                                col00[z+1], col10[z+1], col11[z+1], col01[z+1]);
1227       }
1228     }
1229   } // loop on blocks
1230
1231   return true;
1232 }
1233
1234 //================================================================================
1235 /*!
1236  * \brief Evaluate nb of hexa
1237  */
1238 //================================================================================
1239
1240 bool StdMeshers_HexaFromSkin_3D::Evaluate(SMESH_Mesh &         aMesh,
1241                                           const TopoDS_Shape & aShape,
1242                                           MapShapeNbElems&     aResMap)
1243 {
1244   _Skin skin;
1245   int nbBlocks = skin.findBlocks(aMesh);
1246   if ( nbBlocks == 0 )
1247     return error( skin.error());
1248
1249   bool secondOrder = aMesh.NbFaces( ORDER_QUADRATIC );
1250
1251   int entity = secondOrder ? SMDSEntity_Quad_Hexa : SMDSEntity_Hexa;
1252   vector<int>& nbByType = aResMap[ aMesh.GetSubMesh( aShape )];
1253   if ( entity >= nbByType.size() )
1254     nbByType.resize( SMDSEntity_Last, 0 );
1255
1256   for ( int i = 0; i < nbBlocks; ++i )
1257   {
1258     const _Block& block = skin.getBlock( i );
1259
1260     int nbX = block.getSide(B_BOTTOM).getHoriSize();
1261     int nbY = block.getSide(B_BOTTOM).getVertSize();
1262     int nbZ = block.getSide(B_FRONT ).getVertSize();
1263
1264     int nbHexa  = (nbX-1) * (nbY-1) * (nbZ-1);
1265     int nbNodes = (nbX-2) * (nbY-2) * (nbZ-2);
1266     if ( secondOrder )
1267       nbNodes +=
1268         (nbX-2) * (nbY-2) * (nbZ-1) +
1269         (nbX-2) * (nbY-1) * (nbZ-2) +
1270         (nbX-1) * (nbY-2) * (nbZ-2);
1271
1272
1273     nbByType[ entity ] += nbHexa;
1274     nbByType[ SMDSEntity_Node ] += nbNodes;
1275   }
1276
1277   return true;
1278 }
1279
1280 //================================================================================
1281 /*!
1282  * \brief Abstract method must be defined but does nothing
1283  */
1284 //================================================================================
1285
1286 bool StdMeshers_HexaFromSkin_3D::CheckHypothesis(SMESH_Mesh&, const TopoDS_Shape&,
1287                                                  Hypothesis_Status& aStatus)
1288 {
1289   aStatus = SMESH_Hypothesis::HYP_OK;
1290   return true;
1291 }
1292
1293 //================================================================================
1294 /*!
1295  * \brief Abstract method must be defined but just reports an error as this
1296  *  algo is not intended to work with shapes
1297  */
1298 //================================================================================
1299
1300 bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh&, const TopoDS_Shape&)
1301 {
1302   return error("Algorithm can't work with geometrical shapes");
1303 }