1 // Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE
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.
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 : StdMeshers_HexaFromSkin_3D.cxx
21 // Created : Wed Jan 27 12:28:07 2010
22 // Author : Edward AGAPOV (eap)
24 #include "StdMeshers_HexaFromSkin_3D.hxx"
26 #include "SMDS_VolumeOfNodes.hxx"
27 #include "SMDS_VolumeTool.hxx"
28 #include "SMESH_Block.hxx"
29 #include "SMESH_MesherHelper.hxx"
33 //#include "utilities.h"
35 // Define error message
37 #define BAD_MESH_ERR \
38 error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh.\n" \
39 __FILE__ ":" )<<__LINE__)
41 #define BAD_MESH_ERR \
42 error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh"))
46 #define _DUMP_(msg) //cout << msg << endl
52 enum EBoxSides //!< sides of the block
54 B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_UNDEFINED
56 const char* SBoxSides[] = //!< names of block sides
58 "BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED"
60 enum EQuadEdge //!< edges of quadrangle side
62 Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT, Q_UNDEFINED
66 //================================================================================
68 * \brief return logical coordinates (i.e. min or max) of ends of edge
70 //================================================================================
72 bool getEdgeEnds(EQuadEdge edge, bool& xMax1, bool& yMax1, bool& xMax2, bool& yMax2 )
74 xMax1=0, yMax1=0, xMax2=1, yMax2=1;
77 case Q_BOTTOM: yMax2 = 0; break;
78 case Q_RIGHT: xMax1 = 1; break;
79 case Q_TOP: yMax1 = 1; break;
80 case Q_LEFT: xMax2 = 0; break;
87 //================================================================================
89 * \brief Description of node used to detect corner nodes
91 struct _NodeDescriptor
93 int nbInverseFaces, nbNodesInInverseFaces;
95 _NodeDescriptor(const SMDS_MeshNode* n): nbInverseFaces(0), nbNodesInInverseFaces(0)
99 set<const SMDS_MeshNode*> nodes;
100 SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face );
101 while ( fIt->more() )
103 const SMDS_MeshElement* face = fIt->next();
104 nodes.insert( face->begin_nodes(), face->end_nodes() );
107 nbNodesInInverseFaces = nodes.size();
110 bool operator==(const _NodeDescriptor& other) const
113 nbInverseFaces == other.nbInverseFaces &&
114 nbNodesInInverseFaces == other.nbNodesInInverseFaces;
118 //================================================================================
120 * \brief return true if a node is at block corner
122 * This check is valid for simple cases only
124 //================================================================================
126 bool isCornerNode( const SMDS_MeshNode* n )
128 return n && n->NbInverseElements( SMDSAbs_Face ) % 2;
131 //================================================================================
133 * \brief check element type
135 //================================================================================
137 bool isQuadrangle(const SMDS_MeshElement* e)
139 return ( e && e->NbNodes() == ( e->IsQuadratic() ? 8 : 4 ));
142 //================================================================================
144 * \brief return opposite node of a quadrangle face
146 //================================================================================
148 const SMDS_MeshNode* oppositeNode(const SMDS_MeshElement* quad, int iNode)
150 return quad->GetNode( (iNode+2) % 4 );
153 //================================================================================
155 * \brief Convertor of a pair of integers to a sole index
160 _Indexer( int xSize=0, int ySize=0 ): _xSize(xSize), _ySize(ySize) {}
161 int size() const { return _xSize * _ySize; }
162 int operator()(int x, int y) const { return y * _xSize + x; }
164 //================================================================================
166 * \brief Oriented convertor of a pair of integers to a sole index
168 class _OrientedIndexer : public _Indexer
171 enum OriFlags //!< types of block side orientation
173 REV_X = 1, REV_Y = 2, SWAP_XY = 4, MAX_ORI = REV_X|REV_Y|SWAP_XY
175 _OrientedIndexer( const _Indexer& indexer, const int oriFlags ):
176 _Indexer( indexer._xSize, indexer._ySize ),
177 _xSize (indexer._xSize), _ySize(indexer._ySize),
178 _xRevFun((oriFlags & REV_X) ? & reverse : & lazy),
179 _yRevFun((oriFlags & REV_Y) ? & reverse : & lazy),
180 _swapFun((oriFlags & SWAP_XY ) ? & swap : & lazy)
182 (*_swapFun)( _xSize, _ySize );
184 //!< Return index by XY
185 int operator()(int x, int y) const
187 (*_xRevFun)( x, const_cast<int&>( _xSize ));
188 (*_yRevFun)( y, const_cast<int&>( _ySize ));
190 return _Indexer::operator()(x,y);
192 //!< Return index for a corner
193 int corner(bool xMax, bool yMax) const
195 int x = xMax, y = yMax, size = 2;
196 (*_xRevFun)( x, size );
197 (*_yRevFun)( y, size );
199 return _Indexer::operator()(x ? _Indexer::_xSize-1 : 0 , y ? _Indexer::_ySize-1 : 0);
201 int xSize() const { return _xSize; }
202 int ySize() const { return _ySize; }
207 typedef void (*TFun)(int& x, int& y);
208 TFun _xRevFun, _yRevFun, _swapFun;
210 static void lazy(int&, int&) {}
211 static void reverse(int& x, int& size) { x = size - x - 1; }
212 static void swap(int& x, int& y) { std::swap(x,y); }
214 //================================================================================
216 * \brief Structure corresponding to the meshed side of block
220 vector<const SMDS_MeshNode*> _grid;
222 int _nbBlocksExpected;
225 #ifdef _DEBUG_ // want to get SIGSEGV in case of invalid index
226 #define _grid_access_(pobj, i) pobj->_grid[ ((i) < pobj->_grid.size()) ? i : int(1e100)]
228 #define _grid_access_(pobj, i) pobj->_grid[ i ]
230 //!< Return node at XY
231 const SMDS_MeshNode* getNode(int x, int y) const { return _grid_access_(this, _index( x,y ));}
233 void setNode(int x, int y, const SMDS_MeshNode* n) { _grid_access_(this, _index( x,y )) = n; }
235 SMESH_OrientedLink getEdge(EQuadEdge edge) const
237 bool x1, y1, x2, y2; getEdgeEnds( edge, x1, y1, x2, y2 );
238 return SMESH_OrientedLink( getCornerNode ( x1, y1 ), getCornerNode( x2, y2 ));
240 //!< Return a corner node
241 const SMDS_MeshNode* getCornerNode(bool isXMax, bool isYMax) const
243 return getNode( isXMax ? _index._xSize-1 : 0 , isYMax ? _index._ySize-1 : 0 );
245 const SMDS_MeshElement* getCornerFace(const SMDS_MeshNode* cornerNode) const;
246 //!< True if all blocks this side belongs to have beem found
247 bool isBound() const { return _nbBlocksExpected <= _nbBlocksFound; }
248 //!< Return coordinates of node at XY
249 gp_XYZ getXYZ(int x, int y) const { return SMESH_MeshEditor::TNodeXYZ( getNode( x, y )); }
250 //!< Return gravity center of the four corners and the middle node
255 getXYZ( _index._xSize-1, 0 ) +
256 getXYZ( 0, _index._ySize-1 ) +
257 getXYZ( _index._xSize-1, _index._ySize-1 ) +
258 getXYZ( _index._xSize/2, _index._ySize/2 );
261 //!< Return number of mesh faces
262 int getNbFaces() const { return (_index._xSize-1) * (_index._ySize-1); }
264 //================================================================================
266 * \brief _BlockSide with changed orientation
268 struct _OrientedBlockSide
271 _OrientedIndexer _index;
273 _OrientedBlockSide( _BlockSide* side=0, const int oriFlags=0 ):
274 _side(side), _index(side ? side->_index : _Indexer(), oriFlags ) {}
275 //!< return coordinates by XY
276 gp_XYZ xyz(int x, int y) const
278 return SMESH_MeshEditor::TNodeXYZ( _grid_access_(_side, _index( x, y )) );
280 //!< safely return a node by XY
281 const SMDS_MeshNode* node(int x, int y) const
283 int i = _index( x, y );
284 return ( i < 0 || i >= _side->_grid.size()) ? 0 : _side->_grid[i];
287 SMESH_OrientedLink edge(EQuadEdge edge) const
289 bool x1, y1, x2, y2; getEdgeEnds( edge, x1, y1, x2, y2 );
290 return SMESH_OrientedLink( cornerNode ( x1, y1 ), cornerNode( x2, y2 ));
292 //!< Return a corner node
293 const SMDS_MeshNode* cornerNode(bool isXMax, bool isYMax) const
295 return _grid_access_(_side, _index.corner( isXMax, isYMax ));
297 //!< return its size in nodes
298 int getHoriSize() const { return _index.xSize(); }
299 int getVertSize() const { return _index.ySize(); }
300 //!< True if _side has been initialized
301 operator bool() const { return _side; }
302 //! Direct access to _side
303 const _BlockSide* operator->() const { return _side; }
304 _BlockSide* operator->() { return _side; }
306 //================================================================================
308 * \brief Meshed skin of block
312 _OrientedBlockSide _side[6]; // 6 sides of a sub-block
314 const _OrientedBlockSide& getSide(int i) const { return _side[i]; }
315 bool hasSide( const _OrientedBlockSide& s) const
317 if ( s ) for (int i=0;i<6;++i) if ( _side[i] && _side[i]._side == s._side ) return true;
321 //================================================================================
323 * \brief Skin mesh possibly containing several meshed blocks
329 int findBlocks(SMESH_Mesh& mesh);
330 //!< return i-th block
331 const _Block& getBlock(int i) const { return _blocks[i]; }
332 //!< return error description
333 const SMESH_Comment& error() const { return _error; }
336 bool fillSide( _BlockSide& side,
337 const SMDS_MeshElement* cornerQuad,
338 const SMDS_MeshNode* cornerNode);
339 bool fillRowsUntilCorner(const SMDS_MeshElement* quad,
340 const SMDS_MeshNode* n1,
341 const SMDS_MeshNode* n2,
342 vector<const SMDS_MeshNode*>& verRow1,
343 vector<const SMDS_MeshNode*>& verRow2,
345 _OrientedBlockSide findBlockSide( EBoxSides startBlockSide,
346 EQuadEdge sharedSideEdge1,
347 EQuadEdge sharedSideEdge2);
348 //!< update own data and data of the side bound to block
349 void setSideBoundToBlock( _BlockSide& side )
351 side._nbBlocksFound++;
352 if ( side.isBound() )
354 for ( int e = 0; e < int(Q_UNDEFINED); ++e )
355 _edge2sides[ side.getEdge( (EQuadEdge) e ) ].erase( &side );
358 //!< store reason of error
359 int error(const SMESH_Comment& reason) { _error = reason; return 0; }
361 SMESH_Comment _error;
363 list< _BlockSide > _allSides;
364 vector< _Block > _blocks;
366 //map< const SMDS_MeshNode*, set< _BlockSide* > > _corner2sides;
367 map< SMESH_OrientedLink, set< _BlockSide* > > _edge2sides;
370 //================================================================================
372 * \brief Find and return number of submeshes corresponding to blocks
374 //================================================================================
376 int _Skin::findBlocks(SMESH_Mesh& mesh)
378 SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
380 // Find a node at any block corner
382 SMDS_NodeIteratorPtr nIt = meshDS->nodesIterator(/*idInceasingOrder=*/true);
383 if ( !nIt->more() ) return error("Empty mesh");
385 const SMDS_MeshNode* nCorner = 0;
386 while ( nIt->more() )
388 nCorner = nIt->next();
389 if ( isCornerNode( nCorner ))
397 // --------------------------------------------------------------------
398 // Find all block sides starting from mesh faces sharing the corner node
399 // --------------------------------------------------------------------
401 int nbFacesOnSides = 0;
402 TIDSortedElemSet cornerFaces; // corner faces of found _BlockSide's
403 list< const SMDS_MeshNode* > corners( 1, nCorner );
404 list< const SMDS_MeshNode* >::iterator corner = corners.begin();
405 while ( corner != corners.end() )
407 SMDS_ElemIteratorPtr faceIt = (*corner)->GetInverseElementIterator( SMDSAbs_Face );
408 while ( faceIt->more() )
410 const SMDS_MeshElement* face = faceIt->next();
411 if ( !cornerFaces.insert( face ).second )
412 continue; // already loaded block side
414 if ( !isQuadrangle( face ))
415 return error("Non-quadrangle elements in the input mesh");
417 if ( _allSides.empty() || !_allSides.back()._grid.empty() )
418 _allSides.push_back( _BlockSide() );
420 _BlockSide& side = _allSides.back();
421 if ( !fillSide( side, face, *corner ) )
423 if ( !_error.empty() )
428 for ( int isXMax = 0; isXMax < 2; ++isXMax )
429 for ( int isYMax = 0; isYMax < 2; ++isYMax )
431 const SMDS_MeshNode* nCorner = side.getCornerNode(isXMax,isYMax );
432 corners.push_back( nCorner );
433 cornerFaces.insert( side.getCornerFace( nCorner ));
435 for ( int e = 0; e < int(Q_UNDEFINED); ++e )
436 _edge2sides[ side.getEdge( (EQuadEdge) e ) ].insert( &side );
438 nbFacesOnSides += side.getNbFaces();
443 // find block sides of other domains if any
444 if ( corner == corners.end() && nbFacesOnSides < mesh.NbQuadrangles() )
446 while ( nIt->more() )
448 nCorner = nIt->next();
449 if ( isCornerNode( nCorner ))
450 corner = corners.insert( corner, nCorner );
452 nbFacesOnSides = mesh.NbQuadrangles();
456 if ( _allSides.empty() )
458 if ( _allSides.back()._grid.empty() )
459 _allSides.pop_back();
460 _DUMP_("Nb detected sides "<< _allSides.size());
462 // ---------------------------
463 // Organize sides into blocks
464 // ---------------------------
466 // analyse sharing of sides by blocks
467 int nbBlockSides = 0; // nb of block sides taking into account their sharing
468 list < _BlockSide >::iterator sideIt = _allSides.begin();
469 for ( ; sideIt != _allSides.end(); ++sideIt )
471 _BlockSide& side = *sideIt;
472 bool isSharedSide = true;
473 for ( int e = 0; e < int(Q_UNDEFINED) && isSharedSide; ++e )
474 isSharedSide = _edge2sides[ side.getEdge( (EQuadEdge) e ) ].size() > 2;
475 side._nbBlocksFound = 0;
476 side._nbBlocksExpected = isSharedSide ? 2 : 1;
477 nbBlockSides += side._nbBlocksExpected;
480 // find sides of each block
482 while ( nbBlockSides >= 6 )
484 // get any side not bound to all blocks it belongs to
485 sideIt = _allSides.begin();
486 while ( sideIt != _allSides.end() && sideIt->isBound())
489 // start searching for block sides from the got side
491 if ( _blocks.empty() || _blocks.back()._side[B_FRONT] )
492 _blocks.resize( _blocks.size() + 1 );
494 _Block& block = _blocks.back();
495 block._side[B_FRONT] = &(*sideIt);
496 setSideBoundToBlock( *sideIt );
499 // edges of neighbour sides of B_FRONT corresponding to front's edges
500 EQuadEdge edgeOfFront[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
501 EQuadEdge edgeToFind [4] = { Q_BOTTOM, Q_LEFT, Q_BOTTOM, Q_LEFT };
502 // TODO: first find all sides where no choice of sides is needed,
503 // then perform search with selection, which is then easier
504 for ( int i = Q_BOTTOM; ok && i <= Q_LEFT; ++i )
505 ok = ( block._side[i] = findBlockSide( B_FRONT, edgeOfFront[i], edgeToFind[i]));
507 ok = ( block._side[B_BACK] = findBlockSide( B_TOP, Q_TOP, Q_TOP ));
511 // check if just found block is same as one of previously found
513 for ( int i = 1; i < _blocks.size() && !isSame; ++i )
515 _Block& prevBlock = _blocks[i-1];
517 for ( int j = 0; j < 6 && isSame; ++j )
518 isSame = prevBlock.hasSide( block._side[ j ]);
524 // check block validity
525 int x = block._side[B_BOTTOM].getHoriSize();
526 ok = ( block._side[B_TOP ].getHoriSize() == x &&
527 block._side[B_FRONT ].getHoriSize() == x &&
528 block._side[B_BACK ].getHoriSize() == x );
529 int y = block._side[B_BOTTOM].getVertSize();
531 block._side[B_TOP ].getVertSize() == y &&
532 block._side[B_LEFT ].getHoriSize() == y &&
533 block._side[B_RIGHT ].getHoriSize() == y);
534 int z = block._side[B_FRONT].getVertSize();
536 block._side[B_BACK ].getVertSize() == z &&
537 block._side[B_LEFT ].getVertSize() == z &&
538 block._side[B_RIGHT ].getVertSize() == z);
541 // count the found sides
543 for (int i = 0; i < B_UNDEFINED; ++i )
545 _DUMP_("** Block side "<< SBoxSides[i] <<" "<< block._side[ i ]._side);
546 if ( block._side[ i ] )
548 if ( ok && i != B_FRONT)
550 setSideBoundToBlock( *block._side[ i ]._side );
553 _DUMP_("Corner 0,0 "<< block._side[ i ].cornerNode(0,0));
554 _DUMP_("Corner 1,0 "<< block._side[ i ].cornerNode(1,0));
555 _DUMP_("Corner 1,1 "<< block._side[ i ].cornerNode(1,1));
556 _DUMP_("Corner 0,1 "<< block._side[ i ].cornerNode(0,1));
560 _DUMP_("Not found"<<endl);
564 block._side[0] = block._side[1] = block._side[2] =
565 block._side[3] = block._side[4] = block._side[5] = 0;
569 if ( nbBlocks == 0 && _error.empty() )
575 //================================================================================
577 * \brief Fill block side data starting from its corner quadrangle
579 //================================================================================
581 bool _Skin::fillSide( _BlockSide& side,
582 const SMDS_MeshElement* cornerQuad,
583 const SMDS_MeshNode* nCorner)
585 // Find out size of block side mesured in nodes and by the way find two rows
586 // of nodes in two directions.
589 const SMDS_MeshElement* firstQuad = cornerQuad;
591 // get a node on block edge
592 int iCorner = firstQuad->GetNodeIndex( nCorner );
593 const SMDS_MeshNode* nOnEdge = firstQuad->GetNode( (iCorner+1) % 4);
595 // find out size of block side
596 vector<const SMDS_MeshNode*> horRow1, horRow2, verRow1, verRow2;
597 if ( !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, horRow1, horRow2, true ) ||
598 !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, verRow1, verRow2, false ))
600 nbX = horRow1.size(), nbY = verRow1.size();
603 side._index._xSize = horRow1.size();
604 side._index._ySize = verRow1.size();
605 side._grid.resize( side._index.size(), NULL );
607 for ( x = 0; x < horRow1.size(); ++x )
609 side.setNode( x, 0, horRow1[x] );
610 side.setNode( x, 1, horRow2[x] );
612 for ( y = 0; y < verRow1.size(); ++y )
614 side.setNode( 0, y, verRow1[y] );
615 side.setNode( 1, y, verRow2[y] );
618 // Find the rest nodes
620 y = 1; // y of the row to fill
621 TIDSortedElemSet emptySet, avoidSet;
624 // get next firstQuad in the next row of quadrangles
629 // o---o o o o o <- found nodes
632 int i1down, i2down, i2up;
633 const SMDS_MeshNode* n1down = side.getNode( 0, y-1 );
634 const SMDS_MeshNode* n2down = side.getNode( 1, y-1 );
635 avoidSet.clear(); avoidSet.insert( firstQuad );
636 firstQuad = SMESH_MeshEditor::FindFaceInSet( n1down, n2down, emptySet, avoidSet,
638 if ( !isQuadrangle( firstQuad ))
641 const SMDS_MeshNode* n2up = oppositeNode( firstQuad, i1down );
642 avoidSet.clear(); avoidSet.insert( firstQuad );
644 // find the rest nodes in the y-th row by faces in the row
649 const SMDS_MeshElement* quad = SMESH_MeshEditor::FindFaceInSet( n2up, n2down, emptySet,
650 avoidSet, &i2up, &i2down);
651 if ( !isQuadrangle( quad ))
654 n2up = oppositeNode( quad, i2down );
655 n2down = oppositeNode( quad, i2up );
656 avoidSet.clear(); avoidSet.insert( quad );
658 side.setNode( x, y, n2up );
662 // check side validity
664 side.getCornerFace( side.getCornerNode( 0, 0 )) &&
665 side.getCornerFace( side.getCornerNode( 1, 0 )) &&
666 side.getCornerFace( side.getCornerNode( 0, 1 )) &&
667 side.getCornerFace( side.getCornerNode( 1, 1 ));
672 //================================================================================
674 * \brief Try to find a block side adjacent to the given side by given edge
676 //================================================================================
678 _OrientedBlockSide _Skin::findBlockSide( EBoxSides startBlockSide,
679 EQuadEdge sharedSideEdge1,
680 EQuadEdge sharedSideEdge2)
682 _Block& block = _blocks.back();
683 _OrientedBlockSide& side1 = block._side[ startBlockSide ];
685 // get corner nodes of the given block edge
686 SMESH_OrientedLink edge = side1.edge( sharedSideEdge1 );
687 const SMDS_MeshNode* n1 = edge.node1();
688 const SMDS_MeshNode* n2 = edge.node2();
689 if ( edge._reversed ) swap( n1, n2 );
691 // find all sides sharing both nodes n1 and n2
692 set< _BlockSide* > sidesOnEdge = _edge2sides[ edge ]; // copy a set
694 // exclude loaded sides of block from sidesOnEdge
695 int nbLoadedSides = 0;
696 for (int i = 0; i < B_UNDEFINED; ++i )
698 if ( block._side[ i ] )
701 sidesOnEdge.erase( block._side[ i ]._side );
704 int nbSidesOnEdge = sidesOnEdge.size();
705 _DUMP_("nbSidesOnEdge "<< nbSidesOnEdge << " " << n1->GetID() << "-" << n2->GetID() );
706 if ( nbSidesOnEdge == 0 )
709 _BlockSide* foundSide = 0;
710 if ( nbSidesOnEdge == 1 /*|| nbSidesOnEdge == 2 && nbLoadedSides == 1 */)
712 foundSide = *sidesOnEdge.begin();
716 set< _BlockSide* >::iterator sideIt = sidesOnEdge.begin();
717 gp_XYZ gc(0,0,0); // gravity center of already loaded block sides
718 if ( nbLoadedSides > 1 )
720 // Find the side having more than 2 corners common with already loaded sides
722 set<const SMDS_MeshNode*> loadedCorners;
723 for (int i = 0; i < B_UNDEFINED; ++i )
724 if ( block._side[ i ] )
726 loadedCorners.insert( block._side[ i ]->getCornerNode(0,0));
727 loadedCorners.insert( block._side[ i ]->getCornerNode(1,0));
728 loadedCorners.insert( block._side[ i ]->getCornerNode(0,1));
729 loadedCorners.insert( block._side[ i ]->getCornerNode(1,1));
730 gc += block._side[ i ]._side->getGC();
734 for (; !foundSide && sideIt != sidesOnEdge.end(); ++sideIt )
736 _BlockSide* sideI = *sideIt;
737 int nbCommonCorners =
738 loadedCorners.count( sideI->getCornerNode(0,0)) +
739 loadedCorners.count( sideI->getCornerNode(1,0)) +
740 loadedCorners.count( sideI->getCornerNode(0,1)) +
741 loadedCorners.count( sideI->getCornerNode(1,1));
742 if ( nbCommonCorners > 2 )
749 // Select one of found sides most close to startBlockSide
751 gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()), p2 (n2->X(),n2->Y(),n2->Z());
752 gp_Vec p1p2( p1, p2 );
754 const SMDS_MeshElement* face1 = side1->getCornerFace( n1 );
755 gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(n1)));
756 gp_Vec side1Dir( p1, p1Op );
757 gp_Ax2 pln( p1, p1p2, side1Dir ); // plane with normal p1p2 and X dir side1Dir
758 _DUMP_(" Select adjacent for "<< side1._side << " - side dir ("
759 << side1Dir.X() << ", " << side1Dir.Y() << ", " << side1Dir.Z() << ")" );
761 map < double , _BlockSide* > angleOfSide;
762 for (sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
764 _BlockSide* sideI = *sideIt;
765 const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 );
766 gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(n1)));
767 gp_Vec sideIDir( p1, p1Op );
768 // compute angle of (sideIDir projection to pln) and (X dir of pln)
769 gp_Vec2d sideIDirProj( sideIDir * pln.XDirection(), sideIDir * pln.YDirection());
770 double angle = sideIDirProj.Angle( gp::DX2d() );
771 if ( angle < 0 ) angle += 2 * PI; // angle [0-2*PI]
772 angleOfSide.insert( make_pair( angle, sideI ));
773 _DUMP_(" "<< sideI << " - side dir ("
774 << sideIDir.X() << ", " << sideIDir.Y() << ", " << sideIDir.Z() << ")"
775 << " angle " << angle);
777 if ( nbLoadedSides == 1 )
779 double angF = angleOfSide.begin()->first, angL = angleOfSide.rbegin()->first;
780 if ( angF > PI ) angF = 2*PI - angF;
781 if ( angL > PI ) angL = 2*PI - angL;
782 foundSide = angF < angL ? angleOfSide.begin()->second : angleOfSide.rbegin()->second;
786 gp_Vec gcDir( p1, gc );
787 gp_Vec2d gcDirProj( gcDir * pln.XDirection(), gcDir * pln.YDirection());
788 double gcAngle = gcDirProj.Angle( gp::DX2d() );
789 foundSide = gcAngle < 0 ? angleOfSide.rbegin()->second : angleOfSide.begin()->second;
792 _DUMP_(" selected "<< foundSide );
795 // Orient the found side correctly
797 // corners of found side corresponding to nodes n1 and n2
798 bool xMax1, yMax1, xMax2, yMax2;
799 if ( !getEdgeEnds( sharedSideEdge2, xMax1, yMax1, xMax2, yMax2 ))
800 return error(SMESH_Comment("Internal error at ")<<__FILE__<<":"<<__LINE__),
801 _OrientedBlockSide(0);
803 for ( int ori = 0; ori < _OrientedIndexer::MAX_ORI+1; ++ori )
805 _OrientedBlockSide orientedSide( foundSide, ori );
806 const SMDS_MeshNode* n12 = orientedSide.cornerNode( xMax1, yMax1);
807 const SMDS_MeshNode* n22 = orientedSide.cornerNode( xMax2, yMax2);
808 if ( n1 == n12 && n2 == n22 )
811 error(SMESH_Comment("Failed to orient a block side found by edge ")<<sharedSideEdge1
812 << " of side " << startBlockSide
813 << " of block " << _blocks.size());
817 //================================================================================
819 * \brief: Fill rows (which are actually columns,if !alongN1N2) of nodes starting
820 * from the given quadrangle until another block corner encounters.
821 * n1 and n2 are at bottom of quad, n1 is at block corner.
823 //================================================================================
825 bool _Skin::fillRowsUntilCorner(const SMDS_MeshElement* quad,
826 const SMDS_MeshNode* n1,
827 const SMDS_MeshNode* n2,
828 vector<const SMDS_MeshNode*>& row1,
829 vector<const SMDS_MeshNode*>& row2,
830 const bool alongN1N2 )
832 const SMDS_MeshNode* corner1 = n1;
834 // Store nodes of quad in the rows and find new n1 and n2 to get
835 // the next face so that new n2 is on block edge
836 int i1 = quad->GetNodeIndex( n1 );
837 int i2 = quad->GetNodeIndex( n2 );
838 row1.clear(); row2.clear();
839 row1.push_back( n1 );
842 row1.push_back( n2 );
843 row2.push_back( oppositeNode( quad, i2 ));
844 row2.push_back( n1 = oppositeNode( quad, i1 ));
848 row2.push_back( n2 );
849 row1.push_back( n2 = oppositeNode( quad, i2 ));
850 row2.push_back( n1 = oppositeNode( quad, i1 ));
853 _NodeDescriptor nDesc( row1[1] );
854 if ( nDesc == _NodeDescriptor( row1[0] ))
855 return true; // row of 2 nodes
857 // Find the rest nodes
858 TIDSortedElemSet emptySet, avoidSet;
859 //while ( !isCornerNode( n2 ))
860 while ( nDesc == _NodeDescriptor( n2 ))
862 avoidSet.clear(); avoidSet.insert( quad );
863 quad = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2 );
864 if ( !isQuadrangle( quad ))
867 row1.push_back( n2 = oppositeNode( quad, i1 ));
868 row2.push_back( n1 = oppositeNode( quad, i2 ));
870 return n1 != corner1;
873 //================================================================================
875 * \brief Return a corner face by a corner node
877 //================================================================================
879 const SMDS_MeshElement* _BlockSide::getCornerFace(const SMDS_MeshNode* cornerNode) const
881 int x, y, isXMax, isYMax, found = 0;
882 for ( isXMax = 0; isXMax < 2; ++isXMax )
884 for ( isYMax = 0; isYMax < 2; ++isYMax )
886 x = isXMax ? _index._xSize-1 : 0;
887 y = isYMax ? _index._ySize-1 : 0;
888 found = ( getNode(x,y) == cornerNode );
893 if ( !found ) return 0;
894 int dx = isXMax ? -1 : +1;
895 int dy = isYMax ? -1 : +1;
896 const SMDS_MeshNode* n1 = getNode(x,y);
897 const SMDS_MeshNode* n2 = getNode(x+dx,y);
898 const SMDS_MeshNode* n3 = getNode(x,y+dy);
899 const SMDS_MeshNode* n4 = getNode(x+dx,y+dy);
900 return SMDS_Mesh::FindFace(n1, n2, n3, n4 );
905 //=======================================================================
906 //function : StdMeshers_HexaFromSkin_3D
908 //=======================================================================
910 StdMeshers_HexaFromSkin_3D::StdMeshers_HexaFromSkin_3D(int hypId, int studyId, SMESH_Gen* gen)
911 :SMESH_3D_Algo(hypId, studyId, gen)
913 MESSAGE("StdMeshers_HexaFromSkin_3D::StdMeshers_HexaFromSkin_3D");
914 _name = "HexaFromSkin_3D";
917 StdMeshers_HexaFromSkin_3D::~StdMeshers_HexaFromSkin_3D()
919 MESSAGE("StdMeshers_HexaFromSkin_3D::~StdMeshers_HexaFromSkin_3D");
922 //================================================================================
924 * \brief Main method, which generates hexaheda
926 //================================================================================
928 bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
931 int nbBlocks = skin.findBlocks(aMesh);
933 return error( skin.error());
935 vector< vector< const SMDS_MeshNode* > > columns;
936 int x, xSize, y, ySize, z, zSize;
939 for ( int i = 0; i < nbBlocks; ++i )
941 const _Block& block = skin.getBlock( i );
943 // ------------------------------------------
944 // Fill columns of nodes with existing nodes
945 // ------------------------------------------
947 xSize = block.getSide(B_BOTTOM).getHoriSize();
948 ySize = block.getSide(B_BOTTOM).getVertSize();
949 zSize = block.getSide(B_FRONT ).getVertSize();
950 int X = xSize - 1, Y = ySize - 1, Z = zSize - 1;
951 colIndex = _Indexer( xSize, ySize );
952 columns.resize( colIndex.size() );
954 // fill node columns by front and back box sides
955 for ( x = 0; x < xSize; ++x ) {
956 vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
957 vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
958 column0.resize( zSize );
959 column1.resize( zSize );
960 for ( z = 0; z < zSize; ++z ) {
961 column0[ z ] = block.getSide(B_FRONT).node( x, z );
962 column1[ z ] = block.getSide(B_BACK) .node( x, z );
965 // fill node columns by left and right box sides
966 for ( y = 1; y < ySize-1; ++y ) {
967 vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
968 vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
969 column0.resize( zSize );
970 column1.resize( zSize );
971 for ( z = 0; z < zSize; ++z ) {
972 column0[ z ] = block.getSide(B_LEFT) .node( y, z );
973 column1[ z ] = block.getSide(B_RIGHT).node( y, z );
976 // get nodes from top and bottom box sides
977 for ( x = 1; x < xSize-1; ++x ) {
978 for ( y = 1; y < ySize-1; ++y ) {
979 vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
980 column.resize( zSize );
981 column.front() = block.getSide(B_BOTTOM).node( x, y );
982 column.back() = block.getSide(B_TOP) .node( x, y );
986 // ----------------------------
987 // Add internal nodes of a box
988 // ----------------------------
989 // projection points of internal nodes on box subshapes by which
990 // coordinates of internal nodes are computed
991 vector<gp_XYZ> pointOnShape( SMESH_Block::ID_Shell );
993 // projections on vertices are constant
994 pointOnShape[ SMESH_Block::ID_V000 ] = block.getSide(B_BOTTOM).xyz( 0, 0 );
995 pointOnShape[ SMESH_Block::ID_V100 ] = block.getSide(B_BOTTOM).xyz( X, 0 );
996 pointOnShape[ SMESH_Block::ID_V010 ] = block.getSide(B_BOTTOM).xyz( 0, Y );
997 pointOnShape[ SMESH_Block::ID_V110 ] = block.getSide(B_BOTTOM).xyz( X, Y );
998 pointOnShape[ SMESH_Block::ID_V001 ] = block.getSide(B_TOP).xyz( 0, 0 );
999 pointOnShape[ SMESH_Block::ID_V101 ] = block.getSide(B_TOP).xyz( X, 0 );
1000 pointOnShape[ SMESH_Block::ID_V011 ] = block.getSide(B_TOP).xyz( 0, Y );
1001 pointOnShape[ SMESH_Block::ID_V111 ] = block.getSide(B_TOP).xyz( X, Y );
1003 for ( x = 1; x < xSize-1; ++x )
1005 gp_XYZ params; // normalized parameters of internal node within a unit box
1006 params.SetCoord( 1, x / double(X) );
1007 for ( y = 1; y < ySize-1; ++y )
1009 params.SetCoord( 2, y / double(Y) );
1010 // column to fill during z loop
1011 vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
1012 // projections on horizontal edges
1013 pointOnShape[ SMESH_Block::ID_Ex00 ] = block.getSide(B_BOTTOM).xyz( x, 0 );
1014 pointOnShape[ SMESH_Block::ID_Ex10 ] = block.getSide(B_BOTTOM).xyz( x, Y );
1015 pointOnShape[ SMESH_Block::ID_E0y0 ] = block.getSide(B_BOTTOM).xyz( 0, y );
1016 pointOnShape[ SMESH_Block::ID_E1y0 ] = block.getSide(B_BOTTOM).xyz( X, y );
1017 pointOnShape[ SMESH_Block::ID_Ex01 ] = block.getSide(B_TOP).xyz( x, 0 );
1018 pointOnShape[ SMESH_Block::ID_Ex11 ] = block.getSide(B_TOP).xyz( x, Y );
1019 pointOnShape[ SMESH_Block::ID_E0y1 ] = block.getSide(B_TOP).xyz( 0, y );
1020 pointOnShape[ SMESH_Block::ID_E1y1 ] = block.getSide(B_TOP).xyz( X, y );
1021 // projections on horizontal sides
1022 pointOnShape[ SMESH_Block::ID_Fxy0 ] = block.getSide(B_BOTTOM).xyz( x, y );
1023 pointOnShape[ SMESH_Block::ID_Fxy1 ] = block.getSide(B_TOP) .xyz( x, y );
1024 for ( z = 1; z < zSize-1; ++z ) // z loop
1026 params.SetCoord( 3, z / double(Z) );
1027 // projections on vertical edges
1028 pointOnShape[ SMESH_Block::ID_E00z ] = block.getSide(B_FRONT).xyz( 0, z );
1029 pointOnShape[ SMESH_Block::ID_E10z ] = block.getSide(B_FRONT).xyz( X, z );
1030 pointOnShape[ SMESH_Block::ID_E01z ] = block.getSide(B_BACK).xyz( 0, z );
1031 pointOnShape[ SMESH_Block::ID_E11z ] = block.getSide(B_BACK).xyz( X, z );
1032 // projections on vertical sides
1033 pointOnShape[ SMESH_Block::ID_Fx0z ] = block.getSide(B_FRONT).xyz( x, z );
1034 pointOnShape[ SMESH_Block::ID_Fx1z ] = block.getSide(B_BACK) .xyz( x, z );
1035 pointOnShape[ SMESH_Block::ID_F0yz ] = block.getSide(B_LEFT) .xyz( y, z );
1036 pointOnShape[ SMESH_Block::ID_F1yz ] = block.getSide(B_RIGHT).xyz( y, z );
1038 // compute internal node coordinates
1040 SMESH_Block::ShellPoint( params, pointOnShape, coords );
1041 column[ z ] = aHelper->AddNode( coords.X(), coords.Y(), coords.Z() );
1045 //cout << "----------------------------------------------------------------------"<<endl;
1046 //for ( int id = SMESH_Block::ID_V000; id < SMESH_Block::ID_Shell; ++id)
1048 // gp_XYZ p = pointOnShape[ id ];
1049 // SMESH_Block::DumpShapeID( id,cout)<<" ( "<<p.X()<<", "<<p.Y()<<", "<<p.Z()<<" )"<<endl;
1051 //cout << "Params: ( "<< params.X()<<", "<<params.Y()<<", "<<params.Z()<<" )"<<endl;
1052 //cout << "coords: ( "<< coords.X()<<", "<<coords.Y()<<", "<<coords.Z()<<" )"<<endl;
1061 // find out orientation
1062 const SMDS_MeshNode* n000 = block.getSide(B_BOTTOM).cornerNode( 0, 0 );
1063 const SMDS_MeshNode* n100 = block.getSide(B_BOTTOM).cornerNode( 1, 0 );
1064 const SMDS_MeshNode* n010 = block.getSide(B_BOTTOM).cornerNode( 0, 1 );
1065 const SMDS_MeshNode* n110 = block.getSide(B_BOTTOM).cornerNode( 1, 1 );
1066 const SMDS_MeshNode* n001 = block.getSide(B_TOP).cornerNode( 0, 0 );
1067 const SMDS_MeshNode* n101 = block.getSide(B_TOP).cornerNode( 1, 0 );
1068 const SMDS_MeshNode* n011 = block.getSide(B_TOP).cornerNode( 0, 1 );
1069 const SMDS_MeshNode* n111 = block.getSide(B_TOP).cornerNode( 1, 1 );
1070 SMDS_VolumeOfNodes probeVolume (n000,n010,n110,n100,
1071 n001,n011,n111,n101);
1072 bool isForw = SMDS_VolumeTool( &probeVolume ).IsForward();
1075 for ( x = 0; x < xSize-1; ++x ) {
1076 for ( y = 0; y < ySize-1; ++y ) {
1077 vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
1078 vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
1079 vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
1080 vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
1081 // bottom face normal of a hexa mush point outside the volume
1083 for ( z = 0; z < zSize-1; ++z )
1084 aHelper->AddVolume(col00[z], col01[z], col11[z], col10[z],
1085 col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
1087 for ( z = 0; z < zSize-1; ++z )
1088 aHelper->AddVolume(col00[z], col10[z], col11[z], col01[z],
1089 col00[z+1], col10[z+1], col11[z+1], col01[z+1]);
1097 //================================================================================
1099 * \brief Evaluate nb of hexa
1101 //================================================================================
1103 bool StdMeshers_HexaFromSkin_3D::Evaluate(SMESH_Mesh & aMesh,
1104 const TopoDS_Shape & aShape,
1105 MapShapeNbElems& aResMap)
1108 int nbBlocks = skin.findBlocks(aMesh);
1109 if ( nbBlocks == 0 )
1110 return error( skin.error());
1112 bool secondOrder = aMesh.NbFaces( ORDER_QUADRATIC );
1114 int entity = secondOrder ? SMDSEntity_Quad_Hexa : SMDSEntity_Hexa;
1115 vector<int>& nbByType = aResMap[ aMesh.GetSubMesh( aShape )];
1116 if ( entity >= nbByType.size() )
1117 nbByType.resize( SMDSEntity_Last, 0 );
1119 for ( int i = 0; i < nbBlocks; ++i )
1121 const _Block& block = skin.getBlock( i );
1123 int nbX = block.getSide(B_BOTTOM).getHoriSize();
1124 int nbY = block.getSide(B_BOTTOM).getVertSize();
1125 int nbZ = block.getSide(B_FRONT ).getVertSize();
1127 int nbHexa = (nbX-1) * (nbY-1) * (nbZ-1);
1128 int nbNodes = (nbX-2) * (nbY-2) * (nbZ-2);
1131 (nbX-2) * (nbY-2) * (nbZ-1) +
1132 (nbX-2) * (nbY-1) * (nbZ-2) +
1133 (nbX-1) * (nbY-2) * (nbZ-2);
1136 nbByType[ entity ] += nbHexa;
1137 nbByType[ SMDSEntity_Node ] += nbNodes;
1143 //================================================================================
1145 * \brief Abstract method must be defined but does nothing
1147 //================================================================================
1149 bool StdMeshers_HexaFromSkin_3D::CheckHypothesis(SMESH_Mesh&, const TopoDS_Shape&,
1150 Hypothesis_Status& aStatus)
1152 aStatus = SMESH_Hypothesis::HYP_OK;
1156 //================================================================================
1158 * \brief Abstract method must be defined but just reports an error as this
1159 * algo is not intended to work with shapes
1161 //================================================================================
1163 bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh&, const TopoDS_Shape&)
1165 return error("Algorithm can't work with geometrical shapes");