From a4d5a23480fdb1aba7fe76a0d15774148598e09c Mon Sep 17 00:00:00 2001 From: fps Date: Wed, 18 Jul 2012 09:46:32 +0000 Subject: [PATCH] Livraison CS du 18/07/2012 --- configure.ac | 2 +- .../HEXABLOCKPlugin_FromSkin_3D.cxx | 283 +++++++++++++----- src/HEXABLOCKPlugin/HEXABLOCKPlugin_mesh.cxx | 38 ++- 3 files changed, 229 insertions(+), 94 deletions(-) diff --git a/configure.ac b/configure.ac index 0bc8458..1fec193 100755 --- a/configure.ac +++ b/configure.ac @@ -35,7 +35,7 @@ AM_INIT_AUTOMAKE([-Wno-portability]) XVERSION=`echo $VERSION | awk -F. '{printf("0x%02x%02x%02x",$1,$2,$3)}'` AC_SUBST(XVERSION) -VERSION_DEV=1 +VERSION_DEV=0 AC_SUBST(VERSION_DEV) # set up MODULE_NAME variable for dynamic construction of directories (resources, etc.) diff --git a/src/HEXABLOCKPlugin/HEXABLOCKPlugin_FromSkin_3D.cxx b/src/HEXABLOCKPlugin/HEXABLOCKPlugin_FromSkin_3D.cxx index 43d201b..1dc0e31 100755 --- a/src/HEXABLOCKPlugin/HEXABLOCKPlugin_FromSkin_3D.cxx +++ b/src/HEXABLOCKPlugin/HEXABLOCKPlugin_FromSkin_3D.cxx @@ -1,4 +1,4 @@ -// 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 @@ -31,19 +31,25 @@ #include //#include "utilities.h" +#include -// 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_ @@ -59,10 +65,12 @@ namespace { 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 @@ -194,9 +202,9 @@ namespace 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); } }; //================================================================================ /*! @@ -230,7 +238,7 @@ namespace 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 )); } @@ -344,10 +352,11 @@ namespace vector& verRow1, vector& 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 ) { @@ -510,20 +519,25 @@ namespace 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 ); @@ -663,16 +677,77 @@ namespace 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 ]; @@ -723,45 +798,75 @@ namespace 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::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 ] ) @@ -1061,7 +1166,7 @@ bool SMESH_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHel // ---------------------------- // 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 pointOnShape( SMESH_Block::ID_Shell ); @@ -1133,18 +1238,49 @@ bool SMESH_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHel // 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::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 ) { @@ -1178,7 +1314,7 @@ bool SMESH_HexaFromSkin_3D::Compute(SMESH_Mesh & aMesh, SMESH_MesherHelper* aHel bool SMESH_HexaFromSkin_3D::Compute( SMESH_Mesh & aMesh, SMESH_MesherHelper* aHelper, std::map& volumesOnHexa, std::map vertexNode ) -{ + { if(MYDEBUG) MESSAGE("SMESH_HexaFromSkin_3D::Compute BEGIN"); _Skin skin; int nbBlocks = skin.findBlocks(aMesh); @@ -1213,7 +1349,7 @@ bool SMESH_HexaFromSkin_3D::Compute( SMESH_Mesh & aMesh, SMESH_MesherHelper* aHe 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 ) { @@ -1254,7 +1390,7 @@ bool SMESH_HexaFromSkin_3D::Compute( SMESH_Mesh & aMesh, SMESH_MesherHelper* aHe 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 ) @@ -1304,8 +1440,8 @@ bool SMESH_HexaFromSkin_3D::Compute( SMESH_Mesh & aMesh, SMESH_MesherHelper* aHe //cout << "Params: ( "<< params.X()<<", "<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 ){ @@ -1370,14 +1506,6 @@ bool SMESH_HexaFromSkin_3D::Compute( SMESH_Mesh & aMesh, SMESH_MesherHelper* aHe return true; } - - - - - - - - //================================================================================ /*! * \brief Evaluate nb of hexa @@ -1449,3 +1577,4 @@ bool SMESH_HexaFromSkin_3D::Compute(SMESH_Mesh&, const TopoDS_Shape&) return error("Algorithm can't work with geometrical shapes"); } + diff --git a/src/HEXABLOCKPlugin/HEXABLOCKPlugin_mesh.cxx b/src/HEXABLOCKPlugin/HEXABLOCKPlugin_mesh.cxx index c16bfd1..b2a4f9f 100755 --- a/src/HEXABLOCKPlugin/HEXABLOCKPlugin_mesh.cxx +++ b/src/HEXABLOCKPlugin/HEXABLOCKPlugin_mesh.cxx @@ -1274,6 +1274,8 @@ void SMESH_HexaBlocks::_buildMyCurve( 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 ); @@ -1285,18 +1287,18 @@ void SMESH_HexaBlocks::_buildMyCurve( 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 ->"<"<FirstParameter()); MESSAGE("LastParameter ->"<LastParameter()); MESSAGE("theCurve_length ->"<debut ->"<<(*assoc)->debut ); - MESSAGE("(*assoc)->fin ->"<<(*assoc)->fin ); + MESSAGE("(*assoc)->debut ->"<< ass_debut ); + MESSAGE("(*assoc)->fin ->"<< ass_fin ); MESSAGE("u_start ->"<"<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)"); @@ -1414,12 +1412,19 @@ gp_Pnt SMESH_HexaBlocks::_getPtOnMyCurve( GCPnts_AbscissaPoint discret; if (MYDEBUG){ - MESSAGE("looking for curve: myCurve_u = "<