From b0d030f73dace135d1b2bc520e4b532d7d150a98 Mon Sep 17 00:00:00 2001 From: eap Date: Tue, 24 Dec 2013 07:04:11 +0000 Subject: [PATCH] 22360: EDF SMESH: Body Fitting algorithm: incorporate edges --- idl/SMESH_BasicHypothesis.idl | 15 +- .../StdMeshers_CartesianParameters3D.cxx | 93 +- .../StdMeshers_CartesianParameters3D.hxx | 15 + src/StdMeshers/StdMeshers_Cartesian_3D.cxx | 1137 ++++++++++++++--- .../StdMeshersGUI_CartesianParamCreator.cxx | 20 +- .../StdMeshersGUI_CartesianParamCreator.h | 20 +- src/StdMeshersGUI/StdMeshers_msg_en.ts | 4 + src/StdMeshersGUI/StdMeshers_msg_fr.ts | 44 + .../StdMeshers_CartesianParameters3D_i.cxx | 22 + .../StdMeshers_CartesianParameters3D_i.hxx | 8 + 10 files changed, 1188 insertions(+), 190 deletions(-) diff --git a/idl/SMESH_BasicHypothesis.idl b/idl/SMESH_BasicHypothesis.idl index aa8940da8..500c3fba5 100644 --- a/idl/SMESH_BasicHypothesis.idl +++ b/idl/SMESH_BasicHypothesis.idl @@ -952,7 +952,7 @@ module StdMeshers /*! * Set size threshold. A polyhedral cell got by cutting an initial * hexahedron by geometry boundary is considered small and is removed if - * it's size is \athreshold times less than the size of the initial hexahedron. + * it's size is \a threshold times less than the size of the initial hexahedron. * threshold must be > 1.0 */ void SetSizeThreshold(in double threshold) raises (SALOME::SALOME_Exception); @@ -987,6 +987,13 @@ module StdMeshers void GetGridSpacing(out SMESH::string_array spaceFunctions, out SMESH::double_array internalPoints, in short axis) raises (SALOME::SALOME_Exception); + /*! + * Enables implementation of geometrical edges into the mesh. If this feature + * is disabled, sharp edges of the shape are lost ("smoothed") in the mesh if + * they don't coincide with the grid lines + */ + void SetToAddEdges(in boolean toAdd); + boolean GetToAddEdges(); /*! * \brief Computes node coordinates by spacing functions @@ -994,13 +1001,15 @@ module StdMeshers * \param x1 - upper coordinate * \param spaceFuns - space functions * \param points - internal points - * \param coords - the computed coordinates + * \param axisName - e.g. "X" + * \return the computed coordinates */ SMESH::double_array ComputeCoordinates(in double x0, in double x1, in SMESH::string_array spaceFuns, in SMESH::double_array points, - in string axisName ) raises (SALOME::SALOME_Exception); + in string axisName ) + raises (SALOME::SALOME_Exception); }; /*! diff --git a/src/StdMeshers/StdMeshers_CartesianParameters3D.cxx b/src/StdMeshers/StdMeshers_CartesianParameters3D.cxx index fa4aa909f..369d4e424 100644 --- a/src/StdMeshers/StdMeshers_CartesianParameters3D.cxx +++ b/src/StdMeshers/StdMeshers_CartesianParameters3D.cxx @@ -32,11 +32,12 @@ #include "utilities.h" -#include -#include - #include +#include +#include +#include + using namespace std; //======================================================================= @@ -48,10 +49,23 @@ StdMeshers_CartesianParameters3D::StdMeshers_CartesianParameters3D(int h int studyId, SMESH_Gen * gen) : SMESH_Hypothesis(hypId, studyId, gen), - _sizeThreshold( 4.0 ) // default according to the customer specification + _sizeThreshold( 4.0 ), // default according to the customer specification + _toAddEdges( false ) { _name = "CartesianParameters3D"; // used by "Cartesian_3D" _param_algo_dim = 3; // 3D + + _axisDirs[0] = 1.; + _axisDirs[1] = 0.; + _axisDirs[2] = 0.; + + _axisDirs[3] = 0.; + _axisDirs[4] = 1.; + _axisDirs[5] = 0.; + + _axisDirs[6] = 0.; + _axisDirs[7] = 0.; + _axisDirs[8] = 1.; } @@ -308,6 +322,44 @@ void StdMeshers_CartesianParameters3D::GetCoordinates(std::vector& xNode zNodes = _coords[2]; } +//======================================================================= +//function : SetAxisDirs +//purpose : Sets directions of axes +//======================================================================= + +void StdMeshers_CartesianParameters3D::SetAxisDirs(const double* the9DirComps) + throw ( SALOME_Exception ) +{ + gp_Vec x( the9DirComps[0], + the9DirComps[1], + the9DirComps[2] ); + gp_Vec y( the9DirComps[3], + the9DirComps[4], + the9DirComps[5] ); + gp_Vec z( the9DirComps[6], + the9DirComps[7], + the9DirComps[8] ); + if ( x.Magnitude() < RealSmall() || + y.Magnitude() < RealSmall() || + z.Magnitude() < RealSmall() ) + throw SALOME_Exception("Zero magnitude of axis direction"); + + if ( x.IsParallel( y, M_PI / 180. ) || + x.IsParallel( z, M_PI / 180. ) || + y.IsParallel( z, M_PI / 180. )) + throw SALOME_Exception("Parallel axis directions"); + + bool isChanged = false; + for ( int i = 0; i < 9; ++i ) + { + if ( Abs( _axisDirs[i] - the9DirComps[i] ) > 1e-7 ) + isChanged = true; + _axisDirs[i] = the9DirComps[i]; + } + if ( isChanged ) + NotifySubMeshesHypothesisModification(); +} + //======================================================================= //function : GetGrid //purpose : Return coordinates of node positions along the three axes @@ -332,6 +384,33 @@ double StdMeshers_CartesianParameters3D::GetSizeThreshold() const return _sizeThreshold; } +//======================================================================= +//function : SetToAddEdges +//purpose : Enables implementation of geometrical edges into the mesh. If this feature +// is disabled, sharp edges of the shape are lost ("smoothed") in the mesh if +// they don't coincide with the grid lines +//======================================================================= + +void StdMeshers_CartesianParameters3D::SetToAddEdges(bool toAdd) +{ + if ( _toAddEdges != toAdd ) + { + _toAddEdges = toAdd; + NotifySubMeshesHypothesisModification(); + } +} + +//======================================================================= +//function : GetToAddEdges +//purpose : Returns true if implementation of geometrical edges into the +// mesh is enabled +//======================================================================= + +bool StdMeshers_CartesianParameters3D::GetToAddEdges() const +{ + return _toAddEdges; +} + //======================================================================= //function : IsDefined //purpose : Return true if parameters are well defined @@ -369,6 +448,7 @@ std::ostream & StdMeshers_CartesianParameters3D::SaveTo(std::ostream & save) for ( size_t j = 0; j < _spaceFunctions[i].size(); ++j ) save << _spaceFunctions[i][j] << " "; } + save << _toAddEdges << " "; return save; } @@ -382,7 +462,7 @@ std::istream & StdMeshers_CartesianParameters3D::LoadFrom(std::istream & load) { bool ok; - ok = (load >> _sizeThreshold ); + ok = ( load >> _sizeThreshold ); for ( int ax = 0; ax < 3; ++ax ) { if (ok) @@ -419,6 +499,9 @@ std::istream & StdMeshers_CartesianParameters3D::LoadFrom(std::istream & load) } } } + + load >> _toAddEdges; + return load; } diff --git a/src/StdMeshers/StdMeshers_CartesianParameters3D.hxx b/src/StdMeshers/StdMeshers_CartesianParameters3D.hxx index adb22f00c..98970956b 100644 --- a/src/StdMeshers/StdMeshers_CartesianParameters3D.hxx +++ b/src/StdMeshers/StdMeshers_CartesianParameters3D.hxx @@ -100,6 +100,10 @@ public: std::vector& yNodes, std::vector& zNodes, const Bnd_Box& bndBox) const throw ( SALOME_Exception ); + + void SetAxisDirs(const double* the9DirComps) throw ( SALOME_Exception ); + const double* GetAxisDirs() const { return _axisDirs; } + /*! * Set size threshold. A polyhedral cell got by cutting an initial * hexahedron by geometry boundary is considered small and is removed if @@ -111,6 +115,14 @@ public: */ double GetSizeThreshold() const; + /*! + * \brief Enables implementation of geometrical edges into the mesh. If this feature + * is disabled, sharp edges of the shape are lost ("smoothed") in the mesh if + * they don't coincide with the grid lines + */ + void SetToAddEdges(bool toAdd); + bool GetToAddEdges() const; + /*! * \brief Return true if parameters are well defined */ @@ -138,7 +150,10 @@ public: std::vector _spaceFunctions[3]; std::vector _internalPoints[3]; + double _axisDirs[9]; + double _sizeThreshold; + bool _toAddEdges; }; #endif diff --git a/src/StdMeshers/StdMeshers_Cartesian_3D.cxx b/src/StdMeshers/StdMeshers_Cartesian_3D.cxx index ad13adede..b94506f68 100644 --- a/src/StdMeshers/StdMeshers_Cartesian_3D.cxx +++ b/src/StdMeshers/StdMeshers_Cartesian_3D.cxx @@ -37,6 +37,7 @@ #include #include +#include #include #include #include @@ -44,6 +45,7 @@ #include #include #include +#include #include #include #include @@ -63,7 +65,6 @@ #include #include #include -#include #include #include #include @@ -76,7 +77,7 @@ #include #include -//#undef WITH_TBB +#undef WITH_TBB #ifdef WITH_TBB #include //#include @@ -150,6 +151,8 @@ bool StdMeshers_Cartesian_3D::CheckHypothesis (SMESH_Mesh& aMesh, namespace { + typedef int TGeomID; + //============================================================================= // Definitions of internal utils // -------------------------------------------------------------------------- @@ -161,17 +164,40 @@ namespace }; // -------------------------------------------------------------------------- /*! - * \brief Data of intersection between a GridLine and a TopoDS_Face + * \brief Common data of any intersection between a Grid and a shape */ - struct IntersectionPoint + struct B_IntersectPoint { - double _paramOnLine; - mutable Transition _transition; mutable const SMDS_MeshNode* _node; - mutable size_t _indexOnLine; + mutable vector< TGeomID > _faceIDs; + + B_IntersectPoint(): _node(NULL) {} + void Add( const vector< TGeomID >& fIDs, const SMDS_MeshNode* n=0 ) const; + bool HasCommonFace( const B_IntersectPoint * other ) const; + bool IsOnFace( int faceID ) const; + virtual ~B_IntersectPoint() {} + }; + // -------------------------------------------------------------------------- + /*! + * \brief Data of intersection between a GridLine and a TopoDS_Face + */ + struct F_IntersectPoint : public B_IntersectPoint + { + double _paramOnLine; + mutable Transition _transition; + mutable size_t _indexOnLine; - IntersectionPoint(): _node(0) {} - bool operator< ( const IntersectionPoint& o ) const { return _paramOnLine < o._paramOnLine; } + bool operator< ( const F_IntersectPoint& o ) const { return _paramOnLine < o._paramOnLine; } + }; + // -------------------------------------------------------------------------- + /*! + * \brief Data of intersection between GridPlanes and a TopoDS_EDGE + */ + struct E_IntersectPoint : public B_IntersectPoint + { + gp_Pnt _point; + double _uvw[3]; + TGeomID _shapeID; }; // -------------------------------------------------------------------------- /*! @@ -181,10 +207,23 @@ namespace { gp_Lin _line; double _length; // line length - multiset< IntersectionPoint > _intPoints; + multiset< F_IntersectPoint > _intPoints; void RemoveExcessIntPoints( const double tol ); - bool GetIsOutBefore( multiset< IntersectionPoint >::iterator ip, bool prevIsOut ); + bool GetIsOutBefore( multiset< F_IntersectPoint >::iterator ip, bool prevIsOut ); + }; + // -------------------------------------------------------------------------- + /*! + * \brief Planes of the grid used to find intersections of an EDGE with a hexahedron + */ + struct GridPlanes + { + double _factor; + gp_XYZ _uNorm, _vNorm, _zNorm; + vector< gp_XYZ > _origins; // origin points of all planes in one direction + vector< double > _zProjs; // projections of origins to _zNorm + + gp_XY GetUV( const gp_Pnt& p, const gp_Pnt& origin ); }; // -------------------------------------------------------------------------- /*! @@ -234,11 +273,15 @@ namespace struct Grid { vector< double > _coords[3]; // coordinates of grid nodes + gp_XYZ _axes [3]; // axis directions vector< GridLine > _lines [3]; // in 3 directions double _tol, _minCellSize; - vector< const SMDS_MeshNode* > _nodes; // mesh nodes at grid nodes - vector< bool > _isBndNode; // is mesh node at intersection with geometry + vector< const SMDS_MeshNode* > _nodes; // mesh nodes at grid nodes + vector< const F_IntersectPoint* > _gridIntP; // grid node intersection with geometry + + list< E_IntersectPoint > _edgeIntP; // intersections with EDGEs + TopTools_IndexedMapOfShape _shapes; size_t CellIndex( size_t i, size_t j, size_t k ) const { @@ -257,6 +300,7 @@ namespace void SetCoordinates(const vector& xCoords, const vector& yCoords, const vector& zCoords, + const double* axesDirs, const TopoDS_Shape& shape ); void ComputeNodes(SMESH_MesherHelper& helper); }; @@ -302,10 +346,11 @@ namespace struct FaceGridIntersector { TopoDS_Face _face; + TGeomID _faceID; Grid* _grid; Bnd_Box _bndBox; __IntCurvesFace_Intersector* _surfaceInt; - vector< std::pair< GridLine*, IntersectionPoint > > _intersections; + vector< std::pair< GridLine*, F_IntersectPoint > > _intersections; FaceGridIntersector(): _grid(0), _surfaceInt(0) {} void Intersect(); @@ -314,7 +359,12 @@ namespace void StoreIntersections() { for ( size_t i = 0; i < _intersections.size(); ++i ) - _intersections[i].first->_intPoints.insert( _intersections[i].second ); + { + multiset< F_IntersectPoint >::iterator ip = + _intersections[i].first->_intPoints.insert( _intersections[i].second ); + ip->_faceIDs.reserve( 1 ); + ip->_faceIDs.push_back( _faceID ); + } } const Bnd_Box& GetFaceBndBox() { @@ -352,7 +402,7 @@ namespace gp_Torus _torus; __IntCurvesFace_Intersector* _surfaceInt; - vector< IntersectionPoint > _intPoints; + vector< F_IntersectPoint > _intPoints; void IntersectWithPlane (const GridLine& gridLine); void IntersectWithCylinder(const GridLine& gridLine); @@ -381,20 +431,55 @@ namespace struct _Face; struct _Link; // -------------------------------------------------------------------------------- - struct _Node //!< node either at a hexahedron corner or at GridLine intersection + struct _Node //!< node either at a hexahedron corner or at intersection { - const SMDS_MeshNode* _node; // mesh node at hexahedron corner - const IntersectionPoint* _intPoint; - - _Node(const SMDS_MeshNode* n=0, const IntersectionPoint* ip=0):_node(n), _intPoint(ip) {} - const SMDS_MeshNode* Node() const { return _intPoint ? _intPoint->_node : _node; } - //bool IsCorner() const { return _node; } + const SMDS_MeshNode* _node; // mesh node at hexahedron corner + const B_IntersectPoint* _intPoint; + + _Node(const SMDS_MeshNode* n=0, const B_IntersectPoint* ip=0):_node(n), _intPoint(ip) {} + const SMDS_MeshNode* Node() const + { return _intPoint ? _intPoint->_node : _node; } + const F_IntersectPoint* FaceIntPnt() const + { return static_cast< const F_IntersectPoint* >( _intPoint ); } + const E_IntersectPoint* EdgeIntPnt() const + { return static_cast< const E_IntersectPoint* >( _intPoint ); } + void Add( const E_IntersectPoint* ip ) + { + if ( !_intPoint ) { + _intPoint = ip; + } + else if ( !_intPoint->_node ) { + ip->Add( _intPoint->_faceIDs ); + _intPoint = ip; + } + else { + _intPoint->Add( ip->_faceIDs ); + } + } + bool IsLinked( const B_IntersectPoint* other ) const + { + // node inside a SOLID is considered "linked" with any other node + return ( !other || !_intPoint ) ? true : _intPoint->HasCommonFace( other ); + } + bool IsOnFace( int faceID ) const // returns true if faceID is found + { + return _intPoint ? _intPoint->IsOnFace( faceID ) : false; + } + gp_Pnt Point() const + { + if ( const SMDS_MeshNode* n = Node() ) + return SMESH_TNodeXYZ( n ); + if ( const E_IntersectPoint* eip = + dynamic_cast< const E_IntersectPoint* >( _intPoint )) + return eip->_point; + return gp_Pnt( 1e100, 0, 0 ); + } }; // -------------------------------------------------------------------------------- struct _Link // link connecting two _Node's { _Node* _nodes[2]; - vector< _Node> _intNodes; // _Node's at GridLine intersections + vector< _Node > _intNodes; // _Node's at GridLine intersections vector< _Link > _splits; vector< _Face*> _faces; }; @@ -412,21 +497,50 @@ namespace } _Node* FirstNode() const { return _link->_nodes[ _reverse ]; } _Node* LastNode() const { return _link->_nodes[ !_reverse ]; } + vector< TGeomID > GetNotUsedFace(const set& usedIDs ) const // returns a supporting FACEs + { + vector< TGeomID > faces; + const B_IntersectPoint *ip0, *ip1; + if (( ip0 = _link->_nodes[0]->_intPoint ) && + ( ip1 = _link->_nodes[1]->_intPoint )) + { + for ( size_t i = 0; i < ip0->_faceIDs.size(); ++i ) + if ( ip1->IsOnFace ( ip0->_faceIDs[i] ) && + !usedIDs.count( ip0->_faceIDs[i] ) ) + faces.push_back( ip0->_faceIDs[i] ); + } + return faces; + } + // TGeomID GetNotUsedFace( set& usedIDs ) const // returns a supporting FACE + // { + // const B_IntersectPoint *ip0, *ip1; + // if (( ip0 = _link->_nodes[0]->_intPoint ) && + // ( ip1 = _link->_nodes[1]->_intPoint )) + // { + // for ( size_t i = 0; i < ip0->_faceIDs.size(); ++i ) + // if ( ip1->IsOnFace ( ip0->_faceIDs[i] ) && + // usedIDs.insert( ip0->_faceIDs[i] ).second ) + // return ip0->_faceIDs[i]; + // } + // return -100; + // } }; // -------------------------------------------------------------------------------- struct _Face { - vector< _OrientedLink > _links; - vector< _Link > _polyLinks; // links added to close a polygonal face + vector< _OrientedLink > _links; // links on GridLine's + vector< _Link > _polyLinks; // links added to close a polygonal face + vector< _Node > _edgeNodes; // nodes at intersection with EDGEs }; // -------------------------------------------------------------------------------- struct _volumeDef // holder of nodes of a volume mesh element { - vector< const SMDS_MeshNode* > _nodes; - vector< int > _quantities; + //vector< const SMDS_MeshNode* > _nodes; + vector< _Node* > _nodes; + vector< int > _quantities; typedef boost::shared_ptr<_volumeDef> Ptr; - void set( const vector< const SMDS_MeshNode* >& nodes, - const vector< int > quant = vector< int >() ) + void set( const vector< _Node* >& nodes, + const vector< int >& quant = vector< int >() ) { _nodes = nodes; _quantities = quant; } // static Ptr New( const vector< const SMDS_MeshNode* >& nodes, // const vector< int > quant = vector< int >() ) @@ -440,13 +554,19 @@ namespace // topology of a hexahedron int _nodeShift[8]; - _Node _hexNodes[8]; - _Link _hexLinks[12]; - _Face _hexQuads[6]; + _Node _hexNodes [8]; + _Link _hexLinks [12]; + _Face _hexQuads [6]; // faces resulted from hexahedron intersection vector< _Face > _polygons; + // intresections with EDGEs + vector< const E_IntersectPoint* > _edgeIntPnts; + + // nodes inside the hexahedron (at VERTEXes) + vector< _Node > _vertexNodes; + // computed volume elements //vector< _volumeDef::Ptr > _volumeDefs; _volumeDef _volumeDefs; @@ -459,7 +579,8 @@ namespace public: Hexahedron(const double sizeThreshold, Grid* grid); - int MakeElements(SMESH_MesherHelper& helper); + int MakeElements(SMESH_MesherHelper& helper, + const map< TGeomID, vector< TGeomID > >& edge2faceIDsMap); void ComputeElements(); void Init() { init( _i, _j, _k ); } @@ -467,6 +588,12 @@ namespace Hexahedron(const Hexahedron& other ); void init( size_t i, size_t j, size_t k ); void init( size_t i ); + void addEdges(SMESH_MesherHelper& helper, + vector< Hexahedron* >& intersectedHex, + const map< TGeomID, vector< TGeomID > >& edge2faceIDsMap); + int getEntity( const E_IntersectPoint* ip, int* facets, int& sub ); + void addIntersection( const E_IntersectPoint& ip ); + bool findChain( _Node* n1, _Node* n2, _Face& quad, vector<_Node*>& chainNodes ); int addElements(SMESH_MesherHelper& helper); bool isInHole() const; bool checkPolyhedronSize() const; @@ -474,8 +601,17 @@ namespace bool addTetra(); bool addPenta(); bool addPyra (); + _Node* FindEqualNode( vector< _Node >& nodes, + const E_IntersectPoint* ip, + const double tol2 ) + { + for ( size_t i = 0; i < nodes.size(); ++i ) + if ( nodes[i].Point().SquareDistance( ip->_point ) <= tol2 ) + return & nodes[i]; + return 0; + } }; - + #ifdef WITH_TBB // -------------------------------------------------------------------------- /*! @@ -507,11 +643,26 @@ namespace _faceVec[i].Intersect(); } }; - #endif + //============================================================================= // Implementation of internal utils //============================================================================= + /*! + * \brief adjust \a i to have \a val between values[i] and values[i+1] + */ + inline void locateValue( int & i, double val, const vector& values ) + { + val += values[0]; // input \a val is measured from 0. + if ( i > values.size()-2 ) + i = values.size()-2; + else + while ( i+2 < values.size() && val > values[ i+1 ]) + ++i; + while ( i > 0 && val < values[ i ]) + --i; + } + //============================================================================= /* * Remove coincident intersection points */ @@ -520,7 +671,7 @@ namespace if ( _intPoints.size() < 2 ) return; set< Transition > tranSet; - multiset< IntersectionPoint >::iterator ip1, ip2 = _intPoints.begin(); + multiset< F_IntersectPoint >::iterator ip1, ip2 = _intPoints.begin(); while ( ip2 != _intPoints.end() ) { tranSet.clear(); @@ -529,6 +680,7 @@ namespace { tranSet.insert( ip1->_transition ); tranSet.insert( ip2->_transition ); + ip2->Add( ip1->_faceIDs ); _intPoints.erase( ip1 ); ip1 = ip2++; } @@ -547,7 +699,7 @@ namespace /* * Return "is OUT" state for nodes before the given intersection point */ - bool GridLine::GetIsOutBefore( multiset< IntersectionPoint >::iterator ip, bool prevIsOut ) + bool GridLine::GetIsOutBefore( multiset< F_IntersectPoint >::iterator ip, bool prevIsOut ) { if ( ip->_transition == Trans_IN ) return true; @@ -558,7 +710,7 @@ namespace // singularity point (apex of a cone) if ( _intPoints.size() == 1 || ip == _intPoints.begin() ) return true; - multiset< IntersectionPoint >::iterator ipBef = ip, ipAft = ++ip; + multiset< F_IntersectPoint >::iterator ipBef = ip, ipAft = ++ip; if ( ipAft == _intPoints.end() ) return false; --ipBef; @@ -569,13 +721,65 @@ namespace return prevIsOut; // _transition == Trans_TANGENT } //================================================================================ + /* + * Returns parameters of a point in i-th plane + */ + gp_XY GridPlanes::GetUV( const gp_Pnt& p, const gp_Pnt& origin ) + { + gp_Vec v( origin, p ); + return gp_XY( v.Dot( _uNorm ) * _factor, + v.Dot( _vNorm ) * _factor ); + } + //================================================================================ + /* + * Adds face IDs + */ + void B_IntersectPoint::Add( const vector< TGeomID >& fIDs, + const SMDS_MeshNode* n) const + { + if ( _faceIDs.empty() ) + _faceIDs = fIDs; + else + for ( size_t i = 0; i < fIDs.size(); ++i ) + { + vector< TGeomID >::iterator it = + std::find( _faceIDs.begin(), _faceIDs.end(), fIDs[i] ); + if ( it == _faceIDs.end() ) + _faceIDs.push_back( fIDs[i] ); + } + if ( !_node ) + _node = n; + } + //================================================================================ + /* + * Returns \c true if \a other B_IntersectPoint holds the same face ID + */ + bool B_IntersectPoint::HasCommonFace( const B_IntersectPoint * other ) const + { + if ( other ) + for ( size_t i = 0; i < other->_faceIDs.size(); ++i ) + if ( IsOnFace( other->_faceIDs[i] ) ) + return true; + return false; + } + //================================================================================ + /* + * Returns \c true if \a faceID in in this->_faceIDs + */ + bool B_IntersectPoint::IsOnFace( int faceID ) const // returns true if faceID is found + { + vector< TGeomID >::const_iterator it = + std::find( _faceIDs.begin(), _faceIDs.end(), faceID ); + return ( it != _faceIDs.end() ); + } + //================================================================================ /* * Return an iterator on GridLine's in a given direction */ LineIndexer Grid::GetLineIndexer(size_t iDir) const { const size_t indices[] = { 1,2,0, 0,2,1, 0,1,2 }; - const string s[] = { "X", "Y", "Z" }; + const string s [] = { "X", "Y", "Z" }; LineIndexer li( _coords[0].size(), _coords[1].size(), _coords[2].size(), indices[iDir*3], indices[iDir*3+1], indices[iDir*3+2], s[indices[iDir*3]], s[indices[iDir*3+1]], s[indices[iDir*3+2]]); @@ -588,11 +792,21 @@ namespace void Grid::SetCoordinates(const vector& xCoords, const vector& yCoords, const vector& zCoords, + const double* axesDirs, const TopoDS_Shape& shape) { _coords[0] = xCoords; _coords[1] = yCoords; _coords[2] = zCoords; + _axes[0].SetCoord( axesDirs[0], + axesDirs[1], + axesDirs[2]); + _axes[1].SetCoord( axesDirs[3], + axesDirs[4], + axesDirs[5]); + _axes[2].SetCoord( axesDirs[6], + axesDirs[7], + axesDirs[8]); // compute tolerance _minCellSize = Precision::Infinite(); @@ -649,7 +863,7 @@ namespace const size_t nbGridNodes = _coords[0].size() * _coords[1].size() * _coords[2].size(); vector< bool > isNodeOut( nbGridNodes, false ); _nodes.resize( nbGridNodes, 0 ); - _isBndNode.resize( nbGridNodes, false ); + _gridIntP.resize( nbGridNodes, NULL ); for ( int iDir = 0; iDir < 3; ++iDir ) // loop on 3 line directions { @@ -669,8 +883,8 @@ namespace GridLine& line = _lines[ iDir ][ li.LineIndex() ]; line.RemoveExcessIntPoints( _tol ); - multiset< IntersectionPoint >& intPnts = _lines[ iDir ][ li.LineIndex() ]._intPoints; - multiset< IntersectionPoint >::iterator ip = intPnts.begin(); + multiset< F_IntersectPoint >& intPnts = _lines[ iDir ][ li.LineIndex() ]._intPoints; + multiset< F_IntersectPoint >::iterator ip = intPnts.begin(); bool isOut = true; const double* nodeCoord = & coords[0], *coord0 = nodeCoord, *coordEnd = coord0 + coords.size(); @@ -706,14 +920,18 @@ namespace else { int nodeIndex = nIndex0 + nShift * ( nodeCoord-coord0 ); - if ( ! _nodes[ nodeIndex ] ) + if ( !_nodes[ nodeIndex ] ) { li.SetIndexOnLine( nodeCoord-coord0 ); double xyz[3] = { _coords[0][ li.I() ], _coords[1][ li.J() ], _coords[2][ li.K() ]}; - _nodes[ nodeIndex ] = helper.AddNode( xyz[0], xyz[1], xyz[2] ); - _isBndNode[ nodeIndex ] = true; + _nodes [ nodeIndex ] = helper.AddNode( xyz[0], xyz[1], xyz[2] ); + _gridIntP[ nodeIndex ] = & * ip; } - //ip->_node = _nodes[ nodeIndex ]; + if ( _gridIntP[ nodeIndex ] ) + _gridIntP[ nodeIndex ]->Add( ip->_faceIDs ); + else + _gridIntP[ nodeIndex ] = & * ip; + // ip->_node = _nodes[ nodeIndex ]; -- to differ from ip on links ip->_indexOnLine = nodeCoord-coord0; if ( ++nodeCoord < coordEnd ) nodeParam = *nodeCoord - *coord0; @@ -744,7 +962,7 @@ namespace LineIndexer li = GetLineIndexer( iDir ); for ( ; li.More(); ++li ) { - multiset< IntersectionPoint >& intPnts = _lines[ iDir ][ li.LineIndex() ]._intPoints; + multiset< F_IntersectPoint >& intPnts = _lines[ iDir ][ li.LineIndex() ]._intPoints; if ( intPnts.empty() ) continue; if ( intPnts.size() == 1 ) { @@ -944,7 +1162,7 @@ namespace { if ( !toClassify || UVIsOnFace() ) { - IntersectionPoint p; + F_IntersectPoint p; p._paramOnLine = _w; p._transition = _transition; _intPoints.push_back( p ); @@ -1291,16 +1509,21 @@ namespace _origNodeInd = _grid->NodeIndex( i,j,k ); for ( int iN = 0; iN < 8; ++iN ) { - _hexNodes[iN]._node = _grid->_nodes[ _origNodeInd + _nodeShift[iN] ]; + _hexNodes[iN]._node = _grid->_nodes [ _origNodeInd + _nodeShift[iN] ]; + _hexNodes[iN]._intPoint = _grid->_gridIntP[ _origNodeInd + _nodeShift[iN] ]; _nbCornerNodes += bool( _hexNodes[iN]._node ); - _nbBndNodes += _grid->_isBndNode[ _origNodeInd + _nodeShift[iN] ]; + if ( _hexNodes[iN]._intPoint ) + { + ++_nbBndNodes; + _hexNodes[iN]._intPoint->_node = _hexNodes[iN]._node; + } } _sideLength[0] = _grid->_coords[0][i+1] - _grid->_coords[0][i]; _sideLength[1] = _grid->_coords[1][j+1] - _grid->_coords[1][j]; _sideLength[2] = _grid->_coords[2][k+1] - _grid->_coords[2][k]; - if ( _nbCornerNodes < 8 && _nbIntNodes + _nbCornerNodes > 3) + if ( _nbCornerNodes < 8 && _nbIntNodes + _nbCornerNodes + _edgeIntPnts.size() > 3) { _Link split; // create sub-links (_splits) by splitting links with _intNodes @@ -1324,6 +1547,97 @@ namespace link._splits.push_back( split ); } } + + // create _Node's at intersections with EDGEs + + const double tol2 = _grid->_tol * _grid->_tol; + int facets[3], nbFacets, subEntity; + + for ( size_t iP = 0; iP < _edgeIntPnts.size(); ++iP ) + { + nbFacets = getEntity( _edgeIntPnts[iP], facets, subEntity ); + _Node* equalNode = 0; + switch( nbFacets ) { + case 1: // in a _Face + { + _Face& quad = _hexQuads[ facets[0] - SMESH_Block::ID_FirstF ]; + equalNode = FindEqualNode( quad._edgeNodes, _edgeIntPnts[ iP ], tol2 ); + if ( equalNode ) { + equalNode->Add( _edgeIntPnts[ iP ] ); + } + else { + quad._edgeNodes.push_back( _Node( 0, _edgeIntPnts[ iP ])); + ++_nbIntNodes; + } + break; + } + case 2: // on a _Link + { + _Link& link = _hexLinks[ subEntity - SMESH_Block::ID_FirstE ]; + if ( link._splits.size() > 0 ) + { + equalNode = FindEqualNode( link._intNodes, _edgeIntPnts[ iP ], tol2 ); + if ( equalNode ) + equalNode->Add( _edgeIntPnts[ iP ] ); + } + else + { + for ( int iF = 0; iF < 2; ++iF ) + { + _Face& quad = _hexQuads[ facets[iF] - SMESH_Block::ID_FirstF ]; + equalNode = FindEqualNode( quad._edgeNodes, _edgeIntPnts[ iP ], tol2 ); + if ( equalNode ) { + equalNode->Add( _edgeIntPnts[ iP ] ); + } + else { + quad._edgeNodes.push_back( _Node( 0, _edgeIntPnts[ iP ])); + ++_nbIntNodes; + } + } + } + break; + } + case 3: // at a corner + { + _Node& node = _hexNodes[ subEntity - SMESH_Block::ID_FirstV ]; + if ( node.Node() > 0 ) + { + if ( node._intPoint ) + node._intPoint->Add( _edgeIntPnts[ iP ]->_faceIDs, _edgeIntPnts[ iP ]->_node ); + } + else + { + for ( int iF = 0; iF < 3; ++iF ) + { + _Face& quad = _hexQuads[ facets[iF] - SMESH_Block::ID_FirstF ]; + equalNode = FindEqualNode( quad._edgeNodes, _edgeIntPnts[ iP ], tol2 ); + if ( equalNode ) { + equalNode->Add( _edgeIntPnts[ iP ] ); + } + else { + quad._edgeNodes.push_back( _Node( 0, _edgeIntPnts[ iP ])); + ++_nbIntNodes; + } + } + } + break; + } + default: // inside a hex + { + equalNode = FindEqualNode( _vertexNodes, _edgeIntPnts[ iP ], tol2 ); + if ( equalNode ) { + equalNode->Add( _edgeIntPnts[ iP ] ); + } + else { + _vertexNodes.push_back( _Node( 0, _edgeIntPnts[iP] )); + ++_nbIntNodes; + } + } + } // switch( nbFacets ) + + } // loop on _edgeIntPnts + + //SMESHUtils::FreeVector( _edgeIntPnts ); } } //================================================================================ @@ -1355,51 +1669,69 @@ namespace return; _polygons.clear(); - - vector polyhedraNodes; - vector quantities; + _polygons.reserve( 10 ); // create polygons from quadrangles and get their nodes - vector<_Node*> nodes; - nodes.reserve( _nbCornerNodes + _nbIntNodes ); + vector<_Node*> chainNodes; + //vector<_Node*> nodes; // nodes of a face + //nodes.reserve( _nbCornerNodes + _nbIntNodes ); _Link polyLink; - polyLink._faces.reserve( 1 ); + + bool hasEdgeIntersections = !_vertexNodes.empty(); for ( int iF = 0; iF < 6; ++iF ) // loop on 6 sides of a hexahedron { - const _Face& quad = _hexQuads[ iF ] ; + _Face& quad = _hexQuads[ iF ] ; + + if ( !quad._edgeNodes.empty() ) + hasEdgeIntersections = true; _polygons.resize( _polygons.size() + 1 ); _Face& polygon = _polygons.back(); - polygon._links.clear(); - polygon._polyLinks.clear(); polygon._polyLinks.reserve( 10 ); + polygon._polyLinks.reserve( 20 ); - // add splits of a link to a polygon and collect info on nodes - //int nbIn = 0, nbOut = 0, nbCorners = 0; - nodes.clear(); + // add splits of links to a polygon and add _polyLinks to make + // polygon's boundary closed for ( int iE = 0; iE < 4; ++iE ) // loop on 4 sides of a quadrangle { int nbSpits = quad._links[ iE ].NbResultLinks(); for ( int iS = 0; iS < nbSpits; ++iS ) { _OrientedLink split = quad._links[ iE ].ResultLink( iS ); - _Node* n = split.FirstNode(); + _Node* n1 = split.FirstNode(); + _Node* n2 = split.LastNode(); + if ( !n1->IsLinked( n2->_intPoint )) + { + if ( findChain( n1, n2, quad, chainNodes )) + { + for ( size_t i = 1; i < chainNodes.size(); ++i ) + { + polyLink._nodes[0] = chainNodes[i-1]; + polyLink._nodes[1] = chainNodes[i]; + polygon._polyLinks.push_back( polyLink ); + polygon._links.push_back( _OrientedLink( &polygon._polyLinks.back() )); + } + continue; + } + } if ( !polygon._links.empty() ) { _Node* nPrev = polygon._links.back().LastNode(); - if ( nPrev != n ) + if ( nPrev != n1 ) { - polyLink._nodes[0] = nPrev; - polyLink._nodes[1] = n; - polygon._polyLinks.push_back( polyLink ); - polygon._links.push_back( _OrientedLink( &polygon._polyLinks.back() )); - nodes.push_back( nPrev ); + findChain( nPrev, n1, quad, chainNodes ); + for ( size_t i = 1; i < chainNodes.size(); ++i ) + { + polyLink._nodes[0] = chainNodes[i-1]; + polyLink._nodes[1] = chainNodes[i]; + polygon._polyLinks.push_back( polyLink ); + polygon._links.push_back( _OrientedLink( &polygon._polyLinks.back() )); + } } } polygon._links.push_back( split ); - nodes.push_back( n ); } } if ( polygon._links.size() > 1 ) @@ -1408,19 +1740,21 @@ namespace _Node* n2 = polygon._links.front().FirstNode(); if ( n1 != n2 ) { - polyLink._nodes[0] = n1; - polyLink._nodes[1] = n2; - polygon._polyLinks.push_back( polyLink ); - polygon._links.push_back( _OrientedLink( &polygon._polyLinks.back() )); - nodes.push_back( n1 ); + findChain( n1, n2, quad, chainNodes ); + for ( size_t i = 1; i < chainNodes.size(); ++i ) + { + polyLink._nodes[0] = chainNodes[i-1]; + polyLink._nodes[1] = chainNodes[i]; + polygon._polyLinks.push_back( polyLink ); + polygon._links.push_back( _OrientedLink( &polygon._polyLinks.back() )); + } } // add polygon to its links for ( size_t iL = 0; iL < polygon._links.size(); ++iL ) + { + polygon._links[ iL ]._link->_faces.reserve( 2 ); polygon._links[ iL ]._link->_faces.push_back( &polygon ); - // store polygon nodes - quantities.push_back( nodes.size() ); - for ( size_t i = 0; i < nodes.size(); ++i ) - polyhedraNodes.push_back( nodes[i]->Node() ); + } } else { @@ -1432,6 +1766,7 @@ namespace // find free links vector< _OrientedLink* > freeLinks; + freeLinks.reserve(20); for ( size_t iP = 0; iP < _polygons.size(); ++iP ) { _Face& polygon = _polygons[ iP ]; @@ -1439,58 +1774,171 @@ namespace if ( polygon._links[ iL ]._link->_faces.size() < 2 ) freeLinks.push_back( & polygon._links[ iL ]); } - // make closed chains of free links int nbFreeLinks = freeLinks.size(); if ( 0 < nbFreeLinks && nbFreeLinks < 3 ) return; + + set usedFaceIDs; + // if ( hasEdgeIntersections ) + // { + // // detect FACEs supporting remaining polygons + // map nbUsesOfFace; + // map::iterator f2nb; + // for ( size_t iL = 0; iL < freeLinks.size(); ++iL ) + // if ( const B_IntersectPoint* ip = freeLinks[ iL ]->FirstNode()->_intPoint ) + // for ( size_t i = 0; i < ip->_faceIDs.size(); ++i ) + // { + // f2nb = nbUsesOfFace.insert( make_pair( ip->_faceIDs[i], 0 )).first; + // f2nb->second++; + // } + // for ( f2nb = nbUsesOfFace.begin(); f2nb != nbUsesOfFace.end(); ++f2nb ) + // if ( f2nb->second >= 3 ) + // usefulFaceIDs.insert( usefulFaceIDs.end(), f2nb->first ); + // } + // make closed chains of free links while ( nbFreeLinks > 0 ) { - nodes.clear(); _polygons.resize( _polygons.size() + 1 ); _Face& polygon = _polygons.back(); - polygon._links.clear(); + polygon._polyLinks.reserve( 20 ); - // get a remaining link to start from _OrientedLink* curLink = 0; - for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL ) - if (( curLink = freeLinks[ iL ] )) - freeLinks[ iL ] = 0; - nodes.push_back( curLink->LastNode() ); - polygon._links.push_back( *curLink ); - - // find all links connected to curLink - _Node* curNode = 0; - do + _Node* curNode; + if ( !hasEdgeIntersections ) { - curNode = curLink->FirstNode(); - curLink = 0; + // get a remaining link to start from for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL ) - if ( freeLinks[ iL ] && freeLinks[ iL ]->LastNode() == curNode ) - { - curLink = freeLinks[ iL ]; + if (( curLink = freeLinks[ iL ] )) freeLinks[ iL ] = 0; - nodes.push_back( curNode ); - polygon._links.push_back( *curLink ); + polygon._links.push_back( *curLink ); + --nbFreeLinks; + do + { + // find all links connected to curLink + curNode = curLink->FirstNode(); + curLink = 0; + for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL ) + if ( freeLinks[ iL ] && freeLinks[ iL ]->LastNode() == curNode ) + { + curLink = freeLinks[ iL ]; + freeLinks[ iL ] = 0; + polygon._links.push_back( *curLink ); + --nbFreeLinks; + } + } while ( curLink ); + } + else // there are intersections with EDGEs + { + TGeomID curFace; + // get a remaining link to start from, one lying on minimal + // nb of FACEs + { + map< vector< TGeomID >, int > facesOfLink; + map< vector< TGeomID >, int >::iterator f2l; + for ( size_t iL = 0; iL < freeLinks.size(); ++iL ) + if ( freeLinks[ iL ] ) + { + f2l = facesOfLink.insert + ( make_pair( freeLinks[ iL ]->GetNotUsedFace( usedFaceIDs ), iL )).first; + if ( f2l->first.size() == 1 ) + break; + } + f2l = facesOfLink.begin(); + if ( f2l->first.empty() ) + return; + curFace = f2l->first[0]; + curLink = freeLinks[ f2l->second ]; + freeLinks[ f2l->second ] = 0; + } + usedFaceIDs.insert( curFace ); + polygon._links.push_back( *curLink ); + --nbFreeLinks; + + // find all links bounding a FACE of curLink + do + { + // go forward from curLink + curNode = curLink->LastNode(); + curLink = 0; + for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL ) + if ( freeLinks[ iL ] && + freeLinks[ iL ]->FirstNode() == curNode && + freeLinks[ iL ]->LastNode()->IsOnFace( curFace )) + { + curLink = freeLinks[ iL ]; + freeLinks[ iL ] = 0; + polygon._links.push_back( *curLink ); + --nbFreeLinks; + } + } while ( curLink ); + + std::reverse( polygon._links.begin(), polygon._links.end() ); + + curLink = & polygon._links.back(); + do + { + // go backward from curLink + curNode = curLink->FirstNode(); + curLink = 0; + for ( size_t iL = 0; iL < freeLinks.size() && !curLink; ++iL ) + if ( freeLinks[ iL ] && + freeLinks[ iL ]->LastNode() == curNode && + freeLinks[ iL ]->FirstNode()->IsOnFace( curFace )) + { + curLink = freeLinks[ iL ]; + freeLinks[ iL ] = 0; + polygon._links.push_back( *curLink ); + --nbFreeLinks; + } + } while ( curLink ); + + curNode = polygon._links.back().FirstNode(); + + if ( polygon._links[0].LastNode() != curNode ) + { + if ( !_vertexNodes.empty() ) + { + // add links with _vertexNodes if not already used + for ( size_t iN = 0; iN < _vertexNodes.size(); ++iN ) + if ( _vertexNodes[ iN ].IsOnFace( curFace )) + { + bool used = ( curNode == &_vertexNodes[ iN ] ); + for ( size_t iL = 0; iL < polygon._links.size() && !used; ++iL ) + used = ( &_vertexNodes[ iN ] == polygon._links[ iL ].LastNode() ); + if ( !used ) + { + polyLink._nodes[0] = &_vertexNodes[ iN ]; + polyLink._nodes[1] = curNode; + polygon._polyLinks.push_back( polyLink ); + polygon._links.push_back( _OrientedLink( &polygon._polyLinks.back() )); + freeLinks.push_back( &polygon._links.back() ); + ++nbFreeLinks; + curNode = &_vertexNodes[ iN ]; + } + // TODO: to reorder _vertexNodes within polygon, if there are several ones + } } - } while ( curLink ); + polyLink._nodes[0] = polygon._links[0].LastNode(); + polyLink._nodes[1] = curNode; + polygon._polyLinks.push_back( polyLink ); + polygon._links.push_back( _OrientedLink( &polygon._polyLinks.back() )); + freeLinks.push_back( &polygon._links.back() ); + ++nbFreeLinks; + } - nbFreeLinks -= polygon._links.size(); + } // if there are intersections with EDGEs - if ( curNode != nodes.front() || polygon._links.size() < 3 ) + if ( polygon._links.size() < 3 || + polygon._links[0].LastNode() != polygon._links.back().FirstNode() ) return; // closed polygon not found -> invalid polyhedron - quantities.push_back( nodes.size() ); - for ( size_t i = 0; i < nodes.size(); ++i ) - polyhedraNodes.push_back( nodes[i]->Node() ); - - // add polygon to its links and reverse links - for ( size_t i = 0; i < polygon._links.size(); ++i ) + // add polygon to its links + for ( size_t iL = 0; iL < polygon._links.size(); ++iL ) { - polygon._links[i].Reverse(); - polygon._links[i]._link->_faces.push_back( &polygon ); + polygon._links[ iL ]._link->_faces.reserve( 2 ); + polygon._links[ iL ]._link->_faces.push_back( &polygon ); + polygon._links[ iL ].Reverse(); } - - //const size_t firstPoly = _polygons.size(); - } + } // while ( nbFreeLinks > 0 ) if ( ! checkPolyhedronSize() ) { @@ -1505,20 +1953,32 @@ namespace else if ( nbNodes == 6 && _polygons.size() == 5 ) isClassicElem = addPenta(); else if ( nbNodes == 5 && _polygons.size() == 5 ) isClassicElem = addPyra (); if ( !isClassicElem ) - _volumeDefs.set( polyhedraNodes, quantities ); + { + _volumeDefs._nodes.clear(); + _volumeDefs._quantities.clear(); + + for ( size_t iF = 0; iF < _polygons.size(); ++iF ) + { + const size_t nbLinks = _polygons[ iF ]._links.size(); + _volumeDefs._quantities.push_back( nbLinks ); + for ( size_t iL = 0; iL < nbLinks; ++iL ) + _volumeDefs._nodes.push_back( _polygons[ iF ]._links[ iL ].FirstNode() ); + } + } } //================================================================================ /*! * \brief Create elements in the mesh */ - int Hexahedron::MakeElements(SMESH_MesherHelper& helper) + int Hexahedron::MakeElements(SMESH_MesherHelper& helper, + const map< TGeomID, vector< TGeomID > >& edge2faceIDsMap) { SMESHDS_Mesh* mesh = helper.GetMeshDS(); size_t nbCells[3] = { _grid->_coords[0].size() - 1, _grid->_coords[1].size() - 1, _grid->_coords[2].size() - 1 }; - const size_t nbGridCells = nbCells[0] *nbCells [1] * nbCells[2]; + const size_t nbGridCells = nbCells[0] * nbCells[1] * nbCells[2]; vector< Hexahedron* > intersectedHex( nbGridCells, 0 ); int nbIntHex = 0; @@ -1535,7 +1995,7 @@ namespace for ( ; lineInd.More(); ++lineInd ) { GridLine& line = _grid->_lines[ iDir ][ lineInd.LineIndex() ]; - multiset< IntersectionPoint >::const_iterator ip = line._intPoints.begin(); + multiset< F_IntersectPoint >::const_iterator ip = line._intPoints.begin(); for ( ; ip != line._intPoints.end(); ++ip ) { if ( !ip->_node ) continue; @@ -1567,6 +2027,9 @@ namespace } } + // implement geom edges into the mesh + addEdges( helper, intersectedHex, edge2faceIDsMap ); + // add not split hexadrons to the mesh int nbAdded = 0; vector intHexInd( nbIntHex ); @@ -1636,6 +2099,316 @@ namespace return nbAdded; } + //================================================================================ + /*! + * \brief Implements geom edges into the mesh + */ + void Hexahedron::addEdges(SMESH_MesherHelper& helper, + vector< Hexahedron* >& hexes, + const map< TGeomID, vector< TGeomID > >& edge2faceIDsMap) + { + if ( edge2faceIDsMap.empty() ) return; + + // Prepare planes for intersecting with EDGEs + GridPlanes pln[3]; + { + gp_XYZ origPnt = ( _grid->_coords[0][0] * _grid->_axes[0] + + _grid->_coords[1][0] * _grid->_axes[1] + + _grid->_coords[2][0] * _grid->_axes[2] ); + for ( int iDirZ = 0; iDirZ < 3; ++iDirZ ) // iDirZ gives normal direction to planes + { + GridPlanes& planes = pln[ iDirZ ]; + int iDirX = ( iDirZ + 1 ) % 3; + int iDirY = ( iDirZ + 2 ) % 3; + planes._uNorm = ( _grid->_axes[ iDirY ] ^ _grid->_axes[ iDirZ ] ).Normalized(); + planes._vNorm = ( _grid->_axes[ iDirZ ] ^ _grid->_axes[ iDirX ] ).Normalized(); + planes._zNorm = ( _grid->_axes[ iDirX ] ^ _grid->_axes[ iDirY ] ).Normalized(); + double uvDot = planes._uNorm * planes._vNorm; + planes._factor = sqrt( 1. - uvDot * uvDot ); + planes._origins.resize( _grid->_coords[ iDirZ ].size() ); + planes._zProjs.resize ( _grid->_coords[ iDirZ ].size() ); + planes._origins[0] = origPnt; + planes._zProjs [0] = 0; + const double zFactor = _grid->_axes[ iDirZ ] * planes._zNorm; + const vector< double > & u = _grid->_coords[ iDirZ ]; + for ( int i = 1; i < planes._origins.size(); ++i ) + { + planes._origins[i] = origPnt + _grid->_axes[ iDirZ ] * ( u[i] - u[0] ); + planes._zProjs [i] = zFactor * ( u[i] - u[0] ); + } + } + } + const double deflection = _grid->_minCellSize / 20.; + // int facets[6] = { SMESH_Block::ID_F0yz, SMESH_Block::ID_F1yz, + // SMESH_Block::ID_Fx0z, SMESH_Block::ID_Fx1z, + // SMESH_Block::ID_Fxy0, SMESH_Block::ID_Fxy1 }; + E_IntersectPoint ip; + //ip._faceIDs.reserve(2); + + // Intersect EDGEs with the planes + map< TGeomID, vector< TGeomID > >::const_iterator e2fIt = edge2faceIDsMap.begin(); + for ( ; e2fIt != edge2faceIDsMap.end(); ++e2fIt ) + { + const TGeomID edgeID = e2fIt->first; + const TopoDS_Edge & E = TopoDS::Edge( _grid->_shapes( edgeID )); + BRepAdaptor_Curve curve( E ); + + ip._faceIDs = e2fIt->second; + ip._shapeID = edgeID; + + // discretize the EGDE + GCPnts_UniformDeflection discret( curve, deflection, true ); + if ( !discret.IsDone() || discret.NbPoints() < 2 ) + continue; + + // perform intersection + for ( int iDirZ = 0; iDirZ < 3; ++iDirZ ) + { + GridPlanes& planes = pln[ iDirZ ]; + int iDirX = ( iDirZ + 1 ) % 3; + int iDirY = ( iDirZ + 2 ) % 3; + double xLen = _grid->_coords[ iDirX ].back() - _grid->_coords[ iDirX ][0]; + double yLen = _grid->_coords[ iDirY ].back() - _grid->_coords[ iDirY ][0]; + double zFactor = _grid->_axes[ iDirZ ] * planes._zNorm; + //int* zFacets = facets + iDirZ * 2; + + // locate the 1st point of a segment within the grid + gp_XYZ p1 = discret.Value( 1 ).XYZ(); + double u1 = discret.Parameter( 1 ); + double zProj1 = planes._zNorm * ( p1 - planes._origins[0] ); + gp_Pnt orig = planes._origins[0] + planes._zNorm * zProj1; + gp_XY uv = planes.GetUV( p1, orig ); + int iX1 = int( uv.X() / xLen * ( _grid->_coords[ iDirX ].size() - 1. )); + int iY1 = int( uv.Y() / yLen * ( _grid->_coords[ iDirY ].size() - 1. )); + int iZ1 = int( zProj1 / planes._zProjs.back() * ( planes._zProjs.size() - 1. )); + locateValue( iX1, uv.X(), _grid->_coords[ iDirX ] ); + locateValue( iY1, uv.Y(), _grid->_coords[ iDirY ] ); + locateValue( iZ1, zProj1, planes._zProjs ); + + int ijk[3]; // grid index where a segment intersect a plane + ijk[ iDirX ] = iX1; + ijk[ iDirY ] = iY1; + ijk[ iDirZ ] = iZ1; + ip._uvw[ iDirX ] = uv.X() + _grid->_coords[ iDirX ][0]; + ip._uvw[ iDirY ] = uv.Y() + _grid->_coords[ iDirY ][0]; + ip._uvw[ iDirZ ] = zProj1 / zFactor + _grid->_coords[ iDirZ ][0]; + + // add the 1st vertex point to a hexahedron + size_t hexIndex = _grid->CellIndex( ijk[0], ijk[1], ijk[2] ); + if ( iDirZ == 0 && hexIndex < hexes.size() && hexes[ hexIndex ] ) + { + //ip._shapeID = _grid->_shapes.Add( helper.IthVertex( 0, curve.Edge(),/*CumOri=*/false)); + ip._point = p1; + _grid->_edgeIntP.push_back( ip ); + hexes[ hexIndex ]->addIntersection( _grid->_edgeIntP.back() ); + } + + for ( int iP = 2; iP <= discret.NbPoints(); ++iP ) + { + // locate the 2nd point of a segment within the grid + gp_XYZ p2 = discret.Value( iP ).XYZ(); + double u2 = discret.Parameter( iP ); + double zProj2 = planes._zNorm * ( p2 - planes._origins[0] ); + int iZ2 = iZ1; + locateValue( iZ2, zProj2, planes._zProjs ); + + // treat intersections with planes between 2 points of a segment + int dZ = ( iZ1 <= iZ2 ) ? +1 : -1; + int iZ = iZ1 + ( iZ1 < iZ2 ); + for ( int i = 0, nb = Abs( iZ1 - iZ2 ); i < nb; ++i, iZ += dZ ) + { + double r = (( planes._zProjs[ iZ ] - zProj1 ) / ( zProj2 - zProj1 )); + ip._point = curve.Value( u1 * ( 1 - r ) + u2 * r ); + gp_XY uv = planes.GetUV( ip._point, planes._origins[ iZ ]); + locateValue( ijk[ iDirX ], uv.X(), _grid->_coords[ iDirX ] ); + locateValue( ijk[ iDirY ], uv.Y(), _grid->_coords[ iDirY ] ); + ijk[ iDirZ ] = iZ; + ip._uvw[ iDirX ] = uv.X() + _grid->_coords[ iDirX ][0]; + ip._uvw[ iDirY ] = uv.Y() + _grid->_coords[ iDirY ][0]; + ip._uvw[ iDirZ ] = planes._zProjs[ iZ ] / zFactor + _grid->_coords[ iDirZ ][0]; + + _grid->_edgeIntP.push_back( ip ); + + // add ip to hex "above" the plane + hexIndex = _grid->CellIndex( ijk[0], ijk[1], ijk[2] ); + if ( hexIndex < hexes.size() && hexes[ hexIndex ] ) + hexes[ hexIndex ]->addIntersection( _grid->_edgeIntP.back() ); + + // add ip to hex "below" the plane + ijk[ iDirZ ] = iZ-1; + hexIndex = _grid->CellIndex( ijk[0], ijk[1], ijk[2] ); + if ( hexIndex < hexes.size() && hexes[ hexIndex ] ) + hexes[ hexIndex ]->addIntersection( _grid->_edgeIntP.back() ); + } + iZ1 = iZ2; + p1 = p2; + u1 = u2; + zProj1 = zProj2; + } + // add the 2nd vertex point to a hexahedron + if ( iDirZ == 0 ) + { + orig = planes._origins[0] + planes._zNorm * zProj1; + uv = planes.GetUV( p1, orig ); + locateValue( ijk[ iDirX ], uv.X(), _grid->_coords[ iDirX ] ); + locateValue( ijk[ iDirY ], uv.Y(), _grid->_coords[ iDirY ] ); + ijk[ iDirZ ] = iZ1; + hexIndex = _grid->CellIndex( ijk[0], ijk[1], ijk[2] ); + if ( hexIndex < hexes.size() && hexes[ hexIndex ] ) + { + ip._uvw[ iDirX ] = uv.X() + _grid->_coords[ iDirX ][0]; + ip._uvw[ iDirY ] = uv.Y() + _grid->_coords[ iDirY ][0]; + ip._uvw[ iDirZ ] = zProj1 / zFactor + _grid->_coords[ iDirZ ][0]; + ip._point = p1; + _grid->_edgeIntP.push_back( ip ); + hexes[ hexIndex ]->addIntersection( _grid->_edgeIntP.back() ); + } + } + } // loop on 3 grid directions + } // loop on EDGEs + + // Create nodes at found intersections + // const E_IntersectPoint* eip; + // for ( size_t i = 0; i < hexes.size(); ++i ) + // { + // Hexahedron* h = hexes[i]; + // if ( !h ) continue; + // for ( int iF = 0; iF < 6; ++iF ) + // { + // _Face& quad = h->_hexQuads[ iF ]; + // for ( size_t iP = 0; iP < quad._edgeNodes.size(); ++iP ) + // if ( !quad._edgeNodes[ iP ]._node ) + // if (( eip = quad._edgeNodes[ iP ].EdgeIntPnt() )) + // quad._edgeNodes[ iP ]._intPoint->_node = helper.AddNode( eip->_point.X(), + // eip->_point.Y(), + // eip->_point.Z() ); + // } + // for ( size_t iP = 0; iP < hexes[i]->_vertexNodes.size(); ++iP ) + // if (( eip = h->_vertexNodes[ iP ].EdgeIntPnt() )) + // h->_vertexNodes[ iP ]._intPoint->_node = helper.AddNode( eip->_point.X(), + // eip->_point.Y(), + // eip->_point.Z() ); + // } + } + + //================================================================================ + /*! + * \brief Returns index of a hexahedron sub-entities holding a point + * \param [in] ip - intersection point + * \param [out] facets - 0-3 facets holding a point + * \param [out] sub - index of a vertex or an edge holding a point + * \return int - number of facets holding a point + */ + int Hexahedron::getEntity( const E_IntersectPoint* ip, int* facets, int& sub ) + { + enum { X = 4, Y = 2, Z = 1 }; // == 100, 010, 001 + int nbFacets = 0; + int vertex = 0, egdeMask = 0; + + if ( Abs( _grid->_coords[0][ _i ] - ip->_uvw[0] ) < _grid->_tol ) { + facets[ nbFacets++ ] = SMESH_Block::ID_F0yz; + egdeMask |= X; + } + if ( Abs( _grid->_coords[0][ _i+1 ] - ip->_uvw[0] ) < _grid->_tol ) { + facets[ nbFacets++ ] = SMESH_Block::ID_F1yz; + vertex |= X; + egdeMask |= X; + } + if ( Abs( _grid->_coords[1][ _j ] - ip->_uvw[1] ) < _grid->_tol ) { + facets[ nbFacets++ ] = SMESH_Block::ID_Fx0z; + egdeMask |= Y; + } + if ( Abs( _grid->_coords[1][ _j+1 ] - ip->_uvw[1] ) < _grid->_tol ) { + facets[ nbFacets++ ] = SMESH_Block::ID_Fx1z; + vertex |= Y; + egdeMask |= Y; + } + if ( Abs( _grid->_coords[2][ _k ] - ip->_uvw[2] ) < _grid->_tol ) { + facets[ nbFacets++ ] = SMESH_Block::ID_Fxy0; + egdeMask |= Z; + } + if ( Abs( _grid->_coords[2][ _k+1 ] - ip->_uvw[2] ) < _grid->_tol ) { + facets[ nbFacets++ ] = SMESH_Block::ID_Fxy1; + vertex |= Z; + egdeMask |= Z; + } + if ( nbFacets == 3 ) + { + sub = vertex + SMESH_Block::ID_FirstV; + } + else if ( nbFacets == 2 ) + { + const int edge [3][8] = { + { SMESH_Block::ID_E00z, 0, SMESH_Block::ID_E01z, 0, + SMESH_Block::ID_E10z, 0, SMESH_Block::ID_E11z }, + { SMESH_Block::ID_E0y0, SMESH_Block::ID_E0y1, 0, 0, + SMESH_Block::ID_E1y0, SMESH_Block::ID_E1y1, 0, 0 }, + { SMESH_Block::ID_Ex00, SMESH_Block::ID_Ex01, + SMESH_Block::ID_Ex10, SMESH_Block::ID_Ex11, 0, 0, 0, 0 } + }; + switch ( egdeMask ) { + case X | Y: sub = edge[ 0 ][ vertex ]; break; + case X | Z: sub = edge[ 1 ][ vertex ]; break; + default: sub = edge[ 2 ][ vertex ]; + } + } + else if ( nbFacets == 1 ) + { + sub = facets[0]; + } + else // nbFacets == 0 + { + } + return nbFacets; + } + //================================================================================ + /*! + * \brief Adds intersection with an EDGE, either to a _Face or inside a hex + */ + void Hexahedron::addIntersection( const E_IntersectPoint& ip ) + { + // check if ip is really inside the hex +#ifdef _DEBUG_ + if (( _grid->_coords[0][ _i ] - _grid->_tol > ip._uvw[0] ) || + ( _grid->_coords[0][ _i+1 ] + _grid->_tol < ip._uvw[0] ) || + ( _grid->_coords[1][ _j ] - _grid->_tol > ip._uvw[1] ) || + ( _grid->_coords[1][ _j+1 ] + _grid->_tol < ip._uvw[1] ) || + ( _grid->_coords[2][ _k ] - _grid->_tol > ip._uvw[2] ) || + ( _grid->_coords[2][ _k+1 ] + _grid->_tol < ip._uvw[2] )) + throw SALOME_Exception("ip outside a hex"); +#endif + _edgeIntPnts.push_back( & ip ); + } + //================================================================================ + /*! + * \brief Finds nodes at a path from one node to another via intersections with EDGEs + */ + bool Hexahedron::findChain( _Node* n1, + _Node* n2, + _Face& quad, + vector<_Node*>& chn ) + { + chn.clear(); + chn.push_back( n1 ); + bool found = false; + do + { + found = false; + for ( size_t iP = 0; iP < quad._edgeNodes.size(); ++iP ) + if (( std::find( ++chn.begin(), chn.end(), & quad._edgeNodes[iP]) == chn.end()) && + chn.back()->IsLinked( quad._edgeNodes[ iP ]._intPoint )) + { + chn.push_back( & quad._edgeNodes[ iP ]); + found = true; + break; + } + } while ( found && chn.back() != n2 ); + + if ( chn.back() != n2 ) + chn.push_back( n2 ); + + return chn.size() > 2; + } //================================================================================ /*! * \brief Adds computed elements to the mesh @@ -1646,8 +2419,19 @@ namespace // add elements resulted from hexahedron intersection //for ( size_t i = 0; i < _volumeDefs.size(); ++i ) { - vector< const SMDS_MeshNode* >& nodes = _volumeDefs._nodes; - + vector< const SMDS_MeshNode* > nodes( _volumeDefs._nodes.size() ); + for ( size_t iN = 0; iN < nodes.size(); ++iN ) + if ( !( nodes[iN] = _volumeDefs._nodes[iN]->Node() )) + { + if ( const E_IntersectPoint* eip = _volumeDefs._nodes[iN]->EdgeIntPnt() ) + nodes[iN] = _volumeDefs._nodes[iN]->_intPoint->_node = + helper.AddNode( eip->_point.X(), + eip->_point.Y(), + eip->_point.Z() ); + else + throw SALOME_Exception("Bug: no node at intersection point"); + } + if ( !_volumeDefs._quantities.empty() ) { helper.AddPolyhedralVolume( nodes, _volumeDefs._quantities ); @@ -1680,7 +2464,7 @@ namespace bool Hexahedron::isInHole() const { const int ijk[3] = { _i, _j, _k }; - IntersectionPoint curIntPnt; + F_IntersectPoint curIntPnt; // consider a cell to be in a hole if all links in any direction // comes OUT of geometry @@ -1698,19 +2482,19 @@ namespace { const _Link& link = _hexLinks[ iL + 4*iDir ]; // check transition of the first node of a link - const IntersectionPoint* firstIntPnt = 0; + const F_IntersectPoint* firstIntPnt = 0; if ( link._nodes[0]->Node() ) // 1st node is a hexa corner { curIntPnt._paramOnLine = coords[ ijk[ iDir ]] - coords[0]; const GridLine& line = _grid->_lines[ iDir ][ lineIndex[ iL ]]; - multiset< IntersectionPoint >::const_iterator ip = + multiset< F_IntersectPoint >::const_iterator ip = line._intPoints.upper_bound( curIntPnt ); --ip; firstIntPnt = &(*ip); } else if ( !link._intNodes.empty() ) { - firstIntPnt = link._intNodes[0]._intPoint; + firstIntPnt = link._intNodes[0].FaceIntPnt(); } if ( firstIntPnt ) @@ -1736,10 +2520,10 @@ namespace { const _Face& polygon = _polygons[iP]; gp_XYZ area (0,0,0); - SMESH_TNodeXYZ p1 ( polygon._links[ 0 ].FirstNode()->Node() ); + gp_XYZ p1 = polygon._links[ 0 ].FirstNode()->Point().XYZ(); for ( size_t iL = 0; iL < polygon._links.size(); ++iL ) { - SMESH_TNodeXYZ p2 ( polygon._links[ iL ].LastNode()->Node() ); + gp_XYZ p2 = polygon._links[ iL ].LastNode()->Point().XYZ(); area += p1 ^ p2; p1 = p2; } @@ -1764,12 +2548,12 @@ namespace _polygons[4]._links.size() != 4 || _polygons[5]._links.size() != 4 ) return false; - const SMDS_MeshNode* nodes[8]; + _Node* nodes[8]; int nbN = 0; for ( int iL = 0; iL < 4; ++iL ) { // a base node - nodes[iL] = _polygons[0]._links[iL].FirstNode()->Node(); + nodes[iL] = _polygons[0]._links[iL].FirstNode(); ++nbN; // find a top node above the base node @@ -1781,13 +2565,13 @@ namespace if ( quad->_links[i]._link == link ) { // 1st node of a link opposite to in - nodes[iL+4] = quad->_links[(i+2)%4].FirstNode()->Node(); + nodes[iL+4] = quad->_links[(i+2)%4].FirstNode(); ++nbN; break; } } if ( nbN == 8 ) - _volumeDefs.set( vector< const SMDS_MeshNode* >( nodes, nodes+8 )); + _volumeDefs.set( vector< _Node* >( nodes, nodes+8 )); return nbN == 8; } @@ -1797,10 +2581,10 @@ namespace */ bool Hexahedron::addTetra() { - const SMDS_MeshNode* nodes[4]; - nodes[0] = _polygons[0]._links[0].FirstNode()->Node(); - nodes[1] = _polygons[0]._links[1].FirstNode()->Node(); - nodes[2] = _polygons[0]._links[2].FirstNode()->Node(); + _Node* nodes[4]; + nodes[0] = _polygons[0]._links[0].FirstNode(); + nodes[1] = _polygons[0]._links[1].FirstNode(); + nodes[2] = _polygons[0]._links[2].FirstNode(); _Link* link = _polygons[0]._links[0]._link; ASSERT( link->_faces.size() > 1 ); @@ -1810,8 +2594,8 @@ namespace for ( int i = 0; i < 3; ++i ) if ( tria->_links[i]._link == link ) { - nodes[3] = tria->_links[(i+1)%3].LastNode()->Node(); - _volumeDefs.set( vector< const SMDS_MeshNode* >( nodes, nodes+4 )); + nodes[3] = tria->_links[(i+1)%3].LastNode(); + _volumeDefs.set( vector< _Node* >( nodes, nodes+4 )); return true; } @@ -1831,12 +2615,12 @@ namespace if ( iTri < 0 ) return false; // find nodes - const SMDS_MeshNode* nodes[6]; + _Node* nodes[6]; int nbN = 0; for ( int iL = 0; iL < 3; ++iL ) { // a base node - nodes[iL] = _polygons[ iTri ]._links[iL].FirstNode()->Node(); + nodes[iL] = _polygons[ iTri ]._links[iL].FirstNode(); ++nbN; // find a top node above the base node @@ -1849,13 +2633,13 @@ namespace if ( quad->_links[i]._link == link ) { // 1st node of a link opposite to in - nodes[iL+3] = quad->_links[(i+2)%4].FirstNode()->Node(); + nodes[iL+3] = quad->_links[(i+2)%4].FirstNode(); ++nbN; break; } } if ( nbN == 6 ) - _volumeDefs.set( vector< const SMDS_MeshNode* >( nodes, nodes+6 )); + _volumeDefs.set( vector< _Node* >( nodes, nodes+6 )); return ( nbN == 6 ); } @@ -1873,11 +2657,11 @@ namespace if ( iQuad < 0 ) return false; // find nodes - const SMDS_MeshNode* nodes[5]; - nodes[0] = _polygons[iQuad]._links[0].FirstNode()->Node(); - nodes[1] = _polygons[iQuad]._links[1].FirstNode()->Node(); - nodes[2] = _polygons[iQuad]._links[2].FirstNode()->Node(); - nodes[3] = _polygons[iQuad]._links[3].FirstNode()->Node(); + _Node* nodes[5]; + nodes[0] = _polygons[iQuad]._links[0].FirstNode(); + nodes[1] = _polygons[iQuad]._links[1].FirstNode(); + nodes[2] = _polygons[iQuad]._links[2].FirstNode(); + nodes[3] = _polygons[iQuad]._links[3].FirstNode(); _Link* link = _polygons[iQuad]._links[0]._link; ASSERT( link->_faces.size() > 1 ); @@ -1888,8 +2672,8 @@ namespace for ( int i = 0; i < 3; ++i ) if ( tria->_links[i]._link == link ) { - nodes[4] = tria->_links[(i+1)%3].LastNode()->Node(); - _volumeDefs.set( vector< const SMDS_MeshNode* >( nodes, nodes+5 )); + nodes[4] = tria->_links[(i+1)%3].LastNode(); + _volumeDefs.set( vector< _Node* >( nodes, nodes+5 )); return true; } @@ -1930,25 +2714,37 @@ bool StdMeshers_Cartesian_3D::Compute(SMESH_Mesh & theMesh, { Grid grid; - TopTools_MapOfShape faceMap; - for ( TopExp_Explorer fExp( theShape, TopAbs_FACE ); fExp.More(); fExp.Next() ) - if ( !faceMap.Add( fExp.Current() )) - faceMap.Remove( fExp.Current() ); // remove a face shared by two solids - + vector< TopoDS_Shape > faceVec; + { + TopTools_MapOfShape faceMap; + for ( TopExp_Explorer fExp( theShape, TopAbs_FACE ); fExp.More(); fExp.Next() ) + if ( faceMap.Add( fExp.Current() )) // skip a face shared by two solids + faceVec.push_back( fExp.Current() ); + } Bnd_Box shapeBox; - vector facesItersectors( faceMap.Extent() ); - TopTools_MapIteratorOfMapOfShape faceMppIt( faceMap ); - for ( int i = 0; faceMppIt.More(); faceMppIt.Next(), ++i ) + vector facesItersectors( faceVec.size() ); + map< TGeomID, vector< TGeomID > > edge2faceIDsMap; + TopExp_Explorer eExp; + for ( int i = 0; i < faceVec.size(); ++i ) { - facesItersectors[i]._face = TopoDS::Face( faceMppIt.Key() ); - facesItersectors[i]._grid = &grid; + facesItersectors[i]._face = TopoDS::Face ( faceVec[i] ); + facesItersectors[i]._faceID = grid._shapes.Add( faceVec[i] ); + facesItersectors[i]._grid = &grid; shapeBox.Add( facesItersectors[i].GetFaceBndBox() ); + + if ( _hyp->GetToAddEdges() ) + for ( eExp.Init( faceVec[i], TopAbs_EDGE ); eExp.More(); eExp.Next() ) + { + const TopoDS_Edge& edge = TopoDS::Edge( eExp.Current() ); + if ( !SMESH_Algo::isDegenerated( edge )) + edge2faceIDsMap[ grid._shapes.Add( edge )].push_back( facesItersectors[i]._faceID ); + } } vector xCoords, yCoords, zCoords; _hyp->GetCoordinates( xCoords, yCoords, zCoords, shapeBox ); - grid.SetCoordinates( xCoords, yCoords, zCoords, theShape ); + grid.SetCoordinates( xCoords, yCoords, zCoords, _hyp->GetAxisDirs(), theShape ); // check if the grid encloses the shape if ( !_hyp->IsGridBySpacing(0) || @@ -2015,12 +2811,12 @@ bool StdMeshers_Cartesian_3D::Compute(SMESH_Mesh & theMesh, // create volume elements Hexahedron hex( _hyp->GetSizeThreshold(), &grid ); - int nbAdded = hex.MakeElements( helper ); + int nbAdded = hex.MakeElements( helper, edge2faceIDsMap ); SMESHDS_Mesh* meshDS = theMesh.GetMeshDS(); if ( nbAdded > 0 ) { - // make all SOLIDS computed + // make all SOLIDs computed if ( SMESHDS_SubMesh* sm1 = meshDS->MeshElements( solidExp.Current()) ) { SMDS_ElemIteratorPtr volIt = sm1->GetElements(); @@ -2038,23 +2834,28 @@ bool StdMeshers_Cartesian_3D::Compute(SMESH_Mesh & theMesh, // remove free nodes if ( SMESHDS_SubMesh * smDS = meshDS->MeshElements( helper.GetSubShapeID() )) { - // intersection nodes + TIDSortedNodeSet nodesToRemove; + // get intersection nodes for ( int iDir = 0; iDir < 3; ++iDir ) { vector< GridLine >& lines = grid._lines[ iDir ]; for ( size_t i = 0; i < lines.size(); ++i ) { - multiset< IntersectionPoint >::iterator ip = lines[i]._intPoints.begin(); + multiset< F_IntersectPoint >::iterator ip = lines[i]._intPoints.begin(); for ( ; ip != lines[i]._intPoints.end(); ++ip ) if ( ip->_node && ip->_node->NbInverseElements() == 0 ) - meshDS->RemoveFreeNode( ip->_node, smDS, /*fromGroups=*/false ); + nodesToRemove.insert( nodesToRemove.end(), ip->_node ); } } - // grid nodes + // get grid nodes for ( size_t i = 0; i < grid._nodes.size(); ++i ) - if ( !grid._isBndNode[i] ) // nodes on boundary are already removed - if ( grid._nodes[i] && grid._nodes[i]->NbInverseElements() == 0 ) - meshDS->RemoveFreeNode( grid._nodes[i], smDS, /*fromGroups=*/false ); + if ( grid._nodes[i] && grid._nodes[i]->NbInverseElements() == 0 ) + nodesToRemove.insert( nodesToRemove.end(), grid._nodes[i] ); + + // do remove + TIDSortedNodeSet::iterator n = nodesToRemove.begin(); + for ( ; n != nodesToRemove.end(); ++n ) + meshDS->RemoveFreeNode( *n, smDS, /*fromGroups=*/false ); } return nbAdded; diff --git a/src/StdMeshersGUI/StdMeshersGUI_CartesianParamCreator.cxx b/src/StdMeshersGUI/StdMeshersGUI_CartesianParamCreator.cxx index 937c79066..c2b58de8d 100644 --- a/src/StdMeshersGUI/StdMeshersGUI_CartesianParamCreator.cxx +++ b/src/StdMeshersGUI/StdMeshersGUI_CartesianParamCreator.cxx @@ -42,6 +42,7 @@ #include #include #include +#include #include #include #include @@ -49,14 +50,14 @@ #include #include #include +#include #include #include #include +#include #include #include #include -#include -#include #define SPACING 6 #define MARGIN 11 @@ -162,7 +163,7 @@ namespace StdMeshersGUI myStepSpin->SetValue( myStep = 1. ); // 3) Coodrinates/Spacing group - QFrame* csFrame = new QFrame( this ); + QFrame* csFrame = new QFrame( this ); QVBoxLayout* scLay = new QVBoxLayout( csFrame ); scLay->setMargin( 0 ); scLay->setSpacing( SPACING ); @@ -610,7 +611,12 @@ QFrame* StdMeshersGUI_CartesianParamCreator::buildFrame() argGroupLayout->addWidget( myThreshold, row, 1 ); row++; - // 2) Grid definition + // 2) "Implement edges" + myAddEdges = new QCheckBox( tr("ADD_EDGES"), GroupC1 ); + argGroupLayout->addWidget( myAddEdges, row, 0, 1, 2 ); + row++; + + // 3) Grid definition QTabWidget* tabWdg = new QTabWidget( fr ); myAxisTabs[ 0 ] = new StdMeshersGUI::GridAxisTab( tabWdg, 0 ); myAxisTabs[ 1 ] = new StdMeshersGUI::GridAxisTab( tabWdg, 1 ); @@ -637,6 +643,8 @@ void StdMeshersGUI_CartesianParamCreator::retrieveParams() const else myThreshold->setText( varName ); + myAddEdges->setChecked( h->GetToAddEdges() ); + for ( int ax = 0; ax < 3; ++ax ) { if ( h->IsGridBySpacing( ax )) @@ -653,7 +661,8 @@ void StdMeshersGUI_CartesianParamCreator::retrieveParams() const } } if ( dlg() ) - dlg()->setMinimumSize( dlg()->minimumSizeHint().width(), dlg()->minimumSizeHint().height() ); + dlg()->setMinimumSize( dlg()->minimumSizeHint().width(), + dlg()->minimumSizeHint().height() ); } QString StdMeshersGUI_CartesianParamCreator::storeParams() const @@ -668,6 +677,7 @@ QString StdMeshersGUI_CartesianParamCreator::storeParams() const h->SetVarParameter( myThreshold->text().toLatin1().constData(), "SetSizeThreshold" ); h->SetSizeThreshold( myThreshold->text().toDouble() ); + h->SetToAddEdges( myAddEdges->isChecked() ); for ( int ax = 0; ax < 3; ++ax ) { diff --git a/src/StdMeshersGUI/StdMeshersGUI_CartesianParamCreator.h b/src/StdMeshersGUI/StdMeshersGUI_CartesianParamCreator.h index 5a02d5ca1..1187bb90b 100644 --- a/src/StdMeshersGUI/StdMeshersGUI_CartesianParamCreator.h +++ b/src/StdMeshersGUI/StdMeshersGUI_CartesianParamCreator.h @@ -39,18 +39,19 @@ #include #include -class SMESHGUI_SpinBox; -class QLineEdit; +class QAbstractItemModel; class QButtonGroup; -class QTreeWidgetItem; -class QString; -class QWidget; -class QTreeWidget; +class QCheckBox; +class QLineEdit; class QListWidget; -class QStyleOptionViewItem; -class QModelIndex; -class QAbstractItemModel; class QListWidgetItem; +class QModelIndex; +class QString; +class QStyleOptionViewItem; +class QTreeWidget; +class QTreeWidgetItem; +class QWidget; +class SMESHGUI_SpinBox; namespace StdMeshersGUI { @@ -136,6 +137,7 @@ protected: private: QLineEdit* myName; SMESHGUI_SpinBox* myThreshold; + QCheckBox* myAddEdges; StdMeshersGUI::GridAxisTab* myAxisTabs[3]; }; diff --git a/src/StdMeshersGUI/StdMeshers_msg_en.ts b/src/StdMeshersGUI/StdMeshers_msg_en.ts index c9982405e..a1b6307cc 100644 --- a/src/StdMeshersGUI/StdMeshers_msg_en.ts +++ b/src/StdMeshersGUI/StdMeshers_msg_en.ts @@ -493,6 +493,10 @@ THRESHOLD Threshold + + ADD_EDGES + Implement Edges + AXIS_X Axis X diff --git a/src/StdMeshersGUI/StdMeshers_msg_fr.ts b/src/StdMeshersGUI/StdMeshers_msg_fr.ts index 17fac5a04..6c191f2de 100755 --- a/src/StdMeshersGUI/StdMeshers_msg_fr.ts +++ b/src/StdMeshersGUI/StdMeshers_msg_fr.ts @@ -3,6 +3,31 @@ @default + + SMESH_EDGES_WITH_LAYERS + Edges with layers + + + SMESH_FACES_WITH_LAYERS + Faces with layers +(walls) + + + SMESH_ADAPTIVE1D_TITLE + Hypothesis Construction + + + SMESH_MAX_SIZE + Max size + + + SMESH_MIN_SIZE + Min size + + + SMESH_ADAPTIVE1D_HYPOTHESIS + Adaptive + SMESH_ARITHMETIC_1D_HYPOTHESIS Arithmétique 1D @@ -450,6 +475,10 @@ StdMeshersGUI_CartesianParamCreator + + ADD_EDGES + Implement Edges + THRESHOLD Seuil @@ -486,4 +515,19 @@ Pas + + StdMeshersGUI_StdHypothesisCreator + + TO_IGNORE_EDGES + Edges without layers (inlets and oulets) + + + NOT_TO_IGNORE_EDGES + Edges with layers (walls) + + + TO_IGNORE_EDGES_OR_NOT + Specified edges are + + diff --git a/src/StdMeshers_I/StdMeshers_CartesianParameters3D_i.cxx b/src/StdMeshers_I/StdMeshers_CartesianParameters3D_i.cxx index 5221f2cec..3db734054 100644 --- a/src/StdMeshers_I/StdMeshers_CartesianParameters3D_i.cxx +++ b/src/StdMeshers_I/StdMeshers_CartesianParameters3D_i.cxx @@ -224,6 +224,28 @@ void StdMeshers_CartesianParameters3D_i::GetGridSpacing(SMESH::string_array_out } } +//======================================================================= +//function : SetToAddEdges +//purpose : Enables implementation of geometrical edges into the mesh. +//======================================================================= + +void StdMeshers_CartesianParameters3D_i::SetToAddEdges(CORBA::Boolean toAdd) +{ + GetImpl()->SetToAddEdges( toAdd ); + SMESH::TPythonDump() << _this() << ".SetToAddEdges( " << toAdd << " )"; +} + +//======================================================================= +//function : GetToAddEdges +//purpose : Returns true if implementation of geometrical edges into the +// mesh is enabled +//======================================================================= + +CORBA::Boolean StdMeshers_CartesianParameters3D_i::GetToAddEdges() +{ + return GetImpl()->GetToAddEdges(); +} + //======================================================================= //function : IsGridBySpacing //purpose : Return true if the grid is defined by spacing functions and diff --git a/src/StdMeshers_I/StdMeshers_CartesianParameters3D_i.hxx b/src/StdMeshers_I/StdMeshers_CartesianParameters3D_i.hxx index a1f00f605..abdda05f4 100644 --- a/src/StdMeshers_I/StdMeshers_CartesianParameters3D_i.hxx +++ b/src/StdMeshers_I/StdMeshers_CartesianParameters3D_i.hxx @@ -84,6 +84,14 @@ class STDMESHERS_I_EXPORT StdMeshers_CartesianParameters3D_i: SMESH::double_array_out xInternalPoints, CORBA::Short axis) throw (SALOME::SALOME_Exception); + /*! + * \brief Enables implementation of geometrical edges into the mesh. If this feature + * is disabled, sharp edges of the shape are lost ("smoothed") in the mesh if + * they don't coincide with the grid lines + */ + void SetToAddEdges(CORBA::Boolean toAdd); + CORBA::Boolean GetToAddEdges(); + /*! * \brief Return true if the grid is defined by spacing functions and * not by node coordinates -- 2.39.2