-// Copyright (C) 2009-2012 CEA/DEN, EDF R&D
+// Copyright (C) 2009-2011 CEA/DEN, EDF R&D
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
#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
+#ifdef _MYDEBUG_
+#define _DUMP_(msg) cout << msg << endl
+#else
+#define _DUMP_(msg)
+#endif
#ifdef _DEBUG_
{
B_BOTTOM=0, B_RIGHT, B_TOP, B_LEFT, B_FRONT, B_BACK, NB_BLOCK_SIDES
};
+#ifdef _MYDEBUG_
const char* SBoxSides[] = //!< names of block sides
{
"BOTTOM", "RIGHT", "TOP", "LEFT", "FRONT", "BACK", "UNDEFINED"
};
+#endif
enum EQuadEdge //!< edges of quadrangle side
{
Q_BOTTOM = 0, Q_RIGHT, Q_TOP, Q_LEFT, NB_QUAD_SIDES
typedef void (*TFun)(int& x, int& y);
TFun _xRevFun, _yRevFun, _swapFun;
- static void lazy(int&, int&) {}
+ static void lazy (int&, int&) {}
static void reverse(int& x, int& size) { x = size - x - 1; }
- static void swap(int& x, int& y) { std::swap(x,y); }
+ static void swap (int& x, int& y) { std::swap(x,y); }
};
//================================================================================
/*!
return getNode( isXMax ? _index._xSize-1 : 0 , isYMax ? _index._ySize-1 : 0 );
}
const SMDS_MeshElement* getCornerFace(const SMDS_MeshNode* cornerNode) const;
- //!< True if all blocks this side belongs to have beem found
+ //!< True if all blocks this side belongs to have been found
bool isBound() const { return _nbBlocksExpected <= _nbBlocksFound; }
//!< Return coordinates of node at XY
gp_XYZ getXYZ(int x, int y) const { return SMESH_TNodeXYZ( getNode( x, y )); }
vector<const SMDS_MeshNode*>& verRow1,
vector<const SMDS_MeshNode*>& verRow2,
bool alongN1N2 );
- _OrientedBlockSide findBlockSide( EBoxSides startBlockSide,
- EQuadEdge sharedSideEdge1,
- EQuadEdge sharedSideEdge2,
- bool withGeometricAnalysis);
+ _OrientedBlockSide findBlockSide( EBoxSides startBlockSide,
+ EQuadEdge sharedSideEdge1,
+ EQuadEdge sharedSideEdge2,
+ bool withGeometricAnalysis,
+ set< _BlockSide* >& sidesAround);
//!< update own data and data of the side bound to block
void setSideBoundToBlock( _BlockSide& side )
{
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
+ set< _BlockSide* > sidesAround;
for ( int advAnalys = 0; advAnalys < 2; ++advAnalys )
{
+ // try to find 4 sides adjacent to a FRONT side
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 ( !block._side[i] )
+ ok = block.setSide( i, findBlockSide( B_FRONT, edgeOfFront[i], edgeOfAdj[i],
+ advAnalys, sidesAround));
+ // try to find a BACK side by a TOP one
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 ( !block._side[B_BACK] && block._side[B_TOP] )
+ ok = block.setSide( B_BACK, findBlockSide( B_TOP, Q_TOP, Q_TOP,
+ advAnalys, sidesAround ));
if ( !advAnalys ) ok = true;
}
ok = block.isValid();
if ( ok )
{
- // check if just found block is same as one of previously found
+ // check if just found block is same as one of previously found blocks
bool isSame = false;
for ( int i = 1; i < _blocks.size() && !isSame; ++i )
isSame = ( block._corners == _blocks[i-1]._corners );
return ok;
}
+ //================================================================================
+ /*!
+ * \brief Return true if it's possible to make a loop over corner2Sides starting
+ * from the startSide
+ */
+ //================================================================================
+
+ bool isClosedChainOfSides( _BlockSide* startSide,
+ map< const SMDS_MeshNode*, list< _BlockSide* > > & corner2Sides )
+ {
+ // get start and end nodes
+ const SMDS_MeshNode *n1 = 0, *n2 = 0, *n;
+ for ( int y = 0; y < 2; ++y )
+ for ( int x = 0; x < 2; ++x )
+ {
+ n = startSide->getCornerNode(x,y);
+ if ( !corner2Sides.count( n )) continue;
+ if ( n1 )
+ n2 = n;
+ else
+ n1 = n;
+ }
+ if ( !n2 ) return false;
+
+ map< const SMDS_MeshNode*, list< _BlockSide* > >::iterator
+ c2sides = corner2Sides.find( n1 );
+ if ( c2sides == corner2Sides.end() ) return false;
+
+ int nbChainLinks = 1;
+ n = n1;
+ _BlockSide* prevSide = startSide;
+ while ( n != n2 )
+ {
+ // get the next side sharing n
+ list< _BlockSide* > & sides = c2sides->second;
+ _BlockSide* nextSide = ( sides.back() == prevSide ? sides.front() : sides.back() );
+ if ( nextSide == prevSide ) return false;
+
+ // find the next corner of the nextSide being in corner2Sides
+ n1 = n;
+ n = 0;
+ for ( int y = 0; y < 2 && !n; ++y )
+ for ( int x = 0; x < 2; ++x )
+ {
+ n = nextSide->getCornerNode(x,y);
+ c2sides = corner2Sides.find( n );
+ if ( n == n1 || c2sides == corner2Sides.end() )
+ n = 0;
+ else
+ break;
+ }
+ if ( !n ) return false;
+
+ prevSide = nextSide;
+ nbChainLinks++;
+ }
+
+ return ( n == n2 && nbChainLinks == NB_QUAD_SIDES );
+ }
+
//================================================================================
/*!
* \brief Try to find a block side adjacent to the given side by given edge
*/
//================================================================================
- _OrientedBlockSide _Skin::findBlockSide( EBoxSides startBlockSide,
- EQuadEdge sharedSideEdge1,
- EQuadEdge sharedSideEdge2,
- bool withGeometricAnalysis)
+ _OrientedBlockSide _Skin::findBlockSide( EBoxSides startBlockSide,
+ EQuadEdge sharedSideEdge1,
+ EQuadEdge sharedSideEdge2,
+ bool withGeometricAnalysis,
+ set< _BlockSide* >& sidesAround)
{
_Block& block = _blocks.back();
_OrientedBlockSide& side1 = block._side[ startBlockSide ];
if ( !foundSide )
{
- 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 )
+ if ( !withGeometricAnalysis )
{
- _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 * PI; // angle [0-2*PI]
- angleOfSide.insert( make_pair( angle, sideI ));
- _DUMP_(" "<< sideI << " - side dir ("
- << sideIDir.X() << ", " << sideIDir.Y() << ", " << sideIDir.Z() << ")"
- << " angle " << angle);
+ sidesAround.insert( sidesOnEdge.begin(), sidesOnEdge.end() );
+ return 0;
}
if ( nbLoadedSides == 1 )
{
- double angF = angleOfSide.begin()->first, angL = angleOfSide.rbegin()->first;
- if ( angF > PI ) angF = 2*PI - angF;
- if ( angL > PI ) angL = 2*PI - angL;
- foundSide = angF < angL ? angleOfSide.begin()->second : angleOfSide.rbegin()->second;
+ // Issue 0021529. There are at least 2 sides by each edge and
+ // position of block gravity center is undefined.
+ // Find a side starting from which we can walk around the startBlockSide
+
+ // fill in corner2Sides
+ map< const SMDS_MeshNode*, list< _BlockSide* > > corner2Sides;
+ for ( sideIt = sidesAround.begin(); sideIt != sidesAround.end(); ++sideIt )
+ {
+ _BlockSide* sideI = *sideIt;
+ corner2Sides[ sideI->getCornerNode(0,0) ].push_back( sideI );
+ corner2Sides[ sideI->getCornerNode(1,0) ].push_back( sideI );
+ corner2Sides[ sideI->getCornerNode(0,1) ].push_back( sideI );
+ corner2Sides[ sideI->getCornerNode(1,1) ].push_back( sideI );
+ }
+ // remove corners of startBlockSide from corner2Sides
+ set<const SMDS_MeshNode*>::iterator nIt = block._corners.begin();
+ for ( ; nIt != block._corners.end(); ++nIt )
+ corner2Sides.erase( *nIt );
+
+ // select a side
+ for ( sideIt = sidesOnEdge.begin(); sideIt != sidesOnEdge.end(); ++sideIt )
+ {
+ if ( isClosedChainOfSides( *sideIt, corner2Sides ))
+ {
+ foundSide = *sideIt;
+ break;
+ }
+ }
+ if ( !foundSide )
+ return 0;
}
else
{
+ // 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);
+ }
+
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 ] )
// ----------------------------
// 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 ) {
bool SMESH_HexaFromSkin_3D::Compute( SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper,
std::map<HEXA_NS::Hexa*, SMESH_HexaBlocks::SMESHVolumes>& volumesOnHexa,
std::map<HEXA_NS::Vertex*, SMDS_MeshNode*> vertexNode )
-{
+ {
if(MYDEBUG) MESSAGE("SMESH_HexaFromSkin_3D::Compute BEGIN");
_Skin skin;
int nbBlocks = skin.findBlocks(aMesh);
for ( z = 0; z < zSize; ++z ) {
column0[ z ] = block.getSide(B_FRONT).node( x, z );
column1[ z ] = block.getSide(B_BACK) .node( x, z );
- }
+ }
}
// fill node columns by left and right box sides
for ( y = 1; y < ySize-1; ++y ) {
pointOnShape[ SMESH_Block::ID_V111 ] = block.getSide(B_TOP).xyz( X, Y );
for ( x = 1; x < xSize-1; ++x )
- {
+ {
gp_XYZ params; // normalized parameters of internal node within a unit box
params.SetCoord( 1, x / double(X) );
for ( y = 1; y < ySize-1; ++y )
//cout << "Params: ( "<< params.X()<<", "<<params.Y()<<", "<<params.Z()<<" )"<<endl;
//cout << "coords: ( "<< coords.X()<<", "<<coords.Y()<<", "<<coords.Z()<<" )"<<endl;
#endif
+ }
}
- }
}
// ----------------
// Add hexahedrons
aHelper->AddVolume(col00[z], col10[z], col11[z], col01[z],
col00[z+1], col10[z+1], col11[z+1], col01[z+1]);
volumesOnBlock.push_back( newVolume );
- }
}
}
+ }
// std::cout << "block i = " << i << std::endl;
HEXA_NS::Hexa* currentHexa = _block2Hexa( block, _doc, vertexNode );
if ( currentHexa != NULL ){
return true;
}
-
-
-
-
-
-
-
-
//================================================================================
/*!
* \brief Evaluate nb of hexa
return error("Algorithm can't work with geometrical shapes");
}
+
assoc != associations.end();
++assoc ){
string theBrep = (*assoc)->getBrep();
+ double ass_debut = std::min ((*assoc)->debut, (*assoc)->fin);
+ double ass_fin = std::max ((*assoc)->debut, (*assoc)->fin);
TopoDS_Shape theShape = string2shape( theBrep );
TopoDS_Edge theEdge = TopoDS::Edge( theShape );
double theCurve_length = _edgeLength( theEdge );
Handle(Geom_Curve) testCurve = BRep_Tool::Curve(theEdge, f, l);
theCurve = new BRepAdaptor_Curve( theEdge );
- GCPnts_AbscissaPoint discret_start(*theCurve, theCurve_length*(*assoc)->debut, theCurve->FirstParameter() );
- GCPnts_AbscissaPoint discret_end(*theCurve, theCurve_length*(*assoc)->fin, theCurve->FirstParameter() );
+ GCPnts_AbscissaPoint discret_start (*theCurve,
+ theCurve_length*ass_debut,
+ theCurve->FirstParameter() );
+ GCPnts_AbscissaPoint discret_end (*theCurve,
+ theCurve_length*ass_fin,
+ theCurve->FirstParameter() );
double u_start = discret_start.Parameter();
double u_end = discret_end.Parameter();
ASSERT( discret_start.IsDone() && discret_end.IsDone() );
theCurve_start = theCurve->Value( u_start);
theCurve_end = theCurve->Value( u_end );
-// double u_start = (l-f)*(*assoc)->debut;
-// double u_end = (l-f)*(*assoc)->fin;
-// theCurve_start = theCurve->Value( (l-f)*(*assoc)->debut );
-// theCurve_end = theCurve->Value( (l-f)*(*assoc)->fin );
- theCurve_length = theCurve_length*( (*assoc)->fin - (*assoc)->debut );
+ theCurve_length = theCurve_length*( ass_fin - ass_debut);
if (MYDEBUG){
MESSAGE("testCurve->f ->"<<f);
MESSAGE("FirstParameter ->"<<theCurve->FirstParameter());
MESSAGE("LastParameter ->"<<theCurve->LastParameter());
MESSAGE("theCurve_length ->"<<theCurve_length);
- MESSAGE("(*assoc)->debut ->"<<(*assoc)->debut );
- MESSAGE("(*assoc)->fin ->"<<(*assoc)->fin );
+ MESSAGE("(*assoc)->debut ->"<< ass_debut );
+ MESSAGE("(*assoc)->fin ->"<< ass_fin );
MESSAGE("u_start ->"<<u_start);
MESSAGE("u_end ->"<<u_end);
MESSAGE("myCurve_start( "<<myCurve_start.X()<<", "<<myCurve_start.Y()<<", "<<myCurve_start.Z()<<") ");
}
if ( thePreviousCurve == NULL ){
- // working on first valid association: it can be the first or last curve.
- // using myCurve_start and myCurve_end to check it out.
- // gp_Pnt theCurve_start = theCurve->Value( f + theCurveLength*assoc->debut );
-
// setting myCurve_way and first curve way
if ( myCurve_start.IsEqual(theCurve_start, HEXA_EPSILON) ){
if(MYDEBUG) MESSAGE("myCurve_start.IsEqual(theCurve_start, HEXA_EPSILON)");
GCPnts_AbscissaPoint discret;
if (MYDEBUG){
- MESSAGE("looking for curve: myCurve_u = "<<myCurve_u);
+ MESSAGE("looking for curve: c = "<<myCurve_u);
MESSAGE("looking for curve: curve_start = "<<curve_start);
MESSAGE("looking for curve: curve_end = "<<curve_end);
MESSAGE("looking for curve: curve_lenght = "<<myCurve_lengths[curve]);
+ MESSAGE("looking for curve: curve.size _lenght= "<<myCurve.size());
}
while ( not ( (myCurve_u >= curve_start) and (myCurve_u <= curve_end) ) ) {
+ if (myCurve.size() == 0 )
+ {
+ PutData (myCurve.size());
+ PutData (myCurve_u);
+ }
+
ASSERT( myCurve.size() != 0 );
myCurve.pop_front();
curve = myCurve.front();
curve_end = curve_start + myCurve_lengths[curve];
if (MYDEBUG){
MESSAGE("go next curve: curve_lenght = "<<myCurve_lengths[curve]);
- MESSAGE("go next curve: curve_start = "<<curve_start);
- MESSAGE("go next curve: curve_end = "<<curve_end);
+ MESSAGE("go next curve: curve_start = "<<curve_start);
+ MESSAGE("go next curve: curve_end = "<<curve_end);
+ MESSAGE("go next curve: myCurve_u = "<<myCurve_u);
}
}
myCurve_start = curve_start;
