-// Copyright (C) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2011 CEA/DEN, EDF R&D, OPEN CASCADE
//
-// This library is free software; you can redistribute it and/or
-// modify it under the terms of the GNU Lesser General Public
-// License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
//
-// This library is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-// Lesser General Public License for more details.
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
//
-// You should have received a copy of the GNU Lesser General Public
-// License along with this library; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
-//
-// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
// File : StdMeshers_HexaFromSkin_3D.cxx
// Created : Wed Jan 27 12:28:07 2010
// Author : Edward AGAPOV (eap)
-//
+
#include "StdMeshers_HexaFromSkin_3D.hxx"
#include "SMDS_VolumeOfNodes.hxx"
#include <gp_Ax2.hxx>
//#include "utilities.h"
+#include <limits>
-// Define error message
+// Define error message and _MYDEBUG_ if needed
#ifdef _DEBUG_
#define BAD_MESH_ERR \
error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh.\n" \
__FILE__ ":" )<<__LINE__)
+//#define _MYDEBUG_
#else
#define BAD_MESH_ERR \
error(SMESH_Comment("Can't detect block-wise structure of the input 2D mesh"))
#endif
-// Debug output
-#define _DUMP_(msg) //cout << msg << endl
+// Debug output
+#ifdef _MYDEBUG_
+#define _DUMP_(msg) cout << msg << endl
+#else
+#define _DUMP_(msg)
+#endif
namespace
{
enum EBoxSides //!< sides of the block
{
- B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, B_UNDEFINED
+ B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, NB_BLOCK_SIDES
};
- const char* SBoxSides[] = //!< names of block sides
+#ifdef _MYDEBUG_
+ const char* SBoxSides[] = //!< names of block sides -- needed for DEBUG only
{
"BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED"
};
+#endif
enum EQuadEdge //!< edges of quadrangle side
{
- Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT, Q_UNDEFINED
+ Q_BOTTOM = 0, Q_RIGHT, Q_TOP, Q_LEFT, NB_QUAD_SIDES
};
//================================================================================
/*!
* \brief return true if a node is at block corner
+ *
+ * This check is valid for simple cases only
*/
//================================================================================
bool isCornerNode( const SMDS_MeshNode* n )
{
- return n && n->NbInverseElements( SMDSAbs_Face ) % 2;
+ int nbF = n ? n->NbInverseElements( SMDSAbs_Face ) : 1;
+ if ( nbF % 2 )
+ return true;
+
+ set<const SMDS_MeshNode*> nodesInInverseFaces;
+ SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator(SMDSAbs_Face );
+ while ( fIt->more() )
+ {
+ const SMDS_MeshElement* face = fIt->next();
+ nodesInInverseFaces.insert( face->begin_nodes(), face->end_nodes() );
+ }
+
+ return nodesInInverseFaces.size() != ( 6 + (nbF/2-1)*3 );
}
//================================================================================
bool isQuadrangle(const SMDS_MeshElement* e)
{
- return ( e && e->NbNodes() == ( e->IsQuadratic() ? 8 : 4 ));
+ return ( e && e->NbCornerNodes() == 4 );
}
//================================================================================
int _nbBlocksExpected;
int _nbBlocksFound;
-#ifdef _DEBUG_
-#define _grid_access_(args) _grid.at( args )
+#ifdef _DEBUG_ // want to get SIGSEGV in case of invalid index
+#define _grid_access_(pobj, i) pobj->_grid[ ((i) < pobj->_grid.size()) ? i : int(1e100)]
#else
-#define _grid_access_(args) _grid[ args ]
+#define _grid_access_(pobj, i) pobj->_grid[ i ]
#endif
//!< Return node at XY
- const SMDS_MeshNode* getNode(int x, int y) const { return _grid_access_( _index( x, y )); }
+ const SMDS_MeshNode* getNode(int x, int y) const { return _grid_access_(this, _index( x,y ));}
//!< Set node at XY
- void setNode(int x, int y, const SMDS_MeshNode* n) { _grid_access_( _index( x, y )) = n; }
+ void setNode(int x, int y, const SMDS_MeshNode* n) { _grid_access_(this, _index( x,y )) = n; }
//!< Return an edge
SMESH_OrientedLink getEdge(EQuadEdge edge) const
{
//!< True if all blocks this side belongs to have beem found
bool isBound() const { return _nbBlocksExpected <= _nbBlocksFound; }
//!< Return coordinates of node at XY
- gp_XYZ getXYZ(int x, int y) const { return SMESH_MeshEditor::TNodeXYZ( getNode( x, y )); }
+ gp_XYZ getXYZ(int x, int y) const { return SMESH_TNodeXYZ( getNode( x, y )); }
//!< Return gravity center of the four corners and the middle node
gp_XYZ getGC() const
{
//!< return coordinates by XY
gp_XYZ xyz(int x, int y) const
{
- return SMESH_MeshEditor::TNodeXYZ( _side->_grid_access_( _index( x, y )) );
+ return SMESH_TNodeXYZ( _grid_access_(_side, _index( x, y )) );
}
//!< safely return a node by XY
const SMDS_MeshNode* node(int x, int y) const
//!< Return a corner node
const SMDS_MeshNode* cornerNode(bool isXMax, bool isYMax) const
{
- return _side->_grid_access_( _index.corner( isXMax, isYMax ));
+ return _grid_access_(_side, _index.corner( isXMax, isYMax ));
}
//!< return its size in nodes
int getHoriSize() const { return _index.xSize(); }
*/
struct _Block
{
- _OrientedBlockSide _side[6]; // 6 sides of a sub-block
+ _OrientedBlockSide _side[6]; // 6 sides of a sub-block
+ set<const SMDS_MeshNode*> _corners;
const _OrientedBlockSide& getSide(int i) const { return _side[i]; }
+ bool setSide( int i, const _OrientedBlockSide& s)
+ {
+ if (( _side[i] = s ))
+ {
+ _corners.insert( s.cornerNode(0,0));
+ _corners.insert( s.cornerNode(1,0));
+ _corners.insert( s.cornerNode(0,1));
+ _corners.insert( s.cornerNode(1,1));
+ }
+ return s;
+ }
+ void clear() { for (int i=0;i<6;++i) _side[i]=0; _corners.clear(); }
bool hasSide( const _OrientedBlockSide& s) const
{
if ( s ) for (int i=0;i<6;++i) if ( _side[i] && _side[i]._side == s._side ) return true;
return false;
}
+ int nbSides() const { int n=0; for (int i=0;i<6;++i) if ( _side[i] ) ++n; return n; }
+ bool isValid() const;
};
//================================================================================
/*!
const SMESH_Comment& error() const { return _error; }
private:
- bool fillSide( _BlockSide& side, const SMDS_MeshElement* cornerQuad);
+ bool fillSide( _BlockSide& side,
+ const SMDS_MeshElement* cornerQuad,
+ const SMDS_MeshNode* cornerNode);
bool fillRowsUntilCorner(const SMDS_MeshElement* quad,
const SMDS_MeshNode* n1,
const SMDS_MeshNode* n2,
bool alongN1N2 );
_OrientedBlockSide findBlockSide( EBoxSides startBlockSide,
EQuadEdge sharedSideEdge1,
- EQuadEdge sharedSideEdge2);
+ EQuadEdge sharedSideEdge2,
+ bool withGeometricAnalysis);
//!< update own data and data of the side bound to block
void setSideBoundToBlock( _BlockSide& side )
{
- side._nbBlocksFound++;
- if ( side.isBound() )
- {
- for ( int e = 0; e < int(Q_UNDEFINED); ++e )
+ if ( side._nbBlocksFound++, side.isBound() )
+ for ( int e = 0; e < int(NB_QUAD_SIDES); ++e )
_edge2sides[ side.getEdge( (EQuadEdge) e ) ].erase( &side );
- }
}
//!< store reason of error
int error(const SMESH_Comment& reason) { _error = reason; return 0; }
_allSides.push_back( _BlockSide() );
_BlockSide& side = _allSides.back();
- if ( !fillSide( side, face ) )
+ if ( !fillSide( side, face, *corner ) )
{
if ( !_error.empty() )
return false;
for ( int isYMax = 0; isYMax < 2; ++isYMax )
{
const SMDS_MeshNode* nCorner = side.getCornerNode(isXMax,isYMax );
- if ( !isCornerNode( nCorner ))
- return BAD_MESH_ERR;
- //_corner2sides[ nCorner ].insert( &side );
corners.push_back( nCorner );
cornerFaces.insert( side.getCornerFace( nCorner ));
}
- for ( int e = 0; e < int(Q_UNDEFINED); ++e )
+ for ( int e = 0; e < int(NB_QUAD_SIDES); ++e )
_edge2sides[ side.getEdge( (EQuadEdge) e ) ].insert( &side );
nbFacesOnSides += side.getNbFaces();
return BAD_MESH_ERR;
if ( _allSides.back()._grid.empty() )
_allSides.pop_back();
+ _DUMP_("Nb detected sides "<< _allSides.size());
// ---------------------------
// Organize sides into blocks
// ---------------------------
- // analyse sharing of sides by blocks
- int nbBlockSides = 0; // nb of block sides taking into account their sharing
- list < _BlockSide >::iterator sideIt = _allSides.begin();
- for ( ; sideIt != _allSides.end(); ++sideIt )
+ // analyse sharing of sides by blocks and sort sides by nb of adjacent sides
+ int nbBlockSides = 0; // total nb of block sides taking into account their sharing
+ multimap<int, _BlockSide* > sortedSides;
{
- _BlockSide& side = *sideIt;
- bool isSharedSide = true;
- for ( int e = 0; e < int(Q_UNDEFINED) && isSharedSide; ++e )
- isSharedSide = _edge2sides[ side.getEdge( (EQuadEdge) e ) ].size() > 2;
- side._nbBlocksFound = 0;
- side._nbBlocksExpected = isSharedSide ? 2 : 1;
- nbBlockSides += side._nbBlocksExpected;
+ list < _BlockSide >::iterator sideIt = _allSides.begin();
+ for ( ; sideIt != _allSides.end(); ++sideIt )
+ {
+ _BlockSide& side = *sideIt;
+ bool isSharedSide = true;
+ int nbAdjacent = 0;
+ for ( int e = 0; e < int(NB_QUAD_SIDES) && isSharedSide; ++e )
+ {
+ int nbAdj = _edge2sides[ side.getEdge( (EQuadEdge) e ) ].size();
+ nbAdjacent += nbAdj;
+ isSharedSide = ( nbAdj > 2 );
+ }
+ side._nbBlocksFound = 0;
+ side._nbBlocksExpected = isSharedSide ? 2 : 1;
+ nbBlockSides += side._nbBlocksExpected;
+ sortedSides.insert( make_pair( nbAdjacent, & side ));
+ }
}
// find sides of each block
while ( nbBlockSides >= 6 )
{
// get any side not bound to all blocks it belongs to
- sideIt = _allSides.begin();
- while ( sideIt != _allSides.end() && sideIt->isBound())
- ++sideIt;
+ multimap<int, _BlockSide*>::iterator i_side = sortedSides.begin();
+ while ( i_side != sortedSides.end() && i_side->second->isBound())
+ ++i_side;
// start searching for block sides from the got side
bool ok = true;
_blocks.resize( _blocks.size() + 1 );
_Block& block = _blocks.back();
- block._side[B_FRONT] = &(*sideIt);
- setSideBoundToBlock( *sideIt );
+ block.setSide( B_FRONT, i_side->second );
+ setSideBoundToBlock( *i_side->second );
nbBlockSides--;
-
- // edges of neighbour sides of B_FRONT corresponding to front's edges
- EQuadEdge edgeOfFront[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
- EQuadEdge edgeToFind [4] = { Q_BOTTOM, Q_LEFT, Q_BOTTOM, Q_LEFT };
- for ( int i = Q_BOTTOM; ok && i <= Q_LEFT; ++i )
- ok = ( block._side[i] = findBlockSide( B_FRONT, edgeOfFront[i], edgeToFind[i]));
- if ( ok )
- ok = ( block._side[B_BACK] = findBlockSide( B_TOP, Q_TOP, Q_TOP ));
+ // edges of adjacent sides of B_FRONT corresponding to front's edges
+ EQuadEdge edgeOfFront[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
+ EQuadEdge edgeOfAdj [4] = { Q_BOTTOM, Q_LEFT, Q_BOTTOM, Q_LEFT };
+ // first find all sides detectable w/o advanced analysis,
+ // then repeat the search, which then may pass without advanced analysis
+ for ( int advAnalys = 0; advAnalys < 2; ++advAnalys )
+ {
+ for ( int i = 0; (ok || !advAnalys) && i < NB_QUAD_SIDES; ++i )
+ if ( !block._side[i] ) // try to find 4 sides adjacent to front side
+ ok = block.setSide( i, findBlockSide( B_FRONT, edgeOfFront[i],edgeOfAdj[i],advAnalys));
+ if ( ok || !advAnalys)
+ if ( !block._side[B_BACK] && block._side[B_TOP] ) // try to find back side by top one
+ ok = block.setSide( B_BACK, findBlockSide( B_TOP, Q_TOP, Q_TOP, advAnalys ));
+ if ( !advAnalys ) ok = true;
+ }
+ ok = block.isValid();
if ( ok )
{
// check if just found block is same as one of previously found
bool isSame = false;
for ( int i = 1; i < _blocks.size() && !isSame; ++i )
- {
- _Block& prevBlock = _blocks[i-1];
- isSame = true;
- for ( int j = 0; j < 6 && isSame; ++j )
- isSame = prevBlock.hasSide( block._side[ j ]);
- }
+ isSame = ( block._corners == _blocks[i-1]._corners );
ok = !isSame;
}
// count the found sides
- _DUMP_(endl);
- for (int i = 0; i < B_UNDEFINED; ++i )
+ _DUMP_(endl << "** Block " << _blocks.size() << " valid: " << block.isValid());
+ for (int i = 0; i < NB_BLOCK_SIDES; ++i )
{
- _DUMP_("** Block side "<< SBoxSides[i] <<" "<< block._side[ i ]._side);
+ _DUMP_("\tSide "<< SBoxSides[i] <<" "<< block._side[ i ]._side);
if ( block._side[ i ] )
{
if ( ok && i != B_FRONT)
setSideBoundToBlock( *block._side[ i ]._side );
nbBlockSides--;
}
- _DUMP_("Corner 0,0 "<< block._side[ i ].cornerNode(0,0));
- _DUMP_("Corner 1,0 "<< block._side[ i ].cornerNode(1,0));
- _DUMP_("Corner 1,1 "<< block._side[ i ].cornerNode(1,1));
- _DUMP_("Corner 0,1 "<< block._side[ i ].cornerNode(0,1));
+ _DUMP_("\t corners "<<
+ block._side[ i ].cornerNode(0,0)->GetID() << ", " <<
+ block._side[ i ].cornerNode(1,0)->GetID() << ", " <<
+ block._side[ i ].cornerNode(1,1)->GetID() << ", " <<
+ block._side[ i ].cornerNode(0,1)->GetID() << ", "<<endl);
}
else
{
- _DUMP_("Not found"<<endl);
+ _DUMP_("\t not found"<<endl);
}
}
if ( !ok )
- block._side[0] = block._side[1] = block._side[2] =
- block._side[3] = block._side[4] = block._side[5] = 0;
+ block.clear();
else
nbBlocks++;
}
+ _DUMP_("Nb found blocks "<< nbBlocks <<endl);
+
if ( nbBlocks == 0 && _error.empty() )
return BAD_MESH_ERR;
*/
//================================================================================
- bool _Skin::fillSide( _BlockSide& side, const SMDS_MeshElement* cornerQuad)
+ bool _Skin::fillSide( _BlockSide& side,
+ const SMDS_MeshElement* cornerQuad,
+ const SMDS_MeshNode* nCorner)
{
// Find out size of block side mesured in nodes and by the way find two rows
// of nodes in two directions.
int x, y, nbX, nbY;
const SMDS_MeshElement* firstQuad = cornerQuad;
{
- // get corner node of cornerQuad
- const SMDS_MeshNode* nCorner = 0;
- for ( int i = 0; i < 4 && !nCorner; ++i )
- if ( isCornerNode( firstQuad->GetNode(i)))
- nCorner = firstQuad->GetNode(i);
- if ( !nCorner ) return false;
-
// get a node on block edge
int iCorner = firstQuad->GetNodeIndex( nCorner );
const SMDS_MeshNode* nOnEdge = firstQuad->GetNode( (iCorner+1) % 4);
}
}
- return true;
+ // check side validity
+ bool ok =
+ side.getCornerFace( side.getCornerNode( 0, 0 )) &&
+ side.getCornerFace( side.getCornerNode( 1, 0 )) &&
+ side.getCornerFace( side.getCornerNode( 0, 1 )) &&
+ side.getCornerFace( side.getCornerNode( 1, 1 ));
+
+ return ok;
}
//================================================================================
_OrientedBlockSide _Skin::findBlockSide( EBoxSides startBlockSide,
EQuadEdge sharedSideEdge1,
- EQuadEdge sharedSideEdge2)
+ EQuadEdge sharedSideEdge2,
+ bool withGeometricAnalysis)
{
_Block& block = _blocks.back();
_OrientedBlockSide& side1 = block._side[ startBlockSide ];
set< _BlockSide* > sidesOnEdge = _edge2sides[ edge ]; // copy a set
// exclude loaded sides of block from sidesOnEdge
- int nbLoadedSides = 0;
- for (int i = 0; i < B_UNDEFINED; ++i )
- {
+ for (int i = 0; i < NB_BLOCK_SIDES; ++i )
if ( block._side[ i ] )
- {
- nbLoadedSides++;
sidesOnEdge.erase( block._side[ i ]._side );
- }
- }
+
int nbSidesOnEdge = sidesOnEdge.size();
_DUMP_("nbSidesOnEdge "<< nbSidesOnEdge << " " << n1->GetID() << "-" << n2->GetID() );
if ( nbSidesOnEdge == 0 )
return 0;
_BlockSide* foundSide = 0;
- if ( nbSidesOnEdge == 1 /*|| nbSidesOnEdge == 2 && nbLoadedSides == 1 */)
+ if ( nbSidesOnEdge == 1 )
{
foundSide = *sidesOnEdge.begin();
}
else
{
- // Select one of found sides most close to startBlockSide
-
- // gravity center of already loaded block sides
- gp_XYZ gc(0,0,0);
- for (int i = 0; i < B_UNDEFINED; ++i )
- if ( block._side[ i ] )
- gc += block._side[ i ]._side->getGC();
- gc /= nbLoadedSides;
-
- gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()), p2 (n2->X(),n2->Y(),n2->Z());
- gp_Vec p1p2( p1, p2 );
-
- const SMDS_MeshElement* face1 = side1->getCornerFace( n1 );
- gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(n1)));
- gp_Vec side1Dir( p1, p1Op );
- gp_Ax2 pln( p1, p1p2, side1Dir ); // plane with normal p1p2 and X dir side1Dir
- _DUMP_(" Select adjacent for "<< side1._side << " - side dir ("
- << side1Dir.X() << ", " << side1Dir.Y() << ", " << side1Dir.Z() << ")" );
-
- map < double , _BlockSide* > angleOfSide;
set< _BlockSide* >::iterator sideIt = sidesOnEdge.begin();
- for (; sideIt != sidesOnEdge.end(); ++sideIt )
+ int nbLoadedSides = block.nbSides();
+ if ( nbLoadedSides > 1 )
+ {
+ // Find the side having more than 2 corners common with already loaded sides
+ for (; !foundSide && sideIt != sidesOnEdge.end(); ++sideIt )
+ {
+ _BlockSide* sideI = *sideIt;
+ int nbCommonCorners =
+ block._corners.count( sideI->getCornerNode(0,0)) +
+ block._corners.count( sideI->getCornerNode(1,0)) +
+ block._corners.count( sideI->getCornerNode(0,1)) +
+ block._corners.count( sideI->getCornerNode(1,1));
+ if ( nbCommonCorners > 2 )
+ foundSide = sideI;
+ }
+ }
+
+ if ( !foundSide )
{
- _BlockSide* sideI = *sideIt;
- const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 );
- gp_XYZ p1Op = SMESH_MeshEditor::TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(n1)));
- gp_Vec sideIDir( p1, p1Op );
- // compute angle of (sideIDir projection to pln) and (X dir of pln)
- gp_Vec2d sideIDirProj( sideIDir * pln.XDirection(), sideIDir * pln.YDirection());
- double angle = sideIDirProj.Angle( gp::DX2d() );
- if ( angle < 0 ) angle += 2 * PI;
- angleOfSide.insert( make_pair( angle, sideI ));
- _DUMP_(" "<< sideI << " - side dir ("
- << sideIDir.X() << ", " << sideIDir.Y() << ", " << sideIDir.Z() << ")"
- << " angle " << angle);
+ if ( !withGeometricAnalysis ) return 0;
+
+ // Select one of found sides most close to startBlockSide
+
+ gp_XYZ p1 ( n1->X(),n1->Y(),n1->Z()), p2 (n2->X(),n2->Y(),n2->Z());
+ gp_Vec p1p2( p1, p2 );
+
+ const SMDS_MeshElement* face1 = side1->getCornerFace( n1 );
+ gp_XYZ p1Op = SMESH_TNodeXYZ( oppositeNode( face1, face1->GetNodeIndex(n1)));
+ gp_Vec side1Dir( p1, p1Op );
+ gp_Ax2 pln( p1, p1p2, side1Dir ); // plane with normal p1p2 and X dir side1Dir
+ _DUMP_(" Select adjacent for "<< side1._side << " - side dir ("
+ << side1Dir.X() << ", " << side1Dir.Y() << ", " << side1Dir.Z() << ")" );
+
+ map < double , _BlockSide* > angleOfSide;
+ for (sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
+ {
+ _BlockSide* sideI = *sideIt;
+ const SMDS_MeshElement* faceI = sideI->getCornerFace( n1 );
+ gp_XYZ p1Op = SMESH_TNodeXYZ( oppositeNode( faceI, faceI->GetNodeIndex(n1)));
+ gp_Vec sideIDir( p1, p1Op );
+ // compute angle of (sideIDir projection to pln) and (X dir of pln)
+ gp_Vec2d sideIDirProj( sideIDir * pln.XDirection(), sideIDir * pln.YDirection());
+ double angle = sideIDirProj.Angle( gp::DX2d() );
+ if ( angle < 0 ) angle += 2. * M_PI; // angle [0-2*PI]
+ angleOfSide.insert( make_pair( angle, sideI ));
+ _DUMP_(" "<< sideI << " - side dir ("
+ << sideIDir.X() << ", " << sideIDir.Y() << ", " << sideIDir.Z() << ")"
+ << " angle " << angle);
+ }
+ if ( nbLoadedSides == 1 )
+ {
+ double angF = angleOfSide.begin()->first, angL = angleOfSide.rbegin()->first;
+ if ( angF > M_PI ) angF = 2.*M_PI - angF;
+ if ( angL > M_PI ) angL = 2.*M_PI - angL;
+ foundSide = angF < angL ? angleOfSide.begin()->second : angleOfSide.rbegin()->second;
+ }
+ else
+ {
+ gp_XYZ gc(0,0,0); // gravity center of already loaded block sides
+ for (int i = 0; i < NB_BLOCK_SIDES; ++i )
+ if ( block._side[ i ] )
+ gc += block._side[ i ]._side->getGC();
+ gc /= nbLoadedSides;
+
+ gp_Vec gcDir( p1, gc );
+ gp_Vec2d gcDirProj( gcDir * pln.XDirection(), gcDir * pln.YDirection());
+ double gcAngle = gcDirProj.Angle( gp::DX2d() );
+ foundSide = gcAngle < 0 ? angleOfSide.rbegin()->second : angleOfSide.begin()->second;
+ }
}
- gp_Vec gcDir( p1, gc );
- gp_Vec2d gcDirProj( gcDir * pln.XDirection(), gcDir * pln.YDirection());
- double gcAngle = gcDirProj.Angle( gp::DX2d() );
- foundSide = gcAngle < 0 ? angleOfSide.rbegin()->second : angleOfSide.begin()->second;
_DUMP_(" selected "<< foundSide );
}
row2.push_back( n1 = oppositeNode( quad, i1 ));
}
+ if ( isCornerNode( row1[1] ))
+ return true;
+
// Find the rest nodes
TIDSortedElemSet emptySet, avoidSet;
- while ( !isCornerNode( n2 ))
+ while ( !isCornerNode( n2 ) )
{
avoidSet.clear(); avoidSet.insert( quad );
quad = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, avoidSet, &i1, &i2 );
return SMDS_Mesh::FindFace(n1, n2, n3, n4 );
}
-}
+ //================================================================================
+ /*!
+ * \brief Checks own validity
+ */
+ //================================================================================
+
+ bool _Block::isValid() const
+ {
+ bool ok = ( nbSides() == 6 );
+
+ // check only corners depending on side selection
+ EBoxSides adjacent[4] = { B_BOTTOM, B_RIGHT, B_TOP, B_LEFT };
+ EQuadEdge edgeAdj [4] = { Q_TOP, Q_RIGHT, Q_TOP, Q_RIGHT };
+ EQuadEdge edgeBack[4] = { Q_BOTTOM, Q_RIGHT, Q_TOP, Q_LEFT };
+
+ for ( int i=0; ok && i < NB_QUAD_SIDES; ++i )
+ {
+ SMESH_OrientedLink eBack = _side[ B_BACK ].edge( edgeBack[i] );
+ SMESH_OrientedLink eAdja = _side[ adjacent[i] ].edge( edgeAdj[i] );
+ ok = ( eBack == eAdja );
+ }
+ return ok;
+ }
+
+} // namespace
//=======================================================================
//function : StdMeshers_HexaFromSkin_3D
// ----------------------------
// Add internal nodes of a box
// ----------------------------
- // projection points of internal nodes on box subshapes by which
+ // projection points of internal nodes on box sub-shapes by which
// coordinates of internal nodes are computed
vector<gp_XYZ> pointOnShape( SMESH_Block::ID_Shell );
// Add hexahedrons
// ----------------
- // find out orientation
- const SMDS_MeshNode* n000 = block.getSide(B_BOTTOM).cornerNode( 0, 0 );
- const SMDS_MeshNode* n100 = block.getSide(B_BOTTOM).cornerNode( 1, 0 );
- const SMDS_MeshNode* n010 = block.getSide(B_BOTTOM).cornerNode( 0, 1 );
- const SMDS_MeshNode* n110 = block.getSide(B_BOTTOM).cornerNode( 1, 1 );
- const SMDS_MeshNode* n001 = block.getSide(B_TOP).cornerNode( 0, 0 );
- const SMDS_MeshNode* n101 = block.getSide(B_TOP).cornerNode( 1, 0 );
- const SMDS_MeshNode* n011 = block.getSide(B_TOP).cornerNode( 0, 1 );
- const SMDS_MeshNode* n111 = block.getSide(B_TOP).cornerNode( 1, 1 );
- SMDS_VolumeOfNodes probeVolume (n000,n010,n110,n100,
- n001,n011,n111,n101);
- bool isForw = SMDS_VolumeTool( &probeVolume ).IsForward();
+ // find out orientation by a least distorted hexahedron (issue 0020855);
+ // the last is defined by evaluating sum of face normals of 8 corner hexahedrons
+ double badness = numeric_limits<double>::max();
+ bool isForw = true;
+ for ( int xMax = 0; xMax < 2; ++xMax )
+ for ( int yMax = 0; yMax < 2; ++yMax )
+ for ( int zMax = 0; zMax < 2; ++zMax )
+ {
+ x = xMax ? xSize-1 : 1;
+ y = yMax ? ySize-1 : 1;
+ z = zMax ? zSize-1 : 1;
+ vector< const SMDS_MeshNode* >& col00 = columns[ colIndex( x-1, y-1 )];
+ vector< const SMDS_MeshNode* >& col10 = columns[ colIndex( x , y-1 )];
+ vector< const SMDS_MeshNode* >& col01 = columns[ colIndex( x-1, y )];
+ vector< const SMDS_MeshNode* >& col11 = columns[ colIndex( x , y )];
+
+ const SMDS_MeshNode* n000 = col00[z-1];
+ const SMDS_MeshNode* n100 = col10[z-1];
+ const SMDS_MeshNode* n010 = col01[z-1];
+ const SMDS_MeshNode* n110 = col11[z-1];
+ const SMDS_MeshNode* n001 = col00[z];
+ const SMDS_MeshNode* n101 = col10[z];
+ const SMDS_MeshNode* n011 = col01[z];
+ const SMDS_MeshNode* n111 = col11[z];
+ SMDS_VolumeOfNodes probeVolume (n000,n010,n110,n100,
+ n001,n011,n111,n101);
+ SMDS_VolumeTool volTool( &probeVolume );
+ double Nx=0.,Ny=0.,Nz=0.;
+ for ( int iFace = 0; iFace < volTool.NbFaces(); ++iFace )
+ {
+ double nx,ny,nz;
+ volTool.GetFaceNormal( iFace, nx,ny,nz );
+ Nx += nx;
+ Ny += ny;
+ Nz += nz;
+ }
+ double quality = Nx*Nx + Ny*Ny + Nz*Nz;
+ if ( quality < badness )
+ {
+ badness = quality;
+ isForw = volTool.IsForward();
+ }
+ }
// add elements
for ( x = 0; x < xSize-1; ++x ) {
{
return error("Algorithm can't work with geometrical shapes");
}
-