1 // Copyright (C) 2009-2016 CEA/DEN, EDF R&D
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, or (at your option) any later version.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
20 // File : SMESH_HexaFromSkin_3D.cxx
21 // Created : Wed Jan 27 12:28:07 2010
22 // Author : Edward AGAPOV (eap)
24 #include "HEXABLOCKPlugin_FromSkin_3D.hxx"
26 #include "SMDS_VolumeOfNodes.hxx"
27 #include "SMDS_VolumeTool.hxx"
28 #include "SMESH_Block.hxx"
29 #include "SMESH_MesherHelper.hxx"
30 #include "SMESH_MeshEditor.hxx"
31 #include "SMESH_MeshAlgos.hxx"
35 //#include "utilities.h"
38 // Define error message and _MYDEBUG_ if needed
40 #define BAD_MESH_ERR \
41 error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh.\n" \
42 __FILE__ ":" )<<__LINE__)
45 #define BAD_MESH_ERR \
46 error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh"))
51 #define _DUMP_(msg) cout << msg << endl
58 static int MYDEBUG = HEXA_NS::on_debug ();
60 static int MYDEBUG = 0;
66 enum EBoxSides //!< sides of the block
68 B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, NB_BLOCK_SIDES
71 const char* SBoxSides[] = //!< names of block sides
73 "BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED"
76 enum EQuadEdge //!< edges of quadrangle side
78 Q_BOTTOM = 0, Q_RIGHT, Q_TOP, Q_LEFT, NB_QUAD_SIDES
82 //================================================================================
84 * \brief return logical coordinates (i.e. min or max) of ends of edge
86 //================================================================================
88 bool getEdgeEnds(EQuadEdge edge, bool& xMax1, bool& yMax1, bool& xMax2, bool& yMax2 )
90 xMax1=0, yMax1=0, xMax2=1, yMax2=1;
93 case Q_BOTTOM: yMax2 = 0; break;
94 case Q_RIGHT: xMax1 = 1; break;
95 case Q_TOP: yMax1 = 1; break;
96 case Q_LEFT: xMax2 = 0; break;
103 //================================================================================
105 * \brief return true if a node is at block corner
107 * This check is valid for simple cases only
109 //================================================================================
111 bool isCornerNode( const SMDS_MeshNode* n )
113 int nbF = n ? n->NbInverseElements( SMDSAbs_Face ) : 1;
117 set<const SMDS_MeshNode*> nodesInInverseFaces;
118 SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face );
119 while ( fIt->more() )
121 const SMDS_MeshElement* face = fIt->next();
122 nodesInInverseFaces.insert( face->begin_nodes(), face->end_nodes() );
125 return nodesInInverseFaces.size() != ( 6 + (nbF/2-1)*3 );
128 //================================================================================
130 * \brief check element type
132 //================================================================================
134 bool isQuadrangle(const SMDS_MeshElement* e)
136 return ( e && e->NbCornerNodes() == 4 );
139 //================================================================================
141 * \brief return opposite node of a quadrangle face
143 //================================================================================
145 const SMDS_MeshNode* oppositeNode(const SMDS_MeshElement* quad, int iNode)
147 return quad->GetNode( (iNode+2) % 4 );
150 //================================================================================
152 * \brief Convertor of a pair of integers to a sole index
157 _Indexer( int xSize=0, int ySize=0 ): _xSize(xSize), _ySize(ySize) {}
158 int size() const { return _xSize * _ySize; }
159 int operator()(int x, int y) const { return y * _xSize + x; }
161 //================================================================================
163 * \brief Oriented convertor of a pair of integers to a sole index
165 class _OrientedIndexer : public _Indexer
168 enum OriFlags //!< types of block side orientation
170 REV_X = 1, REV_Y = 2, SWAP_XY = 4, MAX_ORI = REV_X|REV_Y|SWAP_XY
172 _OrientedIndexer( const _Indexer& indexer, const int oriFlags ):
173 _Indexer( indexer._xSize, indexer._ySize ),
174 _xSize (indexer._xSize), _ySize(indexer._ySize),
175 _xRevFun((oriFlags & REV_X) ? & reverse : & lazy),
176 _yRevFun((oriFlags & REV_Y) ? & reverse : & lazy),
177 _swapFun((oriFlags & SWAP_XY ) ? & swap : & lazy)
179 (*_swapFun)( _xSize, _ySize );
181 //!< Return index by XY
182 int operator()(int x, int y) const
184 (*_xRevFun)( x, const_cast<int&>( _xSize ));
185 (*_yRevFun)( y, const_cast<int&>( _ySize ));
187 return _Indexer::operator()(x,y);
189 //!< Return index for a corner
190 int corner(bool xMax, bool yMax) const
192 int x = xMax, y = yMax, size = 2;
193 (*_xRevFun)( x, size );
194 (*_yRevFun)( y, size );
196 return _Indexer::operator()(x ? _Indexer::_xSize-1 : 0 , y ? _Indexer::_ySize-1 : 0);
198 int xSize() const { return _xSize; }
199 int ySize() const { return _ySize; }
204 typedef void (*TFun)(int& x, int& y);
205 TFun _xRevFun, _yRevFun, _swapFun;
207 static void lazy (int&, int&) {}
208 static void reverse(int& x, int& size) { x = size - x - 1; }
209 static void swap (int& x, int& y) { std::swap(x,y); }
211 //================================================================================
213 * \brief Structure corresponding to the meshed side of block
217 vector<const SMDS_MeshNode*> _grid;
219 int _nbBlocksExpected;
222 #ifdef _DEBUG_ // want to get SIGSEGV in case of invalid index
223 #define _grid_access_(pobj, i) pobj->_grid[ ((i) < pobj->_grid.size()) ? i : int(1e100)]
225 #define _grid_access_(pobj, i) pobj->_grid[ i ]
227 //!< Return node at XY
228 const SMDS_MeshNode* getNode(int x, int y) const { return _grid_access_(this, _index( x,y ));}
230 void setNode(int x, int y, const SMDS_MeshNode* n) { _grid_access_(this, _index( x,y )) = n; }
232 SMESH_OrientedLink getEdge(EQuadEdge edge) const
234 bool x1, y1, x2, y2; getEdgeEnds( edge, x1, y1, x2, y2 );
235 return SMESH_OrientedLink( getCornerNode ( x1, y1 ), getCornerNode( x2, y2 ));
237 //!< Return a corner node
238 const SMDS_MeshNode* getCornerNode(bool isXMax, bool isYMax) const
240 return getNode( isXMax ? _index._xSize-1 : 0 , isYMax ? _index._ySize-1 : 0 );
242 const SMDS_MeshElement* getCornerFace(const SMDS_MeshNode* cornerNode) const;
243 //!< True if all blocks this side belongs to have been found
244 bool isBound() const { return _nbBlocksExpected <= _nbBlocksFound; }
245 //!< Return coordinates of node at XY
246 gp_XYZ getXYZ(int x, int y) const { return SMESH_TNodeXYZ( getNode( x, y )); }
247 //!< Return gravity center of the four corners and the middle node
252 getXYZ( _index._xSize-1, 0 ) +
253 getXYZ( 0, _index._ySize-1 ) +
254 getXYZ( _index._xSize-1, _index._ySize-1 ) +
255 getXYZ( _index._xSize/2, _index._ySize/2 );
258 //!< Return number of mesh faces
259 int getNbFaces() const { return (_index._xSize-1) * (_index._ySize-1); }
261 //================================================================================
263 * \brief _BlockSide with changed orientation
265 struct _OrientedBlockSide
268 _OrientedIndexer _index;
270 _OrientedBlockSide( _BlockSide* side=0, const int oriFlags=0 ):
271 _side(side), _index(side ? side->_index : _Indexer(), oriFlags ) {}
272 //!< return coordinates by XY
273 gp_XYZ xyz(int x, int y) const
275 return SMESH_TNodeXYZ( _grid_access_(_side, _index( x, y )) );
277 //!< safely return a node by XY
278 const SMDS_MeshNode* node(int x, int y) const
280 int i = _index( x, y );
281 return ( i < 0 || i >= _side->_grid.size()) ? 0 : _side->_grid[i];
284 SMESH_OrientedLink edge(EQuadEdge edge) const
286 bool x1, y1, x2, y2; getEdgeEnds( edge, x1, y1, x2, y2 );
287 return SMESH_OrientedLink( cornerNode ( x1, y1 ), cornerNode( x2, y2 ));
289 //!< Return a corner node
290 const SMDS_MeshNode* cornerNode(bool isXMax, bool isYMax) const
292 return _grid_access_(_side, _index.corner( isXMax, isYMax ));
294 //!< return its size in nodes
295 int getHoriSize() const { return _index.xSize(); }
296 int getVertSize() const { return _index.ySize(); }
297 //!< True if _side has been initialized
298 operator bool() const { return _side; }
299 //! Direct access to _side
300 const _BlockSide* operator->() const { return _side; }
301 _BlockSide* operator->() { return _side; }
303 //================================================================================
305 * \brief Meshed skin of block
309 _OrientedBlockSide _side[6]; // 6 sides of a sub-block
310 set<const SMDS_MeshNode*> _corners;
312 const _OrientedBlockSide& getSide(int i) const { return _side[i]; }
313 bool setSide( int i, const _OrientedBlockSide& s)
315 if (( _side[i] = s ))
317 _corners.insert( s.cornerNode(0,0));
318 _corners.insert( s.cornerNode(1,0));
319 _corners.insert( s.cornerNode(0,1));
320 _corners.insert( s.cornerNode(1,1));
324 void clear() { for (int i=0;i<6;++i) _side[i]=0; _corners.clear(); }
325 bool hasSide( const _OrientedBlockSide& s) const
327 if ( s ) for (int i=0;i<6;++i) if ( _side[i] && _side[i]._side == s._side ) return true;
330 int nbSides() const { int n=0; for (int i=0;i<6;++i) if ( _side[i] ) ++n; return n; }
331 bool isValid() const;
333 //================================================================================
335 * \brief Skin mesh possibly containing several meshed blocks
341 int findBlocks(SMESH_Mesh& mesh);
342 //!< return i-th block
343 const _Block& getBlock(int i) const { return _blocks[i]; }
344 //!< return error description
345 const SMESH_Comment& error() const { return _error; }
348 bool fillSide( _BlockSide& side,
349 const SMDS_MeshElement* cornerQuad,
350 const SMDS_MeshNode* cornerNode);
351 bool fillRowsUntilCorner(const SMDS_MeshElement* quad,
352 const SMDS_MeshNode* n1,
353 const SMDS_MeshNode* n2,
354 vector<const SMDS_MeshNode*>& verRow1,
355 vector<const SMDS_MeshNode*>& verRow2,
357 _OrientedBlockSide findBlockSide( EBoxSides startBlockSide,
358 EQuadEdge sharedSideEdge1,
359 EQuadEdge sharedSideEdge2,
360 bool withGeometricAnalysis,
361 set< _BlockSide* >& sidesAround);
362 //!< update own data and data of the side bound to block
363 void setSideBoundToBlock( _BlockSide& side )
365 if ( side._nbBlocksFound++, side.isBound() )
366 for ( int e = 0; e < int(NB_QUAD_SIDES); ++e )
367 _edge2sides[ side.getEdge( (EQuadEdge) e ) ].erase( &side );
369 //!< store reason of error
370 int error(const SMESH_Comment& reason) { _error = reason; return 0; }
372 SMESH_Comment _error;
374 list< _BlockSide > _allSides;
375 vector< _Block > _blocks;
377 //map< const SMDS_MeshNode*, set< _BlockSide* > > _corner2sides;
378 map< SMESH_OrientedLink, set< _BlockSide* > > _edge2sides;
381 //================================================================================
383 * \brief Find and return number of submeshes corresponding to blocks
385 //================================================================================
387 int _Skin::findBlocks(SMESH_Mesh& mesh)
389 SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
391 // Find a node at any block corner
393 SMDS_NodeIteratorPtr nIt = meshDS->nodesIterator(/*idInceasingOrder=*/true);
394 if ( !nIt->more() ) return error("Empty mesh");
396 const SMDS_MeshNode* nCorner = 0;
397 while ( nIt->more() )
399 nCorner = nIt->next();
400 if ( isCornerNode( nCorner ))
408 // --------------------------------------------------------------------
409 // Find all block sides starting from mesh faces sharing the corner node
410 // --------------------------------------------------------------------
412 int nbFacesOnSides = 0;
413 TIDSortedElemSet cornerFaces; // corner faces of found _BlockSide's
414 list< const SMDS_MeshNode* > corners( 1, nCorner );
415 list< const SMDS_MeshNode* >::iterator corner = corners.begin();
416 while ( corner != corners.end() )
418 SMDS_ElemIteratorPtr faceIt = (*corner)->GetInverseElementIterator( SMDSAbs_Face );
419 while ( faceIt->more() )
421 const SMDS_MeshElement* face = faceIt->next();
422 if ( !cornerFaces.insert( face ).second )
423 continue; // already loaded block side
425 if ( !isQuadrangle( face ))
426 return error("Non-quadrangle elements in the input mesh");
428 if ( _allSides.empty() || !_allSides.back()._grid.empty() )
429 _allSides.push_back( _BlockSide() );
431 _BlockSide& side = _allSides.back();
432 if ( !fillSide( side, face, *corner ) )
434 if ( !_error.empty() )
439 for ( int isXMax = 0; isXMax < 2; ++isXMax )
440 for ( int isYMax = 0; isYMax < 2; ++isYMax )
442 const SMDS_MeshNode* nCorner = side.getCornerNode(isXMax,isYMax );
443 corners.push_back( nCorner );
444 cornerFaces.insert( side.getCornerFace( nCorner ));
446 for ( int e = 0; e < int(NB_QUAD_SIDES); ++e )
447 _edge2sides[ side.getEdge( (EQuadEdge) e ) ].insert( &side );
449 nbFacesOnSides += side.getNbFaces();
454 // find block sides of other domains if any
455 if ( corner == corners.end() && nbFacesOnSides < mesh.NbQuadrangles() )
457 while ( nIt->more() )
459 nCorner = nIt->next();
460 if ( isCornerNode( nCorner ))
461 corner = corners.insert( corner, nCorner );
463 nbFacesOnSides = mesh.NbQuadrangles();
467 if ( _allSides.empty() )
469 if ( _allSides.back()._grid.empty() )
470 _allSides.pop_back();
471 _DUMP_("Nb detected sides "<< _allSides.size());
473 // ---------------------------
474 // Organize sides into blocks
475 // ---------------------------
477 // analyse sharing of sides by blocks and sort sides by nb of adjacent sides
478 int nbBlockSides = 0; // total nb of block sides taking into account their sharing
479 multimap<int, _BlockSide* > sortedSides;
481 list < _BlockSide >::iterator sideIt = _allSides.begin();
482 for ( ; sideIt != _allSides.end(); ++sideIt )
484 _BlockSide& side = *sideIt;
485 bool isSharedSide = true;
487 for ( int e = 0; e < int(NB_QUAD_SIDES) && isSharedSide; ++e )
489 int nbAdj = _edge2sides[ side.getEdge( (EQuadEdge) e ) ].size();
491 isSharedSide = ( nbAdj > 2 );
493 side._nbBlocksFound = 0;
494 side._nbBlocksExpected = isSharedSide ? 2 : 1;
495 nbBlockSides += side._nbBlocksExpected;
496 sortedSides.insert( make_pair( nbAdjacent, & side ));
500 // find sides of each block
502 while ( nbBlockSides >= 6 )
504 // get any side not bound to all blocks it belongs to
505 multimap<int, _BlockSide*>::iterator i_side = sortedSides.begin();
506 while ( i_side != sortedSides.end() && i_side->second->isBound())
509 // start searching for block sides from the got side
511 if ( _blocks.empty() || _blocks.back()._side[B_FRONT] )
512 _blocks.resize( _blocks.size() + 1 );
514 _Block& block = _blocks.back();
515 block.setSide( B_FRONT, i_side->second );
516 setSideBoundToBlock( *i_side->second );
519 // edges of adjacent sides of B_FRONT corresponding to front's edges
520 EQuadEdge edgeOfFront[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
521 EQuadEdge edgeOfAdj [4] = { Q_BOTTOM, Q_LEFT, Q_BOTTOM, Q_LEFT };
522 // first find all sides detectable w/o advanced analysis,
523 // then repeat the search, which then may pass without advanced analysis
524 set< _BlockSide* > sidesAround;
525 for ( int advAnalys = 0; advAnalys < 2; ++advAnalys )
527 // try to find 4 sides adjacent to a FRONT side
528 for ( int i = 0; (ok || !advAnalys) && i < NB_QUAD_SIDES; ++i )
529 if ( !block._side[i] )
530 ok = block.setSide( i, findBlockSide( B_FRONT, edgeOfFront[i], edgeOfAdj[i],
531 advAnalys, sidesAround));
532 // try to find a BACK side by a TOP one
533 if ( ok || !advAnalys)
534 if ( !block._side[B_BACK] && block._side[B_TOP] )
535 ok = block.setSide( B_BACK, findBlockSide( B_TOP, Q_TOP, Q_TOP,
536 advAnalys, sidesAround ));
537 if ( !advAnalys ) ok = true;
539 ok = block.isValid();
542 // check if just found block is same as one of previously found blocks
544 for ( int i = 1; i < _blocks.size() && !isSame; ++i )
545 isSame = ( block._corners == _blocks[i-1]._corners );
549 // count the found sides
550 _DUMP_(endl << "** Block " << _blocks.size() << " valid: " << block.isValid());
551 for (int i = 0; i < NB_BLOCK_SIDES; ++i )
553 _DUMP_("\tSide "<< SBoxSides[i] <<" "<< block._side[ i ]._side);
554 if ( block._side[ i ] )
556 if ( ok && i != B_FRONT)
558 setSideBoundToBlock( *block._side[ i ]._side );
561 _DUMP_("\t corners "<<
562 block._side[ i ].cornerNode(0,0)->GetID() << ", " <<
563 block._side[ i ].cornerNode(1,0)->GetID() << ", " <<
564 block._side[ i ].cornerNode(1,1)->GetID() << ", " <<
565 block._side[ i ].cornerNode(0,1)->GetID() << ", "<<endl);
569 _DUMP_("\t not found"<<endl);
577 _DUMP_("Nb found blocks "<< nbBlocks <<endl);
579 if ( nbBlocks == 0 && _error.empty() )
585 //================================================================================
587 * \brief Fill block side data starting from its corner quadrangle
589 //================================================================================
591 bool _Skin::fillSide( _BlockSide& side,
592 const SMDS_MeshElement* cornerQuad,
593 const SMDS_MeshNode* nCorner)
595 // Find out size of block side mesured in nodes and by the way find two rows
596 // of nodes in two directions.
599 const SMDS_MeshElement* firstQuad = cornerQuad;
601 // get a node on block edge
602 int iCorner = firstQuad->GetNodeIndex( nCorner );
603 const SMDS_MeshNode* nOnEdge = firstQuad->GetNode( (iCorner+1) % 4);
605 // find out size of block side
606 vector<const SMDS_MeshNode*> horRow1, horRow2, verRow1, verRow2;
607 if ( !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, horRow1, horRow2, true ) ||
608 !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, verRow1, verRow2, false ))
610 nbX = horRow1.size(), nbY = verRow1.size();
613 side._index._xSize = horRow1.size();
614 side._index._ySize = verRow1.size();
615 side._grid.resize( side._index.size(), NULL );
617 for ( x = 0; x < horRow1.size(); ++x )
619 side.setNode( x, 0, horRow1[x] );
620 side.setNode( x, 1, horRow2[x] );
622 for ( y = 0; y < verRow1.size(); ++y )
624 side.setNode( 0, y, verRow1[y] );
625 side.setNode( 1, y, verRow2[y] );
628 // Find the rest nodes
630 y = 1; // y of the row to fill
631 TIDSortedElemSet emptySet, avoidSet;
634 // get next firstQuad in the next row of quadrangles
639 // o---o o o o o <- found nodes
642 int i1down, i2down, i2up;
643 const SMDS_MeshNode* n1down = side.getNode( 0, y-1 );
644 const SMDS_MeshNode* n2down = side.getNode( 1, y-1 );
645 avoidSet.clear(); avoidSet.insert( firstQuad );
646 firstQuad = SMESH_MeshAlgos::FindFaceInSet( n1down, n2down, emptySet, avoidSet,
648 if ( !isQuadrangle( firstQuad ))
651 const SMDS_MeshNode* n2up = oppositeNode( firstQuad, i1down );
652 avoidSet.clear(); avoidSet.insert( firstQuad );
654 // find the rest nodes in the y-th row by faces in the row
659 const SMDS_MeshElement* quad = SMESH_MeshAlgos::FindFaceInSet( n2up, n2down, emptySet,
660 avoidSet, &i2up, &i2down);
661 if ( !isQuadrangle( quad ))
664 n2up = oppositeNode( quad, i2down );
665 n2down = oppositeNode( quad, i2up );
666 avoidSet.clear(); avoidSet.insert( quad );
668 side.setNode( x, y, n2up );
672 // check side validity
674 side.getCornerFace( side.getCornerNode( 0, 0 )) &&
675 side.getCornerFace( side.getCornerNode( 1, 0 )) &&
676 side.getCornerFace( side.getCornerNode( 0, 1 )) &&
677 side.getCornerFace( side.getCornerNode( 1, 1 ));
682 //================================================================================
684 * \brief Return true if it's possible to make a loop over corner2Sides starting
687 //================================================================================
689 bool isClosedChainOfSides( _BlockSide* startSide,
690 map< const SMDS_MeshNode*, list< _BlockSide* > > & corner2Sides )
692 // get start and end nodes
693 const SMDS_MeshNode *n1 = 0, *n2 = 0, *n;
694 for ( int y = 0; y < 2; ++y )
695 for ( int x = 0; x < 2; ++x )
697 n = startSide->getCornerNode(x,y);
698 if ( !corner2Sides.count( n )) continue;
704 if ( !n2 ) return false;
706 map< const SMDS_MeshNode*, list< _BlockSide* > >::iterator
707 c2sides = corner2Sides.find( n1 );
708 if ( c2sides == corner2Sides.end() ) return false;
710 int nbChainLinks = 1;
712 _BlockSide* prevSide = startSide;
715 // get the next side sharing n
716 list< _BlockSide* > & sides = c2sides->second;
717 _BlockSide* nextSide = ( sides.back() == prevSide ? sides.front() : sides.back() );
718 if ( nextSide == prevSide ) return false;
720 // find the next corner of the nextSide being in corner2Sides
723 for ( int y = 0; y < 2 && !n; ++y )
724 for ( int x = 0; x < 2; ++x )
726 n = nextSide->getCornerNode(x,y);
727 c2sides = corner2Sides.find( n );
728 if ( n == n1 || c2sides == corner2Sides.end() )
733 if ( !n ) return false;
739 return ( n == n2 && nbChainLinks == NB_QUAD_SIDES );
742 //================================================================================
744 * \brief Try to find a block side adjacent to the given side by given edge
746 //================================================================================
748 _OrientedBlockSide _Skin::findBlockSide( EBoxSides startBlockSide,
749 EQuadEdge sharedSideEdge1,
750 EQuadEdge sharedSideEdge2,
751 bool withGeometricAnalysis,
752 set< _BlockSide* >& sidesAround)
754 _Block& block = _blocks.back();
755 _OrientedBlockSide& side1 = block._side[ startBlockSide ];
757 // get corner nodes of the given block edge
758 SMESH_OrientedLink edge = side1.edge( sharedSideEdge1 );
759 const SMDS_MeshNode* n1 = edge.node1();
760 const SMDS_MeshNode* n2 = edge.node2();
761 if ( edge._reversed ) swap( n1, n2 );
763 // find all sides sharing both nodes n1 and n2
764 set< _BlockSide* > sidesOnEdge = _edge2sides[ edge ]; // copy a set
766 // exclude loaded sides of block from sidesOnEdge
767 for (int i = 0; i < NB_BLOCK_SIDES; ++i )
768 if ( block._side[ i ] )
769 sidesOnEdge.erase( block._side[ i ]._side );
771 int nbSidesOnEdge = sidesOnEdge.size();
772 _DUMP_("nbSidesOnEdge "<< nbSidesOnEdge << " " << n1->GetID() << "-" << n2->GetID() );
773 if ( nbSidesOnEdge == 0 )
776 _BlockSide* foundSide = 0;
777 if ( nbSidesOnEdge == 1 )
779 foundSide = *sidesOnEdge.begin();
783 set< _BlockSide* >::iterator sideIt = sidesOnEdge.begin();
784 int nbLoadedSides = block.nbSides();
785 if ( nbLoadedSides > 1 )
787 // Find the side having more than 2 corners common with already loaded sides
788 for (; !foundSide && sideIt != sidesOnEdge.end(); ++sideIt )
790 _BlockSide* sideI = *sideIt;
791 int nbCommonCorners =
792 block._corners.count( sideI->getCornerNode(0,0)) +
793 block._corners.count( sideI->getCornerNode(1,0)) +
794 block._corners.count( sideI->getCornerNode(0,1)) +
795 block._corners.count( sideI->getCornerNode(1,1));
796 if ( nbCommonCorners > 2 )
803 if ( !withGeometricAnalysis )
805 sidesAround.insert( sidesOnEdge.begin(), sidesOnEdge.end() );
808 if ( nbLoadedSides == 1 )
810 // Issue 0021529. There are at least 2 sides by each edge and
811 // position of block gravity center is undefined.
812 // Find a side starting from which we can walk around the startBlockSide
814 // fill in corner2Sides
815 map< const SMDS_MeshNode*, list< _BlockSide* > > corner2Sides;
816 for ( sideIt = sidesAround.begin(); sideIt != sidesAround.end(); ++sideIt )
818 _BlockSide* sideI = *sideIt;
819 corner2Sides[ sideI->getCornerNode(0,0) ].push_back( sideI );
820 corner2Sides[ sideI->getCornerNode(1,0) ].push_back( sideI );
821 corner2Sides[ sideI->getCornerNode(0,1) ].push_back( sideI );
822 corner2Sides[ sideI->getCornerNode(1,1) ].push_back( sideI );
824 // remove corners of startBlockSide from corner2Sides
825 set<const SMDS_MeshNode*>::iterator nIt = block._corners.begin();
826 for ( ; nIt != block._corners.end(); ++nIt )
827 corner2Sides.erase( *nIt );
830 for ( sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
832 if ( isClosedChainOfSides( *sideIt, corner2Sides ))
843 // Select one of found sides most close to startBlockSide
845 gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()), p2 (n2->X(),n2->Y(),n2->Z());
846 gp_Vec p1p2( p1, p2 );
848 const SMDS_MeshElement* face1 = side1->getCornerFace( n1 );
849 gp_XYZ p1Op = SMESH_TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(n1)));
850 gp_Vec side1Dir( p1, p1Op );
851 gp_Ax2 pln( p1, p1p2, side1Dir ); // plane with normal p1p2 and X dir side1Dir
852 _DUMP_(" Select adjacent for "<< side1._side << " - side dir ("
853 << side1Dir.X() << ", " << side1Dir.Y() << ", " << side1Dir.Z() << ")" );
855 map < double , _BlockSide* > angleOfSide;
856 for (sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
858 _BlockSide* sideI = *sideIt;
859 const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 );
860 gp_XYZ p1Op = SMESH_TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(n1)));
861 gp_Vec sideIDir( p1, p1Op );
862 // compute angle of (sideIDir projection to pln) and (X dir of pln)
863 gp_Vec2d sideIDirProj( sideIDir * pln.XDirection(), sideIDir * pln.YDirection());
864 double angle = sideIDirProj.Angle( gp::DX2d() );
865 if ( angle < 0 ) angle += 2. * M_PI; // angle [0-2*PI]
866 angleOfSide.insert( make_pair( angle, sideI ));
867 _DUMP_(" "<< sideI << " - side dir ("
868 << sideIDir.X() << ", " << sideIDir.Y() << ", " << sideIDir.Z() << ")"
869 << " angle " << angle);
872 gp_XYZ gc(0,0,0); // gravity center of already loaded block sides
873 for (int i = 0; i < NB_BLOCK_SIDES; ++i )
874 if ( block._side[ i ] )
875 gc += block._side[ i ]._side->getGC();
878 gp_Vec gcDir( p1, gc );
879 gp_Vec2d gcDirProj( gcDir * pln.XDirection(), gcDir * pln.YDirection());
880 double gcAngle = gcDirProj.Angle( gp::DX2d() );
881 foundSide = gcAngle < 0 ? angleOfSide.rbegin()->second : angleOfSide.begin()->second;
884 _DUMP_(" selected "<< foundSide );
887 // Orient the found side correctly
889 // corners of found side corresponding to nodes n1 and n2
890 bool xMax1, yMax1, xMax2, yMax2;
891 if ( !getEdgeEnds( sharedSideEdge2, xMax1, yMax1, xMax2, yMax2 ))
892 return error(SMESH_Comment("Internal error at ")<<__FILE__<<":"<<__LINE__),
893 _OrientedBlockSide(0);
895 for ( int ori = 0; ori < _OrientedIndexer::MAX_ORI+1; ++ori )
897 _OrientedBlockSide orientedSide( foundSide, ori );
898 const SMDS_MeshNode* n12 = orientedSide.cornerNode( xMax1, yMax1);
899 const SMDS_MeshNode* n22 = orientedSide.cornerNode( xMax2, yMax2);
900 if ( n1 == n12 && n2 == n22 )
903 error(SMESH_Comment("Failed to orient a block side found by edge ")<<sharedSideEdge1
904 << " of side " << startBlockSide
905 << " of block " << _blocks.size());
909 //================================================================================
911 * \brief: Fill rows (which are actually columns,if !alongN1N2) of nodes starting
912 * from the given quadrangle until another block corner encounters.
913 * n1 and n2 are at bottom of quad, n1 is at block corner.
915 //================================================================================
917 bool _Skin::fillRowsUntilCorner(const SMDS_MeshElement* quad,
918 const SMDS_MeshNode* n1,
919 const SMDS_MeshNode* n2,
920 vector<const SMDS_MeshNode*>& row1,
921 vector<const SMDS_MeshNode*>& row2,
922 const bool alongN1N2 )
924 const SMDS_MeshNode* corner1 = n1;
926 // Store nodes of quad in the rows and find new n1 and n2 to get
927 // the next face so that new n2 is on block edge
928 int i1 = quad->GetNodeIndex( n1 );
929 int i2 = quad->GetNodeIndex( n2 );
930 row1.clear(); row2.clear();
931 row1.push_back( n1 );
934 row1.push_back( n2 );
935 row2.push_back( oppositeNode( quad, i2 ));
936 row2.push_back( n1 = oppositeNode( quad, i1 ));
940 row2.push_back( n2 );
941 row1.push_back( n2 = oppositeNode( quad, i2 ));
942 row2.push_back( n1 = oppositeNode( quad, i1 ));
945 if ( isCornerNode( row1[1] ))
948 // Find the rest nodes
949 TIDSortedElemSet emptySet, avoidSet;
950 while ( !isCornerNode( n2 ) )
952 avoidSet.clear(); avoidSet.insert( quad );
953 quad = SMESH_MeshAlgos::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2 );
954 if ( !isQuadrangle( quad ))
957 row1.push_back( n2 = oppositeNode( quad, i1 ));
958 row2.push_back( n1 = oppositeNode( quad, i2 ));
960 return n1 != corner1;
963 //================================================================================
965 * \brief Return a corner face by a corner node
967 //================================================================================
969 const SMDS_MeshElement* _BlockSide::getCornerFace(const SMDS_MeshNode* cornerNode) const
971 int x, y, isXMax, isYMax, found = 0;
972 for ( isXMax = 0; isXMax < 2; ++isXMax )
974 for ( isYMax = 0; isYMax < 2; ++isYMax )
976 x = isXMax ? _index._xSize-1 : 0;
977 y = isYMax ? _index._ySize-1 : 0;
978 found = ( getNode(x,y) == cornerNode );
983 if ( !found ) return 0;
984 int dx = isXMax ? -1 : +1;
985 int dy = isYMax ? -1 : +1;
986 const SMDS_MeshNode* n1 = getNode(x,y);
987 const SMDS_MeshNode* n2 = getNode(x+dx,y);
988 const SMDS_MeshNode* n3 = getNode(x,y+dy);
989 const SMDS_MeshNode* n4 = getNode(x+dx,y+dy);
990 return SMDS_Mesh::FindFace(n1, n2, n3, n4 );
993 //================================================================================
995 * \brief Checks own validity
997 //================================================================================
999 bool _Block::isValid() const
1001 bool ok = ( nbSides() == 6 );
1003 // check only corners depending on side selection
1004 EBoxSides adjacent[4] = { B_BOTTOM, B_RIGHT, B_TOP, B_LEFT };
1005 EQuadEdge edgeAdj [4] = { Q_TOP, Q_RIGHT, Q_TOP, Q_RIGHT };
1006 EQuadEdge edgeBack[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
1008 for ( int i=0; ok && i < NB_QUAD_SIDES; ++i )
1010 SMESH_OrientedLink eBack = _side[ B_BACK ].edge( edgeBack[i] );
1011 SMESH_OrientedLink eAdja = _side[ adjacent[i] ].edge( edgeAdj[i] );
1012 ok = ( eBack == eAdja );
1023 HEXA_NS::Hexa* _block2Hexa( const _Block& block,
1024 HEXA_NS::Document* doc,
1025 std::map<HEXA_NS::Vertex*, SMDS_MeshNode*> vertexNode )
1027 const SMDS_MeshNode* n000 = block.getSide(B_BOTTOM).cornerNode( 0, 0 );
1028 const SMDS_MeshNode* n100 = block.getSide(B_BOTTOM).cornerNode( 1, 0 );
1029 const SMDS_MeshNode* n010 = block.getSide(B_BOTTOM).cornerNode( 0, 1 );
1030 const SMDS_MeshNode* n110 = block.getSide(B_BOTTOM).cornerNode( 1, 1 );
1031 const SMDS_MeshNode* n001 = block.getSide(B_TOP).cornerNode( 0, 0 );
1032 const SMDS_MeshNode* n101 = block.getSide(B_TOP).cornerNode( 1, 0 );
1033 const SMDS_MeshNode* n011 = block.getSide(B_TOP).cornerNode( 0, 1 );
1034 const SMDS_MeshNode* n111 = block.getSide(B_TOP).cornerNode( 1, 1 );
1036 list<const SMDS_MeshNode* > nodeFromBlock;
1037 nodeFromBlock.push_back(n000);
1038 nodeFromBlock.push_back(n100);
1039 nodeFromBlock.push_back(n010);
1040 nodeFromBlock.push_back(n110);
1041 nodeFromBlock.push_back(n001);
1042 nodeFromBlock.push_back(n101);
1043 nodeFromBlock.push_back(n011);
1044 nodeFromBlock.push_back(n111);
1045 nodeFromBlock.sort();
1048 HEXA_NS::Hexa* hexa = NULL;
1049 int nHexa = doc->countUsedHexa();
1050 for (int j=0; j <nHexa; ++j ){
1051 hexa = doc->getUsedHexa(j);
1052 list<const SMDS_MeshNode* > nodeFromHexa;
1053 int nVx = hexa->countVertex();
1054 for ( int i=0; i <nVx; ++i ){
1055 HEXA_NS::Vertex* v = hexa->getVertex(i);
1056 const SMDS_MeshNode* n = vertexNode[v];
1057 nodeFromHexa.push_back(n);
1059 nodeFromHexa.sort();
1061 if ( nodeFromBlock == nodeFromHexa ){
1062 // std::cout << "OK block match hexa "<< hexa <<" id = "<<hexa->getId()<<std::endl;
1065 std::cout << "KO block does not match hexa "<< hexa <<" id = "<<hexa->getId()<<std::endl;
1078 //=======================================================================
1079 //function : SMESH_HexaFromSkin_3D
1081 //=======================================================================
1083 // SMESH_HexaFromSkin_3D::SMESH_HexaFromSkin_3D(int hypId, SMESH_Gen* gen)
1084 // :SMESH_3D_Algo(hypId, gen)
1086 // MESSAGE("SMESH_HexaFromSkin_3D::SMESH_HexaFromSkin_3D");
1087 // _name = "HexaFromSkin_3D";
1090 SMESH_HexaFromSkin_3D::SMESH_HexaFromSkin_3D(int hypId, SMESH_Gen* gen, HEXA_NS::Document* doc)
1091 :SMESH_3D_Algo(hypId, gen),
1094 if(MYDEBUG) MESSAGE("SMESH_HexaFromSkin_3D::SMESH_HexaFromSkin_3D");
1095 _name = "HexaFromSkin_3D";
1099 SMESH_HexaFromSkin_3D::~SMESH_HexaFromSkin_3D()
1101 if(MYDEBUG) MESSAGE("SMESH_HexaFromSkin_3D::~SMESH_HexaFromSkin_3D");
1104 //================================================================================
1106 * \brief Main method, which generates hexaheda
1108 //================================================================================
1110 bool SMESH_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
1113 int nbBlocks = skin.findBlocks(aMesh);
1114 if ( nbBlocks == 0 )
1115 return error( skin.error());
1117 vector< vector< const SMDS_MeshNode* > > columns;
1118 int x, xSize, y, ySize, z, zSize;
1121 for ( int i = 0; i < nbBlocks; ++i )
1123 const _Block& block = skin.getBlock( i );
1125 // ------------------------------------------
1126 // Fill columns of nodes with existing nodes
1127 // ------------------------------------------
1129 xSize = block.getSide(B_BOTTOM).getHoriSize();
1130 ySize = block.getSide(B_BOTTOM).getVertSize();
1131 zSize = block.getSide(B_FRONT ).getVertSize();
1132 int X = xSize - 1, Y = ySize - 1, Z = zSize - 1;
1133 colIndex = _Indexer( xSize, ySize );
1134 columns.resize( colIndex.size() );
1136 // fill node columns by front and back box sides
1137 for ( x = 0; x < xSize; ++x ) {
1138 vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
1139 vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
1140 column0.resize( zSize );
1141 column1.resize( zSize );
1142 for ( z = 0; z < zSize; ++z ) {
1143 column0[ z ] = block.getSide(B_FRONT).node( x, z );
1144 column1[ z ] = block.getSide(B_BACK) .node( x, z );
1147 // fill node columns by left and right box sides
1148 for ( y = 1; y < ySize-1; ++y ) {
1149 vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
1150 vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
1151 column0.resize( zSize );
1152 column1.resize( zSize );
1153 for ( z = 0; z < zSize; ++z ) {
1154 column0[ z ] = block.getSide(B_LEFT) .node( y, z );
1155 column1[ z ] = block.getSide(B_RIGHT).node( y, z );
1158 // get nodes from top and bottom box sides
1159 for ( x = 1; x < xSize-1; ++x ) {
1160 for ( y = 1; y < ySize-1; ++y ) {
1161 vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
1162 column.resize( zSize );
1163 column.front() = block.getSide(B_BOTTOM).node( x, y );
1164 column.back() = block.getSide(B_TOP) .node( x, y );
1168 // ----------------------------
1169 // Add internal nodes of a box
1170 // ----------------------------
1171 // projection points of internal nodes on box sub-shapes by which
1172 // coordinates of internal nodes are computed
1173 vector<gp_XYZ> pointOnShape( SMESH_Block::ID_Shell );
1175 // projections on vertices are constant
1176 pointOnShape[ SMESH_Block::ID_V000 ] = block.getSide(B_BOTTOM).xyz( 0, 0 );
1177 pointOnShape[ SMESH_Block::ID_V100 ] = block.getSide(B_BOTTOM).xyz( X, 0 );
1178 pointOnShape[ SMESH_Block::ID_V010 ] = block.getSide(B_BOTTOM).xyz( 0, Y );
1179 pointOnShape[ SMESH_Block::ID_V110 ] = block.getSide(B_BOTTOM).xyz( X, Y );
1180 pointOnShape[ SMESH_Block::ID_V001 ] = block.getSide(B_TOP).xyz( 0, 0 );
1181 pointOnShape[ SMESH_Block::ID_V101 ] = block.getSide(B_TOP).xyz( X, 0 );
1182 pointOnShape[ SMESH_Block::ID_V011 ] = block.getSide(B_TOP).xyz( 0, Y );
1183 pointOnShape[ SMESH_Block::ID_V111 ] = block.getSide(B_TOP).xyz( X, Y );
1185 for ( x = 1; x < xSize-1; ++x )
1187 gp_XYZ params; // normalized parameters of internal node within a unit box
1188 params.SetCoord( 1, x / double(X) );
1189 for ( y = 1; y < ySize-1; ++y )
1191 params.SetCoord( 2, y / double(Y) );
1192 // column to fill during z loop
1193 vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
1194 // projections on horizontal edges
1195 pointOnShape[ SMESH_Block::ID_Ex00 ] = block.getSide(B_BOTTOM).xyz( x, 0 );
1196 pointOnShape[ SMESH_Block::ID_Ex10 ] = block.getSide(B_BOTTOM).xyz( x, Y );
1197 pointOnShape[ SMESH_Block::ID_E0y0 ] = block.getSide(B_BOTTOM).xyz( 0, y );
1198 pointOnShape[ SMESH_Block::ID_E1y0 ] = block.getSide(B_BOTTOM).xyz( X, y );
1199 pointOnShape[ SMESH_Block::ID_Ex01 ] = block.getSide(B_TOP).xyz( x, 0 );
1200 pointOnShape[ SMESH_Block::ID_Ex11 ] = block.getSide(B_TOP).xyz( x, Y );
1201 pointOnShape[ SMESH_Block::ID_E0y1 ] = block.getSide(B_TOP).xyz( 0, y );
1202 pointOnShape[ SMESH_Block::ID_E1y1 ] = block.getSide(B_TOP).xyz( X, y );
1203 // projections on horizontal sides
1204 pointOnShape[ SMESH_Block::ID_Fxy0 ] = block.getSide(B_BOTTOM).xyz( x, y );
1205 pointOnShape[ SMESH_Block::ID_Fxy1 ] = block.getSide(B_TOP) .xyz( x, y );
1206 for ( z = 1; z < zSize-1; ++z ) // z loop
1208 params.SetCoord( 3, z / double(Z) );
1209 // projections on vertical edges
1210 pointOnShape[ SMESH_Block::ID_E00z ] = block.getSide(B_FRONT).xyz( 0, z );
1211 pointOnShape[ SMESH_Block::ID_E10z ] = block.getSide(B_FRONT).xyz( X, z );
1212 pointOnShape[ SMESH_Block::ID_E01z ] = block.getSide(B_BACK).xyz( 0, z );
1213 pointOnShape[ SMESH_Block::ID_E11z ] = block.getSide(B_BACK).xyz( X, z );
1214 // projections on vertical sides
1215 pointOnShape[ SMESH_Block::ID_Fx0z ] = block.getSide(B_FRONT).xyz( x, z );
1216 pointOnShape[ SMESH_Block::ID_Fx1z ] = block.getSide(B_BACK) .xyz( x, z );
1217 pointOnShape[ SMESH_Block::ID_F0yz ] = block.getSide(B_LEFT) .xyz( y, z );
1218 pointOnShape[ SMESH_Block::ID_F1yz ] = block.getSide(B_RIGHT).xyz( y, z );
1220 // compute internal node coordinates
1222 SMESH_Block::ShellPoint( params, pointOnShape, coords );
1223 column[ z ] = aHelper->AddNode( coords.X(), coords.Y(), coords.Z() );
1227 //cout << "----------------------------------------------------------------------"<<endl;
1228 //for ( int id = SMESH_Block::ID_V000; id < SMESH_Block::ID_Shell; ++id)
1230 // gp_XYZ p = pointOnShape[ id ];
1231 // SMESH_Block::DumpShapeID( id,cout)<<" ( "<<p.X()<<", "<<p.Y()<<", "<<p.Z()<<" )"<<endl;
1233 //cout << "Params: ( "<< params.X()<<", "<<params.Y()<<", "<<params.Z()<<" )"<<endl;
1234 //cout << "coords: ( "<< coords.X()<<", "<<coords.Y()<<", "<<coords.Z()<<" )"<<endl;
1243 // find out orientation by a least distorted hexahedron (issue 0020855);
1244 // the last is defined by evaluating sum of face normals of 8 corner hexahedrons
1245 double badness = numeric_limits<double>::max();
1247 for ( int xMax = 0; xMax < 2; ++xMax )
1248 for ( int yMax = 0; yMax < 2; ++yMax )
1249 for ( int zMax = 0; zMax < 2; ++zMax )
1251 x = xMax ? xSize-1 : 1;
1252 y = yMax ? ySize-1 : 1;
1253 z = zMax ? zSize-1 : 1;
1254 vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x-1, y-1 )];
1255 vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x , y-1 )];
1256 vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x-1, y )];
1257 vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x , y )];
1259 const SMDS_MeshNode* n000 = col00[z-1];
1260 const SMDS_MeshNode* n100 = col10[z-1];
1261 const SMDS_MeshNode* n010 = col01[z-1];
1262 const SMDS_MeshNode* n110 = col11[z-1];
1263 const SMDS_MeshNode* n001 = col00[z];
1264 const SMDS_MeshNode* n101 = col10[z];
1265 const SMDS_MeshNode* n011 = col01[z];
1266 const SMDS_MeshNode* n111 = col11[z];
1267 SMDS_VolumeOfNodes probeVolume (n000,n010,n110,n100,
1268 n001,n011,n111,n101);
1269 SMDS_VolumeTool volTool( &probeVolume );
1270 double Nx=0.,Ny=0.,Nz=0.;
1271 for ( int iFace = 0; iFace < volTool.NbFaces(); ++iFace )
1274 volTool.GetFaceNormal( iFace, nx,ny,nz );
1279 double quality = Nx*Nx + Ny*Ny + Nz*Nz;
1280 if ( quality < badness )
1283 isForw = volTool.IsForward();
1288 for ( x = 0; x < xSize-1; ++x ) {
1289 for ( y = 0; y < ySize-1; ++y ) {
1290 vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
1291 vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
1292 vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
1293 vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
1294 // bottom face normal of a hexa mush point outside the volume
1296 for ( z = 0; z < zSize-1; ++z )
1297 aHelper->AddVolume(col00[z], col01[z], col11[z], col10[z],
1298 col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
1300 for ( z = 0; z < zSize-1; ++z )
1301 aHelper->AddVolume(col00[z], col10[z], col11[z], col01[z],
1302 col00[z+1], col10[z+1], col11[z+1], col01[z+1]);
1311 //================================================================================
1313 * \brief Main method, which generates hexaheda
1315 //================================================================================
1316 bool SMESH_HexaFromSkin_3D::Compute( SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper,
1317 std::map<HEXA_NS::Hexa*, SMESH_HexaBlocks::SMESHVolumes>& volumesOnHexa,
1318 std::map<HEXA_NS::Vertex*, SMDS_MeshNode*> vertexNode )
1320 if(MYDEBUG) MESSAGE("SMESH_HexaFromSkin_3D::Compute BEGIN");
1322 int nbBlocks = skin.findBlocks(aMesh);
1323 if ( nbBlocks == 0 )
1324 return error( skin.error());
1326 vector< vector< const SMDS_MeshNode* > > columns;
1327 int x, xSize, y, ySize, z, zSize;
1330 for ( int i = 0; i < nbBlocks; ++i )
1332 const _Block& block = skin.getBlock( i );
1334 // ------------------------------------------
1335 // Fill columns of nodes with existing nodes
1336 // ------------------------------------------
1338 xSize = block.getSide(B_BOTTOM).getHoriSize();
1339 ySize = block.getSide(B_BOTTOM).getVertSize();
1340 zSize = block.getSide(B_FRONT ).getVertSize();
1341 int X = xSize - 1, Y = ySize - 1, Z = zSize - 1;
1342 colIndex = _Indexer( xSize, ySize );
1343 columns.resize( colIndex.size() );
1345 // fill node columns by front and back box sides
1346 for ( x = 0; x < xSize; ++x ) {
1347 vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
1348 vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
1349 column0.resize( zSize );
1350 column1.resize( zSize );
1351 for ( z = 0; z < zSize; ++z ) {
1352 column0[ z ] = block.getSide(B_FRONT).node( x, z );
1353 column1[ z ] = block.getSide(B_BACK) .node( x, z );
1356 // fill node columns by left and right box sides
1357 for ( y = 1; y < ySize-1; ++y ) {
1358 vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
1359 vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
1360 column0.resize( zSize );
1361 column1.resize( zSize );
1362 for ( z = 0; z < zSize; ++z ) {
1363 column0[ z ] = block.getSide(B_LEFT) .node( y, z );
1364 column1[ z ] = block.getSide(B_RIGHT).node( y, z );
1367 // get nodes from top and bottom box sides
1368 for ( x = 1; x < xSize-1; ++x ) {
1369 for ( y = 1; y < ySize-1; ++y ) {
1370 vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
1371 column.resize( zSize );
1372 column.front() = block.getSide(B_BOTTOM).node( x, y );
1373 column.back() = block.getSide(B_TOP) .node( x, y );
1377 // ----------------------------
1378 // Add internal nodes of a box
1379 // ----------------------------
1380 // projection points of internal nodes on box subshapes by which
1381 // coordinates of internal nodes are computed
1382 vector<gp_XYZ> pointOnShape( SMESH_Block::ID_Shell );
1384 // projections on vertices are constant
1385 pointOnShape[ SMESH_Block::ID_V000 ] = block.getSide(B_BOTTOM).xyz( 0, 0 );
1386 pointOnShape[ SMESH_Block::ID_V100 ] = block.getSide(B_BOTTOM).xyz( X, 0 );
1387 pointOnShape[ SMESH_Block::ID_V010 ] = block.getSide(B_BOTTOM).xyz( 0, Y );
1388 pointOnShape[ SMESH_Block::ID_V110 ] = block.getSide(B_BOTTOM).xyz( X, Y );
1389 pointOnShape[ SMESH_Block::ID_V001 ] = block.getSide(B_TOP).xyz( 0, 0 );
1390 pointOnShape[ SMESH_Block::ID_V101 ] = block.getSide(B_TOP).xyz( X, 0 );
1391 pointOnShape[ SMESH_Block::ID_V011 ] = block.getSide(B_TOP).xyz( 0, Y );
1392 pointOnShape[ SMESH_Block::ID_V111 ] = block.getSide(B_TOP).xyz( X, Y );
1394 for ( x = 1; x < xSize-1; ++x )
1396 gp_XYZ params; // normalized parameters of internal node within a unit box
1397 params.SetCoord( 1, x / double(X) );
1398 for ( y = 1; y < ySize-1; ++y )
1400 params.SetCoord( 2, y / double(Y) );
1401 // column to fill during z loop
1402 vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
1403 // projections on horizontal edges
1404 pointOnShape[ SMESH_Block::ID_Ex00 ] = block.getSide(B_BOTTOM).xyz( x, 0 );
1405 pointOnShape[ SMESH_Block::ID_Ex10 ] = block.getSide(B_BOTTOM).xyz( x, Y );
1406 pointOnShape[ SMESH_Block::ID_E0y0 ] = block.getSide(B_BOTTOM).xyz( 0, y );
1407 pointOnShape[ SMESH_Block::ID_E1y0 ] = block.getSide(B_BOTTOM).xyz( X, y );
1408 pointOnShape[ SMESH_Block::ID_Ex01 ] = block.getSide(B_TOP).xyz( x, 0 );
1409 pointOnShape[ SMESH_Block::ID_Ex11 ] = block.getSide(B_TOP).xyz( x, Y );
1410 pointOnShape[ SMESH_Block::ID_E0y1 ] = block.getSide(B_TOP).xyz( 0, y );
1411 pointOnShape[ SMESH_Block::ID_E1y1 ] = block.getSide(B_TOP).xyz( X, y );
1412 // projections on horizontal sides
1413 pointOnShape[ SMESH_Block::ID_Fxy0 ] = block.getSide(B_BOTTOM).xyz( x, y );
1414 pointOnShape[ SMESH_Block::ID_Fxy1 ] = block.getSide(B_TOP) .xyz( x, y );
1415 for ( z = 1; z < zSize-1; ++z ) // z loop
1417 params.SetCoord( 3, z / double(Z) );
1418 // projections on vertical edges
1419 pointOnShape[ SMESH_Block::ID_E00z ] = block.getSide(B_FRONT).xyz( 0, z );
1420 pointOnShape[ SMESH_Block::ID_E10z ] = block.getSide(B_FRONT).xyz( X, z );
1421 pointOnShape[ SMESH_Block::ID_E01z ] = block.getSide(B_BACK).xyz( 0, z );
1422 pointOnShape[ SMESH_Block::ID_E11z ] = block.getSide(B_BACK).xyz( X, z );
1423 // projections on vertical sides
1424 pointOnShape[ SMESH_Block::ID_Fx0z ] = block.getSide(B_FRONT).xyz( x, z );
1425 pointOnShape[ SMESH_Block::ID_Fx1z ] = block.getSide(B_BACK) .xyz( x, z );
1426 pointOnShape[ SMESH_Block::ID_F0yz ] = block.getSide(B_LEFT) .xyz( y, z );
1427 pointOnShape[ SMESH_Block::ID_F1yz ] = block.getSide(B_RIGHT).xyz( y, z );
1429 // compute internal node coordinates
1431 SMESH_Block::ShellPoint( params, pointOnShape, coords );
1432 column[ z ] = aHelper->AddNode( coords.X(), coords.Y(), coords.Z() );
1436 //cout << "----------------------------------------------------------------------"<<endl;
1437 //for ( int id = SMESH_Block::ID_V000; id < SMESH_Block::ID_Shell; ++id)
1439 // gp_XYZ p = pointOnShape[ id ];
1440 // SMESH_Block::DumpShapeID( id,cout)<<" ( "<<p.X()<<", "<<p.Y()<<", "<<p.Z()<<" )"<<endl;
1442 //cout << "Params: ( "<< params.X()<<", "<<params.Y()<<", "<<params.Z()<<" )"<<endl;
1443 //cout << "coords: ( "<< coords.X()<<", "<<coords.Y()<<", "<<coords.Z()<<" )"<<endl;
1452 // find out orientation
1453 const SMDS_MeshNode* n000 = block.getSide(B_BOTTOM).cornerNode( 0, 0 );
1454 const SMDS_MeshNode* n100 = block.getSide(B_BOTTOM).cornerNode( 1, 0 );
1455 const SMDS_MeshNode* n010 = block.getSide(B_BOTTOM).cornerNode( 0, 1 );
1456 const SMDS_MeshNode* n110 = block.getSide(B_BOTTOM).cornerNode( 1, 1 );
1457 const SMDS_MeshNode* n001 = block.getSide(B_TOP).cornerNode( 0, 0 );
1458 const SMDS_MeshNode* n101 = block.getSide(B_TOP).cornerNode( 1, 0 );
1459 const SMDS_MeshNode* n011 = block.getSide(B_TOP).cornerNode( 0, 1 );
1460 const SMDS_MeshNode* n111 = block.getSide(B_TOP).cornerNode( 1, 1 );
1461 SMDS_VolumeOfNodes probeVolume (n000,n010,n110,n100,
1462 n001,n011,n111,n101);
1463 bool isForw = SMDS_VolumeTool( &probeVolume ).IsForward();
1466 SMESH_HexaBlocks::SMESHVolumes volumesOnBlock; //Groups creation
1468 for ( x = 0; x < xSize-1; ++x ) {
1469 for ( y = 0; y < ySize-1; ++y ) {
1470 vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
1471 vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
1472 vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
1473 vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
1474 // bottom face normal of a hexa mush point outside the volume
1476 for ( z = 0; z < zSize-1; ++z ){
1477 SMDS_MeshVolume* newVolume =
1478 aHelper->AddVolume(col00[z], col01[z], col11[z], col10[z],
1479 col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
1480 volumesOnBlock.push_back( newVolume );
1483 for ( z = 0; z < zSize-1; ++z ){
1484 SMDS_MeshVolume* newVolume =
1485 aHelper->AddVolume(col00[z], col10[z], col11[z], col01[z],
1486 col00[z+1], col10[z+1], col11[z+1], col01[z+1]);
1487 volumesOnBlock.push_back( newVolume );
1491 // std::cout << "block i = " << i << std::endl;
1492 HEXA_NS::Hexa* currentHexa = _block2Hexa( block, _doc, vertexNode );
1493 if ( currentHexa != NULL ){
1494 // std::cout<<"===== found ->"<<currentHexa<<" for block "<<i<<std::endl;
1495 if ( volumesOnHexa.count(currentHexa)==0 ) {
1496 // std::cout<<"added!! "<<std::endl;
1497 volumesOnHexa[currentHexa] = volumesOnBlock;
1499 // std::cout<<"already !! "<<std::endl;
1502 // std::cout<<"===== not found ->"<<currentHexa<<" for block "<<i<<std::endl;
1507 if(MYDEBUG) MESSAGE("SMESH_HexaFromSkin_3D::Compute END");
1511 //================================================================================
1513 * \brief Evaluate nb of hexa
1515 //================================================================================
1517 bool SMESH_HexaFromSkin_3D::Evaluate(SMESH_Mesh & aMesh,
1518 const TopoDS_Shape & aShape,
1519 MapShapeNbElems& aResMap)
1522 int nbBlocks = skin.findBlocks(aMesh);
1523 if ( nbBlocks == 0 )
1524 return error( skin.error());
1526 bool secondOrder = aMesh.NbFaces( ORDER_QUADRATIC );
1528 int entity = secondOrder ? SMDSEntity_Quad_Hexa : SMDSEntity_Hexa;
1529 vector<int>& nbByType = aResMap[ aMesh.GetSubMesh( aShape )];
1530 if ( entity >= nbByType.size() )
1531 nbByType.resize( SMDSEntity_Last, 0 );
1533 for ( int i = 0; i < nbBlocks; ++i )
1535 const _Block& block = skin.getBlock( i );
1537 int nbX = block.getSide(B_BOTTOM).getHoriSize();
1538 int nbY = block.getSide(B_BOTTOM).getVertSize();
1539 int nbZ = block.getSide(B_FRONT ).getVertSize();
1541 int nbHexa = (nbX-1) * (nbY-1) * (nbZ-1);
1542 int nbNodes = (nbX-2) * (nbY-2) * (nbZ-2);
1545 (nbX-2) * (nbY-2) * (nbZ-1) +
1546 (nbX-2) * (nbY-1) * (nbZ-2) +
1547 (nbX-1) * (nbY-2) * (nbZ-2);
1550 nbByType[ entity ] += nbHexa;
1551 nbByType[ SMDSEntity_Node ] += nbNodes;
1557 //================================================================================
1559 * \brief Abstract method must be defined but does nothing
1561 //================================================================================
1563 bool SMESH_HexaFromSkin_3D::CheckHypothesis(SMESH_Mesh&, const TopoDS_Shape&,
1564 Hypothesis_Status& aStatus)
1566 aStatus = SMESH_Hypothesis::HYP_OK;
1570 //================================================================================
1572 * \brief Abstract method must be defined but just reports an error as this
1573 * algo is not intended to work with shapes
1575 //================================================================================
1577 bool SMESH_HexaFromSkin_3D::Compute(SMESH_Mesh&, const TopoDS_Shape&)
1579 return error("Algorithm can't work with geometrical shapes");