1 // Copyright (C) 2007-2020 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, 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 : 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 "SMESHDS_Mesh.hxx"
29 #include "SMESH_Block.hxx"
30 #include "SMESH_Indexer.hxx"
31 #include "SMESH_Mesh.hxx"
32 #include "SMESH_MeshAlgos.hxx"
33 #include "SMESH_MesherHelper.hxx"
41 // Define error message and _MYDEBUG_ if needed
43 #define BAD_MESH_ERR \
44 error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh.\n" \
45 __FILE__ ":" )<<__LINE__)
48 #define BAD_MESH_ERR \
49 error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh"))
55 #define _DUMP_(msg) cout << msg << endl
63 // typedefs for struct's moved to SMESHUtils
64 typedef SMESH_Indexer _Indexer;
65 typedef SMESH_OrientedIndexer _OrientedIndexer;
67 enum EBoxSides //!< sides of the block
69 B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, NB_BLOCK_SIDES
72 const char* SBoxSides[] = //!< names of block sides -- needed for DEBUG only
74 "BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED"
77 enum EQuadEdge //!< edges of quadrangle side
79 Q_BOTTOM = 0, Q_RIGHT, Q_TOP, Q_LEFT, NB_QUAD_SIDES
83 //================================================================================
85 * \brief return logical coordinates (i.e. min or max) of ends of edge
87 //================================================================================
89 bool getEdgeEnds(EQuadEdge edge, bool& xMax1, bool& yMax1, bool& xMax2, bool& yMax2 )
91 xMax1=0, yMax1=0, xMax2=1, yMax2=1;
94 case Q_BOTTOM: yMax2 = 0; break;
95 case Q_RIGHT: xMax1 = 1; break;
96 case Q_TOP: yMax1 = 1; break;
97 case Q_LEFT: xMax2 = 0; break;
104 //================================================================================
106 * \brief return true if a node is at block corner
108 * This check is valid for simple cases only
110 //================================================================================
112 bool isCornerNode( const SMDS_MeshNode* n )
114 int nbF = n ? n->NbInverseElements( SMDSAbs_Face ) : 1;
118 set<const SMDS_MeshNode*> nodesInInverseFaces;
119 SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator( SMDSAbs_Face );
120 while ( fIt->more() )
122 const SMDS_MeshElement* face = fIt->next();
123 nodesInInverseFaces.insert( face->begin_nodes(), face->end_nodes() );
126 return (int)nodesInInverseFaces.size() != ( 6 + (nbF/2-1)*3 );
129 //================================================================================
131 * \brief check element type
133 //================================================================================
135 bool isQuadrangle(const SMDS_MeshElement* e)
137 return ( e && e->NbCornerNodes() == 4 );
140 //================================================================================
142 * \brief return opposite node of a quadrangle face
144 //================================================================================
146 const SMDS_MeshNode* oppositeNode(const SMDS_MeshElement* quad, int iNode)
148 return quad->GetNode( (iNode+2) % 4 );
151 //================================================================================
153 * \brief Structure corresponding to the meshed side of block
157 vector<const SMDS_MeshNode*> _grid;
159 int _nbBlocksExpected;
162 #ifdef _DEBUG_ // want to get SIGSEGV in case of invalid index
163 #define _grid_access_(pobj, i) pobj->_grid[ ((i) < pobj->_grid.size()) ? i : int(1e100)]
165 #define _grid_access_(pobj, i) pobj->_grid[ i ]
167 //!< Return node at XY
168 const SMDS_MeshNode* getNode(int x, int y) const { return _grid_access_(this, _index( x,y ));}
170 void setNode(int x, int y, const SMDS_MeshNode* n) { _grid_access_(this, _index( x,y )) = n; }
172 SMESH_OrientedLink getEdge(EQuadEdge edge) const
174 bool x1, y1, x2, y2; getEdgeEnds( edge, x1, y1, x2, y2 );
175 return SMESH_OrientedLink( getCornerNode ( x1, y1 ), getCornerNode( x2, y2 ));
177 //!< Return a corner node
178 const SMDS_MeshNode* getCornerNode(bool isXMax, bool isYMax) const
180 return getNode( isXMax ? _index._xSize-1 : 0 , isYMax ? _index._ySize-1 : 0 );
182 const SMDS_MeshElement* getCornerFace(const SMDS_MeshNode* cornerNode) const;
183 //!< True if all blocks this side belongs to have been found
184 bool isBound() const { return _nbBlocksExpected <= _nbBlocksFound; }
185 //!< Return coordinates of node at XY
186 gp_XYZ getXYZ(int x, int y) const { return SMESH_TNodeXYZ( getNode( x, y )); }
187 //!< Return gravity center of the four corners and the middle node
192 getXYZ( _index._xSize-1, 0 ) +
193 getXYZ( 0, _index._ySize-1 ) +
194 getXYZ( _index._xSize-1, _index._ySize-1 ) +
195 getXYZ( _index._xSize/2, _index._ySize/2 );
198 //!< Return number of mesh faces
199 int getNbFaces() const { return (_index._xSize-1) * (_index._ySize-1); }
201 //================================================================================
203 * \brief _BlockSide with changed orientation
205 struct _OrientedBlockSide
208 _OrientedIndexer _index;
210 _OrientedBlockSide( _BlockSide* side=0, const int oriFlags=0 ):
211 _side(side), _index(side ? side->_index : _Indexer(), oriFlags ) {}
212 //!< return coordinates by XY
213 gp_XYZ xyz(int x, int y) const
215 return SMESH_TNodeXYZ( _grid_access_(_side, _index( x, y )) );
217 //!< safely return a node by XY
218 const SMDS_MeshNode* node(int x, int y) const
220 size_t i = _index( x, y );
221 return ( i >= _side->_grid.size() ) ? 0 : _side->_grid[i];
224 SMESH_OrientedLink edge(EQuadEdge edge) const
226 bool x1, y1, x2, y2; getEdgeEnds( edge, x1, y1, x2, y2 );
227 return SMESH_OrientedLink( cornerNode ( x1, y1 ), cornerNode( x2, y2 ));
229 //!< Return a corner node
230 const SMDS_MeshNode* cornerNode(bool isXMax, bool isYMax) const
232 return _grid_access_(_side, _index.corner( isXMax, isYMax ));
234 //!< return its size in nodes
235 int getHoriSize() const { return _index.xSize(); }
236 int getVertSize() const { return _index.ySize(); }
237 //!< True if _side has been initialized
238 operator bool() const { return _side; }
239 //! Direct access to _side
240 const _BlockSide* operator->() const { return _side; }
241 _BlockSide* operator->() { return _side; }
243 //================================================================================
245 * \brief Meshed skin of block
249 _OrientedBlockSide _side[6]; // 6 sides of a sub-block
250 set<const SMDS_MeshNode*> _corners;
252 const _OrientedBlockSide& getSide(int i) const { return _side[i]; }
253 bool setSide( int i, const _OrientedBlockSide& s)
255 if (( _side[i] = s ))
257 _corners.insert( s.cornerNode(0,0));
258 _corners.insert( s.cornerNode(1,0));
259 _corners.insert( s.cornerNode(0,1));
260 _corners.insert( s.cornerNode(1,1));
264 void clear() { for (int i=0;i<6;++i) _side[i]=0; _corners.clear(); }
265 bool hasSide( const _OrientedBlockSide& s) const
267 if ( s ) for (int i=0;i<6;++i) if ( _side[i] && _side[i]._side == s._side ) return true;
270 int nbSides() const { int n=0; for (int i=0;i<6;++i) if ( _side[i] ) ++n; return n; }
271 bool isValid() const;
273 //================================================================================
275 * \brief Skin mesh possibly containing several meshed blocks
281 int findBlocks(SMESH_Mesh& mesh);
282 //!< return i-th block
283 const _Block& getBlock(int i) const { return _blocks[i]; }
284 //!< return error description
285 const SMESH_Comment& error() const { return _error; }
288 bool fillSide( _BlockSide& side,
289 const SMDS_MeshElement* cornerQuad,
290 const SMDS_MeshNode* cornerNode);
291 bool fillRowsUntilCorner(const SMDS_MeshElement* quad,
292 const SMDS_MeshNode* n1,
293 const SMDS_MeshNode* n2,
294 vector<const SMDS_MeshNode*>& verRow1,
295 vector<const SMDS_MeshNode*>& verRow2,
297 _OrientedBlockSide findBlockSide( EBoxSides startBlockSide,
298 EQuadEdge sharedSideEdge1,
299 EQuadEdge sharedSideEdge2,
300 bool withGeometricAnalysis,
301 set< _BlockSide* >& sidesAround);
302 //!< update own data and data of the side bound to block
303 void setSideBoundToBlock( _BlockSide& side )
305 if ( side._nbBlocksFound++, side.isBound() )
306 for ( int e = 0; e < int(NB_QUAD_SIDES); ++e )
307 _edge2sides[ side.getEdge( (EQuadEdge) e ) ].erase( &side );
309 //!< store reason of error
310 int error(const SMESH_Comment& reason) { _error = reason; return 0; }
312 SMESH_Comment _error;
314 list< _BlockSide > _allSides;
315 vector< _Block > _blocks;
317 //map< const SMDS_MeshNode*, set< _BlockSide* > > _corner2sides;
318 map< SMESH_OrientedLink, set< _BlockSide* > > _edge2sides;
321 //================================================================================
323 * \brief Find blocks and return their number
325 //================================================================================
327 int _Skin::findBlocks(SMESH_Mesh& mesh)
329 SMESHDS_Mesh* meshDS = mesh.GetMeshDS();
331 // Find a node at any block corner
333 SMDS_NodeIteratorPtr nIt = meshDS->nodesIterator();
334 if ( !nIt->more() ) return error("Empty mesh");
336 const SMDS_MeshNode* nCorner = 0;
337 while ( nIt->more() )
339 nCorner = nIt->next();
340 if ( isCornerNode( nCorner ))
348 // --------------------------------------------------------------------
349 // Find all block sides starting from mesh faces sharing the corner node
350 // --------------------------------------------------------------------
352 int nbFacesOnSides = 0;
353 TIDSortedElemSet cornerFaces; // corner faces of found _BlockSide's
354 list< const SMDS_MeshNode* > corners( 1, nCorner );
355 list< const SMDS_MeshNode* >::iterator corner = corners.begin();
356 while ( corner != corners.end() )
358 SMDS_ElemIteratorPtr faceIt = (*corner)->GetInverseElementIterator( SMDSAbs_Face );
359 while ( faceIt->more() )
361 const SMDS_MeshElement* face = faceIt->next();
362 if ( !cornerFaces.insert( face ).second )
363 continue; // already loaded block side
365 if ( !isQuadrangle( face ))
366 return error("Non-quadrangle elements in the input mesh");
368 if ( _allSides.empty() || !_allSides.back()._grid.empty() )
369 _allSides.push_back( _BlockSide() );
371 _BlockSide& side = _allSides.back();
372 if ( !fillSide( side, face, *corner ))
374 if ( !_error.empty() )
379 for ( int isXMax = 0; isXMax < 2; ++isXMax )
380 for ( int isYMax = 0; isYMax < 2; ++isYMax )
382 const SMDS_MeshNode* nCorner = side.getCornerNode(isXMax,isYMax );
383 corners.push_back( nCorner );
384 cornerFaces.insert( side.getCornerFace( nCorner ));
386 for ( int e = 0; e < int(NB_QUAD_SIDES); ++e )
387 _edge2sides[ side.getEdge( (EQuadEdge) e ) ].insert( &side );
389 nbFacesOnSides += side.getNbFaces();
394 // find block sides of other domains if any
395 if ( corner == corners.end() && nbFacesOnSides < mesh.NbQuadrangles() )
397 while ( nIt->more() )
399 nCorner = nIt->next();
400 if ( isCornerNode( nCorner ))
401 corner = corners.insert( corner, nCorner );
403 nbFacesOnSides = mesh.NbQuadrangles();
407 if ( _allSides.empty() )
409 if ( _allSides.back()._grid.empty() )
410 _allSides.pop_back();
411 _DUMP_("Nb detected sides "<< _allSides.size());
413 // ---------------------------
414 // Organize sides into blocks
415 // ---------------------------
417 // analyse sharing of sides by blocks and sort sides by nb of adjacent sides
418 int nbBlockSides = 0; // total nb of block sides taking into account their sharing
419 multimap<int, _BlockSide* > sortedSides;
421 list < _BlockSide >::iterator sideIt = _allSides.begin();
422 for ( ; sideIt != _allSides.end(); ++sideIt )
424 _BlockSide& side = *sideIt;
425 bool isSharedSide = true;
427 for ( int e = 0; e < int(NB_QUAD_SIDES) && isSharedSide; ++e )
429 int nbAdj = _edge2sides[ side.getEdge( (EQuadEdge) e ) ].size();
431 isSharedSide = ( nbAdj > 2 );
433 side._nbBlocksFound = 0;
434 side._nbBlocksExpected = isSharedSide ? 2 : 1;
435 nbBlockSides += side._nbBlocksExpected;
436 sortedSides.insert( make_pair( nbAdjacent, & side ));
440 // find sides of each block
442 while ( nbBlockSides >= 6 )
444 // get any side not bound to all blocks it belongs to
445 multimap<int, _BlockSide*>::iterator i_side = sortedSides.begin();
446 while ( i_side != sortedSides.end() && i_side->second->isBound())
449 // start searching for block sides from the got side
451 if ( _blocks.empty() || _blocks.back()._side[B_FRONT] )
452 _blocks.resize( _blocks.size() + 1 );
454 _Block& block = _blocks.back();
455 block.setSide( B_FRONT, i_side->second );
456 setSideBoundToBlock( *i_side->second );
459 // edges of adjacent sides of B_FRONT corresponding to front's edges
460 EQuadEdge edgeOfFront[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
461 EQuadEdge edgeOfAdj [4] = { Q_BOTTOM, Q_LEFT, Q_BOTTOM, Q_LEFT };
462 // first find all sides detectable w/o advanced analysis,
463 // then repeat the search, which then may pass without advanced analysis
464 set< _BlockSide* > sidesAround;
465 for ( int advAnalys = 0; advAnalys < 2; ++advAnalys )
467 // try to find 4 sides adjacent to a FRONT side
468 for ( int i = 0; (ok || !advAnalys) && i < NB_QUAD_SIDES; ++i )
469 if ( !block._side[i] )
470 ok = block.setSide( i, findBlockSide( B_FRONT, edgeOfFront[i], edgeOfAdj[i],
471 advAnalys, sidesAround));
472 // try to find a BACK side by a TOP one
473 if ( ok || !advAnalys )
474 if ( !block._side[B_BACK] && block._side[B_TOP] )
475 ok = block.setSide( B_BACK, findBlockSide( B_TOP, Q_TOP, Q_TOP,
476 advAnalys, sidesAround ));
477 if ( !advAnalys ) ok = true;
479 ok = block.isValid();
482 // check if just found block is same as one of previously found blocks
484 for ( size_t i = 1; i < _blocks.size() && !isSame; ++i )
485 isSame = ( block._corners == _blocks[i-1]._corners );
489 // count the found sides
490 _DUMP_(endl << "** Block " << _blocks.size() << " valid: " << block.isValid());
491 for (int i = 0; i < NB_BLOCK_SIDES; ++i )
493 _DUMP_("\tSide "<< SBoxSides[i] <<" "<< block._side[ i ]._side);
494 if ( block._side[ i ] )
496 if ( ok && i != B_FRONT)
498 setSideBoundToBlock( *block._side[ i ]._side );
501 _DUMP_("\t corners "<<
502 block._side[ i ].cornerNode(0,0)->GetID() << ", " <<
503 block._side[ i ].cornerNode(1,0)->GetID() << ", " <<
504 block._side[ i ].cornerNode(1,1)->GetID() << ", " <<
505 block._side[ i ].cornerNode(0,1)->GetID() << ", "<<endl);
509 _DUMP_("\t not found"<<endl);
517 _DUMP_("Nb found blocks "<< nbBlocks <<endl);
519 if ( nbBlocks == 0 && _error.empty() )
525 //================================================================================
527 * \brief Fill block side data starting from its corner quadrangle
529 //================================================================================
531 bool _Skin::fillSide( _BlockSide& side,
532 const SMDS_MeshElement* cornerQuad,
533 const SMDS_MeshNode* nCorner)
535 // Find out size of block side measured in nodes and by the way find two rows
536 // of nodes in two directions.
539 const SMDS_MeshElement* firstQuad = cornerQuad;
541 // get a node on block edge
542 int iCorner = firstQuad->GetNodeIndex( nCorner );
543 const SMDS_MeshNode* nOnEdge = firstQuad->GetNode( (iCorner+1) % 4);
545 // find out size of block side
546 vector<const SMDS_MeshNode*> horRow1, horRow2, verRow1, verRow2;
547 if ( !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, horRow1, horRow2, true ) ||
548 !fillRowsUntilCorner( firstQuad, nCorner, nOnEdge, verRow1, verRow2, false ))
550 nbX = horRow1.size(), nbY = verRow1.size();
553 side._index._xSize = horRow1.size();
554 side._index._ySize = verRow1.size();
555 side._grid.resize( side._index.size(), NULL );
557 for ( x = 0; x < nbX; ++x )
559 side.setNode( x, 0, horRow1[x] );
560 side.setNode( x, 1, horRow2[x] );
562 for ( y = 0; y < nbY; ++y )
564 side.setNode( 0, y, verRow1[y] );
565 side.setNode( 1, y, verRow2[y] );
568 // Find the rest nodes
570 y = 1; // y of the row to fill
571 TIDSortedElemSet emptySet, avoidSet;
574 // get next firstQuad in the next row of quadrangles
579 // o---o o o o o <- found nodes
582 int i1down, i2down, i2up;
583 const SMDS_MeshNode* n1down = side.getNode( 0, y-1 );
584 const SMDS_MeshNode* n2down = side.getNode( 1, y-1 );
585 avoidSet.clear(); avoidSet.insert( firstQuad );
586 firstQuad = SMESH_MeshAlgos::FindFaceInSet( n1down, n2down, emptySet, avoidSet,
588 if ( !isQuadrangle( firstQuad ))
591 const SMDS_MeshNode* n2up = oppositeNode( firstQuad, i1down );
592 avoidSet.clear(); avoidSet.insert( firstQuad );
594 // find the rest nodes in the y-th row by faces in the row
599 const SMDS_MeshElement* quad = SMESH_MeshAlgos::FindFaceInSet( n2up, n2down, emptySet,
600 avoidSet, &i2up, &i2down);
601 if ( !isQuadrangle( quad ))
604 n2up = oppositeNode( quad, i2down );
605 n2down = oppositeNode( quad, i2up );
606 avoidSet.clear(); avoidSet.insert( quad );
608 side.setNode( x, y, n2up );
612 // check side validity
614 side.getCornerFace( side.getCornerNode( 0, 0 )) &&
615 side.getCornerFace( side.getCornerNode( 1, 0 )) &&
616 side.getCornerFace( side.getCornerNode( 0, 1 )) &&
617 side.getCornerFace( side.getCornerNode( 1, 1 ));
622 //================================================================================
624 * \brief Return true if it's possible to make a loop over corner2Sides starting
627 //================================================================================
629 bool isClosedChainOfSides( _BlockSide* startSide,
630 map< const SMDS_MeshNode*, list< _BlockSide* > > & corner2Sides )
632 // get start and end nodes
633 const SMDS_MeshNode *n1 = 0, *n2 = 0, *n;
634 for ( int y = 0; y < 2; ++y )
635 for ( int x = 0; x < 2; ++x )
637 n = startSide->getCornerNode(x,y);
638 if ( !corner2Sides.count( n )) continue;
644 if ( !n2 ) return false;
646 map< const SMDS_MeshNode*, list< _BlockSide* > >::iterator
647 c2sides = corner2Sides.find( n1 );
648 if ( c2sides == corner2Sides.end() ) return false;
650 int nbChainLinks = 1;
652 _BlockSide* prevSide = startSide;
655 // get the next side sharing n
656 list< _BlockSide* > & sides = c2sides->second;
657 _BlockSide* nextSide = ( sides.back() == prevSide ? sides.front() : sides.back() );
658 if ( nextSide == prevSide ) return false;
660 // find the next corner of the nextSide being in corner2Sides
663 for ( int y = 0; y < 2 && !n; ++y )
664 for ( int x = 0; x < 2; ++x )
666 n = nextSide->getCornerNode(x,y);
667 c2sides = corner2Sides.find( n );
668 if ( n == n1 || c2sides == corner2Sides.end() )
673 if ( !n ) return false;
677 if ( ++nbChainLinks > NB_QUAD_SIDES )
681 return ( n == n2 && nbChainLinks == NB_QUAD_SIDES );
684 //================================================================================
686 * \brief Try to find a block side adjacent to the given side by given edge
688 //================================================================================
690 _OrientedBlockSide _Skin::findBlockSide( EBoxSides startBlockSide,
691 EQuadEdge sharedSideEdge1,
692 EQuadEdge sharedSideEdge2,
693 bool withGeometricAnalysis,
694 set< _BlockSide* >& sidesAround)
696 _Block& block = _blocks.back();
697 _OrientedBlockSide& side1 = block._side[ startBlockSide ];
699 // get corner nodes of the given block edge
700 SMESH_OrientedLink edge = side1.edge( sharedSideEdge1 );
701 const SMDS_MeshNode* n1 = edge.node1();
702 const SMDS_MeshNode* n2 = edge.node2();
703 if ( edge._reversed ) swap( n1, n2 );
705 // find all sides sharing both nodes n1 and n2
706 set< _BlockSide* > sidesOnEdge = _edge2sides[ edge ]; // copy a set
708 // exclude loaded sides of block from sidesOnEdge
709 for (int i = 0; i < NB_BLOCK_SIDES; ++i )
710 if ( block._side[ i ] )
711 sidesOnEdge.erase( block._side[ i ]._side );
713 int nbSidesOnEdge = sidesOnEdge.size();
714 _DUMP_("nbSidesOnEdge "<< nbSidesOnEdge << " " << n1->GetID() << "-" << n2->GetID() );
715 if ( nbSidesOnEdge == 0 )
718 _BlockSide* foundSide = 0;
719 if ( nbSidesOnEdge == 1 )
721 foundSide = *sidesOnEdge.begin();
725 set< _BlockSide* >::iterator sideIt = sidesOnEdge.begin();
726 int nbLoadedSides = block.nbSides();
727 if ( nbLoadedSides > 1 )
729 // Find the side having more than 2 corners common with already loaded sides
730 for (; !foundSide && sideIt != sidesOnEdge.end(); ++sideIt )
732 _BlockSide* sideI = *sideIt;
733 int nbCommonCorners =
734 block._corners.count( sideI->getCornerNode(0,0)) +
735 block._corners.count( sideI->getCornerNode(1,0)) +
736 block._corners.count( sideI->getCornerNode(0,1)) +
737 block._corners.count( sideI->getCornerNode(1,1));
738 if ( nbCommonCorners > 2 )
745 if ( !withGeometricAnalysis )
747 sidesAround.insert( sidesOnEdge.begin(), sidesOnEdge.end() );
750 if ( nbLoadedSides == 1 )
752 // Issue 0021529. There are at least 2 sides by each edge and
753 // position of block gravity center is undefined.
754 // Find a side starting from which we can walk around the startBlockSide
756 // fill in corner2Sides
757 map< const SMDS_MeshNode*, list< _BlockSide* > > corner2Sides;
758 for ( sideIt = sidesAround.begin(); sideIt != sidesAround.end(); ++sideIt )
760 _BlockSide* sideI = *sideIt;
761 corner2Sides[ sideI->getCornerNode(0,0) ].push_back( sideI );
762 corner2Sides[ sideI->getCornerNode(1,0) ].push_back( sideI );
763 corner2Sides[ sideI->getCornerNode(0,1) ].push_back( sideI );
764 corner2Sides[ sideI->getCornerNode(1,1) ].push_back( sideI );
766 // remove corners of startBlockSide from corner2Sides
767 set<const SMDS_MeshNode*>::iterator nIt = block._corners.begin();
768 for ( ; nIt != block._corners.end(); ++nIt )
769 corner2Sides.erase( *nIt );
772 for ( sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
774 if ( isClosedChainOfSides( *sideIt, corner2Sides ))
785 // Select one of found sides most close to startBlockSide
787 gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()), p2 (n2->X(),n2->Y(),n2->Z());
788 gp_Vec p1p2( p1, p2 );
790 const SMDS_MeshElement* face1 = side1->getCornerFace( n1 );
791 gp_XYZ p1Op = SMESH_TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(n1)));
792 gp_Vec side1Dir( p1, p1Op );
793 gp_Ax2 pln( p1, p1p2, side1Dir ); // plane with normal p1p2 and X dir side1Dir
794 _DUMP_(" Select adjacent for "<< side1._side << " - side dir ("
795 << side1Dir.X() << ", " << side1Dir.Y() << ", " << side1Dir.Z() << ")" );
797 map < double , _BlockSide* > angleOfSide;
798 for (sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
800 _BlockSide* sideI = *sideIt;
801 const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 );
802 gp_XYZ p1Op = SMESH_TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(n1)));
803 gp_Vec sideIDir( p1, p1Op );
804 // compute angle of (sideIDir projection to pln) and (X dir of pln)
805 gp_Vec2d sideIDirProj( sideIDir * pln.XDirection(), sideIDir * pln.YDirection());
806 double angle = sideIDirProj.Angle( gp::DX2d() );
807 if ( angle < 0 ) angle += 2. * M_PI; // angle [0-2*PI]
808 angleOfSide.insert( make_pair( angle, sideI ));
809 _DUMP_(" "<< sideI << " - side dir ("
810 << sideIDir.X() << ", " << sideIDir.Y() << ", " << sideIDir.Z() << ")"
811 << " angle " << angle);
814 gp_XYZ gc(0,0,0); // gravity center of already loaded block sides
815 for (int i = 0; i < NB_BLOCK_SIDES; ++i )
816 if ( block._side[ i ] )
817 gc += block._side[ i ]._side->getGC();
820 gp_Vec gcDir( p1, gc );
821 gp_Vec2d gcDirProj( gcDir * pln.XDirection(), gcDir * pln.YDirection());
822 double gcAngle = gcDirProj.Angle( gp::DX2d() );
823 foundSide = gcAngle < 0 ? angleOfSide.rbegin()->second : angleOfSide.begin()->second;
826 _DUMP_(" selected "<< foundSide );
829 // Orient the found side correctly
831 // corners of found side corresponding to nodes n1 and n2
832 bool xMax1, yMax1, xMax2, yMax2;
833 if ( !getEdgeEnds( sharedSideEdge2, xMax1, yMax1, xMax2, yMax2 ))
834 return error(SMESH_Comment("Internal error at ")<<__FILE__<<":"<<__LINE__),
835 _OrientedBlockSide(0);
837 for ( int ori = 0; ori < _OrientedIndexer::MAX_ORI+1; ++ori )
839 _OrientedBlockSide orientedSide( foundSide, ori );
840 const SMDS_MeshNode* n12 = orientedSide.cornerNode( xMax1, yMax1);
841 const SMDS_MeshNode* n22 = orientedSide.cornerNode( xMax2, yMax2);
842 if ( n1 == n12 && n2 == n22 )
845 error(SMESH_Comment("Failed to orient a block side found by edge ")<<sharedSideEdge1
846 << " of side " << startBlockSide
847 << " of block " << _blocks.size());
851 //================================================================================
853 * \brief: Fill rows (which are actually columns,if !alongN1N2) of nodes starting
854 * from the given quadrangle until another block corner encounters.
855 * n1 and n2 are at bottom of quad, n1 is at block corner.
857 //================================================================================
859 bool _Skin::fillRowsUntilCorner(const SMDS_MeshElement* quad,
860 const SMDS_MeshNode* n1,
861 const SMDS_MeshNode* n2,
862 vector<const SMDS_MeshNode*>& row1,
863 vector<const SMDS_MeshNode*>& row2,
864 const bool alongN1N2 )
866 const SMDS_MeshNode* corner1 = n1;
868 // Store nodes of quad in the rows and find new n1 and n2 to get
869 // the next face so that new n2 is on block edge
870 int i1 = quad->GetNodeIndex( n1 );
871 int i2 = quad->GetNodeIndex( n2 );
872 row1.clear(); row2.clear();
873 row1.push_back( n1 );
876 row1.push_back( n2 );
877 row2.push_back( oppositeNode( quad, i2 ));
878 row2.push_back( n1 = oppositeNode( quad, i1 ));
882 row2.push_back( n2 );
883 row1.push_back( n2 = oppositeNode( quad, i2 ));
884 row2.push_back( n1 = oppositeNode( quad, i1 ));
887 if ( isCornerNode( row1[1] ))
890 // Find the rest nodes
891 TIDSortedElemSet emptySet, avoidSet;
892 while ( !isCornerNode( n2 ) )
894 avoidSet.clear(); avoidSet.insert( quad );
895 quad = SMESH_MeshAlgos::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2 );
896 if ( !isQuadrangle( quad ))
899 row1.push_back( n2 = oppositeNode( quad, i1 ));
900 row2.push_back( n1 = oppositeNode( quad, i2 ));
902 return n1 != corner1;
905 //================================================================================
907 * \brief Return a corner face by a corner node
909 //================================================================================
911 const SMDS_MeshElement* _BlockSide::getCornerFace(const SMDS_MeshNode* cornerNode) const
913 int x, y, isXMax, isYMax, found = 0;
914 for ( isXMax = 0; isXMax < 2; ++isXMax )
916 for ( isYMax = 0; isYMax < 2; ++isYMax )
918 x = isXMax ? _index._xSize-1 : 0;
919 y = isYMax ? _index._ySize-1 : 0;
920 found = ( getNode(x,y) == cornerNode );
925 if ( !found ) return 0;
926 int dx = isXMax ? -1 : +1;
927 int dy = isYMax ? -1 : +1;
928 const SMDS_MeshNode* n1 = getNode(x,y);
929 const SMDS_MeshNode* n2 = getNode(x+dx,y);
930 const SMDS_MeshNode* n3 = getNode(x,y+dy);
931 const SMDS_MeshNode* n4 = getNode(x+dx,y+dy);
932 return SMDS_Mesh::FindFace(n1, n2, n3, n4 );
935 //================================================================================
937 * \brief Checks own validity
939 //================================================================================
941 bool _Block::isValid() const
943 bool ok = ( nbSides() == 6 );
945 // check only corners depending on side selection
946 EBoxSides adjacent[4] = { B_BOTTOM, B_RIGHT, B_TOP, B_LEFT };
947 EQuadEdge edgeAdj [4] = { Q_TOP, Q_RIGHT, Q_TOP, Q_RIGHT };
948 EQuadEdge edgeBack[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
950 for ( int i=0; ok && i < NB_QUAD_SIDES; ++i )
952 SMESH_OrientedLink eBack = _side[ B_BACK ].edge( edgeBack[i] );
953 SMESH_OrientedLink eAdja = _side[ adjacent[i] ].edge( edgeAdj[i] );
954 ok = ( eBack == eAdja );
956 ok = ok && ( _side[ B_BOTTOM ]._index.size() == _side[ B_TOP ]._index.size() &&
957 _side[ B_RIGHT ]._index.size() == _side[ B_LEFT ]._index.size() &&
958 _side[ B_FRONT ]._index.size() == _side[ B_BACK ]._index.size() );
964 //=======================================================================
965 //function : StdMeshers_HexaFromSkin_3D
967 //=======================================================================
969 StdMeshers_HexaFromSkin_3D::StdMeshers_HexaFromSkin_3D(int hypId, SMESH_Gen* gen)
970 :SMESH_3D_Algo(hypId, gen)
972 _name = "HexaFromSkin_3D";
975 StdMeshers_HexaFromSkin_3D::~StdMeshers_HexaFromSkin_3D()
979 //================================================================================
981 * \brief Main method, which generates hexaheda
983 //================================================================================
985 bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper)
988 int nbBlocks = skin.findBlocks(aMesh);
990 return error( skin.error());
992 vector< vector< const SMDS_MeshNode* > > columns;
993 int x, xSize, y, ySize, z, zSize;
996 for ( int i = 0; i < nbBlocks; ++i )
998 const _Block& block = skin.getBlock( i );
1000 // ------------------------------------------
1001 // Fill columns of nodes with existing nodes
1002 // ------------------------------------------
1004 xSize = block.getSide(B_BOTTOM).getHoriSize();
1005 ySize = block.getSide(B_BOTTOM).getVertSize();
1006 zSize = block.getSide(B_FRONT ).getVertSize();
1007 int X = xSize - 1, Y = ySize - 1, Z = zSize - 1;
1008 colIndex = _Indexer( xSize, ySize );
1009 columns.resize( colIndex.size() );
1011 // fill node columns by front and back box sides
1012 for ( x = 0; x < xSize; ++x ) {
1013 vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( x, 0 )];
1014 vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( x, Y )];
1015 column0.resize( zSize );
1016 column1.resize( zSize );
1017 for ( z = 0; z < zSize; ++z ) {
1018 column0[ z ] = block.getSide(B_FRONT).node( x, z );
1019 column1[ z ] = block.getSide(B_BACK) .node( x, z );
1022 // fill node columns by left and right box sides
1023 for ( y = 1; y < ySize-1; ++y ) {
1024 vector< const SMDS_MeshNode* >& column0 = columns[ colIndex( 0, y )];
1025 vector< const SMDS_MeshNode* >& column1 = columns[ colIndex( X, y )];
1026 column0.resize( zSize );
1027 column1.resize( zSize );
1028 for ( z = 0; z < zSize; ++z ) {
1029 column0[ z ] = block.getSide(B_LEFT) .node( y, z );
1030 column1[ z ] = block.getSide(B_RIGHT).node( y, z );
1033 // get nodes from top and bottom box sides
1034 for ( x = 1; x < xSize-1; ++x ) {
1035 for ( y = 1; y < ySize-1; ++y ) {
1036 vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
1037 column.resize( zSize );
1038 column.front() = block.getSide(B_BOTTOM).node( x, y );
1039 column.back() = block.getSide(B_TOP) .node( x, y );
1043 // ----------------------------
1044 // Add internal nodes of a box
1045 // ----------------------------
1046 // projection points of internal nodes on box sub-shapes by which
1047 // coordinates of internal nodes are computed
1048 vector<gp_XYZ> pointOnShape( SMESH_Block::ID_Shell );
1050 // projections on vertices are constant
1051 pointOnShape[ SMESH_Block::ID_V000 ] = block.getSide(B_BOTTOM).xyz( 0, 0 );
1052 pointOnShape[ SMESH_Block::ID_V100 ] = block.getSide(B_BOTTOM).xyz( X, 0 );
1053 pointOnShape[ SMESH_Block::ID_V010 ] = block.getSide(B_BOTTOM).xyz( 0, Y );
1054 pointOnShape[ SMESH_Block::ID_V110 ] = block.getSide(B_BOTTOM).xyz( X, Y );
1055 pointOnShape[ SMESH_Block::ID_V001 ] = block.getSide(B_TOP).xyz( 0, 0 );
1056 pointOnShape[ SMESH_Block::ID_V101 ] = block.getSide(B_TOP).xyz( X, 0 );
1057 pointOnShape[ SMESH_Block::ID_V011 ] = block.getSide(B_TOP).xyz( 0, Y );
1058 pointOnShape[ SMESH_Block::ID_V111 ] = block.getSide(B_TOP).xyz( X, Y );
1060 for ( x = 1; x < xSize-1; ++x )
1062 gp_XYZ params; // normalized parameters of internal node within a unit box
1063 params.SetCoord( 1, x / double(X) );
1064 for ( y = 1; y < ySize-1; ++y )
1066 params.SetCoord( 2, y / double(Y) );
1067 // column to fill during z loop
1068 vector< const SMDS_MeshNode* >& column = columns[ colIndex( x, y )];
1069 // projections on horizontal edges
1070 pointOnShape[ SMESH_Block::ID_Ex00 ] = block.getSide(B_BOTTOM).xyz( x, 0 );
1071 pointOnShape[ SMESH_Block::ID_Ex10 ] = block.getSide(B_BOTTOM).xyz( x, Y );
1072 pointOnShape[ SMESH_Block::ID_E0y0 ] = block.getSide(B_BOTTOM).xyz( 0, y );
1073 pointOnShape[ SMESH_Block::ID_E1y0 ] = block.getSide(B_BOTTOM).xyz( X, y );
1074 pointOnShape[ SMESH_Block::ID_Ex01 ] = block.getSide(B_TOP).xyz( x, 0 );
1075 pointOnShape[ SMESH_Block::ID_Ex11 ] = block.getSide(B_TOP).xyz( x, Y );
1076 pointOnShape[ SMESH_Block::ID_E0y1 ] = block.getSide(B_TOP).xyz( 0, y );
1077 pointOnShape[ SMESH_Block::ID_E1y1 ] = block.getSide(B_TOP).xyz( X, y );
1078 // projections on horizontal sides
1079 pointOnShape[ SMESH_Block::ID_Fxy0 ] = block.getSide(B_BOTTOM).xyz( x, y );
1080 pointOnShape[ SMESH_Block::ID_Fxy1 ] = block.getSide(B_TOP) .xyz( x, y );
1081 for ( z = 1; z < zSize-1; ++z ) // z loop
1083 params.SetCoord( 3, z / double(Z) );
1084 // projections on vertical edges
1085 pointOnShape[ SMESH_Block::ID_E00z ] = block.getSide(B_FRONT).xyz( 0, z );
1086 pointOnShape[ SMESH_Block::ID_E10z ] = block.getSide(B_FRONT).xyz( X, z );
1087 pointOnShape[ SMESH_Block::ID_E01z ] = block.getSide(B_BACK).xyz( 0, z );
1088 pointOnShape[ SMESH_Block::ID_E11z ] = block.getSide(B_BACK).xyz( X, z );
1089 // projections on vertical sides
1090 pointOnShape[ SMESH_Block::ID_Fx0z ] = block.getSide(B_FRONT).xyz( x, z );
1091 pointOnShape[ SMESH_Block::ID_Fx1z ] = block.getSide(B_BACK) .xyz( x, z );
1092 pointOnShape[ SMESH_Block::ID_F0yz ] = block.getSide(B_LEFT) .xyz( y, z );
1093 pointOnShape[ SMESH_Block::ID_F1yz ] = block.getSide(B_RIGHT).xyz( y, z );
1095 // compute internal node coordinates
1097 SMESH_Block::ShellPoint( params, pointOnShape, coords );
1098 column[ z ] = aHelper->AddNode( coords.X(), coords.Y(), coords.Z() );
1102 //cout << "----------------------------------------------------------------------"<<endl;
1103 //for ( int id = SMESH_Block::ID_V000; id < SMESH_Block::ID_Shell; ++id)
1105 // gp_XYZ p = pointOnShape[ id ];
1106 // SMESH_Block::DumpShapeID( id,cout)<<" ( "<<p.X()<<", "<<p.Y()<<", "<<p.Z()<<" )"<<endl;
1108 //cout << "Params: ( "<< params.X()<<", "<<params.Y()<<", "<<params.Z()<<" )"<<endl;
1109 //cout << "coords: ( "<< coords.X()<<", "<<coords.Y()<<", "<<coords.Z()<<" )"<<endl;
1118 // find out orientation by a least distorted hexahedron (issue 0020855);
1119 // the last is defined by evaluating sum of face normals of 8 corner hexahedrons
1120 double badness = numeric_limits<double>::max();
1122 for ( int xMax = 0; xMax < 2; ++xMax )
1123 for ( int yMax = 0; yMax < 2; ++yMax )
1124 for ( int zMax = 0; zMax < 2; ++zMax )
1126 x = xMax ? xSize-1 : 1;
1127 y = yMax ? ySize-1 : 1;
1128 z = zMax ? zSize-1 : 1;
1129 vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x-1, y-1 )];
1130 vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x , y-1 )];
1131 vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x-1, y )];
1132 vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x , y )];
1134 const SMDS_MeshNode* n000 = col00[z-1];
1135 const SMDS_MeshNode* n100 = col10[z-1];
1136 const SMDS_MeshNode* n010 = col01[z-1];
1137 const SMDS_MeshNode* n110 = col11[z-1];
1138 const SMDS_MeshNode* n001 = col00[z];
1139 const SMDS_MeshNode* n101 = col10[z];
1140 const SMDS_MeshNode* n011 = col01[z];
1141 const SMDS_MeshNode* n111 = col11[z];
1142 SMDS_VolumeOfNodes probeVolume (n000,n010,n110,n100,
1143 n001,n011,n111,n101);
1144 SMDS_VolumeTool volTool( &probeVolume );
1145 double Nx=0.,Ny=0.,Nz=0.;
1146 for ( int iFace = 0; iFace < volTool.NbFaces(); ++iFace )
1149 volTool.GetFaceNormal( iFace, nx,ny,nz );
1154 double quality = Nx*Nx + Ny*Ny + Nz*Nz;
1155 if ( quality < badness )
1158 isForw = volTool.IsForward();
1163 for ( x = 0; x < xSize-1; ++x ) {
1164 for ( y = 0; y < ySize-1; ++y ) {
1165 vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x, y )];
1166 vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x+1, y )];
1167 vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x, y+1 )];
1168 vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x+1, y+1 )];
1169 // bottom face normal of a hexa mush point outside the volume
1171 for ( z = 0; z < zSize-1; ++z )
1172 aHelper->AddVolume(col00[z], col01[z], col11[z], col10[z],
1173 col00[z+1], col01[z+1], col11[z+1], col10[z+1]);
1175 for ( z = 0; z < zSize-1; ++z )
1176 aHelper->AddVolume(col00[z], col10[z], col11[z], col01[z],
1177 col00[z+1], col10[z+1], col11[z+1], col01[z+1]);
1185 //================================================================================
1187 * \brief Evaluate nb of hexa
1189 //================================================================================
1191 bool StdMeshers_HexaFromSkin_3D::Evaluate(SMESH_Mesh & aMesh,
1192 const TopoDS_Shape & aShape,
1193 MapShapeNbElems& aResMap)
1196 int nbBlocks = skin.findBlocks(aMesh);
1197 if ( nbBlocks == 0 )
1198 return error( skin.error());
1200 bool secondOrder = aMesh.NbFaces( ORDER_QUADRATIC );
1202 int entity = secondOrder ? SMDSEntity_Quad_Hexa : SMDSEntity_Hexa;
1203 vector<int>& nbByType = aResMap[ aMesh.GetSubMesh( aShape )];
1204 if ( entity >= (int) nbByType.size() )
1205 nbByType.resize( SMDSEntity_Last, 0 );
1207 for ( int i = 0; i < nbBlocks; ++i )
1209 const _Block& block = skin.getBlock( i );
1211 int nbX = block.getSide(B_BOTTOM).getHoriSize();
1212 int nbY = block.getSide(B_BOTTOM).getVertSize();
1213 int nbZ = block.getSide(B_FRONT ).getVertSize();
1215 int nbHexa = (nbX-1) * (nbY-1) * (nbZ-1);
1216 int nbNodes = (nbX-2) * (nbY-2) * (nbZ-2);
1219 (nbX-2) * (nbY-2) * (nbZ-1) +
1220 (nbX-2) * (nbY-1) * (nbZ-2) +
1221 (nbX-1) * (nbY-2) * (nbZ-2);
1224 nbByType[ entity ] += nbHexa;
1225 nbByType[ SMDSEntity_Node ] += nbNodes;
1231 //================================================================================
1233 * \brief Abstract method must be defined but does nothing
1235 //================================================================================
1237 bool StdMeshers_HexaFromSkin_3D::CheckHypothesis(SMESH_Mesh&, const TopoDS_Shape&,
1238 Hypothesis_Status& aStatus)
1240 aStatus = SMESH_Hypothesis::HYP_OK;
1244 //================================================================================
1246 * \brief Abstract method must be defined but just reports an error as this
1247 * algo is not intended to work with shapes
1249 //================================================================================
1251 bool StdMeshers_HexaFromSkin_3D::Compute(SMESH_Mesh&, const TopoDS_Shape&)
1253 return error("Algorithm can't work with geometrical shapes");