1 // Copyright (C) 2009-2023 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
39 #define BAD_MESH_ERR \
40 error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh.\n" \
41 __FILE__ ":" )<<__LINE__)
44 #define _DUMP_(msg) cout << msg << endl
49 enum EBoxSides //!< sides of the block
51 B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, NB_BLOCK_SIDES
54 const char* SBoxSides[] = //!< names of block sides
56 "BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED"
59 enum EQuadEdge //!< edges of quadrangle side
61 Q_BOTTOM = 0, Q_RIGHT, Q_TOP, Q_LEFT, NB_QUAD_SIDES
65 //================================================================================
67 * \brief return logical coordinates (i.e. min or max) of ends of edge
69 //================================================================================
71 bool getEdgeEnds(EQuadEdge edge, bool& xMax1, bool& yMax1, bool& xMax2, bool& yMax2 )
73 xMax1=0, yMax1=0, xMax2=1, yMax2=1;
76 case Q_BOTTOM: yMax2 = 0; break;
77 case Q_RIGHT: xMax1 = 1; break;
78 case Q_TOP: yMax1 = 1; break;
79 case Q_LEFT: xMax2 = 0; break;
86 //================================================================================
88 * \brief return true if a node is at block corner
90 * This check is valid for simple cases only
92 //================================================================================
94 bool isCornerNode( const SMDS_MeshNode* n )
96 int nbF = n ? n->NbInverseElements( SMDSAbs_Face ) : 1;
100 std::set<const SMDS_MeshNode*> nodesInInverseFaces;
101 SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face );
102 while ( fIt->more() )
104 const SMDS_MeshElement* face = fIt->next();
105 nodesInInverseFaces.insert( face->begin_nodes(), face->end_nodes() );
108 return nodesInInverseFaces.size() != ( 6 + (nbF/2-1)*3 );
111 //================================================================================
113 * \brief check element type
115 //================================================================================
117 bool isQuadrangle(const SMDS_MeshElement* e)
119 return ( e && e->NbCornerNodes() == 4 );
122 //================================================================================
124 * \brief return opposite node of a quadrangle face
126 //================================================================================
128 const SMDS_MeshNode* oppositeNode(const SMDS_MeshElement* quad, int iNode)
130 return quad->GetNode( (iNode+2) % 4 );
133 //================================================================================
135 * \brief Convertor of a pair of integers to a sole index
140 _Indexer( int xSize=0, int ySize=0 ): _xSize(xSize), _ySize(ySize) {}
141 int size() const { return _xSize * _ySize; }
142 int operator()(int x, int y) const { return y * _xSize + x; }
144 //================================================================================
146 * \brief Oriented convertor of a pair of integers to a sole index
148 class _OrientedIndexer : public _Indexer
151 enum OriFlags //!< types of block side orientation
153 REV_X = 1, REV_Y = 2, SWAP_XY = 4, MAX_ORI = REV_X|REV_Y|SWAP_XY
155 _OrientedIndexer( const _Indexer& indexer, const int oriFlags ):
156 _Indexer( indexer._xSize, indexer._ySize ),
157 _xSize (indexer._xSize), _ySize(indexer._ySize),
158 _xRevFun((oriFlags & REV_X) ? & reverse : & lazy),
159 _yRevFun((oriFlags & REV_Y) ? & reverse : & lazy),
160 _swapFun((oriFlags & SWAP_XY ) ? & swap : & lazy)
162 (*_swapFun)( _xSize, _ySize );
164 //!< Return index by XY
165 int operator()(int x, int y) const
167 (*_xRevFun)( x, const_cast<int&>( _xSize ));
168 (*_yRevFun)( y, const_cast<int&>( _ySize ));
170 return _Indexer::operator()(x,y);
172 //!< Return index for a corner
173 int corner(bool xMax, bool yMax) const
175 int x = xMax, y = yMax, size = 2;
176 (*_xRevFun)( x, size );
177 (*_yRevFun)( y, size );
179 return _Indexer::operator()(x ? _Indexer::_xSize-1 : 0 , y ? _Indexer::_ySize-1 : 0);
181 int xSize() const { return _xSize; }
182 int ySize() const { return _ySize; }
187 typedef void (*TFun)(int& x, int& y);
188 TFun _xRevFun, _yRevFun, _swapFun;
190 static void lazy (int&, int&) {}
191 static void reverse(int& x, int& size) { x = size - x - 1; }
192 static void swap (int& x, int& y) { std::swap(x,y); }
194 //================================================================================
196 * \brief Structure corresponding to the meshed side of block
200 std::vector<const SMDS_MeshNode*> _grid;
202 int _nbBlocksExpected;
205 #ifdef _DEBUG_ // want to get SIGSEGV in case of invalid index
206 #define _grid_access_(pobj, i) pobj->_grid[ ((i) < pobj->_grid.size()) ? i : int(1e100)]
208 #define _grid_access_(pobj, i) pobj->_grid[ i ]
210 //!< Return node at XY
211 const SMDS_MeshNode* getNode(int x, int y) const { return _grid_access_(this, _index( x,y ));}
213 void setNode(int x, int y, const SMDS_MeshNode* n) { _grid_access_(this, _index( x,y )) = n; }
215 SMESH_OrientedLink getEdge(EQuadEdge edge) const
217 bool x1, y1, x2, y2; getEdgeEnds( edge, x1, y1, x2, y2 );
218 return SMESH_OrientedLink( getCornerNode ( x1, y1 ), getCornerNode( x2, y2 ));
220 //!< Return a corner node
221 const SMDS_MeshNode* getCornerNode(bool isXMax, bool isYMax) const
223 return getNode( isXMax ? _index._xSize-1 : 0 , isYMax ? _index._ySize-1 : 0 );
225 const SMDS_MeshElement* getCornerFace(const SMDS_MeshNode* cornerNode) const;
226 //!< True if all blocks this side belongs to have been found
227 bool isBound() const { return _nbBlocksExpected <= _nbBlocksFound; }
228 //!< Return coordinates of node at XY
229 gp_XYZ getXYZ(int x, int y) const { return SMESH_TNodeXYZ( getNode( x, y )); }
230 //!< Return gravity center of the four corners and the middle node
235 getXYZ( _index._xSize-1, 0 ) +
236 getXYZ( 0, _index._ySize-1 ) +
237 getXYZ( _index._xSize-1, _index._ySize-1 ) +
238 getXYZ( _index._xSize/2, _index._ySize/2 );
241 //!< Return number of mesh faces
242 int getNbFaces() const { return (_index._xSize-1) * (_index._ySize-1); }
244 //================================================================================
246 * \brief _BlockSide with changed orientation
248 struct _OrientedBlockSide
251 _OrientedIndexer _index;
253 _OrientedBlockSide( _BlockSide* side=0, const int oriFlags=0 ):
254 _side(side), _index(side ? side->_index : _Indexer(), oriFlags ) {}
255 //!< return coordinates by XY
256 gp_XYZ xyz(int x, int y) const
258 return SMESH_TNodeXYZ( _grid_access_(_side, _index( x, y )) );
260 //!< safely return a node by XY
261 const SMDS_MeshNode* node(int x, int y) const
263 int i = _index( x, y );
264 return ( i < 0 || i >= _side->_grid.size()) ? 0 : _side->_grid[i];
267 SMESH_OrientedLink edge(EQuadEdge edge) const
269 bool x1, y1, x2, y2; getEdgeEnds( edge, x1, y1, x2, y2 );
270 return SMESH_OrientedLink( cornerNode ( x1, y1 ), cornerNode( x2, y2 ));
272 //!< Return a corner node
273 const SMDS_MeshNode* cornerNode(bool isXMax, bool isYMax) const
275 return _grid_access_(_side, _index.corner( isXMax, isYMax ));
277 //!< return its size in nodes
278 int getHoriSize() const { return _index.xSize(); }
279 int getVertSize() const { return _index.ySize(); }
280 //!< True if _side has been initialized
281 operator bool() const { return _side; }
282 //! Direct access to _side
283 const _BlockSide* operator->() const { return _side; }
284 _BlockSide* operator->() { return _side; }
286 //================================================================================
288 * \brief Meshed skin of block
292 _OrientedBlockSide _side[6]; // 6 sides of a sub-block
293 std::set<const SMDS_MeshNode*> _corners;
295 const _OrientedBlockSide& getSide(int i) const { return _side[i]; }
296 bool setSide( int i, const _OrientedBlockSide& s)
298 if (( _side[i] = s ))
300 _corners.insert( s.cornerNode(0,0));
301 _corners.insert( s.cornerNode(1,0));
302 _corners.insert( s.cornerNode(0,1));
303 _corners.insert( s.cornerNode(1,1));
307 void clear() { for (int i=0;i<6;++i) _side[i]=0; _corners.clear(); }
308 bool hasSide( const _OrientedBlockSide& s) const
310 if ( s ) for (int i=0;i<6;++i) if ( _side[i] && _side[i]._side == s._side ) return true;
313 int nbSides() const { int n=0; for (int i=0;i<6;++i) if ( _side[i] ) ++n; return n; }
314 bool isValid() const;
316 //================================================================================
318 * \brief Skin mesh possibly containing several meshed blocks
324 int findBlocks(SMESH_Mesh& mesh);
325 //!< return i-th block
326 const _Block& getBlock(int i) const { return _blocks[i]; }
327 //!< return error description
328 const SMESH_Comment& error() const { return _error; }
331 bool fillSide( _BlockSide& side,
332 const SMDS_MeshElement* cornerQuad,
333 const SMDS_MeshNode* cornerNode);
334 bool fillRowsUntilCorner(const SMDS_MeshElement* quad,
335 const SMDS_MeshNode* n1,
336 const SMDS_MeshNode* n2,
337 std::vector<const SMDS_MeshNode*>& verRow1,
338 std::vector<const SMDS_MeshNode*>& verRow2,
340 _OrientedBlockSide findBlockSide( EBoxSides startBlockSide,
341 EQuadEdge sharedSideEdge1,
342 EQuadEdge sharedSideEdge2,
343 bool withGeometricAnalysis,
344 std::set< _BlockSide* >& sidesAround);
345 //!< update own data and data of the side bound to block
346 void setSideBoundToBlock( _BlockSide& side )
348 if ( side._nbBlocksFound++, side.isBound() )
349 for ( int e = 0; e < int(NB_QUAD_SIDES); ++e )
350 _edge2sides[ side.getEdge( (EQuadEdge) e ) ].erase( &side );
352 //!< store reason of error
353 int error(const SMESH_Comment& reason) { _error = reason; return 0; }
355 SMESH_Comment _error;
357 std::list< _BlockSide > _allSides;
358 std::vector< _Block > _blocks;
360 //map< const SMDS_MeshNode*, set< _BlockSide* > > _corner2sides;
361 std::map< SMESH_OrientedLink, std::set< _BlockSide* > > _edge2sides;
364 //================================================================================
366 * \brief Find and return number of submeshes corresponding to blocks
368 //================================================================================
370 int _Skin::findBlocks(SMESH_Mesh& mesh)
372 SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
374 // Find a node at any block corner
376 SMDS_NodeIteratorPtr nIt = meshDS->nodesIterator();
377 if ( !nIt->more() ) return error("Empty mesh");
379 const SMDS_MeshNode* nCorner = 0;
380 while ( nIt->more() )
382 nCorner = nIt->next();
383 if ( isCornerNode( nCorner ))
391 // --------------------------------------------------------------------
392 // Find all block sides starting from mesh faces sharing the corner node
393 // --------------------------------------------------------------------
395 smIdType nbFacesOnSides = 0;
396 TIDSortedElemSet cornerFaces; // corner faces of found _BlockSide's
397 std::list< const SMDS_MeshNode* > corners( 1, nCorner );
398 std::list< const SMDS_MeshNode* >::iterator corner = corners.begin();
399 while ( corner != corners.end() )
401 SMDS_ElemIteratorPtr faceIt = (*corner)->GetInverseElementIterator( SMDSAbs_Face );
402 while ( faceIt->more() )
404 const SMDS_MeshElement* face = faceIt->next();
405 if ( !cornerFaces.insert( face ).second )
406 continue; // already loaded block side
408 if ( !isQuadrangle( face ))
409 return error("Non-quadrangle elements in the input mesh");
411 if ( _allSides.empty() || !_allSides.back()._grid.empty() )
412 _allSides.push_back( _BlockSide() );
414 _BlockSide& side = _allSides.back();
415 if ( !fillSide( side, face, *corner ) )
417 if ( !_error.empty() )
422 for ( int isXMax = 0; isXMax < 2; ++isXMax )
423 for ( int isYMax = 0; isYMax < 2; ++isYMax )
425 const SMDS_MeshNode* nCorner = side.getCornerNode(isXMax,isYMax );
426 corners.push_back( nCorner );
427 cornerFaces.insert( side.getCornerFace( nCorner ));
429 for ( int e = 0; e < int(NB_QUAD_SIDES); ++e )
430 _edge2sides[ side.getEdge( (EQuadEdge) e ) ].insert( &side );
432 nbFacesOnSides += side.getNbFaces();
437 // find block sides of other domains if any
438 if ( corner == corners.end() && nbFacesOnSides < mesh.NbQuadrangles() )
440 while ( nIt->more() )
442 nCorner = nIt->next();
443 if ( isCornerNode( nCorner ))
444 corner = corners.insert( corner, nCorner );
446 nbFacesOnSides = mesh.NbQuadrangles();
450 if ( _allSides.empty() )
452 if ( _allSides.back()._grid.empty() )
453 _allSides.pop_back();
454 _DUMP_("Nb detected sides "<< _allSides.size());
456 // ---------------------------
457 // Organize sides into blocks
458 // ---------------------------
460 // analyse sharing of sides by blocks and sort sides by nb of adjacent sides
461 int nbBlockSides = 0; // total nb of block sides taking into account their sharing
462 std::multimap<int, _BlockSide* > sortedSides;
464 std::list < _BlockSide >::iterator sideIt = _allSides.begin();
465 for ( ; sideIt != _allSides.end(); ++sideIt )
467 _BlockSide& side = *sideIt;
468 bool isSharedSide = true;
470 for ( int e = 0; e < int(NB_QUAD_SIDES) && isSharedSide; ++e )
472 int nbAdj = _edge2sides[ side.getEdge( (EQuadEdge) e ) ].size();
474 isSharedSide = ( nbAdj > 2 );
476 side._nbBlocksFound = 0;
477 side._nbBlocksExpected = isSharedSide ? 2 : 1;
478 nbBlockSides += side._nbBlocksExpected;
479 sortedSides.insert( std::make_pair( nbAdjacent, & side ));
483 // find sides of each block
485 while ( nbBlockSides >= 6 )
487 // get any side not bound to all blocks it belongs to
488 std::multimap<int, _BlockSide*>::iterator i_side = sortedSides.begin();
489 while ( i_side != sortedSides.end() && i_side->second->isBound())
492 // start searching for block sides from the got side
494 if ( _blocks.empty() || _blocks.back()._side[B_FRONT] )
495 _blocks.resize( _blocks.size() + 1 );
497 _Block& block = _blocks.back();
498 block.setSide( B_FRONT, i_side->second );
499 setSideBoundToBlock( *i_side->second );
502 // edges of adjacent sides of B_FRONT corresponding to front's edges
503 EQuadEdge edgeOfFront[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
504 EQuadEdge edgeOfAdj [4] = { Q_BOTTOM, Q_LEFT, Q_BOTTOM, Q_LEFT };
505 // first find all sides detectable w/o advanced analysis,
506 // then repeat the search, which then may pass without advanced analysis
507 std::set< _BlockSide* > sidesAround;
508 for ( int advAnalys = 0; advAnalys < 2; ++advAnalys )
510 // try to find 4 sides adjacent to a FRONT side
511 for ( int i = 0; (ok || !advAnalys) && i < NB_QUAD_SIDES; ++i )
512 if ( !block._side[i] )
513 ok = block.setSide( i, findBlockSide( B_FRONT, edgeOfFront[i], edgeOfAdj[i],
514 advAnalys, sidesAround));
515 // try to find a BACK side by a TOP one
516 if ( ok || !advAnalys)
517 if ( !block._side[B_BACK] && block._side[B_TOP] )
518 ok = block.setSide( B_BACK, findBlockSide( B_TOP, Q_TOP, Q_TOP,
519 advAnalys, sidesAround ));
520 if ( !advAnalys ) ok = true;
522 ok = block.isValid();
525 // check if just found block is same as one of previously found blocks
527 for ( int i = 1; i < _blocks.size() && !isSame; ++i )
528 isSame = ( block._corners == _blocks[i-1]._corners );
532 // count the found sides
533 _DUMP_(endl << "** Block " << _blocks.size() << " valid: " << block.isValid());
534 for (int i = 0; i < NB_BLOCK_SIDES; ++i )
536 _DUMP_("\tSide "<< SBoxSides[i] <<" "<< block._side[ i ]._side);
537 if ( block._side[ i ] )
539 if ( ok && i != B_FRONT)
541 setSideBoundToBlock( *block._side[ i ]._side );
544 _DUMP_("\t corners "<<
545 block._side[ i ].cornerNode(0,0)->GetID() << ", " <<
546 block._side[ i ].cornerNode(1,0)->GetID() << ", " <<
547 block._side[ i ].cornerNode(1,1)->GetID() << ", " <<
548 block._side[ i ].cornerNode(0,1)->GetID() << ", "<<endl);
552 _DUMP_("\t not found"<<endl);
560 _DUMP_("Nb found blocks "<< nbBlocks <<endl);
562 if ( nbBlocks == 0 && _error.empty() )
568 //================================================================================
570 * \brief Fill block side data starting from its corner quadrangle
572 //================================================================================
574 bool _Skin::fillSide( _BlockSide& side,
575 const SMDS_MeshElement* cornerQuad,
576 const SMDS_MeshNode* nCorner)
578 // Find out size of block side mesured in nodes and by the way find two rows
579 // of nodes in two directions.
582 const SMDS_MeshElement* firstQuad = cornerQuad;
584 // get a node on block edge
585 int iCorner = firstQuad->GetNodeIndex( nCorner );
586 const SMDS_MeshNode* nOnEdge = firstQuad->GetNode( (iCorner+1) % 4);
588 // find out size of block side
589 std::vector<const SMDS_MeshNode*> horRow1, horRow2, verRow1, verRow2;
590 if ( !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, horRow1, horRow2, true ) ||
591 !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, verRow1, verRow2, false ))
593 nbX = horRow1.size(), nbY = verRow1.size();
596 side._index._xSize = horRow1.size();
597 side._index._ySize = verRow1.size();
598 side._grid.resize( side._index.size(), NULL );
600 for ( x = 0; x < horRow1.size(); ++x )
602 side.setNode( x, 0, horRow1[x] );
603 side.setNode( x, 1, horRow2[x] );
605 for ( y = 0; y < verRow1.size(); ++y )
607 side.setNode( 0, y, verRow1[y] );
608 side.setNode( 1, y, verRow2[y] );
611 // Find the rest nodes
613 y = 1; // y of the row to fill
614 TIDSortedElemSet emptySet, avoidSet;
617 // get next firstQuad in the next row of quadrangles
622 // o---o o o o o <- found nodes
625 int i1down, i2down, i2up;
626 const SMDS_MeshNode* n1down = side.getNode( 0, y-1 );
627 const SMDS_MeshNode* n2down = side.getNode( 1, y-1 );
628 avoidSet.clear(); avoidSet.insert( firstQuad );
629 firstQuad = SMESH_MeshAlgos::FindFaceInSet( n1down, n2down, emptySet, avoidSet,
631 if ( !isQuadrangle( firstQuad ))
634 const SMDS_MeshNode* n2up = oppositeNode( firstQuad, i1down );
635 avoidSet.clear(); avoidSet.insert( firstQuad );
637 // find the rest nodes in the y-th row by faces in the row
642 const SMDS_MeshElement* quad = SMESH_MeshAlgos::FindFaceInSet( n2up, n2down, emptySet,
643 avoidSet, &i2up, &i2down);
644 if ( !isQuadrangle( quad ))
647 n2up = oppositeNode( quad, i2down );
648 n2down = oppositeNode( quad, i2up );
649 avoidSet.clear(); avoidSet.insert( quad );
651 side.setNode( x, y, n2up );
655 // check side validity
657 side.getCornerFace( side.getCornerNode( 0, 0 )) &&
658 side.getCornerFace( side.getCornerNode( 1, 0 )) &&
659 side.getCornerFace( side.getCornerNode( 0, 1 )) &&
660 side.getCornerFace( side.getCornerNode( 1, 1 ));
665 //================================================================================
667 * \brief Return true if it's possible to make a loop over corner2Sides starting
670 //================================================================================
672 bool isClosedChainOfSides( _BlockSide* startSide,
673 std::map< const SMDS_MeshNode*, std::list< _BlockSide* > > & corner2Sides )
675 // get start and end nodes
676 const SMDS_MeshNode *n1 = 0, *n2 = 0, *n;
677 for ( int y = 0; y < 2; ++y )
678 for ( int x = 0; x < 2; ++x )
680 n = startSide->getCornerNode(x,y);
681 if ( !corner2Sides.count( n )) continue;
687 if ( !n2 ) return false;
689 std::map< const SMDS_MeshNode*, std::list< _BlockSide* > >::iterator
690 c2sides = corner2Sides.find( n1 );
691 if ( c2sides == corner2Sides.end() ) return false;
693 int nbChainLinks = 1;
695 _BlockSide* prevSide = startSide;
698 // get the next side sharing n
699 std::list< _BlockSide* > & sides = c2sides->second;
700 _BlockSide* nextSide = ( sides.back() == prevSide ? sides.front() : sides.back() );
701 if ( nextSide == prevSide ) return false;
703 // find the next corner of the nextSide being in corner2Sides
706 for ( int y = 0; y < 2 && !n; ++y )
707 for ( int x = 0; x < 2; ++x )
709 n = nextSide->getCornerNode(x,y);
710 c2sides = corner2Sides.find( n );
711 if ( n == n1 || c2sides == corner2Sides.end() )
716 if ( !n ) return false;
722 return ( n == n2 && nbChainLinks == NB_QUAD_SIDES );
725 //================================================================================
727 * \brief Try to find a block side adjacent to the given side by given edge
729 //================================================================================
731 _OrientedBlockSide _Skin::findBlockSide( EBoxSides startBlockSide,
732 EQuadEdge sharedSideEdge1,
733 EQuadEdge sharedSideEdge2,
734 bool withGeometricAnalysis,
735 std::set< _BlockSide* >& sidesAround)
737 _Block& block = _blocks.back();
738 _OrientedBlockSide& side1 = block._side[ startBlockSide ];
740 // get corner nodes of the given block edge
741 SMESH_OrientedLink edge = side1.edge( sharedSideEdge1 );
742 const SMDS_MeshNode* n1 = edge.node1();
743 const SMDS_MeshNode* n2 = edge.node2();
744 if ( edge._reversed ) std::swap( n1, n2 );
746 // find all sides sharing both nodes n1 and n2
747 std::set< _BlockSide* > sidesOnEdge = _edge2sides[ edge ]; // copy a set
749 // exclude loaded sides of block from sidesOnEdge
750 for (int i = 0; i < NB_BLOCK_SIDES; ++i )
751 if ( block._side[ i ] )
752 sidesOnEdge.erase( block._side[ i ]._side );
754 int nbSidesOnEdge = sidesOnEdge.size();
755 _DUMP_("nbSidesOnEdge "<< nbSidesOnEdge << " " << n1->GetID() << "-" << n2->GetID() );
756 if ( nbSidesOnEdge == 0 )
759 _BlockSide* foundSide = 0;
760 if ( nbSidesOnEdge == 1 )
762 foundSide = *sidesOnEdge.begin();
766 std::set< _BlockSide* >::iterator sideIt = sidesOnEdge.begin();
767 int nbLoadedSides = block.nbSides();
768 if ( nbLoadedSides > 1 )
770 // Find the side having more than 2 corners common with already loaded sides
771 for (; !foundSide && sideIt != sidesOnEdge.end(); ++sideIt )
773 _BlockSide* sideI = *sideIt;
774 int nbCommonCorners =
775 block._corners.count( sideI->getCornerNode(0,0)) +
776 block._corners.count( sideI->getCornerNode(1,0)) +
777 block._corners.count( sideI->getCornerNode(0,1)) +
778 block._corners.count( sideI->getCornerNode(1,1));
779 if ( nbCommonCorners > 2 )
786 if ( !withGeometricAnalysis )
788 sidesAround.insert( sidesOnEdge.begin(), sidesOnEdge.end() );
791 if ( nbLoadedSides == 1 )
793 // Issue 0021529. There are at least 2 sides by each edge and
794 // position of block gravity center is undefined.
795 // Find a side starting from which we can walk around the startBlockSide
797 // fill in corner2Sides
798 std::map< const SMDS_MeshNode*, std::list< _BlockSide* > > corner2Sides;
799 for ( sideIt = sidesAround.begin(); sideIt != sidesAround.end(); ++sideIt )
801 _BlockSide* sideI = *sideIt;
802 corner2Sides[ sideI->getCornerNode(0,0) ].push_back( sideI );
803 corner2Sides[ sideI->getCornerNode(1,0) ].push_back( sideI );
804 corner2Sides[ sideI->getCornerNode(0,1) ].push_back( sideI );
805 corner2Sides[ sideI->getCornerNode(1,1) ].push_back( sideI );
807 // remove corners of startBlockSide from corner2Sides
808 std::set<const SMDS_MeshNode*>::iterator nIt = block._corners.begin();
809 for ( ; nIt != block._corners.end(); ++nIt )
810 corner2Sides.erase( *nIt );
813 for ( sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
815 if ( isClosedChainOfSides( *sideIt, corner2Sides ))
826 // Select one of found sides most close to startBlockSide
828 gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()), p2 (n2->X(),n2->Y(),n2->Z());
829 gp_Vec p1p2( p1, p2 );
831 const SMDS_MeshElement* face1 = side1->getCornerFace( n1 );
832 gp_XYZ p1Op = SMESH_TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(n1)));
833 gp_Vec side1Dir( p1, p1Op );
834 gp_Ax2 pln( p1, p1p2, side1Dir ); // plane with normal p1p2 and X dir side1Dir
835 _DUMP_(" Select adjacent for "<< side1._side << " - side dir ("
836 << side1Dir.X() << ", " << side1Dir.Y() << ", " << side1Dir.Z() << ")" );
838 std::map < double , _BlockSide* > angleOfSide;
839 for (sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
841 _BlockSide* sideI = *sideIt;
842 const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 );
843 gp_XYZ p1Op = SMESH_TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(n1)));
844 gp_Vec sideIDir( p1, p1Op );
845 // compute angle of (sideIDir projection to pln) and (X dir of pln)
846 gp_Vec2d sideIDirProj( sideIDir * pln.XDirection(), sideIDir * pln.YDirection());
847 double angle = sideIDirProj.Angle( gp::DX2d() );
848 if ( angle < 0 ) angle += 2. * M_PI; // angle [0-2*PI]
849 angleOfSide.insert( std::make_pair( angle, sideI ));
850 _DUMP_(" "<< sideI << " - side dir ("
851 << sideIDir.X() << ", " << sideIDir.Y() << ", " << sideIDir.Z() << ")"
852 << " angle " << angle);
855 gp_XYZ gc(0,0,0); // gravity center of already loaded block sides
856 for (int i = 0; i < NB_BLOCK_SIDES; ++i )
857 if ( block._side[ i ] )
858 gc += block._side[ i ]._side->getGC();
861 gp_Vec gcDir( p1, gc );
862 gp_Vec2d gcDirProj( gcDir * pln.XDirection(), gcDir * pln.YDirection());
863 double gcAngle = gcDirProj.Angle( gp::DX2d() );
864 foundSide = gcAngle < 0 ? angleOfSide.rbegin()->second : angleOfSide.begin()->second;
867 _DUMP_(" selected "<< foundSide );
870 // Orient the found side correctly
872 // corners of found side corresponding to nodes n1 and n2
873 bool xMax1, yMax1, xMax2, yMax2;
874 if ( !getEdgeEnds( sharedSideEdge2, xMax1, yMax1, xMax2, yMax2 ))
875 return error(SMESH_Comment("Internal error at ")<<__FILE__<<":"<<__LINE__),
876 _OrientedBlockSide(0);
878 for ( int ori = 0; ori < _OrientedIndexer::MAX_ORI+1; ++ori )
880 _OrientedBlockSide orientedSide( foundSide, ori );
881 const SMDS_MeshNode* n12 = orientedSide.cornerNode( xMax1, yMax1);
882 const SMDS_MeshNode* n22 = orientedSide.cornerNode( xMax2, yMax2);
883 if ( n1 == n12 && n2 == n22 )
886 error(SMESH_Comment("Failed to orient a block side found by edge ")<<sharedSideEdge1
887 << " of side " << startBlockSide
888 << " of block " << _blocks.size());
892 //================================================================================
894 * \brief: Fill rows (which are actually columns,if !alongN1N2) of nodes starting
895 * from the given quadrangle until another block corner encounters.
896 * n1 and n2 are at bottom of quad, n1 is at block corner.
898 //================================================================================
900 bool _Skin::fillRowsUntilCorner(const SMDS_MeshElement* quad,
901 const SMDS_MeshNode* n1,
902 const SMDS_MeshNode* n2,
903 std::vector<const SMDS_MeshNode*>& row1,
904 std::vector<const SMDS_MeshNode*>& row2,
905 const bool alongN1N2 )
907 const SMDS_MeshNode* corner1 = n1;
909 // Store nodes of quad in the rows and find new n1 and n2 to get
910 // the next face so that new n2 is on block edge
911 int i1 = quad->GetNodeIndex( n1 );
912 int i2 = quad->GetNodeIndex( n2 );
913 row1.clear(); row2.clear();
914 row1.push_back( n1 );
917 row1.push_back( n2 );
918 row2.push_back( oppositeNode( quad, i2 ));
919 row2.push_back( n1 = oppositeNode( quad, i1 ));
923 row2.push_back( n2 );
924 row1.push_back( n2 = oppositeNode( quad, i2 ));
925 row2.push_back( n1 = oppositeNode( quad, i1 ));
928 if ( isCornerNode( row1[1] ))
931 // Find the rest nodes
932 TIDSortedElemSet emptySet, avoidSet;
933 while ( !isCornerNode( n2 ) )
935 avoidSet.clear(); avoidSet.insert( quad );
936 quad = SMESH_MeshAlgos::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2 );
937 if ( !isQuadrangle( quad ))
940 row1.push_back( n2 = oppositeNode( quad, i1 ));
941 row2.push_back( n1 = oppositeNode( quad, i2 ));
943 return n1 != corner1;
946 //================================================================================
948 * \brief Return a corner face by a corner node
950 //================================================================================
952 const SMDS_MeshElement* _BlockSide::getCornerFace(const SMDS_MeshNode* cornerNode) const
954 int x, y, isXMax, isYMax, found = 0;
955 for ( isXMax = 0; isXMax < 2; ++isXMax )
957 for ( isYMax = 0; isYMax < 2; ++isYMax )
959 x = isXMax ? _index._xSize-1 : 0;
960 y = isYMax ? _index._ySize-1 : 0;
961 found = ( getNode(x,y) == cornerNode );
966 if ( !found ) return 0;
967 int dx = isXMax ? -1 : +1;
968 int dy = isYMax ? -1 : +1;
969 const SMDS_MeshNode* n1 = getNode(x,y);
970 const SMDS_MeshNode* n2 = getNode(x+dx,y);
971 const SMDS_MeshNode* n3 = getNode(x,y+dy);
972 const SMDS_MeshNode* n4 = getNode(x+dx,y+dy);
973 return SMDS_Mesh::FindFace(n1, n2, n3, n4 );
976 //================================================================================
978 * \brief Checks own validity
980 //================================================================================
982 bool _Block::isValid() const
984 bool ok = ( nbSides() == 6 );
986 // check only corners depending on side selection
987 EBoxSides adjacent[4] = { B_BOTTOM, B_RIGHT, B_TOP, B_LEFT };
988 EQuadEdge edgeAdj [4] = { Q_TOP, Q_RIGHT, Q_TOP, Q_RIGHT };
989 EQuadEdge edgeBack[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
991 for ( int i=0; ok && i < NB_QUAD_SIDES; ++i )
993 SMESH_OrientedLink eBack = _side[ B_BACK ].edge( edgeBack[i] );
994 SMESH_OrientedLink eAdja = _side[ adjacent[i] ].edge( edgeAdj[i] );
995 ok = ( eBack == eAdja );
1006 HEXA_NS::Hexa* _block2Hexa( const _Block& block,
1007 HEXA_NS::Document* doc,
1008 std::map<HEXA_NS::Vertex*, SMDS_MeshNode*> vertexNode )
1010 const SMDS_MeshNode* n000 = block.getSide(B_BOTTOM).cornerNode( 0, 0 );
1011 const SMDS_MeshNode* n100 = block.getSide(B_BOTTOM).cornerNode( 1, 0 );
1012 const SMDS_MeshNode* n010 = block.getSide(B_BOTTOM).cornerNode( 0, 1 );
1013 const SMDS_MeshNode* n110 = block.getSide(B_BOTTOM).cornerNode( 1, 1 );
1014 const SMDS_MeshNode* n001 = block.getSide(B_TOP).cornerNode( 0, 0 );
1015 const SMDS_MeshNode* n101 = block.getSide(B_TOP).cornerNode( 1, 0 );
1016 const SMDS_MeshNode* n011 = block.getSide(B_TOP).cornerNode( 0, 1 );
1017 const SMDS_MeshNode* n111 = block.getSide(B_TOP).cornerNode( 1, 1 );
1019 std::list<const SMDS_MeshNode* > nodeFromBlock;
1020 nodeFromBlock.push_back(n000);
1021 nodeFromBlock.push_back(n100);
1022 nodeFromBlock.push_back(n010);
1023 nodeFromBlock.push_back(n110);
1024 nodeFromBlock.push_back(n001);
1025 nodeFromBlock.push_back(n101);
1026 nodeFromBlock.push_back(n011);
1027 nodeFromBlock.push_back(n111);
1028 nodeFromBlock.sort();
1031 HEXA_NS::Hexa* hexa = NULL;
1032 int nHexa = doc->countUsedHexa();
1033 for (int j=0; j <nHexa; ++j ){
1034 hexa = doc->getUsedHexa(j);
1035 std::list<const SMDS_MeshNode* > nodeFromHexa;
1036 int nVx = hexa->countVertex();
1037 for ( int i=0; i <nVx; ++i ){
1038 HEXA_NS::Vertex* v = hexa->getVertex(i);
1039 const SMDS_MeshNode* n = vertexNode[v];
1040 nodeFromHexa.push_back(n);
1042 nodeFromHexa.sort();
1044 if ( nodeFromBlock == nodeFromHexa ){
1045 // std::cout << "OK block match hexa "<< hexa <<" id = "<<hexa->getId()<<std::endl;
1048 std::cout << "KO block does not match hexa "<< hexa <<" id = "<<hexa->getId()<<std::endl;
1061 //=======================================================================
1062 //function : SMESH_HexaFromSkin_3D
1064 //=======================================================================
1066 // SMESH_HexaFromSkin_3D::SMESH_HexaFromSkin_3D(int hypId, SMESH_Gen* gen)
1067 // :SMESH_3D_Algo(hypId, gen)
1069 // MESSAGE("SMESH_HexaFromSkin_3D::SMESH_HexaFromSkin_3D");
1070 // _name = "HexaFromSkin_3D";
1073 SMESH_HexaFromSkin_3D::SMESH_HexaFromSkin_3D(int hypId, SMESH_Gen* gen, HEXA_NS::Document* doc)
1074 :SMESH_3D_Algo(hypId, gen),
1077 MESSAGE("SMESH_HexaFromSkin_3D::SMESH_HexaFromSkin_3D");
1078 _name = "HexaFromSkin_3D";
1082 SMESH_HexaFromSkin_3D::~SMESH_HexaFromSkin_3D()
1084 MESSAGE("SMESH_HexaFromSkin_3D::~SMESH_HexaFromSkin_3D");
1087 //================================================================================
1089 * \brief Main method, which generates hexaheda
1091 //================================================================================
1093 bool SMESH_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
1096 int nbBlocks = skin.findBlocks(aMesh);
1097 if ( nbBlocks == 0 )
1098 return error( skin.error());
1100 std::vector< std::vector< const SMDS_MeshNode* > > columns;
1101 int x, xSize, y, ySize, z, zSize;
1104 for ( int i = 0; i < nbBlocks; ++i )
1106 const _Block& block = skin.getBlock( i );
1108 // ------------------------------------------
1109 // Fill columns of nodes with existing nodes
1110 // ------------------------------------------
1112 xSize = block.getSide(B_BOTTOM).getHoriSize();
1113 ySize = block.getSide(B_BOTTOM).getVertSize();
1114 zSize = block.getSide(B_FRONT ).getVertSize();
1115 int X = xSize - 1, Y = ySize - 1, Z = zSize - 1;
1116 colIndex = _Indexer( xSize, ySize );
1117 columns.resize( colIndex.size() );
1119 // fill node columns by front and back box sides
1120 for ( x = 0; x < xSize; ++x ) {
1121 std::vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
1122 std::vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
1123 column0.resize( zSize );
1124 column1.resize( zSize );
1125 for ( z = 0; z < zSize; ++z ) {
1126 column0[ z ] = block.getSide(B_FRONT).node( x, z );
1127 column1[ z ] = block.getSide(B_BACK) .node( x, z );
1130 // fill node columns by left and right box sides
1131 for ( y = 1; y < ySize-1; ++y ) {
1132 std::vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
1133 std::vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
1134 column0.resize( zSize );
1135 column1.resize( zSize );
1136 for ( z = 0; z < zSize; ++z ) {
1137 column0[ z ] = block.getSide(B_LEFT) .node( y, z );
1138 column1[ z ] = block.getSide(B_RIGHT).node( y, z );
1141 // get nodes from top and bottom box sides
1142 for ( x = 1; x < xSize-1; ++x ) {
1143 for ( y = 1; y < ySize-1; ++y ) {
1144 std::vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
1145 column.resize( zSize );
1146 column.front() = block.getSide(B_BOTTOM).node( x, y );
1147 column.back() = block.getSide(B_TOP) .node( x, y );
1151 // ----------------------------
1152 // Add internal nodes of a box
1153 // ----------------------------
1154 // projection points of internal nodes on box sub-shapes by which
1155 // coordinates of internal nodes are computed
1156 std::vector<gp_XYZ> pointOnShape( SMESH_Block::ID_Shell );
1158 // projections on vertices are constant
1159 pointOnShape[ SMESH_Block::ID_V000 ] = block.getSide(B_BOTTOM).xyz( 0, 0 );
1160 pointOnShape[ SMESH_Block::ID_V100 ] = block.getSide(B_BOTTOM).xyz( X, 0 );
1161 pointOnShape[ SMESH_Block::ID_V010 ] = block.getSide(B_BOTTOM).xyz( 0, Y );
1162 pointOnShape[ SMESH_Block::ID_V110 ] = block.getSide(B_BOTTOM).xyz( X, Y );
1163 pointOnShape[ SMESH_Block::ID_V001 ] = block.getSide(B_TOP).xyz( 0, 0 );
1164 pointOnShape[ SMESH_Block::ID_V101 ] = block.getSide(B_TOP).xyz( X, 0 );
1165 pointOnShape[ SMESH_Block::ID_V011 ] = block.getSide(B_TOP).xyz( 0, Y );
1166 pointOnShape[ SMESH_Block::ID_V111 ] = block.getSide(B_TOP).xyz( X, Y );
1168 for ( x = 1; x < xSize-1; ++x )
1170 gp_XYZ params; // normalized parameters of internal node within a unit box
1171 params.SetCoord( 1, x / double(X) );
1172 for ( y = 1; y < ySize-1; ++y )
1174 params.SetCoord( 2, y / double(Y) );
1175 // column to fill during z loop
1176 std::vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
1177 // projections on horizontal edges
1178 pointOnShape[ SMESH_Block::ID_Ex00 ] = block.getSide(B_BOTTOM).xyz( x, 0 );
1179 pointOnShape[ SMESH_Block::ID_Ex10 ] = block.getSide(B_BOTTOM).xyz( x, Y );
1180 pointOnShape[ SMESH_Block::ID_E0y0 ] = block.getSide(B_BOTTOM).xyz( 0, y );
1181 pointOnShape[ SMESH_Block::ID_E1y0 ] = block.getSide(B_BOTTOM).xyz( X, y );
1182 pointOnShape[ SMESH_Block::ID_Ex01 ] = block.getSide(B_TOP).xyz( x, 0 );
1183 pointOnShape[ SMESH_Block::ID_Ex11 ] = block.getSide(B_TOP).xyz( x, Y );
1184 pointOnShape[ SMESH_Block::ID_E0y1 ] = block.getSide(B_TOP).xyz( 0, y );
1185 pointOnShape[ SMESH_Block::ID_E1y1 ] = block.getSide(B_TOP).xyz( X, y );
1186 // projections on horizontal sides
1187 pointOnShape[ SMESH_Block::ID_Fxy0 ] = block.getSide(B_BOTTOM).xyz( x, y );
1188 pointOnShape[ SMESH_Block::ID_Fxy1 ] = block.getSide(B_TOP) .xyz( x, y );
1189 for ( z = 1; z < zSize-1; ++z ) // z loop
1191 params.SetCoord( 3, z / double(Z) );
1192 // projections on vertical edges
1193 pointOnShape[ SMESH_Block::ID_E00z ] = block.getSide(B_FRONT).xyz( 0, z );
1194 pointOnShape[ SMESH_Block::ID_E10z ] = block.getSide(B_FRONT).xyz( X, z );
1195 pointOnShape[ SMESH_Block::ID_E01z ] = block.getSide(B_BACK).xyz( 0, z );
1196 pointOnShape[ SMESH_Block::ID_E11z ] = block.getSide(B_BACK).xyz( X, z );
1197 // projections on vertical sides
1198 pointOnShape[ SMESH_Block::ID_Fx0z ] = block.getSide(B_FRONT).xyz( x, z );
1199 pointOnShape[ SMESH_Block::ID_Fx1z ] = block.getSide(B_BACK) .xyz( x, z );
1200 pointOnShape[ SMESH_Block::ID_F0yz ] = block.getSide(B_LEFT) .xyz( y, z );
1201 pointOnShape[ SMESH_Block::ID_F1yz ] = block.getSide(B_RIGHT).xyz( y, z );
1203 // compute internal node coordinates
1205 SMESH_Block::ShellPoint( params, pointOnShape, coords );
1206 column[ z ] = aHelper->AddNode( coords.X(), coords.Y(), coords.Z() );
1210 //cout << "----------------------------------------------------------------------"<<endl;
1211 //for ( int id = SMESH_Block::ID_V000; id < SMESH_Block::ID_Shell; ++id)
1213 // gp_XYZ p = pointOnShape[ id ];
1214 // SMESH_Block::DumpShapeID( id,cout)<<" ( "<<p.X()<<", "<<p.Y()<<", "<<p.Z()<<" )"<<endl;
1216 //cout << "Params: ( "<< params.X()<<", "<<params.Y()<<", "<<params.Z()<<" )"<<endl;
1217 //cout << "coords: ( "<< coords.X()<<", "<<coords.Y()<<", "<<coords.Z()<<" )"<<endl;
1226 // find out orientation by a least distorted hexahedron (issue 0020855);
1227 // the last is defined by evaluating sum of face normals of 8 corner hexahedrons
1228 double badness = std::numeric_limits<double>::max();
1230 for ( int xMax = 0; xMax < 2; ++xMax )
1231 for ( int yMax = 0; yMax < 2; ++yMax )
1232 for ( int zMax = 0; zMax < 2; ++zMax )
1234 x = xMax ? xSize-1 : 1;
1235 y = yMax ? ySize-1 : 1;
1236 z = zMax ? zSize-1 : 1;
1237 std::vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x-1, y-1 )];
1238 std::vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x , y-1 )];
1239 std::vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x-1, y )];
1240 std::vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x , y )];
1242 const SMDS_MeshNode* n000 = col00[z-1];
1243 const SMDS_MeshNode* n100 = col10[z-1];
1244 const SMDS_MeshNode* n010 = col01[z-1];
1245 const SMDS_MeshNode* n110 = col11[z-1];
1246 const SMDS_MeshNode* n001 = col00[z];
1247 const SMDS_MeshNode* n101 = col10[z];
1248 const SMDS_MeshNode* n011 = col01[z];
1249 const SMDS_MeshNode* n111 = col11[z];
1250 SMDS_VolumeOfNodes probeVolume (n000,n010,n110,n100,
1251 n001,n011,n111,n101);
1252 SMDS_VolumeTool volTool( &probeVolume );
1253 double Nx=0.,Ny=0.,Nz=0.;
1254 for ( int iFace = 0; iFace < volTool.NbFaces(); ++iFace )
1257 volTool.GetFaceNormal( iFace, nx,ny,nz );
1262 double quality = Nx*Nx + Ny*Ny + Nz*Nz;
1263 if ( quality < badness )
1266 isForw = volTool.IsForward();
1271 for ( x = 0; x < xSize-1; ++x ) {
1272 for ( y = 0; y < ySize-1; ++y ) {
1273 std::vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
1274 std::vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
1275 std::vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
1276 std::vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
1277 // bottom face normal of a hexa mush point outside the volume
1279 for ( z = 0; z < zSize-1; ++z )
1280 aHelper->AddVolume(col00[z], col01[z], col11[z], col10[z],
1281 col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
1283 for ( z = 0; z < zSize-1; ++z )
1284 aHelper->AddVolume(col00[z], col10[z], col11[z], col01[z],
1285 col00[z+1], col10[z+1], col11[z+1], col01[z+1]);
1294 //================================================================================
1296 * \brief Main method, which generates hexaheda
1298 //================================================================================
1299 bool SMESH_HexaFromSkin_3D::Compute( SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper,
1300 std::map<HEXA_NS::Hexa*, SMESH_HexaBlocks::SMESHVolumes>& volumesOnHexa,
1301 std::map<HEXA_NS::Vertex*, SMDS_MeshNode*> vertexNode )
1303 MESSAGE("SMESH_HexaFromSkin_3D::Compute BEGIN");
1305 int nbBlocks = skin.findBlocks(aMesh);
1306 if ( nbBlocks == 0 )
1307 return error( skin.error());
1309 std::vector< std::vector< const SMDS_MeshNode* > > columns;
1310 int x, xSize, y, ySize, z, zSize;
1313 for ( int i = 0; i < nbBlocks; ++i )
1315 const _Block& block = skin.getBlock( i );
1317 // ------------------------------------------
1318 // Fill columns of nodes with existing nodes
1319 // ------------------------------------------
1321 xSize = block.getSide(B_BOTTOM).getHoriSize();
1322 ySize = block.getSide(B_BOTTOM).getVertSize();
1323 zSize = block.getSide(B_FRONT ).getVertSize();
1324 int X = xSize - 1, Y = ySize - 1, Z = zSize - 1;
1325 colIndex = _Indexer( xSize, ySize );
1326 columns.resize( colIndex.size() );
1328 // fill node columns by front and back box sides
1329 for ( x = 0; x < xSize; ++x ) {
1330 std::vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
1331 std::vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
1332 column0.resize( zSize );
1333 column1.resize( zSize );
1334 for ( z = 0; z < zSize; ++z ) {
1335 column0[ z ] = block.getSide(B_FRONT).node( x, z );
1336 column1[ z ] = block.getSide(B_BACK) .node( x, z );
1339 // fill node columns by left and right box sides
1340 for ( y = 1; y < ySize-1; ++y ) {
1341 std::vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
1342 std::vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
1343 column0.resize( zSize );
1344 column1.resize( zSize );
1345 for ( z = 0; z < zSize; ++z ) {
1346 column0[ z ] = block.getSide(B_LEFT) .node( y, z );
1347 column1[ z ] = block.getSide(B_RIGHT).node( y, z );
1350 // get nodes from top and bottom box sides
1351 for ( x = 1; x < xSize-1; ++x ) {
1352 for ( y = 1; y < ySize-1; ++y ) {
1353 std::vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
1354 column.resize( zSize );
1355 column.front() = block.getSide(B_BOTTOM).node( x, y );
1356 column.back() = block.getSide(B_TOP) .node( x, y );
1360 // ----------------------------
1361 // Add internal nodes of a box
1362 // ----------------------------
1363 // projection points of internal nodes on box subshapes by which
1364 // coordinates of internal nodes are computed
1365 std::vector<gp_XYZ> pointOnShape( SMESH_Block::ID_Shell );
1367 // projections on vertices are constant
1368 pointOnShape[ SMESH_Block::ID_V000 ] = block.getSide(B_BOTTOM).xyz( 0, 0 );
1369 pointOnShape[ SMESH_Block::ID_V100 ] = block.getSide(B_BOTTOM).xyz( X, 0 );
1370 pointOnShape[ SMESH_Block::ID_V010 ] = block.getSide(B_BOTTOM).xyz( 0, Y );
1371 pointOnShape[ SMESH_Block::ID_V110 ] = block.getSide(B_BOTTOM).xyz( X, Y );
1372 pointOnShape[ SMESH_Block::ID_V001 ] = block.getSide(B_TOP).xyz( 0, 0 );
1373 pointOnShape[ SMESH_Block::ID_V101 ] = block.getSide(B_TOP).xyz( X, 0 );
1374 pointOnShape[ SMESH_Block::ID_V011 ] = block.getSide(B_TOP).xyz( 0, Y );
1375 pointOnShape[ SMESH_Block::ID_V111 ] = block.getSide(B_TOP).xyz( X, Y );
1377 for ( x = 1; x < xSize-1; ++x )
1379 gp_XYZ params; // normalized parameters of internal node within a unit box
1380 params.SetCoord( 1, x / double(X) );
1381 for ( y = 1; y < ySize-1; ++y )
1383 params.SetCoord( 2, y / double(Y) );
1384 // column to fill during z loop
1385 std::vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
1386 // projections on horizontal edges
1387 pointOnShape[ SMESH_Block::ID_Ex00 ] = block.getSide(B_BOTTOM).xyz( x, 0 );
1388 pointOnShape[ SMESH_Block::ID_Ex10 ] = block.getSide(B_BOTTOM).xyz( x, Y );
1389 pointOnShape[ SMESH_Block::ID_E0y0 ] = block.getSide(B_BOTTOM).xyz( 0, y );
1390 pointOnShape[ SMESH_Block::ID_E1y0 ] = block.getSide(B_BOTTOM).xyz( X, y );
1391 pointOnShape[ SMESH_Block::ID_Ex01 ] = block.getSide(B_TOP).xyz( x, 0 );
1392 pointOnShape[ SMESH_Block::ID_Ex11 ] = block.getSide(B_TOP).xyz( x, Y );
1393 pointOnShape[ SMESH_Block::ID_E0y1 ] = block.getSide(B_TOP).xyz( 0, y );
1394 pointOnShape[ SMESH_Block::ID_E1y1 ] = block.getSide(B_TOP).xyz( X, y );
1395 // projections on horizontal sides
1396 pointOnShape[ SMESH_Block::ID_Fxy0 ] = block.getSide(B_BOTTOM).xyz( x, y );
1397 pointOnShape[ SMESH_Block::ID_Fxy1 ] = block.getSide(B_TOP) .xyz( x, y );
1398 for ( z = 1; z < zSize-1; ++z ) // z loop
1400 params.SetCoord( 3, z / double(Z) );
1401 // projections on vertical edges
1402 pointOnShape[ SMESH_Block::ID_E00z ] = block.getSide(B_FRONT).xyz( 0, z );
1403 pointOnShape[ SMESH_Block::ID_E10z ] = block.getSide(B_FRONT).xyz( X, z );
1404 pointOnShape[ SMESH_Block::ID_E01z ] = block.getSide(B_BACK).xyz( 0, z );
1405 pointOnShape[ SMESH_Block::ID_E11z ] = block.getSide(B_BACK).xyz( X, z );
1406 // projections on vertical sides
1407 pointOnShape[ SMESH_Block::ID_Fx0z ] = block.getSide(B_FRONT).xyz( x, z );
1408 pointOnShape[ SMESH_Block::ID_Fx1z ] = block.getSide(B_BACK) .xyz( x, z );
1409 pointOnShape[ SMESH_Block::ID_F0yz ] = block.getSide(B_LEFT) .xyz( y, z );
1410 pointOnShape[ SMESH_Block::ID_F1yz ] = block.getSide(B_RIGHT).xyz( y, z );
1412 // compute internal node coordinates
1414 SMESH_Block::ShellPoint( params, pointOnShape, coords );
1415 column[ z ] = aHelper->AddNode( coords.X(), coords.Y(), coords.Z() );
1419 //cout << "----------------------------------------------------------------------"<<endl;
1420 //for ( int id = SMESH_Block::ID_V000; id < SMESH_Block::ID_Shell; ++id)
1422 // gp_XYZ p = pointOnShape[ id ];
1423 // SMESH_Block::DumpShapeID( id,cout)<<" ( "<<p.X()<<", "<<p.Y()<<", "<<p.Z()<<" )"<<endl;
1425 //cout << "Params: ( "<< params.X()<<", "<<params.Y()<<", "<<params.Z()<<" )"<<endl;
1426 //cout << "coords: ( "<< coords.X()<<", "<<coords.Y()<<", "<<coords.Z()<<" )"<<endl;
1435 // find out orientation
1436 const SMDS_MeshNode* n000 = block.getSide(B_BOTTOM).cornerNode( 0, 0 );
1437 const SMDS_MeshNode* n100 = block.getSide(B_BOTTOM).cornerNode( 1, 0 );
1438 const SMDS_MeshNode* n010 = block.getSide(B_BOTTOM).cornerNode( 0, 1 );
1439 const SMDS_MeshNode* n110 = block.getSide(B_BOTTOM).cornerNode( 1, 1 );
1440 const SMDS_MeshNode* n001 = block.getSide(B_TOP).cornerNode( 0, 0 );
1441 const SMDS_MeshNode* n101 = block.getSide(B_TOP).cornerNode( 1, 0 );
1442 const SMDS_MeshNode* n011 = block.getSide(B_TOP).cornerNode( 0, 1 );
1443 const SMDS_MeshNode* n111 = block.getSide(B_TOP).cornerNode( 1, 1 );
1444 SMDS_VolumeOfNodes probeVolume (n000,n010,n110,n100,
1445 n001,n011,n111,n101);
1446 bool isForw = SMDS_VolumeTool( &probeVolume ).IsForward();
1449 SMESH_HexaBlocks::SMESHVolumes volumesOnBlock; //Groups creation
1451 for ( x = 0; x < xSize-1; ++x ) {
1452 for ( y = 0; y < ySize-1; ++y ) {
1453 std::vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
1454 std::vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
1455 std::vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
1456 std::vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
1457 // bottom face normal of a hexa mush point outside the volume
1459 for ( z = 0; z < zSize-1; ++z ){
1460 SMDS_MeshVolume* newVolume =
1461 aHelper->AddVolume(col00[z], col01[z], col11[z], col10[z],
1462 col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
1463 volumesOnBlock.push_back( newVolume );
1466 for ( z = 0; z < zSize-1; ++z ){
1467 SMDS_MeshVolume* newVolume =
1468 aHelper->AddVolume(col00[z], col10[z], col11[z], col01[z],
1469 col00[z+1], col10[z+1], col11[z+1], col01[z+1]);
1470 volumesOnBlock.push_back( newVolume );
1474 // std::cout << "block i = " << i << std::endl;
1475 HEXA_NS::Hexa* currentHexa = _block2Hexa( block, _doc, vertexNode );
1476 if ( currentHexa != NULL ){
1477 // std::cout<<"===== found ->"<<currentHexa<<" for block "<<i<<std::endl;
1478 if ( volumesOnHexa.count(currentHexa)==0 ) {
1479 // std::cout<<"added!! "<<std::endl;
1480 volumesOnHexa[currentHexa] = volumesOnBlock;
1482 // std::cout<<"already !! "<<std::endl;
1485 // std::cout<<"===== not found ->"<<currentHexa<<" for block "<<i<<std::endl;
1490 MESSAGE("SMESH_HexaFromSkin_3D::Compute END");
1494 //================================================================================
1496 * \brief Evaluate nb of hexa
1498 //================================================================================
1500 bool SMESH_HexaFromSkin_3D::Evaluate(SMESH_Mesh & aMesh,
1501 const TopoDS_Shape & aShape,
1502 MapShapeNbElems& aResMap)
1505 int nbBlocks = skin.findBlocks(aMesh);
1506 if ( nbBlocks == 0 )
1507 return error( skin.error());
1509 bool secondOrder = aMesh.NbFaces( ORDER_QUADRATIC );
1511 int entity = secondOrder ? SMDSEntity_Quad_Hexa : SMDSEntity_Hexa;
1512 std::vector<smIdType>& nbByType = aResMap[ aMesh.GetSubMesh( aShape )];
1513 if ( entity >= nbByType.size() )
1514 nbByType.resize( SMDSEntity_Last, 0 );
1516 for ( int i = 0; i < nbBlocks; ++i )
1518 const _Block& block = skin.getBlock( i );
1520 int nbX = block.getSide(B_BOTTOM).getHoriSize();
1521 int nbY = block.getSide(B_BOTTOM).getVertSize();
1522 int nbZ = block.getSide(B_FRONT ).getVertSize();
1524 int nbHexa = (nbX-1) * (nbY-1) * (nbZ-1);
1525 int nbNodes = (nbX-2) * (nbY-2) * (nbZ-2);
1528 (nbX-2) * (nbY-2) * (nbZ-1) +
1529 (nbX-2) * (nbY-1) * (nbZ-2) +
1530 (nbX-1) * (nbY-2) * (nbZ-2);
1533 nbByType[ entity ] += nbHexa;
1534 nbByType[ SMDSEntity_Node ] += nbNodes;
1540 //================================================================================
1542 * \brief Abstract method must be defined but does nothing
1544 //================================================================================
1546 bool SMESH_HexaFromSkin_3D::CheckHypothesis(SMESH_Mesh&, const TopoDS_Shape&,
1547 Hypothesis_Status& aStatus)
1549 aStatus = SMESH_Hypothesis::HYP_OK;
1553 //================================================================================
1555 * \brief Abstract method must be defined but just reports an error as this
1556 * algo is not intended to work with shapes
1558 //================================================================================
1560 bool SMESH_HexaFromSkin_3D::Compute(SMESH_Mesh&, const TopoDS_Shape&)
1562 return error("Algorithm can't work with geometrical shapes");