X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ViscousLayers.cxx;h=ee90bca54f66b1e62276af0889e5cefd0c43129e;hp=cd43e8f27a1571fdb2f437f521ec11a2754efe77;hb=c65cc4a05360d20c47bc98bc2a3f6464d4f1a327;hpb=c3335563bccf2506186ad293cad2d76e6cca63ee diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index cd43e8f27..ee90bca54 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2021 CEA/DEN, EDF R&D, OPEN CASCADE // // This library is free software; you can redistribute it and/or // modify it under the terms of the GNU Lesser General Public @@ -23,15 +23,18 @@ #include "StdMeshers_ViscousLayers.hxx" +#include "ObjectPool.hxx" #include "SMDS_EdgePosition.hxx" #include "SMDS_FaceOfNodes.hxx" #include "SMDS_FacePosition.hxx" #include "SMDS_MeshNode.hxx" +#include "SMDS_PolygonalFaceOfNodes.hxx" #include "SMDS_SetIterator.hxx" #include "SMESHDS_Group.hxx" #include "SMESHDS_Hypothesis.hxx" #include "SMESHDS_Mesh.hxx" #include "SMESH_Algo.hxx" +#include "SMESH_Block.hxx" #include "SMESH_ComputeError.hxx" #include "SMESH_ControlsDef.hxx" #include "SMESH_Gen.hxx" @@ -39,18 +42,20 @@ #include "SMESH_HypoFilter.hxx" #include "SMESH_Mesh.hxx" #include "SMESH_MeshAlgos.hxx" +#include "SMESH_MeshEditor.hxx" #include "SMESH_MesherHelper.hxx" #include "SMESH_ProxyMesh.hxx" #include "SMESH_subMesh.hxx" #include "SMESH_subMeshEventListener.hxx" #include "StdMeshers_FaceSide.hxx" +#include "StdMeshers_ProjectionUtils.hxx" +#include "StdMeshers_Quadrangle_2D.hxx" #include "StdMeshers_ViscousLayers2D.hxx" #include #include #include #include -//#include #include #include #include @@ -93,6 +98,7 @@ #include #include #include +#include #ifdef _DEBUG_ #define __myDEBUG @@ -114,7 +120,7 @@ namespace VISCOUS_3D enum UIndex { U_TGT = 1, U_SRC, LEN_TGT }; const double theMinSmoothCosin = 0.1; - const double theSmoothThickToElemSizeRatio = 0.3; + const double theSmoothThickToElemSizeRatio = 0.6; const double theMinSmoothTriaAngle = 30; const double theMinSmoothQuadAngle = 45; @@ -200,8 +206,8 @@ namespace VISCOUS_3D virtual void ProcessEvent(const int event, const int eventType, SMESH_subMesh* subMesh, - SMESH_subMeshEventListenerData* data, - const SMESH_Hypothesis* hyp) + SMESH_subMeshEventListenerData* /*data*/, + const SMESH_Hypothesis* /*hyp*/) { if (( SMESH_subMesh::COMPUTE_EVENT == eventType ) && ( SMESH_subMesh::CHECK_COMPUTE_STATE != event && @@ -369,22 +375,7 @@ namespace VISCOUS_3D double _h2lenRatio; // avgNormProj / (2*avgDist) gp_Pnt2d _uv; // UV used in putOnOffsetSurface() public: - static _Curvature* New( double avgNormProj, double avgDist ) - { - _Curvature* c = 0; - if ( fabs( avgNormProj / avgDist ) > 1./200 ) - { - c = new _Curvature; - c->_r = avgDist * avgDist / avgNormProj; - c->_k = avgDist * avgDist / c->_r / c->_r; - //c->_k = avgNormProj / c->_r; - c->_k *= ( c->_r < 0 ? 1/1.1 : 1.1 ); // not to be too restrictive - c->_h2lenRatio = avgNormProj / ( avgDist + avgDist ); - - c->_uv.SetCoord( 0., 0. ); - } - return c; - } + static _Curvature* New( double avgNormProj, double avgDist ); double lenDelta(double len) const { return _k * ( _r + len ); } double lenDeltaByDist(double dist) const { return dist * _h2lenRatio; } }; @@ -394,6 +385,7 @@ namespace VISCOUS_3D struct _LayerEdge; struct _EdgesOnShape; struct _Smoother1D; + struct _Mapper2D; typedef map< const SMDS_MeshNode*, _LayerEdge*, TIDCompare > TNode2Edge; //-------------------------------------------------------------------------------- @@ -409,7 +401,7 @@ namespace VISCOUS_3D gp_XYZ _normal; // to boundary of solid vector _pos; // points computed during inflation - double _len; // length achived with the last inflation step + double _len; // length achieved with the last inflation step double _maxLen; // maximal possible length double _cosin; // of angle (_normal ^ surface) double _minAngle; // of _simplices @@ -427,30 +419,36 @@ namespace VISCOUS_3D enum EFlags { TO_SMOOTH = 0x0000001, MOVED = 0x0000002, // set by _neibors[i]->SetNewLength() - SMOOTHED = 0x0000004, // set by this->Smooth() + SMOOTHED = 0x0000004, // set by _LayerEdge::Smooth() DIFFICULT = 0x0000008, // near concave VERTEX ON_CONCAVE_FACE = 0x0000010, BLOCKED = 0x0000020, // not to inflate any more INTERSECTED = 0x0000040, // close intersection with a face found NORMAL_UPDATED = 0x0000080, - MARKED = 0x0000100, // local usage - MULTI_NORMAL = 0x0000200, // a normal is invisible by some of surrounding faces - NEAR_BOUNDARY = 0x0000400, // is near FACE boundary forcing smooth - SMOOTHED_C1 = 0x0000800, // is on _eosC1 - DISTORTED = 0x0001000, // was bad before smoothing - RISKY_SWOL = 0x0002000, // SWOL is parallel to a source FACE - SHRUNK = 0x0004000, // target node reached a tgt position while shrink() - UNUSED_FLAG = 0x0100000 + UPD_NORMAL_CONV = 0x0000100, // to update normal on boundary of concave FACE + MARKED = 0x0000200, // local usage + MULTI_NORMAL = 0x0000400, // a normal is invisible by some of surrounding faces + NEAR_BOUNDARY = 0x0000800, // is near FACE boundary forcing smooth + SMOOTHED_C1 = 0x0001000, // is on _eosC1 + DISTORTED = 0x0002000, // was bad before smoothing + RISKY_SWOL = 0x0004000, // SWOL is parallel to a source FACE + SHRUNK = 0x0008000, // target node reached a tgt position while shrink() + UNUSED_FLAG = 0x0100000 // to add user flags after }; bool Is ( int flag ) const { return _flags & flag; } void Set ( int flag ) { _flags |= flag; } void Unset( int flag ) { _flags &= ~flag; } + std::string DumpFlags() const; // debug void SetNewLength( double len, _EdgesOnShape& eos, SMESH_MesherHelper& helper ); bool SetNewLength2d( Handle(Geom_Surface)& surface, const TopoDS_Face& F, _EdgesOnShape& eos, SMESH_MesherHelper& helper ); + bool UpdatePositionOnSWOL( SMDS_MeshNode* n, + double tol, + _EdgesOnShape& eos, + SMESH_MesherHelper& helper ); void SetDataByNeighbors( const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, const _EdgesOnShape& eos, @@ -496,11 +494,14 @@ namespace VISCOUS_3D const gp_XYZ& PrevPos() const { return _pos[ _pos.size() - 2 ]; } gp_XYZ PrevCheckPos( _EdgesOnShape* eos=0 ) const; gp_Ax1 LastSegment(double& segLen, _EdgesOnShape& eos) const; - gp_XY LastUV( const TopoDS_Face& F, _EdgesOnShape& eos ) const; + gp_XY LastUV( const TopoDS_Face& F, _EdgesOnShape& eos, int which=-1 ) const; bool IsOnEdge() const { return _2neibors; } + bool IsOnFace() const { return ( _nodes[0]->GetPosition()->GetDim() == 2 ); } + int BaseShapeDim() const { return _nodes[0]->GetPosition()->GetDim(); } gp_XYZ Copy( _LayerEdge& other, _EdgesOnShape& eos, SMESH_MesherHelper& helper ); - void SetCosin( double cosin ); + double SetCosin( double cosin ); void SetNormal( const gp_XYZ& n ) { _normal = n; } + void SetMaxLen( double l ) { _maxLen = l; } int NbSteps() const { return _pos.size() - 1; } // nb inlation steps bool IsNeiborOnEdge( const _LayerEdge* edge ) const; void SetSmooLen( double len ) { // set _len at which smoothing is needed @@ -575,6 +576,7 @@ namespace VISCOUS_3D gp_XYZ* _plnNorm; _2NearEdges() { _edges[0]=_edges[1]=0; _plnNorm = 0; } + ~_2NearEdges(){ delete _plnNorm; } const SMDS_MeshNode* tgtNode(bool is2nd) { return _edges[is2nd] ? _edges[is2nd]->_nodes.back() : 0; } @@ -615,12 +617,30 @@ namespace VISCOUS_3D _thickness = Max( _thickness, hyp->GetTotalThickness() ); _stretchFactor += hyp->GetStretchFactor(); _method = hyp->GetMethod(); + if ( _groupName.empty() ) + _groupName = hyp->GetGroupName(); } } double GetTotalThickness() const { return _thickness; /*_nbHyps ? _thickness / _nbHyps : 0;*/ } double GetStretchFactor() const { return _nbHyps ? _stretchFactor / _nbHyps : 0; } int GetNumberLayers() const { return _nbLayers; } int GetMethod() const { return _method; } + bool ToCreateGroup() const { return !_groupName.empty(); } + const std::string& GetGroupName() const { return _groupName; } + + double Get1stLayerThickness( double realThickness = 0.) const + { + const double T = ( realThickness > 0 ) ? realThickness : GetTotalThickness(); + const double f = GetStretchFactor(); + const int N = GetNumberLayers(); + const double fPowN = pow( f, N ); + double h0; + if ( fPowN - 1 <= numeric_limits::min() ) + h0 = T / N; + else + h0 = T * ( f - 1 )/( fPowN - 1 ); + return h0; + } bool UseSurfaceNormal() const { return _method == StdMeshers_ViscousLayers::SURF_OFFSET_SMOOTH; } @@ -629,9 +649,19 @@ namespace VISCOUS_3D bool IsOffsetMethod() const { return _method == StdMeshers_ViscousLayers::FACE_OFFSET; } + bool operator==( const AverageHyp& other ) const + { + return ( _nbLayers == other._nbLayers && + _method == other._method && + Equals( GetTotalThickness(), other.GetTotalThickness() ) && + Equals( GetStretchFactor(), other.GetStretchFactor() )); + } + static bool Equals( double v1, double v2 ) { return Abs( v1 - v2 ) < 0.01 * ( v1 + v2 ); } + private: - int _nbLayers, _nbHyps, _method; - double _thickness, _stretchFactor; + int _nbLayers, _nbHyps, _method; + double _thickness, _stretchFactor; + std::string _groupName; }; //-------------------------------------------------------------------------------- @@ -655,14 +685,19 @@ namespace VISCOUS_3D vector< _EdgesOnShape* > _eosConcaVer; // edges at concave VERTEXes of a FACE vector< _EdgesOnShape* > _eosC1; // to smooth together several C1 continues shapes - vector< gp_XYZ > _faceNormals; // if _shape is FACE + typedef std::unordered_map< const SMDS_MeshElement*, gp_XYZ > TFace2NormMap; + TFace2NormMap _faceNormals; // if _shape is FACE vector< _EdgesOnShape* > _faceEOS; // to get _faceNormals of adjacent FACEs Handle(ShapeAnalysis_Surface) _offsetSurf; _LayerEdge* _edgeForOffset; + double _offsetValue; + _Mapper2D* _mapper2D; _SolidData* _data; // parent SOLID + _LayerEdge* operator[](size_t i) const { return (_LayerEdge*) _edges[i]; } + size_t size() const { return _edges.size(); } TopAbs_ShapeEnum ShapeType() const { return _shape.IsNull() ? TopAbs_SHAPE : _shape.ShapeType(); } TopAbs_ShapeEnum SWOLType() const @@ -671,13 +706,17 @@ namespace VISCOUS_3D { return std::find( _eosC1.begin(), _eosC1.end(), other ) != _eosC1.end(); } bool GetNormal( const SMDS_MeshElement* face, gp_Vec& norm ); _SolidData& GetData() const { return *_data; } + char ShapeTypeLetter() const + { switch ( ShapeType() ) { case TopAbs_FACE: return 'F'; case TopAbs_EDGE: return 'E'; + case TopAbs_VERTEX: return 'V'; default: return 'S'; }} - _EdgesOnShape(): _shapeID(-1), _subMesh(0), _toSmooth(false), _edgeSmoother(0) {} + _EdgesOnShape(): _shapeID(-1), _subMesh(0), _toSmooth(false), _edgeSmoother(0), _mapper2D(0) {} + ~_EdgesOnShape(); }; //-------------------------------------------------------------------------------- /*! - * \brief Convex FACE whose radius of curvature is less than the thickness of + * \brief Convex FACE whose radius of curvature is less than the thickness of * layers. It is used to detect distortion of prisms based on a convex * FACE and to update normals to enable further increasing the thickness */ @@ -691,7 +730,14 @@ namespace VISCOUS_3D // map a sub-shape to _SolidData::_edgesOnShape map< TGeomID, _EdgesOnShape* > _subIdToEOS; + bool _isTooCurved; bool _normalsFixed; + bool _normalsFixedOnBorders; // used in putOnOffsetSurface() + + double GetMaxCurvature( _SolidData& data, + _EdgesOnShape& eof, + BRepLProp_SLProps& surfProp, + SMESH_MesherHelper& helper); bool GetCenterOfCurvature( _LayerEdge* ledge, BRepLProp_SLProps& surfProp, @@ -725,6 +771,7 @@ namespace VISCOUS_3D TopTools_MapOfShape _before; // SOLIDs to be computed before _solid TGeomID _index; // SOLID id _MeshOfSolid* _proxyMesh; + bool _done; list< THyp > _hyps; list< TopoDS_Shape > _hypShapes; map< TGeomID, THyp > _face2hyp; // filled if _hyps.size() > 1 @@ -741,7 +788,7 @@ namespace VISCOUS_3D // _LayerEdge's with underlying shapes vector< _EdgesOnShape > _edgesOnShape; - // key: an id of shape (EDGE or VERTEX) shared by a FACE with + // key: an ID of shape (EDGE or VERTEX) shared by a FACE with // layers and a FACE w/o layers // value: the shape (FACE or EDGE) to shrink mesh on. // _LayerEdge's basing on nodes on key shape are inflated along the value shape @@ -750,14 +797,12 @@ namespace VISCOUS_3D // Convex FACEs whose radius of curvature is less than the thickness of layers map< TGeomID, _ConvexFace > _convexFaces; - // shapes (EDGEs and VERTEXes) srink from which is forbidden due to collisions with + // shapes (EDGEs and VERTEXes) shrink from which is forbidden due to collisions with // the adjacent SOLID set< TGeomID > _noShrinkShapes; int _nbShapesToSmooth; - //map< TGeomID,Handle(Geom_Curve)> _edge2curve; - vector< _CollisionEdges > _collisionEdges; set< TGeomID > _concaveFaces; @@ -770,8 +815,8 @@ namespace VISCOUS_3D _SolidData(const TopoDS_Shape& s=TopoDS_Shape(), _MeshOfSolid* m=0) - :_solid(s), _proxyMesh(m), _helper(0) {} - ~_SolidData(); + :_solid(s), _proxyMesh(m), _done(false),_helper(0) {} + ~_SolidData() { delete _helper; _helper = 0; } void SortOnEdge( const TopoDS_Edge& E, vector< _LayerEdge* >& edges); void Sort2NeiborsOnEdge( vector< _LayerEdge* >& edges ); @@ -835,6 +880,8 @@ namespace VISCOUS_3D void Append( const gp_Pnt& center, _LayerEdge* ledge ) { + if ( ledge->Is( _LayerEdge::MULTI_NORMAL )) + return; if ( _curvaCenters.size() > 0 ) _segLength2.push_back( center.SquareDistance( _curvaCenters.back() )); _curvaCenters.push_back( center ); @@ -869,6 +916,8 @@ namespace VISCOUS_3D const gp_XY& uvToFix, const double refSign ); }; + struct PyDump; + struct Periodicity; //-------------------------------------------------------------------------------- /*! * \brief Builder of viscous layers @@ -892,13 +941,15 @@ namespace VISCOUS_3D private: - bool findSolidsWithLayers(); + bool findSolidsWithLayers(const bool checkFaceMesh=true); bool setBefore( _SolidData& solidBefore, _SolidData& solidAfter ); bool findFacesWithLayers(const bool onlyWith=false); + void findPeriodicFaces(); void getIgnoreFaces(const TopoDS_Shape& solid, const StdMeshers_ViscousLayers* hyp, const TopoDS_Shape& hypShape, set& ignoreFaces); + int makeEdgesOnShape(); bool makeLayer(_SolidData& data); void setShapeData( _EdgesOnShape& eos, SMESH_subMesh* sm, _SolidData& data ); bool setEdgeData( _LayerEdge& edge, _EdgesOnShape& eos, @@ -941,12 +992,15 @@ namespace VISCOUS_3D void makeOffsetSurface( _EdgesOnShape& eos, SMESH_MesherHelper& ); void putOnOffsetSurface( _EdgesOnShape& eos, int infStep, vector< _EdgesOnShape* >& eosC1, - int smooStep=0, bool moveAll=false ); + int smooStep=0, int moveAll=false ); void findCollisionEdges( _SolidData& data, SMESH_MesherHelper& helper ); + void findEdgesToUpdateNormalNearConvexFace( _ConvexFace & convFace, + _SolidData& data, + SMESH_MesherHelper& helper ); void limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelper& helper ); void limitMaxLenByCurvature( _LayerEdge* e1, _LayerEdge* e2, _EdgesOnShape& eos1, _EdgesOnShape& eos2, - SMESH_MesherHelper& helper ); + const bool isSmoothable ); bool updateNormals( _SolidData& data, SMESH_MesherHelper& helper, int stepNb, double stepSize ); bool updateNormalsOfConvexFaces( _SolidData& data, SMESH_MesherHelper& helper, @@ -978,14 +1032,16 @@ namespace VISCOUS_3D // debug void makeGroupOfLE(); - SMESH_Mesh* _mesh; - SMESH_ComputeErrorPtr _error; + SMESH_Mesh* _mesh; + SMESH_ComputeErrorPtr _error; - vector< _SolidData > _sdVec; - TopTools_IndexedMapOfShape _solids; // to find _SolidData by a solid - TopTools_MapOfShape _shrinkedFaces; + vector< _SolidData > _sdVec; + TopTools_IndexedMapOfShape _solids; // to find _SolidData by a solid + TopTools_MapOfShape _shrunkFaces; + std::unique_ptr _periodicity; - int _tmpFaceID; + int _tmpFaceID; + PyDump* _pyDump; }; //-------------------------------------------------------------------------------- /*! @@ -1033,6 +1089,7 @@ namespace VISCOUS_3D size_t _iSeg[2]; // index of segment where extreme tgt node is projected _EdgesOnShape& _eos; double _curveLen; // length of the EDGE + std::pair _eToSmooth[2]; // indices of _LayerEdge's in _eos static Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge& E, _EdgesOnShape& eos, @@ -1046,61 +1103,71 @@ namespace VISCOUS_3D bool Perform(_SolidData& data, Handle(ShapeAnalysis_Surface)& surface, const TopoDS_Face& F, - SMESH_MesherHelper& helper ) - { - if ( _leParams.empty() || ( !isAnalytic() && _offPoints.empty() )) - prepare( data ); + SMESH_MesherHelper& helper ); - if ( isAnalytic() ) - return smoothAnalyticEdge( data, surface, F, helper ); - else - return smoothComplexEdge ( data, surface, F, helper ); - } void prepare(_SolidData& data ); + void findEdgesToSmooth(); + + bool isToSmooth( int iE ); + bool smoothAnalyticEdge( _SolidData& data, Handle(ShapeAnalysis_Surface)& surface, const TopoDS_Face& F, SMESH_MesherHelper& helper); - - bool smoothComplexEdge( _SolidData& data, + bool smoothComplexEdge( _SolidData& data, Handle(ShapeAnalysis_Surface)& surface, const TopoDS_Face& F, SMESH_MesherHelper& helper); - gp_XYZ getNormalNormal( const gp_XYZ & normal, const gp_XYZ& edgeDir); - _LayerEdge* getLEdgeOnV( bool is2nd ) { return _eos._edges[ is2nd ? _eos._edges.size()-1 : 0 ]->_2neibors->_edges[ is2nd ]; } bool isAnalytic() const { return !_anaCurve.IsNull(); } + + void offPointsToPython() const; // debug + }; + + //-------------------------------------------------------------------------------- + /*! + * \brief Compute positions of nodes of 2D structured mesh using TFI + */ + class _Mapper2D + { + FaceQuadStruct _quadPoints; + + UVPtStruct& uvPnt( size_t i, size_t j ) { return _quadPoints.UVPt( i, j ); } + + public: + _Mapper2D( const TParam2ColumnMap & param2ColumnMap, const TNode2Edge& n2eMap ); + bool ComputeNodePositions(); }; + //-------------------------------------------------------------------------------- /*! * \brief Class of temporary mesh face. * We can't use SMDS_FaceOfNodes since it's impossible to set it's ID which is - * needed because SMESH_ElementSearcher internaly uses set of elements sorted by ID + * needed because SMESH_ElementSearcher internally uses set of elements sorted by ID */ - struct _TmpMeshFace : public SMDS_MeshElement + struct _TmpMeshFace : public SMDS_PolygonalFaceOfNodes { - vector _nn; + const SMDS_MeshElement* _srcFace; + _TmpMeshFace( const vector& nodes, - int id, int faceID=-1, int idInFace=-1): - SMDS_MeshElement(id), _nn(nodes) { setShapeId(faceID); setIdInShape(idInFace); } - virtual const SMDS_MeshNode* GetNode(const int ind) const { return _nn[ind]; } - virtual SMDSAbs_ElementType GetType() const { return SMDSAbs_Face; } - virtual vtkIdType GetVtkType() const { return -1; } - virtual SMDSAbs_EntityType GetEntityType() const { return SMDSEntity_Last; } - virtual SMDSAbs_GeometryType GetGeomType() const - { return _nn.size() == 3 ? SMDSGeom_TRIANGLE : SMDSGeom_QUADRANGLE; } - virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType) const - { return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));} + int ID, + int faceID=-1, + const SMDS_MeshElement* srcFace=0 ): + SMDS_PolygonalFaceOfNodes(nodes), _srcFace( srcFace ) { setID( ID ); setShapeID( faceID ); } + virtual SMDSAbs_EntityType GetEntityType() const + { return _srcFace ? _srcFace->GetEntityType() : SMDSEntity_Quadrangle; } + virtual SMDSAbs_GeometryType GetGeomType() const + { return _srcFace ? _srcFace->GetGeomType() : SMDSGeom_QUADRANGLE; } }; //-------------------------------------------------------------------------------- /*! - * \brief Class of temporary mesh face storing _LayerEdge it's based on + * \brief Class of temporary mesh quadrangle face storing _LayerEdge it's based on */ struct _TmpMeshFaceOnEdge : public _TmpMeshFace { @@ -1108,17 +1175,21 @@ namespace VISCOUS_3D _TmpMeshFaceOnEdge( _LayerEdge* le1, _LayerEdge* le2, int ID ): _TmpMeshFace( vector(4), ID ), _le1(le1), _le2(le2) { - _nn[0]=_le1->_nodes[0]; - _nn[1]=_le1->_nodes.back(); - _nn[2]=_le2->_nodes.back(); - _nn[3]=_le2->_nodes[0]; + myNodes[0]=_le1->_nodes[0]; + myNodes[1]=_le1->_nodes.back(); + myNodes[2]=_le2->_nodes.back(); + myNodes[3]=_le2->_nodes[0]; + } + const SMDS_MeshNode* n( size_t i ) const + { + return myNodes[ i ]; } gp_XYZ GetDir() const // return average direction of _LayerEdge's, normal to EDGE { - SMESH_TNodeXYZ p0s( _nn[0] ); - SMESH_TNodeXYZ p0t( _nn[1] ); - SMESH_TNodeXYZ p1t( _nn[2] ); - SMESH_TNodeXYZ p1s( _nn[3] ); + SMESH_TNodeXYZ p0s( myNodes[0] ); + SMESH_TNodeXYZ p0t( myNodes[1] ); + SMESH_TNodeXYZ p1t( myNodes[2] ); + SMESH_TNodeXYZ p1s( myNodes[3] ); gp_XYZ v0 = p0t - p0s; gp_XYZ v1 = p1t - p1s; gp_XYZ v01 = p1s - p0s; @@ -1129,10 +1200,10 @@ namespace VISCOUS_3D } gp_XYZ GetDir(_LayerEdge* le1, _LayerEdge* le2) // return average direction of _LayerEdge's { - _nn[0]=le1->_nodes[0]; - _nn[1]=le1->_nodes.back(); - _nn[2]=le2->_nodes.back(); - _nn[3]=le2->_nodes[0]; + myNodes[0]=le1->_nodes[0]; + myNodes[1]=le1->_nodes.back(); + myNodes[2]=le2->_nodes.back(); + myNodes[3]=le2->_nodes[0]; return GetDir(); } }; @@ -1190,6 +1261,27 @@ namespace VISCOUS_3D ( dot * dot ) / l1 / l2 >= ( cos * cos )); } + class _Factory + { + ObjectPool< _LayerEdge > _edgePool; + ObjectPool< _Curvature > _curvaturePool; + ObjectPool< _2NearEdges > _nearEdgesPool; + + static _Factory* & me() + { + static _Factory* theFactory = 0; + return theFactory; + } + public: + + _Factory() { me() = this; } + ~_Factory() { me() = 0; } + + static _LayerEdge* NewLayerEdge() { return me()->_edgePool.getNew(); } + static _Curvature * NewCurvature() { return me()->_curvaturePool.getNew(); } + static _2NearEdges* NewNearEdges() { return me()->_nearEdgesPool.getNew(); } + }; + } // namespace VISCOUS_3D @@ -1197,10 +1289,11 @@ namespace VISCOUS_3D //================================================================================ // StdMeshers_ViscousLayers hypothesis // -StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH_Gen* gen) - :SMESH_Hypothesis(hypId, studyId, gen), +StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, SMESH_Gen* gen) + :SMESH_Hypothesis(hypId, gen), _isToIgnoreShapes(1), _nbLayers(1), _thickness(1), _stretchFactor(1), - _method( SURF_OFFSET_SMOOTH ) + _method( SURF_OFFSET_SMOOTH ), + _groupName("") { _name = StdMeshers_ViscousLayers::GetHypType(); _param_algo_dim = -3; // auxiliary hyp used by 3D algos @@ -1232,6 +1325,15 @@ void StdMeshers_ViscousLayers::SetMethod( ExtrusionMethod method ) if ( _method != method ) _method = method, NotifySubMeshesHypothesisModification(); } // -------------------------------------------------------------------------------- +void StdMeshers_ViscousLayers::SetGroupName(const std::string& name) +{ + if ( _groupName != name ) + { + _groupName = name; + if ( !_groupName.empty() ) + NotifySubMeshesHypothesisModification(); + } +} // -------------------------------------------------------------------------------- SMESH_ProxyMesh::Ptr StdMeshers_ViscousLayers::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theShape, @@ -1286,6 +1388,9 @@ std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save) save << " " << _shapeIds[i]; save << " " << !_isToIgnoreShapes; // negate to keep the behavior in old studies. save << " " << _method; + save << " " << _groupName.size(); + if ( !_groupName.empty() ) + save << " " << _groupName; return save; } // -------------------------------------------------------------------------------- std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load) @@ -1298,14 +1403,21 @@ std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load) _isToIgnoreShapes = !shapeToTreat; if ( load >> method ) _method = (ExtrusionMethod) method; + int nameSize = 0; + if ( load >> nameSize && nameSize > 0 ) + { + _groupName.resize( nameSize ); + load.get( _groupName[0] ); // remove a white-space + load.getline( &_groupName[0], nameSize + 1 ); + } } else { _isToIgnoreShapes = true; // old behavior } return load; } // -------------------------------------------------------------------------------- -bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh* theMesh, - const TopoDS_Shape& theShape) +bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh* /*theMesh*/, + const TopoDS_Shape& /*theShape*/) { // TODO return false; @@ -1331,22 +1443,77 @@ bool StdMeshers_ViscousLayers::IsShapeWithLayers(int shapeIndex) const ( std::find( _shapeIds.begin(), _shapeIds.end(), shapeIndex ) != _shapeIds.end() ); return IsToIgnoreShapes() ? !isIn : isIn; } + +// -------------------------------------------------------------------------------- +SMDS_MeshGroup* StdMeshers_ViscousLayers::CreateGroup( const std::string& theName, + SMESH_Mesh& theMesh, + SMDSAbs_ElementType theType) +{ + SMESH_Group* group = 0; + SMDS_MeshGroup* groupDS = 0; + + if ( theName.empty() ) + return groupDS; + + if ( SMESH_Mesh::GroupIteratorPtr grIt = theMesh.GetGroups() ) + while( grIt->more() && !group ) + { + group = grIt->next(); + if ( !group || + group->GetGroupDS()->GetType() != theType || + group->GetName() != theName || + !dynamic_cast< SMESHDS_Group* >( group->GetGroupDS() )) + group = 0; + } + if ( !group ) + group = theMesh.AddGroup( theType, theName.c_str() ); + + groupDS = & dynamic_cast< SMESHDS_Group* >( group->GetGroupDS() )->SMDSGroup(); + + return groupDS; +} + // END StdMeshers_ViscousLayers hypothesis //================================================================================ namespace VISCOUS_3D { - gp_XYZ getEdgeDir( const TopoDS_Edge& E, const TopoDS_Vertex& fromV ) + gp_XYZ getEdgeDir( const TopoDS_Edge& E, const TopoDS_Vertex& fromV, + const double h0, bool* isRegularEdge = nullptr ) { gp_Vec dir; double f,l; Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l ); if ( c.IsNull() ) return gp_XYZ( Precision::Infinite(), 1e100, 1e100 ); - gp_Pnt p = BRep_Tool::Pnt( fromV ); - double distF = p.SquareDistance( c->Value( f )); - double distL = p.SquareDistance( c->Value( l )); + gp_Pnt p = BRep_Tool::Pnt( fromV ); + gp_Pnt pf = c->Value( f ), pl = c->Value( l ); + double distF = p.SquareDistance( pf ); + double distL = p.SquareDistance( pl ); c->D1(( distF < distL ? f : l), p, dir ); if ( distL < distF ) dir.Reverse(); + bool isDifficult = false; + if ( dir.SquareMagnitude() < h0 * h0 ) // check dir orientation + { + gp_Pnt& pClose = distF < distL ? pf : pl; + gp_Pnt& pFar = distF < distL ? pl : pf; + gp_Pnt pMid = 0.9 * pClose.XYZ() + 0.1 * pFar.XYZ(); + gp_Vec vMid( p, pMid ); + double dot = vMid * dir; + double cos2 = dot * dot / dir.SquareMagnitude() / vMid.SquareMagnitude(); + if ( cos2 < 0.7 * 0.7 || dot < 0 ) // large angle between dir and vMid + { + double uClose = distF < distL ? f : l; + double uFar = distF < distL ? l : f; + double r = h0 / SMESH_Algo::EdgeLength( E ); + double uMid = ( 1 - r ) * uClose + r * uFar; + pMid = c->Value( uMid ); + dir = gp_Vec( p, pMid ); + isDifficult = true; + } + } + if ( isRegularEdge ) + *isRegularEdge = !isDifficult; + return dir.XYZ(); } //-------------------------------------------------------------------------------- @@ -1363,8 +1530,8 @@ namespace VISCOUS_3D } //-------------------------------------------------------------------------------- gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV, - const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok, - double* cosin=0); + const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok/*, + double* cosin=0*/); //-------------------------------------------------------------------------------- gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Edge& fromE, const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok) @@ -1405,7 +1572,7 @@ namespace VISCOUS_3D //-------------------------------------------------------------------------------- gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV, const SMDS_MeshNode* node, SMESH_MesherHelper& helper, - bool& ok, double* cosin) + bool& ok/*, double* cosin*/) { TopoDS_Face faceFrw = F; faceFrw.Orientation( TopAbs_FORWARD ); @@ -1437,7 +1604,7 @@ namespace VISCOUS_3D ok = true; for ( size_t i = 0; i < nbEdges && ok; ++i ) { - edgeDir[i] = getEdgeDir( edges[i], fromV ); + edgeDir[i] = getEdgeDir( edges[i], fromV, 0.1 * SMESH_Algo::EdgeLength( edges[i] )); double size2 = edgeDir[i].SquareModulus(); if (( ok = size2 > numeric_limits::min() )) edgeDir[i] /= sqrt( size2 ); @@ -1457,15 +1624,15 @@ namespace VISCOUS_3D if ( angle < 0 ) dir.Reverse(); } - if ( cosin ) { - double angle = gp_Vec( edgeDir[0] ).Angle( dir ); - *cosin = Cos( angle ); - } + // if ( cosin ) { + // double angle = faceNormal.Angle( dir ); + // *cosin = Cos( angle ); + // } } else if ( nbEdges == 1 ) { dir = getFaceDir( faceFrw, edges[ edges[0].IsNull() ], node, helper, ok ); - if ( cosin ) *cosin = 1.; + //if ( cosin ) *cosin = 1.; } else { @@ -1569,7 +1736,7 @@ namespace VISCOUS_3D //================================================================================ /*! - * \brief Computes mimimal distance of face in-FACE nodes from an EDGE + * \brief Computes minimal distance of face in-FACE nodes from an EDGE * \param [in] face - the mesh face to treat * \param [in] nodeOnEdge - a node on the EDGE * \param [out] faceSize - the computed distance @@ -1595,8 +1762,8 @@ namespace VISCOUS_3D // look for two neighbor not in-FACE nodes of face for ( int i = 0; i < 2; ++i ) { - if ( nNext[i]->GetPosition()->GetDim() != 2 && - nNext[i]->GetID() < nodeOnEdge->GetID() ) + if (( nNext[i]->GetPosition()->GetDim() != 2 ) && + ( nodeOnEdge->GetPosition()->GetDim() == 0 || nNext[i]->GetID() < nodeOnEdge->GetID() )) { // look for an in-FACE node for ( int iN = 0; iN < nbN; ++iN ) @@ -1659,21 +1826,22 @@ namespace VISCOUS_3D //-------------------------------------------------------------------------------- // DEBUG. Dump intermediate node positions into a python script - // HOWTO use: run python commands written in a console to see - // construction steps of viscous layers + // HOWTO use: run python commands written in a console and defined in /tmp/viscous.py + // to see construction steps of viscous layers #ifdef __myDEBUG - ofstream* py; - int theNbPyFunc; - struct PyDump { + ostream* py; + int theNbPyFunc; + struct PyDump + { PyDump(SMESH_Mesh& m) { int tag = 3 + m.GetId(); const char* fname = "/tmp/viscous.py"; - cout << "execfile('"< ostream & operator<<( const T &anything ) { return *this ; } + }; + void Pause() { py = &_mystream; } + void Resume() { py = _pyStream; } + MyStream _mystream; + ostream* _pyStream; }; #define dumpFunction(f) { _dumpFunction(f, __LINE__);} #define dumpMove(n) { _dumpMove(n, __LINE__);} @@ -1709,7 +1885,7 @@ namespace VISCOUS_3D #else - struct PyDump { PyDump(SMESH_Mesh&) {} void Finish() {} }; + struct PyDump { PyDump(SMESH_Mesh&) {} void Finish() {} void Pause() {} void Resume() {} }; #define dumpFunction(f) f #define dumpMove(n) #define dumpMoveComm(n,txt) @@ -1842,6 +2018,8 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh& theMesh, { _mesh = & theMesh; + _Factory factory; + // check if proxy mesh already computed TopExp_Explorer exp( theShape, TopAbs_SOLID ); if ( !exp.More() ) @@ -1850,29 +2028,42 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh& theMesh, if ( _ViscousListener::GetSolidMesh( _mesh, exp.Current(), /*toCreate=*/false)) return SMESH_ComputeErrorPtr(); // everything already computed - PyDump debugDump( theMesh ); - - // TODO: ignore already computed SOLIDs + // TODO: ignore already computed SOLIDs if ( !findSolidsWithLayers()) return _error; if ( !findFacesWithLayers() ) return _error; + if ( !makeEdgesOnShape() ) + return _error; + + findPeriodicFaces(); + + PyDump debugDump( theMesh ); + _pyDump = &debugDump; + + for ( size_t i = 0; i < _sdVec.size(); ++i ) { size_t iSD = 0; for ( iSD = 0; iSD < _sdVec.size(); ++iSD ) // find next SOLID to compute if ( _sdVec[iSD]._before.IsEmpty() && - _sdVec[iSD]._n2eMap.empty() ) + !_sdVec[iSD]._solid.IsNull() && + !_sdVec[iSD]._done ) break; + if ( iSD == _sdVec.size() ) + break; // all done if ( ! makeLayer(_sdVec[iSD]) ) // create _LayerEdge's return _error; - if ( _sdVec[iSD]._n2eMap.size() == 0 ) + if ( _sdVec[iSD]._n2eMap.size() == 0 ) // no layers in a SOLID + { + _sdVec[iSD]._solid.Nullify(); continue; - + } + if ( ! inflate(_sdVec[iSD]) ) // increase length of _LayerEdge's return _error; @@ -1884,6 +2075,8 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh& theMesh, addBoundaryElements(_sdVec[iSD]); // create quadrangles on prism bare sides + _sdVec[iSD]._done = true; + const TopoDS_Shape& solid = _sdVec[iSD]._solid; for ( iSD = 0; iSD < _sdVec.size(); ++iSD ) _sdVec[iSD]._before.Remove( solid ); @@ -1910,7 +2103,7 @@ SMESH_ComputeErrorPtr _ViscousBuilder::CheckHypotheses( SMESH_Mesh& mesh return SMESH_ComputeErrorPtr(); // everything already computed - findSolidsWithLayers(); + findSolidsWithLayers( /*checkFaceMesh=*/false ); bool ok = findFacesWithLayers( true ); // remove _MeshOfSolid's of _SolidData's @@ -1929,7 +2122,7 @@ SMESH_ComputeErrorPtr _ViscousBuilder::CheckHypotheses( SMESH_Mesh& mesh */ //================================================================================ -bool _ViscousBuilder::findSolidsWithLayers() +bool _ViscousBuilder::findSolidsWithLayers(const bool checkFaceMesh) { // get all solids TopTools_IndexedMapOfShape allSolids; @@ -1939,13 +2132,28 @@ bool _ViscousBuilder::findSolidsWithLayers() SMESH_HypoFilter filter; for ( int i = 1; i <= allSolids.Extent(); ++i ) { - // find StdMeshers_ViscousLayers hyp assigned to the i-th solid SMESH_subMesh* sm = _mesh->GetSubMesh( allSolids(i) ); if ( sm->GetSubMeshDS() && sm->GetSubMeshDS()->NbElements() > 0 ) continue; // solid is already meshed + // TODO: check if algo is hidden SMESH_Algo* algo = sm->GetAlgo(); if ( !algo ) continue; - // TODO: check if algo is hidden + // check if all FACEs are meshed, which can be false if Compute() a sub-shape + if ( checkFaceMesh ) + { + bool facesMeshed = true; + SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(false,true); + while ( smIt->more() && facesMeshed ) + { + SMESH_subMesh * faceSM = smIt->next(); + if ( faceSM->GetSubShape().ShapeType() != TopAbs_FACE ) + break; + facesMeshed = faceSM->IsMeshComputed(); + } + if ( !facesMeshed ) + continue; + } + // find StdMeshers_ViscousLayers hyp assigned to the i-th solid const list & allHyps = algo->GetUsedHypothesis(*_mesh, allSolids(i), /*ignoreAuxiliary=*/false); _SolidData* soData = 0; @@ -2110,7 +2318,7 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith) for ( ; exp.More(); exp.Next() ) { const TopoDS_Face& face = TopoDS::Face( exp.Current() ); - const TGeomID faceID = getMeshDS()->ShapeToIndex( face ); + const TGeomID faceID = getMeshDS()->ShapeToIndex( face ); if ( //!sdVec[i]._ignoreFaceIds.count( faceID ) && helper.NbAncestors( face, *_mesh, TopAbs_SOLID ) > 1 && helper.IsReversedSubMesh( face )) @@ -2120,7 +2328,7 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith) } } - // Find faces to shrink mesh on (solution 2 in issue 0020832); + // Find FACEs to shrink mesh on (solution 2 in issue 0020832): fill in _shrinkShape2Shape TopTools_IndexedMapOfShape shapes; std::string structAlgoName = "Hexa_3D"; for ( size_t i = 0; i < _sdVec.size(); ++i ) @@ -2130,127 +2338,32 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith) for ( int iE = 1; iE <= shapes.Extent(); ++iE ) { const TopoDS_Shape& edge = shapes(iE); - // find 2 faces sharing an edge + // find 2 FACEs sharing an EDGE TopoDS_Shape FF[2]; - PShapeIteratorPtr fIt = helper.GetAncestors(edge, *_mesh, TopAbs_FACE); + PShapeIteratorPtr fIt = helper.GetAncestors(edge, *_mesh, TopAbs_FACE, &_sdVec[i]._solid); while ( fIt->more()) { const TopoDS_Shape* f = fIt->next(); - if ( helper.IsSubShape( *f, _sdVec[i]._solid)) - FF[ int( !FF[0].IsNull()) ] = *f; + FF[ int( !FF[0].IsNull()) ] = *f; } if( FF[1].IsNull() ) continue; // seam edge can be shared by 1 FACE only + // check presence of layers on them int ignore[2]; for ( int j = 0; j < 2; ++j ) ignore[j] = _sdVec[i]._ignoreFaceIds.count( getMeshDS()->ShapeToIndex( FF[j] )); if ( ignore[0] == ignore[1] ) continue; // nothing interesting - TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ]; - // check presence of layers on fWOL within an adjacent SOLID - bool collision = false; - PShapeIteratorPtr sIt = helper.GetAncestors( fWOL, *_mesh, TopAbs_SOLID ); - while ( const TopoDS_Shape* solid = sIt->next() ) - if ( !solid->IsSame( _sdVec[i]._solid )) - { - int iSolid = _solids.FindIndex( *solid ); - int iFace = getMeshDS()->ShapeToIndex( fWOL ); - if ( iSolid > 0 && !_sdVec[ iSolid-1 ]._ignoreFaceIds.count( iFace )) - { - // check if solid's mesh is unstructured and then try to set it - // to be computed after the i-th solid - SMESH_Algo* algo = _mesh->GetSubMesh( *solid )->GetAlgo(); - bool isStructured = ( algo->GetName() == structAlgoName ); - if ( isStructured || !setBefore( _sdVec[ i ], _sdVec[ iSolid-1 ] )) - collision = true; // don't shrink fWOL - break; - } - } - // add edge to maps + TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ]; // FACE w/o layers + + // add EDGE to maps if ( !fWOL.IsNull()) { TGeomID edgeInd = getMeshDS()->ShapeToIndex( edge ); _sdVec[i]._shrinkShape2Shape.insert( make_pair( edgeInd, fWOL )); - if ( collision ) - { - // _shrinkShape2Shape will be used to temporary inflate _LayerEdge's based - // on the edge but shrink won't be performed - _sdVec[i]._noShrinkShapes.insert( edgeInd ); - for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() ) - _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( vIt.Value() )); - } } } } - // Exclude from _shrinkShape2Shape FACE's that can't be shrinked since - // the algo of the SOLID sharing the FACE does not support it - set< string > notSupportAlgos; notSupportAlgos.insert( structAlgoName ); - for ( size_t i = 0; i < _sdVec.size(); ++i ) - { - map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin(); - for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f ) - { - const TopoDS_Shape& fWOL = e2f->second; - const TGeomID edgeID = e2f->first; - bool notShrinkFace = false; - PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID); - while ( soIt->more() ) - { - const TopoDS_Shape* solid = soIt->next(); - if ( _sdVec[i]._solid.IsSame( *solid )) continue; - SMESH_Algo* algo = _mesh->GetSubMesh( *solid )->GetAlgo(); - if ( !algo || !notSupportAlgos.count( algo->GetName() )) continue; - notShrinkFace = true; - size_t iSolid = 0; - for ( ; iSolid < _sdVec.size(); ++iSolid ) - { - if ( _sdVec[iSolid]._solid.IsSame( *solid ) ) { - if ( _sdVec[iSolid]._shrinkShape2Shape.count( edgeID )) - notShrinkFace = false; - break; - } - } - if ( notShrinkFace ) - { - _sdVec[i]._noShrinkShapes.insert( edgeID ); - - // add VERTEXes of the edge in _noShrinkShapes - TopoDS_Shape edge = getMeshDS()->IndexToShape( edgeID ); - for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() ) - _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( vIt.Value() )); - - // check if there is a collision with to-shrink-from EDGEs in iSolid - if ( iSolid == _sdVec.size() ) - continue; // no VL in the solid - shapes.Clear(); - TopExp::MapShapes( fWOL, TopAbs_EDGE, shapes); - for ( int iE = 1; iE <= shapes.Extent(); ++iE ) - { - const TopoDS_Edge& E = TopoDS::Edge( shapes( iE )); - const TGeomID eID = getMeshDS()->ShapeToIndex( E ); - if ( eID == edgeID || - !_sdVec[iSolid]._shrinkShape2Shape.count( eID ) || - _sdVec[i]._noShrinkShapes.count( eID )) - continue; - for ( int is1st = 0; is1st < 2; ++is1st ) - { - TopoDS_Vertex V = helper.IthVertex( is1st, E ); - if ( _sdVec[i]._noShrinkShapes.count( getMeshDS()->ShapeToIndex( V ) )) - { - // _sdVec[i]._noShrinkShapes.insert( eID ); - // V = helper.IthVertex( !is1st, E ); - // _sdVec[i]._noShrinkShapes.insert( getMeshDS()->ShapeToIndex( V )); - //iE = 0; // re-start the loop on EDGEs of fWOL - return error("No way to make a conformal mesh with " - "the given set of faces with layers", _sdVec[i]._index); - } - } - } - } - - } // while ( soIt->more() ) - } // loop on _sdVec[i]._shrinkShape2Shape - } // loop on _sdVec to fill in _SolidData::_noShrinkShapes // Find the SHAPE along which to inflate _LayerEdge based on VERTEX @@ -2264,18 +2377,14 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith) // find faces WOL sharing the vertex vector< TopoDS_Shape > facesWOL; size_t totalNbFaces = 0; - PShapeIteratorPtr fIt = helper.GetAncestors(vertex, *_mesh, TopAbs_FACE); + PShapeIteratorPtr fIt = helper.GetAncestors(vertex, *_mesh, TopAbs_FACE, &_sdVec[i]._solid ); while ( fIt->more()) { const TopoDS_Shape* f = fIt->next(); - if ( helper.IsSubShape( *f, _sdVec[i]._solid ) ) - { - totalNbFaces++; - const int fID = getMeshDS()->ShapeToIndex( *f ); - if ( _sdVec[i]._ignoreFaceIds.count ( fID ) /*&& - !_sdVec[i]._noShrinkShapes.count( fID )*/) - facesWOL.push_back( *f ); - } + totalNbFaces++; + const int fID = getMeshDS()->ShapeToIndex( *f ); + if ( _sdVec[i]._ignoreFaceIds.count ( fID ) /*&& !_sdVec[i]._noShrinkShapes.count( fID )*/) + facesWOL.push_back( *f ); } if ( facesWOL.size() == totalNbFaces || facesWOL.empty() ) continue; // no layers at this vertex or no WOL @@ -2325,7 +2434,154 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith) } } - // add FACEs of other SOLIDs to _ignoreFaceIds + // Add to _noShrinkShapes sub-shapes of FACE's that can't be shrunk since + // the algo of the SOLID sharing the FACE does not support it or for other reasons + set< string > notSupportAlgos; notSupportAlgos.insert( structAlgoName ); + for ( size_t i = 0; i < _sdVec.size(); ++i ) + { + map< TGeomID, TopoDS_Shape >::iterator e2f = _sdVec[i]._shrinkShape2Shape.begin(); + for ( ; e2f != _sdVec[i]._shrinkShape2Shape.end(); ++e2f ) + { + const TopoDS_Shape& fWOL = e2f->second; + const TGeomID edgeID = e2f->first; + TGeomID faceID = getMeshDS()->ShapeToIndex( fWOL ); + TopoDS_Shape edge = getMeshDS()->IndexToShape( edgeID ); + if ( edge.ShapeType() != TopAbs_EDGE ) + continue; // shrink shape is VERTEX + + TopoDS_Shape solid; + PShapeIteratorPtr soIt = helper.GetAncestors(fWOL, *_mesh, TopAbs_SOLID); + while ( soIt->more() && solid.IsNull() ) + { + const TopoDS_Shape* so = soIt->next(); + if ( !so->IsSame( _sdVec[i]._solid )) + solid = *so; + } + if ( solid.IsNull() ) + continue; + + bool noShrinkE = false; + SMESH_Algo* algo = _mesh->GetSubMesh( solid )->GetAlgo(); + bool isStructured = ( algo && algo->GetName() == structAlgoName ); + size_t iSolid = _solids.FindIndex( solid ) - 1; + if ( iSolid < _sdVec.size() && _sdVec[ iSolid ]._ignoreFaceIds.count( faceID )) + { + // the adjacent SOLID has NO layers on fWOL; + // shrink allowed if + // - there are layers on the EDGE in the adjacent SOLID + // - there are NO layers in the adjacent SOLID && algo is unstructured and computed later + bool hasWLAdj = (_sdVec[iSolid]._shrinkShape2Shape.count( edgeID )); + bool shrinkAllowed = (( hasWLAdj ) || + ( !isStructured && setBefore( _sdVec[ i ], _sdVec[ iSolid ] ))); + noShrinkE = !shrinkAllowed; + } + else if ( iSolid < _sdVec.size() ) + { + // the adjacent SOLID has layers on fWOL; + // check if SOLID's mesh is unstructured and then try to set it + // to be computed after the i-th solid + if ( isStructured || !setBefore( _sdVec[ i ], _sdVec[ iSolid ] )) + noShrinkE = true; // don't shrink fWOL + } + else + { + // the adjacent SOLID has NO layers at all + noShrinkE = isStructured; + } + + if ( noShrinkE ) + { + _sdVec[i]._noShrinkShapes.insert( edgeID ); + + // check if there is a collision with to-shrink-from EDGEs in iSolid + // if ( iSolid < _sdVec.size() ) + // { + // shapes.Clear(); + // TopExp::MapShapes( fWOL, TopAbs_EDGE, shapes); + // for ( int iE = 1; iE <= shapes.Extent(); ++iE ) + // { + // const TopoDS_Edge& E = TopoDS::Edge( shapes( iE )); + // const TGeomID eID = getMeshDS()->ShapeToIndex( E ); + // if ( eID == edgeID || + // !_sdVec[iSolid]._shrinkShape2Shape.count( eID ) || + // _sdVec[i]._noShrinkShapes.count( eID )) + // continue; + // for ( int is1st = 0; is1st < 2; ++is1st ) + // { + // TopoDS_Vertex V = helper.IthVertex( is1st, E ); + // if ( _sdVec[i]._noShrinkShapes.count( getMeshDS()->ShapeToIndex( V ) )) + // { + // return error("No way to make a conformal mesh with " + // "the given set of faces with layers", _sdVec[i]._index); + // } + // } + // } + // } + } + + // add VERTEXes of the edge in _noShrinkShapes, which is necessary if + // _shrinkShape2Shape is different in the adjacent SOLID + for ( TopoDS_Iterator vIt( edge ); vIt.More(); vIt.Next() ) + { + TGeomID vID = getMeshDS()->ShapeToIndex( vIt.Value() ); + bool noShrinkV = false, noShrinkIfAdjMeshed = false; + + if ( iSolid < _sdVec.size() ) + { + if ( _sdVec[ iSolid ]._ignoreFaceIds.count( faceID )) + { + map< TGeomID, TopoDS_Shape >::iterator i2S, i2SAdj; + i2S = _sdVec[i ]._shrinkShape2Shape.find( vID ); + i2SAdj = _sdVec[iSolid]._shrinkShape2Shape.find( vID ); + if ( i2SAdj == _sdVec[iSolid]._shrinkShape2Shape.end() ) + noShrinkV = (( isStructured ) || + ( noShrinkIfAdjMeshed = i2S->second.ShapeType() == TopAbs_EDGE )); + else + noShrinkV = ( ! i2S->second.IsSame( i2SAdj->second )); + } + else + { + noShrinkV = noShrinkE; + } + } + else + { + // the adjacent SOLID has NO layers at all + if ( isStructured ) + { + noShrinkV = true; + } + else + { + noShrinkV = noShrinkIfAdjMeshed = + ( _sdVec[i]._shrinkShape2Shape[ vID ].ShapeType() == TopAbs_EDGE ); + } + } + + if ( noShrinkV && noShrinkIfAdjMeshed ) + { + // noShrinkV if FACEs in the adjacent SOLID are meshed + PShapeIteratorPtr fIt = helper.GetAncestors( _sdVec[i]._shrinkShape2Shape[ vID ], + *_mesh, TopAbs_FACE, &solid ); + while ( fIt->more() ) + { + const TopoDS_Shape* f = fIt->next(); + if ( !f->IsSame( fWOL )) + { + noShrinkV = ! _mesh->GetSubMesh( *f )->IsEmpty(); + break; + } + } + } + if ( noShrinkV ) + _sdVec[i]._noShrinkShapes.insert( vID ); + } + + } // loop on _sdVec[i]._shrinkShape2Shape + } // loop on _sdVec to fill in _SolidData::_noShrinkShapes + + + // add FACEs of other SOLIDs to _ignoreFaceIds for ( size_t i = 0; i < _sdVec.size(); ++i ) { shapes.Clear(); @@ -2403,17 +2659,6 @@ void _ViscousBuilder::getIgnoreFaces(const TopoDS_Shape& solid, bool _ViscousBuilder::makeLayer(_SolidData& data) { - // get all sub-shapes to make layers on - set subIds, faceIds; - subIds = data._noShrinkShapes; - TopExp_Explorer exp( data._solid, TopAbs_FACE ); - for ( ; exp.More(); exp.Next() ) - { - SMESH_subMesh* fSubM = _mesh->GetSubMesh( exp.Current() ); - if ( ! data._ignoreFaceIds.count( fSubM->GetId() )) - faceIds.insert( fSubM->GetId() ); - } - // make a map to find new nodes on sub-shapes shared with other SOLID map< TGeomID, TNode2Edge* >::iterator s2ne; map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin(); @@ -2435,8 +2680,11 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) // Create temporary faces and _LayerEdge's + debugMsg( "######################" ); dumpFunction(SMESH_Comment("makeLayers_")<& edgesByGeom = data._edgesOnShape; + data._stepSize = Precision::Infinite(); data._stepSizeNodes[0] = 0; @@ -2447,34 +2695,19 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) vector< const SMDS_MeshNode*> newNodes; // of a mesh face TNode2Edge::iterator n2e2; - // collect _LayerEdge's of shapes they are based on - vector< _EdgesOnShape >& edgesByGeom = data._edgesOnShape; - const int nbShapes = getMeshDS()->MaxShapeIndex(); - edgesByGeom.resize( nbShapes+1 ); - - // set data of _EdgesOnShape's - if ( SMESH_subMesh* sm = _mesh->GetSubMesh( data._solid )) - { - SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/false); - while ( smIt->more() ) - { - sm = smIt->next(); - if ( sm->GetSubShape().ShapeType() == TopAbs_FACE && - !faceIds.count( sm->GetId() )) - continue; - setShapeData( edgesByGeom[ sm->GetId() ], sm, data ); - } - } // make _LayerEdge's - for ( set::iterator id = faceIds.begin(); id != faceIds.end(); ++id ) + for ( TopExp_Explorer exp( data._solid, TopAbs_FACE ); exp.More(); exp.Next() ) { - const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( *id )); - SMESH_subMesh* sm = _mesh->GetSubMesh( F ); + const TopoDS_Face& F = TopoDS::Face( exp.Current() ); + SMESH_subMesh* sm = _mesh->GetSubMesh( F ); + const TGeomID id = sm->GetId(); + if ( edgesByGeom[ id ]._shape.IsNull() ) + continue; // no layers SMESH_ProxyMesh::SubMesh* proxySub = data._proxyMesh->getFaceSubM( F, /*create=*/true); SMESHDS_SubMesh* smDS = sm->GetSubMeshDS(); - if ( !smDS ) return error(SMESH_Comment("Not meshed face ") << *id, data._index ); + if ( !smDS ) return error(SMESH_Comment("Not meshed face ") << id, data._index ); SMDS_ElemIteratorPtr eIt = smDS->GetElements(); while ( eIt->more() ) @@ -2508,11 +2741,11 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) nbDegenNodes++; } } - TNode2Edge::iterator n2e = data._n2eMap.insert( make_pair( n, (_LayerEdge*)0 )).first; + TNode2Edge::iterator n2e = data._n2eMap.insert({ n, nullptr }).first; if ( !(*n2e).second ) { // add a _LayerEdge - _LayerEdge* edge = new _LayerEdge(); + _LayerEdge* edge = _Factory::NewLayerEdge(); edge->_nodes.push_back( n ); n2e->second = edge; edgesByGeom[ shapeID ]._edges.push_back( edge ); @@ -2542,13 +2775,12 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) if ( !setEdgeData( *edge, edgesByGeom[ shapeID ], helper, data )) return false; - if ( edge->_nodes.size() < 2 ) - edge->Block( data ); - //data._noShrinkShapes.insert( shapeID ); + if ( edge->_nodes.size() < 2 && !noShrink ) + edge->Block( data ); // a sole node is moved only if noShrink } dumpMove(edge->_nodes.back()); - if ( edge->_cosin > faceMaxCosin ) + if ( edge->_cosin > faceMaxCosin && edge->_nodes.size() > 1 ) { faceMaxCosin = edge->_cosin; maxCosinEdge = edge; @@ -2564,7 +2796,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) // create a temporary face const SMDS_MeshElement* newFace = - new _TmpMeshFace( newNodes, --_tmpFaceID, face->getshapeId(), face->getIdInShape() ); + new _TmpMeshFace( newNodes, --_tmpFaceID, face->GetShapeID(), face ); proxySub->AddElement( newFace ); // compute inflation step size by min size of element on a convex surface @@ -2598,6 +2830,42 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) if (( n2e = data._n2eMap.find( *node )) != data._n2eMap.end() ) edge->_neibors.push_back( n2e->second ); } + + // Fix uv of nodes on periodic FACEs (bos #20643) + + if ( eos.ShapeType() != TopAbs_EDGE || + eos.SWOLType() != TopAbs_FACE || + eos.size() == 0 ) + continue; + + const TopoDS_Face& F = TopoDS::Face( eos._sWOL ); + SMESH_MesherHelper faceHelper( *_mesh ); + faceHelper.SetSubShape( F ); + faceHelper.ToFixNodeParameters( true ); + if ( faceHelper.GetPeriodicIndex() == 0 ) + continue; + + SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( F ); + if ( !smDS || smDS->GetNodes() == 0 ) + continue; + + bool toCheck = true; + const double tol = 2 * helper.MaxTolerance( F ); + for ( SMDS_NodeIteratorPtr nIt = smDS->GetNodes(); nIt->more(); ) + { + const SMDS_MeshNode* node = nIt->next(); + gp_XY uvNew( Precision::Infinite(), 0 ); + if ( toCheck ) + { + toCheck = false; + gp_XY uv = faceHelper.GetNodeUV( F, node ); + if ( ! faceHelper.CheckNodeUV( F, node, uvNew, tol, /*force=*/true )) + break; // projection on F failed + if (( uv - uvNew ).Modulus() < Precision::Confusion() ) + break; // current uv is OK + } + faceHelper.CheckNodeUV( F, node, uvNew, tol, /*force=*/true ); + } } data._epsilon = 1e-7; @@ -2748,8 +3016,6 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) { SMESH_MesherHelper helper( *_mesh ); - const int nbTestPnt = 5; // on a FACE sub-shape - BRepLProp_SLProps surfProp( 2, 1e-6 ); data._convexFaces.clear(); @@ -2761,58 +3027,27 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) continue; TopoDS_Face F = TopoDS::Face( eof._shape ); - SMESH_subMesh * sm = eof._subMesh; const TGeomID faceID = eof._shapeID; BRepAdaptor_Surface surface( F, false ); surfProp.SetSurface( surface ); - bool isTooCurved = false; - _ConvexFace cnvFace; - const double oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. ); - SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/true); - while ( smIt->more() ) - { - sm = smIt->next(); - const TGeomID subID = sm->GetId(); - // find _LayerEdge's of a sub-shape - _EdgesOnShape* eos; - if (( eos = data.GetShapeEdges( subID ))) - cnvFace._subIdToEOS.insert( make_pair( subID, eos )); - else - continue; - // check concavity and curvature and limit data._stepSize - const double minCurvature = - 1. / ( eos->_hyp.GetTotalThickness() * ( 1 + theThickToIntersection )); - size_t iStep = Max( 1, eos->_edges.size() / nbTestPnt ); - for ( size_t i = 0; i < eos->_edges.size(); i += iStep ) - { - gp_XY uv = helper.GetNodeUV( F, eos->_edges[ i ]->_nodes[0] ); - surfProp.SetParameters( uv.X(), uv.Y() ); - if ( !surfProp.IsCurvatureDefined() ) - continue; - if ( surfProp.MaxCurvature() * oriFactor > minCurvature ) - { - limitStepSize( data, 0.9 / surfProp.MaxCurvature() * oriFactor ); - isTooCurved = true; - } - if ( surfProp.MinCurvature() * oriFactor > minCurvature ) - { - limitStepSize( data, 0.9 / surfProp.MinCurvature() * oriFactor ); - isTooCurved = true; - } - } - } // loop on sub-shapes of the FACE + cnvFace._face = F; + cnvFace._normalsFixed = false; + cnvFace._isTooCurved = false; - if ( !isTooCurved ) continue; + double maxCurvature = cnvFace.GetMaxCurvature( data, eof, surfProp, helper ); + if ( maxCurvature > 0 ) + { + limitStepSize( data, 0.9 / maxCurvature ); + findEdgesToUpdateNormalNearConvexFace( cnvFace, data, helper ); + } + if ( !cnvFace._isTooCurved ) continue; _ConvexFace & convFace = data._convexFaces.insert( make_pair( faceID, cnvFace )).first->second; - convFace._face = F; - convFace._normalsFixed = false; - // skip a closed surface (data._convexFaces is useful anyway) bool isClosedF = false; helper.SetSubShape( F ); @@ -2825,6 +3060,7 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) if ( isClosedF ) { // limit _LayerEdge::_maxLen on the FACE + const double oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. ); const double minCurvature = 1. / ( eof._hyp.GetTotalThickness() * ( 1 + theThickToIntersection )); map< TGeomID, _EdgesOnShape* >::iterator id2eos = cnvFace._subIdToEOS.find( faceID ); @@ -2836,14 +3072,13 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) _LayerEdge* ledge = eos._edges[ i ]; gp_XY uv = helper.GetNodeUV( F, ledge->_nodes[0] ); surfProp.SetParameters( uv.X(), uv.Y() ); - if ( !surfProp.IsCurvatureDefined() ) - continue; - - if ( surfProp.MaxCurvature() * oriFactor > minCurvature ) - ledge->_maxLen = Min( ledge->_maxLen, 1. / surfProp.MaxCurvature() * oriFactor ); - - if ( surfProp.MinCurvature() * oriFactor > minCurvature ) - ledge->_maxLen = Min( ledge->_maxLen, 1. / surfProp.MinCurvature() * oriFactor ); + if ( surfProp.IsCurvatureDefined() ) + { + double curvature = Max( surfProp.MaxCurvature() * oriFactor, + surfProp.MinCurvature() * oriFactor ); + if ( curvature > minCurvature ) + ledge->SetMaxLen( Min( ledge->_maxLen, 1. / curvature )); + } } } continue; @@ -2863,7 +3098,12 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) for ( size_t j = 0; j < ledge->_simplices.size(); ++j ) if ( ledge->_simplices[j]._nNext->GetPosition()->GetDim() < 2 ) { - convFace._simplexTestEdges.push_back( ledge ); + // do not select _LayerEdge's neighboring sharp EDGEs + bool sharpNbr = false; + for ( size_t iN = 0; iN < ledge->_neibors.size() && !sharpNbr; ++iN ) + sharpNbr = ( ledge->_neibors[iN]->_cosin > theMinSmoothCosin ); + if ( !sharpNbr ) + convFace._simplexTestEdges.push_back( ledge ); break; } } @@ -2913,22 +3153,11 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) // define allowed thickness computeGeomSize( data ); // compute data._geomSize and _LayerEdge::_maxLen - data._maxThickness = 0; - data._minThickness = 1e100; - list< const StdMeshers_ViscousLayers* >::iterator hyp = data._hyps.begin(); - for ( ; hyp != data._hyps.end(); ++hyp ) - { - data._maxThickness = Max( data._maxThickness, (*hyp)->GetTotalThickness() ); - data._minThickness = Min( data._minThickness, (*hyp)->GetTotalThickness() ); - } - //const double tgtThick = /*Min( 0.5 * data._geomSize, */data._maxThickness; // Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's - // boundry inclined to the shape at a sharp angle + // boundary inclined to the shape at a sharp angle - //list< TGeomID > shapesToSmooth; TopTools_MapOfShape edgesOfSmooFaces; - SMESH_MesherHelper helper( *_mesh ); bool ok = true; @@ -2943,32 +3172,34 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) continue; double tgtThick = eos._hyp.GetTotalThickness(); - TopExp_Explorer eExp( edgesByGeom[iS]._shape, TopAbs_EDGE ); - for ( ; eExp.More() && !eos._toSmooth; eExp.Next() ) + SMESH_subMeshIteratorPtr subIt = eos._subMesh->getDependsOnIterator(/*includeSelf=*/false ); + while ( subIt->more() && !eos._toSmooth ) { - TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() ); - vector<_LayerEdge*>& eE = edgesByGeom[ iE ]._edges; - if ( eE.empty() ) continue; + TGeomID iSub = subIt->next()->GetId(); + const vector<_LayerEdge*>& eSub = edgesByGeom[ iSub ]._edges; + if ( eSub.empty() ) continue; double faceSize; - for ( size_t i = 0; i < eE.size() && !eos._toSmooth; ++i ) - if ( eE[i]->_cosin > theMinSmoothCosin ) + for ( size_t i = 0; i < eSub.size() && !eos._toSmooth; ++i ) + if ( eSub[i]->_cosin > theMinSmoothCosin ) { - SMDS_ElemIteratorPtr fIt = eE[i]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face); + SMDS_ElemIteratorPtr fIt = eSub[i]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face); while ( fIt->more() && !eos._toSmooth ) { const SMDS_MeshElement* face = fIt->next(); if ( face->getshapeId() == eos._shapeID && - getDistFromEdge( face, eE[i]->_nodes[0], faceSize )) + getDistFromEdge( face, eSub[i]->_nodes[0], faceSize )) { - eos._toSmooth = needSmoothing( eE[i]->_cosin, tgtThick, faceSize ); + eos._toSmooth = needSmoothing( eSub[i]->_cosin, + tgtThick * eSub[i]->_lenFactor, + faceSize); } } } } if ( eos._toSmooth ) { - for ( eExp.ReInit(); eExp.More(); eExp.Next() ) + for ( TopExp_Explorer eExp( edgesByGeom[iS]._shape, TopAbs_EDGE ); eExp.More(); eExp.Next() ) edgesOfSmooFaces.Add( eExp.Current() ); data.PrepareEdgesToSmoothOnFace( &edgesByGeom[iS], /*substituteSrcNodes=*/false ); @@ -2988,13 +3219,13 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) if ( SMESH_Algo::isDegenerated( E ) || !edgesOfSmooFaces.Contains( E )) continue; - double tgtThick = eos._hyp.GetTotalThickness(); + double tgtThick = eos._hyp.GetTotalThickness(), h0 = eos._hyp.Get1stLayerThickness(); for ( TopoDS_Iterator vIt( E ); vIt.More() && !eos._toSmooth; vIt.Next() ) { TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() ); vector<_LayerEdge*>& eV = edgesByGeom[ iV ]._edges; if ( eV.empty() || eV[0]->Is( _LayerEdge::MULTI_NORMAL )) continue; - gp_Vec eDir = getEdgeDir( E, TopoDS::Vertex( vIt.Value() )); + gp_Vec eDir = getEdgeDir( E, TopoDS::Vertex( vIt.Value() ), h0 ); double angle = eDir.Angle( eV[0]->_normal ); double cosin = Cos( angle ); double cosinAbs = Abs( cosin ); @@ -3012,16 +3243,16 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) if ( endSeg->getshapeId() == (int) iS ) { double segLen = - SMESH_TNodeXYZ( endSeg->GetNode(0) ).Distance( endSeg->GetNode(1 )); - eos._toSmooth = needSmoothing( cosinAbs, tgtThick, segLen ); + SMESH_TNodeXYZ( endSeg->GetNode( 0 )).Distance( endSeg->GetNode( 1 )); + eos._toSmooth = needSmoothing( cosinAbs, tgtThick * eV[0]->_lenFactor, segLen ); } } if ( eos._toSmooth ) { eos._edgeSmoother = new _Smoother1D( curve, eos ); - for ( size_t i = 0; i < eos._edges.size(); ++i ) - eos._edges[i]->Set( _LayerEdge::TO_SMOOTH ); + // for ( size_t i = 0; i < eos._edges.size(); ++i ) + // eos._edges[i]->Set( _LayerEdge::TO_SMOOTH ); } } } @@ -3037,11 +3268,12 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) if ( !eos._hyp.ToSmooth() ) for ( size_t i = 0; i < eos._edges.size(); ++i ) - eos._edges[i]->SetCosin( 0 ); + //eos._edges[i]->SetCosin( 0 ); // keep _cosin to use in limitMaxLenByCurvature() + eos._edges[i]->_lenFactor = 1; } - // Fill _eosC1 to make that C1 FACEs and EGDEs between them to be smoothed as a whole + // Fill _eosC1 to make that C1 FACEs and EDGEs between them to be smoothed as a whole TopTools_MapOfShape c1VV; @@ -3101,19 +3333,21 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) { _EdgesOnShape* eof = data.GetShapeEdges( *face ); if ( !eof ) continue; // other solid - if ( !eos.HasC1( eoe )) - { - eos._eosC1.push_back( eoe ); - eoe->_toSmooth = false; - data.PrepareEdgesToSmoothOnFace( eoe, /*substituteSrcNodes=*/false ); - } - if ( eos._shapeID != eof->_shapeID && !eos.HasC1( eof )) + if ( eos._shapeID == eof->_shapeID ) continue; + if ( !eos.HasC1( eof )) { + // check the FACEs eos._eosC1.push_back( eof ); eof->_toSmooth = false; data.PrepareEdgesToSmoothOnFace( eof, /*substituteSrcNodes=*/false ); smQueue.push_back( eof->_subMesh ); } + if ( !eos.HasC1( eoe )) + { + eos._eosC1.push_back( eoe ); + eoe->_toSmooth = false; + data.PrepareEdgesToSmoothOnFace( eoe, /*substituteSrcNodes=*/false ); + } } } } @@ -3169,7 +3403,7 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) { _EdgesOnShape* eoe = data.GetShapeEdges( *e ); if ( !eoe ) continue; // other solid - gp_XYZ eDir = getEdgeDir( TopoDS::Edge( *e ), V ); + gp_XYZ eDir = getEdgeDir( TopoDS::Edge( *e ), V, eoe->_hyp.Get1stLayerThickness() ); if ( !Precision::IsInfinite( eDir.X() )) dirOfEdges.push_back( make_pair( eoe, eDir.Normalized() )); } @@ -3184,12 +3418,22 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) if ( isC1 ) { double maxEdgeLen = 3 * Min( eov._edges[0]->_maxLen, eov._hyp.GetTotalThickness() ); - double eLen1 = SMESH_Algo::EdgeLength( TopoDS::Edge( dirOfEdges[i].first->_shape )); - double eLen2 = SMESH_Algo::EdgeLength( TopoDS::Edge( dirOfEdges[j].first->_shape )); - if ( eLen1 < maxEdgeLen ) eov._eosC1.push_back( dirOfEdges[i].first ); - if ( eLen2 < maxEdgeLen ) eov._eosC1.push_back( dirOfEdges[j].first ); - dirOfEdges[i].first = 0; - dirOfEdges[j].first = 0; + for ( int isJ = 0; isJ < 2; ++isJ ) // loop on [i,j] + { + size_t k = isJ ? j : i; + const TopoDS_Edge& e = TopoDS::Edge( dirOfEdges[k].first->_shape ); + double eLen = SMESH_Algo::EdgeLength( e ); + if ( eLen < maxEdgeLen ) + { + TopoDS_Shape oppV = SMESH_MesherHelper::IthVertex( 0, e ); + if ( oppV.IsSame( V )) + oppV = SMESH_MesherHelper::IthVertex( 1, e ); + _EdgesOnShape* eovOpp = data.GetShapeEdges( oppV ); + if ( dirOfEdges[k].second * eovOpp->_edges[0]->_normal < 0 ) + eov._eosC1.push_back( dirOfEdges[k].first ); + } + dirOfEdges[k].first = 0; + } } } } // fill _eosC1 of VERTEXes @@ -3199,6 +3443,53 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) return ok; } +//================================================================================ +/*! + * \brief Set up _SolidData::_edgesOnShape + */ +//================================================================================ + +int _ViscousBuilder::makeEdgesOnShape() +{ + const int nbShapes = getMeshDS()->MaxShapeIndex(); + int nbSolidsWL = 0; + + for ( size_t i = 0; i < _sdVec.size(); ++i ) + { + _SolidData& data = _sdVec[ i ]; + vector< _EdgesOnShape >& edgesByGeom = data._edgesOnShape; + edgesByGeom.resize( nbShapes+1 ); + + // set data of _EdgesOnShape's + int nbShapesWL = 0; + if ( SMESH_subMesh* sm = _mesh->GetSubMesh( data._solid )) + { + SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/false); + while ( smIt->more() ) + { + sm = smIt->next(); + if ( sm->GetSubShape().ShapeType() == TopAbs_FACE && + data._ignoreFaceIds.count( sm->GetId() )) + continue; + + setShapeData( edgesByGeom[ sm->GetId() ], sm, data ); + + nbShapesWL += ( sm->GetSubShape().ShapeType() == TopAbs_FACE ); + } + } + if ( nbShapesWL == 0 ) // no shapes with layers in a SOLID + { + data._done = true; + SMESHUtils::FreeVector( edgesByGeom ); + } + else + { + ++nbSolidsWL; + } + } + return nbSolidsWL; +} + //================================================================================ /*! * \brief initialize data of _EdgesOnShape @@ -3222,6 +3513,7 @@ void _ViscousBuilder::setShapeData( _EdgesOnShape& eos, eos._shape.Orientation( helper.GetSubShapeOri( data._solid, eos._shape )); eos._toSmooth = false; eos._data = &data; + eos._mapper2D = nullptr; // set _SWOL map< TGeomID, TopoDS_Shape >::const_iterator s2s = @@ -3269,19 +3561,19 @@ void _ViscousBuilder::setShapeData( _EdgesOnShape& eos, if ( eos.ShapeType() == TopAbs_FACE ) // get normals to elements on a FACE { SMESHDS_SubMesh* smDS = sm->GetSubMeshDS(); - eos._faceNormals.resize( smDS->NbElements() ); + if ( !smDS ) return; + eos._faceNormals.reserve( smDS->NbElements() ); + double oriFactor = helper.IsReversedSubMesh( TopoDS::Face( eos._shape )) ? 1.: -1.; SMDS_ElemIteratorPtr eIt = smDS->GetElements(); - for ( int iF = 0; eIt->more(); ++iF ) + for ( ; eIt->more(); ) { const SMDS_MeshElement* face = eIt->next(); - if ( !SMESH_MeshAlgos::FaceNormal( face, eos._faceNormals[iF], /*normalized=*/true )) - eos._faceNormals[iF].SetCoord( 0,0,0 ); + gp_XYZ& norm = eos._faceNormals[face]; + if ( !SMESH_MeshAlgos::FaceNormal( face, norm, /*normalized=*/true )) + norm.SetCoord( 0,0,0 ); + norm *= oriFactor; } - - if ( !helper.IsReversedSubMesh( TopoDS::Face( eos._shape ))) - for ( size_t iF = 0; iF < eos._faceNormals.size(); ++iF ) - eos._faceNormals[iF].Reverse(); } else // find EOS of adjacent FACEs { @@ -3307,7 +3599,7 @@ void _ViscousBuilder::setShapeData( _EdgesOnShape& eos, bool _EdgesOnShape::GetNormal( const SMDS_MeshElement* face, gp_Vec& norm ) { bool ok = false; - const _EdgesOnShape* eos = 0; + _EdgesOnShape* eos = 0; if ( face->getshapeId() == _shapeID ) { @@ -3321,9 +3613,9 @@ bool _EdgesOnShape::GetNormal( const SMDS_MeshElement* face, gp_Vec& norm ) } if (( eos ) && - ( ok = ( face->getIdInShape() < (int) eos->_faceNormals.size() ))) + ( ok = ( eos->_faceNormals.count( face ) ))) { - norm = eos->_faceNormals[ face->getIdInShape() ]; + norm = eos->_faceNormals[ face ]; } else if ( !eos ) { @@ -3333,6 +3625,17 @@ bool _EdgesOnShape::GetNormal( const SMDS_MeshElement* face, gp_Vec& norm ) return ok; } +//================================================================================ +/*! + * \brief EdgesOnShape destructor + */ +//================================================================================ + +_EdgesOnShape::~_EdgesOnShape() +{ + delete _edgeSmoother; + delete _mapper2D; +} //================================================================================ /*! @@ -3347,12 +3650,13 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, { const SMDS_MeshNode* node = edge._nodes[0]; // source node - edge._len = 0; - edge._maxLen = Precision::Infinite(); - edge._minAngle = 0; - edge._2neibors = 0; - edge._curvature = 0; - edge._flags = 0; + edge._len = 0; + edge._maxLen = Precision::Infinite(); + edge._minAngle = 0; + edge._2neibors = 0; + edge._curvature = 0; + edge._flags = 0; + edge._smooFunction = 0; // -------------------------- // Compute _normal and _cosin @@ -3400,24 +3704,26 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, } // find _normal + bool fromVonF = false; if ( useGeometry ) { - bool fromVonF = ( eos.ShapeType() == TopAbs_VERTEX && - eos.SWOLType() == TopAbs_FACE && - totalNbFaces > 1 ); + fromVonF = ( eos.ShapeType() == TopAbs_VERTEX && + eos.SWOLType() == TopAbs_FACE && + totalNbFaces > 1 ); if ( onShrinkShape && !fromVonF ) // one of faces the node is on has no layers { if ( eos.SWOLType() == TopAbs_EDGE ) { // inflate from VERTEX along EDGE - edge._normal = getEdgeDir( TopoDS::Edge( eos._sWOL ), TopoDS::Vertex( eos._shape )); + edge._normal = getEdgeDir( TopoDS::Edge( eos._sWOL ), TopoDS::Vertex( eos._shape ), + eos._hyp.Get1stLayerThickness(), &eos._isRegularSWOL ); } else if ( eos.ShapeType() == TopAbs_VERTEX ) { // inflate from VERTEX along FACE edge._normal = getFaceDir( TopoDS::Face( eos._sWOL ), TopoDS::Vertex( eos._shape ), - node, helper, normOK, &edge._cosin); + node, helper, normOK/*, &edge._cosin*/); } else { @@ -3506,25 +3812,37 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, break; } case TopAbs_VERTEX: { - //if ( eos.SWOLType() != TopAbs_FACE ) // else _cosin is set by getFaceDir() - { - TopoDS_Vertex V = TopoDS::Vertex( eos._shape ); - gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK ); - double angle = inFaceDir.Angle( edge._normal ); // [0,PI] - edge._cosin = Cos( angle ); - if ( totalNbFaces > 2 || helper.IsSeamShape( node->getshapeId() )) - for ( int iF = totalNbFaces-2; iF >=0; --iF ) - { - F = face2Norm[ iF ].first; - inFaceDir = getFaceDir( F, V, node, helper, normOK=true ); - if ( normOK ) { - double angle = inFaceDir.Angle( edge._normal ); - double cosin = Cos( angle ); - if ( Abs( cosin ) > Abs( edge._cosin )) - edge._cosin = cosin; + TopoDS_Vertex V = TopoDS::Vertex( eos._shape ); + gp_Vec inFaceDir = getFaceDir( face2Norm[0].first , V, node, helper, normOK ); + double angle = inFaceDir.Angle( edge._normal ); // [0,PI] + edge._cosin = Cos( angle ); + if ( fromVonF ) + totalNbFaces--; + if ( totalNbFaces > 1 || helper.IsSeamShape( node->getshapeId() )) + for ( int iF = 1; iF < totalNbFaces; ++iF ) + { + F = face2Norm[ iF ].first; + inFaceDir = getFaceDir( F, V, node, helper, normOK=true ); + if ( normOK ) { + if ( onShrinkShape ) + { + gp_Vec faceNorm = getFaceNormal( node, F, helper, normOK ); + if ( !normOK ) continue; + if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED ) + faceNorm.Reverse(); + angle = 0.5 * M_PI - faceNorm.Angle( edge._normal ); + if ( inFaceDir * edge._normal < 0 ) + angle = M_PI - angle; + } + else + { + angle = inFaceDir.Angle( edge._normal ); } + double cosin = Cos( angle ); + if ( Abs( cosin ) > Abs( edge._cosin )) + edge._cosin = cosin; } - } + } break; } default: @@ -3543,19 +3861,32 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, getMeshDS()->RemoveFreeNode( edge._nodes.back(), 0, /*fromGroups=*/false ); edge._nodes.resize( 1 ); edge._normal.SetCoord( 0,0,0 ); - edge._maxLen = 0; + edge.SetMaxLen( 0 ); } // Set the rest data // -------------------- - edge.SetCosin( edge._cosin ); // to update edge._lenFactor + double realLenFactor = edge.SetCosin( edge._cosin ); // to update edge._lenFactor + // if ( realLenFactor > 3 ) + // { + // edge._cosin = 1; + // if ( onShrinkShape ) + // { + // edge.Set( _LayerEdge::RISKY_SWOL ); + // edge._lenFactor = 2; + // } + // else + // { + // edge._lenFactor = 1; + // } + // } if ( onShrinkShape ) { const SMDS_MeshNode* tgtNode = edge._nodes.back(); if ( SMESHDS_SubMesh* sm = getMeshDS()->MeshElements( data._solid )) - sm->RemoveNode( tgtNode , /*isNodeDeleted=*/false ); + sm->RemoveNode( tgtNode ); // set initial position which is parameters on _sWOL in this case if ( eos.SWOLType() == TopAbs_EDGE ) @@ -3573,7 +3904,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, getMeshDS()->SetNodeOnFace( tgtNode, TopoDS::Face( eos._sWOL ), uv.X(), uv.Y() ); } - if ( edge._nodes.size() > 1 ) + //if ( edge._nodes.size() > 1 ) -- allow RISKY_SWOL on noShrink shape { // check if an angle between a FACE with layers and SWOL is sharp, // else the edge should not inflate @@ -3588,10 +3919,14 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, geomNorm.Reverse(); // inside the SOLID if ( geomNorm * edge._normal < -0.001 ) { - getMeshDS()->RemoveFreeNode( tgtNode, 0, /*fromGroups=*/false ); - edge._nodes.resize( 1 ); + if ( edge._nodes.size() > 1 ) + { + getMeshDS()->RemoveFreeNode( tgtNode, 0, /*fromGroups=*/false ); + edge._nodes.resize( 1 ); + } } - else if ( edge._lenFactor > 3 ) + else if ( realLenFactor > 3 ) /// -- moved to SetCosin() + //else if ( edge._lenFactor > 3 ) { edge._lenFactor = 2; edge.Set( _LayerEdge::RISKY_SWOL ); @@ -3619,7 +3954,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, if ( eos.ShapeType() == TopAbs_EDGE /*|| ( onShrinkShape && posType == SMDS_TOP_VERTEX && fabs( edge._cosin ) < 1e-10 )*/) { - edge._2neibors = new _2NearEdges; + edge._2neibors = _Factory::NewNearEdges(); // target nodes instead of source ones will be set later } @@ -3680,7 +4015,7 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node, isOK = false; return p.XYZ(); } - Quantity_Parameter U,V; + Standard_Real U,V; projector.LowerDistanceParameters(U,V); uv.SetCoord( U,V ); } @@ -3758,7 +4093,7 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node, bool _ViscousBuilder::getFaceNormalAtSingularity( const gp_XY& uv, const TopoDS_Face& face, - SMESH_MesherHelper& helper, + SMESH_MesherHelper& /*helper*/, gp_Dir& normal ) { BRepAdaptor_Surface surface( face ); @@ -3958,7 +4293,7 @@ void _OffsetPlane::ComputeIntersectionLine( _OffsetPlane& pln, // parallel planes - intersection is an offset of the common EDGE gp_Pnt p = BRep_Tool::Pnt( V ); linePos = 0.5 * (( p.XYZ() + n1 ) + ( p.XYZ() + n2 )); - lineDir = getEdgeDir( E, V ); + lineDir = getEdgeDir( E, V, 0.1 * SMESH_Algo::EdgeLength( E )); } else { @@ -4035,7 +4370,7 @@ gp_XYZ _OffsetPlane::GetCommonPoint(bool& isFound, //================================================================================ /*! - * \brief Find 2 neigbor nodes of a node on EDGE + * \brief Find 2 neighbor nodes of a node on EDGE */ //================================================================================ @@ -4080,7 +4415,35 @@ bool _ViscousBuilder::findNeiborsOnEdge(const _LayerEdge* edge, //================================================================================ /*! - * \brief Set _curvature and _2neibors->_plnNorm by 2 neigbor nodes residing the same EDGE + * \brief Create _Curvature + */ +//================================================================================ + +_Curvature* _Curvature::New( double avgNormProj, double avgDist ) +{ + // double _r; // radius + // double _k; // factor to correct node smoothed position + // double _h2lenRatio; // avgNormProj / (2*avgDist) + // gp_Pnt2d _uv; // UV used in putOnOffsetSurface() + + _Curvature* c = 0; + if ( fabs( avgNormProj / avgDist ) > 1./200 ) + { + c = _Factory::NewCurvature(); + c->_r = avgDist * avgDist / avgNormProj; + c->_k = avgDist * avgDist / c->_r / c->_r; + //c->_k = avgNormProj / c->_r; + c->_k *= ( c->_r < 0 ? 1/1.1 : 1.1 ); // not to be too restrictive + c->_h2lenRatio = avgNormProj / ( avgDist + avgDist ); + + c->_uv.SetCoord( 0., 0. ); + } + return c; +} + +//================================================================================ +/*! + * \brief Set _curvature and _2neibors->_plnNorm by 2 neighbor nodes residing the same EDGE */ //================================================================================ @@ -4091,6 +4454,8 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1, { if ( eos.ShapeType() != TopAbs_EDGE ) return; + if ( _curvature && Is( SMOOTHED_C1 )) + return; gp_XYZ pos = SMESH_TNodeXYZ( _nodes[0] ); gp_XYZ vec1 = pos - SMESH_TNodeXYZ( n1 ); @@ -4103,7 +4468,6 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1, _2neibors->_wgt[1] = 1 - vec2.Modulus() / sumLen; double avgNormProj = 0.5 * ( _normal * vec1 + _normal * vec2 ); double avgLen = 0.5 * ( vec1.Modulus() + vec2.Modulus() ); - if ( _curvature ) delete _curvature; _curvature = _Curvature::New( avgNormProj, avgLen ); // if ( _curvature ) // debugMsg( _nodes[0]->GetID() @@ -4133,7 +4497,7 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1, //================================================================================ /*! * \brief Copy data from a _LayerEdge of other SOLID and based on the same node; - * this and other _LayerEdge's are inflated along a FACE or an EDGE + * this and the other _LayerEdge are inflated along a FACE or an EDGE */ //================================================================================ @@ -4147,8 +4511,11 @@ gp_XYZ _LayerEdge::Copy( _LayerEdge& other, _lenFactor = other._lenFactor; _cosin = other._cosin; _2neibors = other._2neibors; - _curvature = 0; std::swap( _curvature, other._curvature ); - _2neibors = 0; std::swap( _2neibors, other._2neibors ); + _curvature = other._curvature; + _2neibors = other._2neibors; + _maxLen = Precision::Infinite();//other._maxLen; + _flags = 0; + _smooFunction = 0; gp_XYZ lastPos( 0,0,0 ); if ( eos.SWOLType() == TopAbs_EDGE ) @@ -4177,12 +4544,23 @@ gp_XYZ _LayerEdge::Copy( _LayerEdge& other, */ //================================================================================ -void _LayerEdge::SetCosin( double cosin ) +double _LayerEdge::SetCosin( double cosin ) { _cosin = cosin; cosin = Abs( _cosin ); - //_lenFactor = ( cosin < 1.-1e-12 ) ? Min( 2., 1./sqrt(1-cosin*cosin )) : 1.0; - _lenFactor = ( cosin < 1.-1e-12 ) ? 1./sqrt(1-cosin*cosin ) : 1.0; + //_lenFactor = ( cosin < 1.-1e-12 ) ? 1./sqrt(1-cosin*cosin ) : 1.0; + double realLenFactor; + if ( cosin < 1.-1e-12 ) + { + _lenFactor = realLenFactor = 1./sqrt(1-cosin*cosin ); + } + else + { + _lenFactor = 1; + realLenFactor = Precision::Infinite(); + } + + return realLenFactor; } //================================================================================ @@ -4257,7 +4635,7 @@ void _Simplex::SortSimplices(vector<_Simplex>& simplices) //================================================================================ /*! - * \brief DEBUG. Create groups contating temorary data of _LayerEdge's + * \brief DEBUG. Create groups containing temporary data of _LayerEdge's */ //================================================================================ @@ -4344,7 +4722,7 @@ void _ViscousBuilder::computeGeomSize( _SolidData& data ) _EdgesOnShape& eos = data._edgesOnShape[ iS ]; if ( eos._edges.empty() ) continue; - // get neighbor faces intersection with which should not be considered since + // get neighbor faces, intersection with which should not be considered since // collisions are avoided by means of smoothing set< TGeomID > neighborFaces; if ( eos._hyp.ToSmooth() ) @@ -4363,17 +4741,89 @@ void _ViscousBuilder::computeGeomSize( _SolidData& data ) double thinkness = eos._hyp.GetTotalThickness(); for ( size_t i = 0; i < eos._edges.size(); ++i ) { - if ( eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue; - eos._edges[i]->_maxLen = thinkness; + if ( eos._edges[i]->_nodes.size() < 2 ) continue; + eos._edges[i]->SetMaxLen( thinkness ); eos._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon, eos, &face ); if ( intersecDist > 0 && face ) { data._geomSize = Min( data._geomSize, intersecDist ); if ( !neighborFaces.count( face->getshapeId() )) - eos._edges[i]->_maxLen = Min( thinkness, intersecDist / ( face->GetID() < 0 ? 3. : 2. )); + eos[i]->SetMaxLen( Min( thinkness, intersecDist / ( face->GetID() < 0 ? 3. : 2. ))); + } + } + } + + data._maxThickness = 0; + data._minThickness = 1e100; + list< const StdMeshers_ViscousLayers* >::iterator hyp = data._hyps.begin(); + for ( ; hyp != data._hyps.end(); ++hyp ) + { + data._maxThickness = Max( data._maxThickness, (*hyp)->GetTotalThickness() ); + data._minThickness = Min( data._minThickness, (*hyp)->GetTotalThickness() ); + } + + // Limit inflation step size by geometry size found by intersecting + // normals of _LayerEdge's with mesh faces + if ( data._stepSize > 0.3 * data._geomSize ) + limitStepSize( data, 0.3 * data._geomSize ); + + if ( data._stepSize > data._minThickness ) + limitStepSize( data, data._minThickness ); + + + // ------------------------------------------------------------------------- + // Detect _LayerEdge which can't intersect with opposite or neighbor layer, + // so no need in detecting intersection at each inflation step + // ------------------------------------------------------------------------- + + int nbSteps = data._maxThickness / data._stepSize; + if ( nbSteps < 3 || nbSteps * data._n2eMap.size() < 100000 ) + return; + + vector< const SMDS_MeshElement* > closeFaces; + int nbDetected = 0; + + for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS ) + { + _EdgesOnShape& eos = data._edgesOnShape[ iS ]; + if ( eos._edges.empty() || eos.ShapeType() != TopAbs_FACE ) + continue; + + for ( size_t i = 0; i < eos.size(); ++i ) + { + SMESH_NodeXYZ p( eos[i]->_nodes[0] ); + double radius = data._maxThickness + 2 * eos[i]->_maxLen; + closeFaces.clear(); + searcher->GetElementsInSphere( p, radius, SMDSAbs_Face, closeFaces ); + + bool toIgnore = true; + for ( size_t iF = 0; iF < closeFaces.size() && toIgnore; ++iF ) + if ( !( toIgnore = ( closeFaces[ iF ]->getshapeId() == eos._shapeID || + data._ignoreFaceIds.count( closeFaces[ iF ]->getshapeId() )))) + { + // check if a _LayerEdge will inflate in a direction opposite to a direction + // toward a close face + bool allBehind = true; + for ( int iN = 0; iN < closeFaces[ iF ]->NbCornerNodes() && allBehind; ++iN ) + { + SMESH_NodeXYZ pi( closeFaces[ iF ]->GetNode( iN )); + allBehind = (( pi - p ) * eos[i]->_normal < 0.1 * data._stepSize ); + } + toIgnore = allBehind; + } + + + if ( toIgnore ) // no need to detect intersection + { + eos[i]->Set( _LayerEdge::INTERSECTED ); + ++nbDetected; } } } + + debugMsg( "Nb LE to intersect " << data._n2eMap.size()-nbDetected << ", ignore " << nbDetected ); + + return; } //================================================================================ @@ -4386,24 +4836,20 @@ bool _ViscousBuilder::inflate(_SolidData& data) { SMESH_MesherHelper helper( *_mesh ); - // Limit inflation step size by geometry size found by itersecting - // normals of _LayerEdge's with mesh faces - if ( data._stepSize > 0.3 * data._geomSize ) - limitStepSize( data, 0.3 * data._geomSize ); - const double tgtThick = data._maxThickness; - if ( data._stepSize > data._minThickness ) - limitStepSize( data, data._minThickness ); if ( data._stepSize < 1. ) data._epsilon = data._stepSize * 1e-7; debugMsg( "-- geomSize = " << data._geomSize << ", stepSize = " << data._stepSize ); + _pyDump->Pause(); findCollisionEdges( data, helper ); limitMaxLenByCurvature( data, helper ); + _pyDump->Resume(); + // limit length of _LayerEdge's around MULTI_NORMAL _LayerEdge's for ( size_t i = 0; i < data._edgesOnShape.size(); ++i ) if ( data._edgesOnShape[i].ShapeType() == TopAbs_VERTEX && @@ -4483,7 +4929,10 @@ bool _ViscousBuilder::inflate(_SolidData& data) const double shapeTgtThick = eos._hyp.GetTotalThickness(); for ( size_t i = 0; i < eos._edges.size(); ++i ) { - avgThick += Min( 1., eos._edges[i]->_len / shapeTgtThick ); + if ( eos._edges[i]->_nodes.size() > 1 ) + avgThick += Min( 1., eos._edges[i]->_len / shapeTgtThick ); + else + avgThick += shapeTgtThick; nbActiveEdges += ( ! eos._edges[i]->Is( _LayerEdge::BLOCKED )); } } @@ -4505,6 +4954,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) break; } #endif + // new step size limitStepSize( data, 0.25 * distToIntersection ); if ( data._stepSizeNodes[0] ) @@ -4663,24 +5113,41 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, for ( size_t iEOS = 0; iEOS < eosC1.size(); ++iEOS ) { isConcaveFace[ iEOS ] = data._concaveFaces.count( eosC1[ iEOS ]->_shapeID ); - vector< _LayerEdge* > & edges = eosC1[ iEOS ]->_edges; - for ( size_t i = 0; i < edges.size(); ++i ) - if ( edges[i]->Is( _LayerEdge::MOVED ) || - edges[i]->Is( _LayerEdge::NEAR_BOUNDARY )) - movedEdges.push_back( edges[i] ); + if ( eosC1[ iEOS ]->_mapper2D ) + { + // compute node position by boundary node position in structured mesh + dumpFunction(SMESH_Comment("map2dS")<_mapper2D->ComputeNodePositions(); + + for ( _LayerEdge* le : eosC1[ iEOS ]->_edges ) + le->_pos.back() = SMESH_NodeXYZ( le->_nodes.back() ); + + dumpFunctionEnd(); + } + else + { + for ( _LayerEdge* le : eosC1[ iEOS ]->_edges ) + if ( le->Is( _LayerEdge::MOVED ) || + le->Is( _LayerEdge::NEAR_BOUNDARY )) + movedEdges.push_back( le ); + } makeOffsetSurface( *eosC1[ iEOS ], helper ); } int step = 0, stepLimit = 5, nbBad = 0; while (( ++step <= stepLimit ) || improved ) { - dumpFunction(SMESH_Comment("smooth")<_mapper2D ) + continue; vector< _LayerEdge* > & edges = eosC1[ iEOS ]->_edges; for ( size_t i = 0; i < edges.size(); ++i ) { @@ -4737,7 +5209,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, if ( nbBad == oldBadNb && nbBad > 0 && - step < stepLimit ) // smooth w/o chech of validity + step < stepLimit ) // smooth w/o check of validity { dumpFunctionEnd(); dumpFunction(SMESH_Comment("smoothWoCheck")<_eosConcaVer.empty() || nbBad > 0 ) - putOnOffsetSurface( *eosC1[ iEOS ], infStep, eosC1 ); + putOnOffsetSurface( *eosC1[ iEOS ], -infStep, eosC1 ); } //if ( !badEdges.empty() ) @@ -4826,6 +5298,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, _LayerEdge* edge = eos._edges[i]; if ( edge->_nodes.size() < 2 ) continue; SMESH_TNodeXYZ tgtXYZ = edge->_nodes.back(); + //SMESH_TNodeXYZ prevXYZ = edge->_nodes[0]; gp_XYZ prevXYZ = edge->PrevCheckPos( &eos ); //const gp_XYZ& prevXYZ = edge->PrevPos(); for ( size_t j = 0; j < edge->_simplices.size(); ++j ) @@ -4865,6 +5338,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, const SMDS_MeshElement* intFace = 0; const SMDS_MeshElement* closestFace = 0; _LayerEdge* le = 0; + bool is1stBlocked = true; // dbg for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS ) { _EdgesOnShape& eos = data._edgesOnShape[ iS ]; @@ -4927,7 +5401,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, // ignore intersection of a _LayerEdge based on a _ConvexFace with a face // lying on this _ConvexFace if ( _ConvexFace* convFace = data.GetConvexFace( intFace->getshapeId() )) - if ( convFace->_subIdToEOS.count ( eos._shapeID )) + if ( convFace->_isTooCurved && convFace->_subIdToEOS.count ( eos._shapeID )) continue; // ignore intersection of a _LayerEdge based on a FACE with an element on this FACE @@ -4936,18 +5410,19 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, continue; // ignore intersection with intFace of an adjacent FACE - if ( dist > 0 ) + if ( dist > 0.01 * eos._edges[i]->_len ) { bool toIgnore = false; - if ( eos._edges[i]->Is( _LayerEdge::TO_SMOOTH )) + if ( eos._toSmooth ) { const TopoDS_Shape& S = getMeshDS()->IndexToShape( intFace->getshapeId() ); if ( !S.IsNull() && S.ShapeType() == TopAbs_FACE ) { - TopExp_Explorer edge( eos._shape, TopAbs_EDGE ); - for ( ; !toIgnore && edge.More(); edge.Next() ) - // is adjacent - has a common EDGE - toIgnore = ( helper.IsSubShape( edge.Current(), S )); + TopExp_Explorer sub( eos._shape, + eos.ShapeType() == TopAbs_FACE ? TopAbs_EDGE : TopAbs_VERTEX ); + for ( ; !toIgnore && sub.More(); sub.Next() ) + // is adjacent - has a common EDGE or VERTEX + toIgnore = ( helper.IsSubShape( sub.Current(), S )); if ( toIgnore ) // check angle between normals { @@ -4974,17 +5449,19 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, if ( toBlockInfaltion && dist < ( eos._edges[i]->_len * theThickToIntersection )) { + if ( is1stBlocked ) { is1stBlocked = false; // debug + dumpFunction(SMESH_Comment("blockIntersected") <Set( _LayerEdge::INTERSECTED ); // not to intersect eos._edges[i]->Block( data ); // not to inflate - if ( _EdgesOnShape* eof = data.GetShapeEdges( intFace->getshapeId() )) + //if ( _EdgesOnShape* eof = data.GetShapeEdges( intFace->getshapeId() )) { // block _LayerEdge's, on top of which intFace is if ( const _TmpMeshFace* f = dynamic_cast< const _TmpMeshFace*>( intFace )) { - const SMDS_MeshElement* srcFace = - eof->_subMesh->GetSubMeshDS()->GetElement( f->getIdInShape() ); - SMDS_ElemIteratorPtr nIt = srcFace->nodesIterator(); + const SMDS_MeshElement* srcFace = f->_srcFace; + SMDS_ElemIteratorPtr nIt = srcFace->nodesIterator(); while ( nIt->more() ) { const SMDS_MeshNode* srcNode = static_cast( nIt->next() ); @@ -5007,11 +5484,16 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, } // loop on eos._edges } // loop on data._edgesOnShape + if ( !is1stBlocked ) + { + dumpFunctionEnd(); + } + if ( closestFace && le ) { #ifdef __myDEBUG SMDS_MeshElement::iterator nIt = closestFace->begin_nodes(); - cout << "Shortest distance: _LayerEdge nodes: tgt " << le->_nodes.back()->GetID() + cout << "#Shortest distance: _LayerEdge nodes: tgt " << le->_nodes.back()->GetID() << " src " << le->_nodes[0]->GetID()<< ", intersection with face (" << (*nIt++)->GetID()<<" "<< (*nIt++)->GetID()<<" "<< (*nIt++)->GetID() << ") distance = " << distToIntersection<< endl; @@ -5199,13 +5681,14 @@ void _ViscousBuilder::makeOffsetSurface( _EdgesOnShape& eos, SMESH_MesherHelper& // find offset gp_Pnt tgtP = SMESH_TNodeXYZ( eos._edgeForOffset->_nodes.back() ); /*gp_Pnt2d uv=*/baseSurface->ValueOfUV( tgtP, Precision::Confusion() ); - double offset = baseSurface->Gap(); + eos._offsetValue = baseSurface->Gap(); eos._offsetSurf.Nullify(); try { - BRepOffsetAPI_MakeOffsetShape offsetMaker( eos._shape, -offset, Precision::Confusion() ); + BRepOffsetAPI_MakeOffsetShape offsetMaker; + offsetMaker.PerformByJoin( eos._shape, -eos._offsetValue, Precision::Confusion() ); if ( !offsetMaker.IsDone() ) return; TopExp_Explorer fExp( offsetMaker.Shape(), TopAbs_FACE ); @@ -5217,7 +5700,7 @@ void _ViscousBuilder::makeOffsetSurface( _EdgesOnShape& eos, SMESH_MesherHelper& eos._offsetSurf = new ShapeAnalysis_Surface( surf ); } - catch ( Standard_Failure ) + catch ( Standard_Failure& ) { } } @@ -5232,7 +5715,7 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, int infStep, vector< _EdgesOnShape* >& eosC1, int smooStep, - bool moveAll ) + int moveAll ) { _EdgesOnShape * eof = & eos; if ( eos.ShapeType() != TopAbs_FACE ) // eos is a boundary of C1 FACE, look for the FACE eos @@ -5257,24 +5740,41 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, return; double preci = BRep_Tool::Tolerance( TopoDS::Face( eof->_shape )), vol; + bool neighborHasRiskySWOL = false; for ( size_t i = 0; i < eos._edges.size(); ++i ) { _LayerEdge* edge = eos._edges[i]; edge->Unset( _LayerEdge::MARKED ); if ( edge->Is( _LayerEdge::BLOCKED ) || !edge->_curvature ) continue; - if ( !moveAll && !edge->Is( _LayerEdge::MOVED )) + if ( moveAll == _LayerEdge::UPD_NORMAL_CONV ) + { + if ( !edge->Is( _LayerEdge::UPD_NORMAL_CONV )) + continue; + } + else if ( moveAll == _LayerEdge::RISKY_SWOL ) + { + if ( !edge->Is( _LayerEdge::RISKY_SWOL ) || + edge->_cosin < 0 ) continue; + } + else if ( !moveAll && !edge->Is( _LayerEdge::MOVED )) + continue; int nbBlockedAround = 0; for ( size_t iN = 0; iN < edge->_neibors.size(); ++iN ) + { nbBlockedAround += edge->_neibors[iN]->Is( _LayerEdge::BLOCKED ); + if ( edge->_neibors[iN]->Is( _LayerEdge::RISKY_SWOL ) && + edge->_neibors[iN]->_cosin > 0 ) + neighborHasRiskySWOL = true; + } if ( nbBlockedAround > 1 ) continue; gp_Pnt tgtP = SMESH_TNodeXYZ( edge->_nodes.back() ); gp_Pnt2d uv = eof->_offsetSurf->NextValueOfUV( edge->_curvature->_uv, tgtP, preci ); - if ( eof->_offsetSurf->Gap() > edge->_len ) continue; // NextValueOfUV() bug + if ( eof->_offsetSurf->Gap() > edge->_len ) continue; // NextValueOfUV() bug edge->_curvature->_uv = uv; if ( eof->_offsetSurf->Gap() < 10 * preci ) continue; // same pos @@ -5293,9 +5793,44 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, edge->_pos.back() = newP; edge->Set( _LayerEdge::MARKED ); + if ( moveAll == _LayerEdge::UPD_NORMAL_CONV ) + { + edge->_normal = ( newP - prevP ).Normalized(); + } + // if ( edge->_len < eof->_offsetValue ) + // edge->_len = eof->_offsetValue; + + if ( !eos._sWOL.IsNull() ) // RISKY_SWOL + { + double change = eof->_offsetSurf->Gap() / eof->_offsetValue; + if (( newP - tgtP.XYZ() ) * edge->_normal < 0 ) + change = 1 - change; + else + change = 1 + change; + gp_XYZ shitfVec = tgtP.XYZ() - SMESH_NodeXYZ( edge->_nodes[0] ); + gp_XYZ newShiftVec = shitfVec * change; + double shift = edge->_normal * shitfVec; + double newShift = edge->_normal * newShiftVec; + newP = tgtP.XYZ() + edge->_normal * ( newShift - shift ); + + uv = eof->_offsetSurf->NextValueOfUV( edge->_curvature->_uv, newP, preci ); + if ( eof->_offsetSurf->Gap() < edge->_len ) + { + edge->_curvature->_uv = uv; + newP = eof->_offsetSurf->Value( uv ).XYZ(); + } + n->setXYZ( newP.X(), newP.Y(), newP.Z()); + if ( !edge->UpdatePositionOnSWOL( n, /*tol=*/10 * edge->_len / ( edge->NbSteps() + 1 ), + eos, eos.GetData().GetHelper() )) + { + debugMsg("UpdatePositionOnSWOL fails in putOnOffsetSurface()" ); + } + } } } + + #ifdef _DEBUG_ // dumpMove() for debug size_t i = 0; @@ -5304,16 +5839,59 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, break; if ( i < eos._edges.size() ) { - dumpFunction(SMESH_Comment("putOnOffsetSurface_F") << eos._shapeID - << "_InfStep" << infStep << "_" << smooStep ); + dumpFunction(SMESH_Comment("putOnOffsetSurface_") << eos.ShapeTypeLetter() << eos._shapeID + << "_InfStep" << infStep << "_" << Abs( smooStep )); for ( ; i < eos._edges.size(); ++i ) { - if ( eos._edges[i]->Is( _LayerEdge::MARKED )) + if ( eos._edges[i]->Is( _LayerEdge::MARKED )) { dumpMove( eos._edges[i]->_nodes.back() ); + } } dumpFunctionEnd(); } #endif + + _ConvexFace* cnvFace; + if ( moveAll != _LayerEdge::UPD_NORMAL_CONV && + eos.ShapeType() == TopAbs_FACE && + (cnvFace = eos.GetData().GetConvexFace( eos._shapeID )) && + !cnvFace->_normalsFixedOnBorders ) + { + // put on the surface nodes built on FACE boundaries + SMESH_subMeshIteratorPtr smIt = eos._subMesh->getDependsOnIterator(/*includeSelf=*/false); + while ( smIt->more() ) + { + SMESH_subMesh* sm = smIt->next(); + _EdgesOnShape* subEOS = eos.GetData().GetShapeEdges( sm->GetId() ); + if ( !subEOS->_sWOL.IsNull() ) continue; + if ( std::find( eosC1.begin(), eosC1.end(), subEOS ) != eosC1.end() ) continue; + + putOnOffsetSurface( *subEOS, infStep, eosC1, smooStep, _LayerEdge::UPD_NORMAL_CONV ); + } + cnvFace->_normalsFixedOnBorders = true; + } + + + // bos #20643 + // negative smooStep means "final step", where we don't treat RISKY_SWOL edges + // as edges based on FACE are a bit late comparing with them + if ( smooStep >= 0 && + neighborHasRiskySWOL && + moveAll != _LayerEdge::RISKY_SWOL && + eos.ShapeType() == TopAbs_FACE ) + { + // put on the surface nodes built on FACE boundaries + SMESH_subMeshIteratorPtr smIt = eos._subMesh->getDependsOnIterator(/*includeSelf=*/false); + while ( smIt->more() ) + { + SMESH_subMesh* sm = smIt->next(); + _EdgesOnShape* subEOS = eos.GetData().GetShapeEdges( sm->GetId() ); + if ( subEOS->_sWOL.IsNull() ) continue; + if ( std::find( eosC1.begin(), eosC1.end(), subEOS ) != eosC1.end() ) continue; + + putOnOffsetSurface( *subEOS, infStep, eosC1, smooStep, _LayerEdge::RISKY_SWOL ); + } + } } //================================================================================ @@ -5413,94 +5991,240 @@ Handle(Geom_Curve) _Smoother1D::CurveForSmooth( const TopoDS_Edge& E, //================================================================================ /*! - * \brief smooth _LayerEdge's on a staight EDGE or circular EDGE + * \brief Smooth edges on EDGE */ //================================================================================ -bool _Smoother1D::smoothAnalyticEdge( _SolidData& data, - Handle(ShapeAnalysis_Surface)& surface, - const TopoDS_Face& F, - SMESH_MesherHelper& helper) +bool _Smoother1D::Perform(_SolidData& data, + Handle(ShapeAnalysis_Surface)& surface, + const TopoDS_Face& F, + SMESH_MesherHelper& helper ) { - if ( !isAnalytic() ) return false; + if ( _leParams.empty() || ( !isAnalytic() && _offPoints.empty() )) + prepare( data ); + + findEdgesToSmooth(); + if ( isAnalytic() ) + return smoothAnalyticEdge( data, surface, F, helper ); + else + return smoothComplexEdge ( data, surface, F, helper ); +} - const size_t iFrom = 0, iTo = _eos._edges.size(); +//================================================================================ +/*! + * \brief Find edges to smooth + */ +//================================================================================ - if ( _anaCurve->IsKind( STANDARD_TYPE( Geom_Line ))) +void _Smoother1D::findEdgesToSmooth() +{ + _LayerEdge* leOnV[2] = { getLEdgeOnV(0), getLEdgeOnV(1) }; + for ( int iEnd = 0; iEnd < 2; ++iEnd ) + if ( leOnV[iEnd]->Is( _LayerEdge::NORMAL_UPDATED )) + _leOnV[iEnd]._cosin = Abs( _edgeDir[iEnd].Normalized() * leOnV[iEnd]->_normal ); + + _eToSmooth[0].first = _eToSmooth[0].second = 0; + + for ( size_t i = 0; i < _eos.size(); ++i ) { - if ( F.IsNull() ) // 3D + if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) { - SMESH_TNodeXYZ p0 ( _eos._edges[iFrom]->_2neibors->tgtNode(0) ); - SMESH_TNodeXYZ p1 ( _eos._edges[iTo-1]->_2neibors->tgtNode(1) ); - SMESH_TNodeXYZ pSrc0( _eos._edges[iFrom]->_2neibors->srcNode(0) ); - SMESH_TNodeXYZ pSrc1( _eos._edges[iTo-1]->_2neibors->srcNode(1) ); - gp_XYZ newPos, lineDir = pSrc1 - pSrc0; - _LayerEdge* vLE0 = _eos._edges[iFrom]->_2neibors->_edges[0]; - _LayerEdge* vLE1 = _eos._edges[iTo-1]->_2neibors->_edges[1]; - bool shiftOnly = ( vLE0->Is( _LayerEdge::NORMAL_UPDATED ) || - vLE0->Is( _LayerEdge::BLOCKED ) || - vLE1->Is( _LayerEdge::NORMAL_UPDATED ) || - vLE1->Is( _LayerEdge::BLOCKED )); - for ( size_t i = iFrom; i < iTo; ++i ) - { - _LayerEdge* edge = _eos._edges[i]; - SMDS_MeshNode* tgtNode = const_cast( edge->_nodes.back() ); - newPos = p0 * ( 1. - _leParams[i] ) + p1 * _leParams[i]; - - if ( shiftOnly || edge->Is( _LayerEdge::NORMAL_UPDATED )) - { - gp_XYZ curPos = SMESH_TNodeXYZ ( tgtNode ); - double shift = ( lineDir * ( newPos - pSrc0 ) - - lineDir * ( curPos - pSrc0 )); - newPos = curPos + lineDir * shift / lineDir.SquareModulus(); - } - if ( edge->Is( _LayerEdge::BLOCKED )) - { - SMESH_TNodeXYZ pSrc( edge->_nodes[0] ); - double curThick = pSrc.SquareDistance( tgtNode ); - double newThink = ( pSrc - newPos ).SquareModulus(); - if ( newThink > curThick ) - continue; - } - edge->_pos.back() = newPos; - tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() ); - dumpMove( tgtNode ); - } + if ( needSmoothing( _leOnV[0]._cosin, + _eos[i]->_len * leOnV[0]->_lenFactor, _curveLen * _leParams[i] ) || + isToSmooth( i ) + ) + _eos[i]->Set( _LayerEdge::TO_SMOOTH ); + else + break; } - else // 2D - { - _LayerEdge* e0 = getLEdgeOnV( 0 ); - _LayerEdge* e1 = getLEdgeOnV( 1 ); - gp_XY uv0 = e0->LastUV( F, *data.GetShapeEdges( e0 )); - gp_XY uv1 = e1->LastUV( F, *data.GetShapeEdges( e1 )); - if ( e0->_nodes.back() == e1->_nodes.back() ) // closed edge - { - int iPeriodic = helper.GetPeriodicIndex(); - if ( iPeriodic == 1 || iPeriodic == 2 ) - { - uv1.SetCoord( iPeriodic, helper.GetOtherParam( uv1.Coord( iPeriodic ))); - if ( uv0.Coord( iPeriodic ) > uv1.Coord( iPeriodic )) - std::swap( uv0, uv1 ); - } - } - const gp_XY rangeUV = uv1 - uv0; - for ( size_t i = iFrom; i < iTo; ++i ) - { - if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue; - gp_XY newUV = uv0 + _leParams[i] * rangeUV; - _eos._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 ); + _eToSmooth[0].second = i+1; + } - gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() ); - SMDS_MeshNode* tgtNode = const_cast( _eos._edges[i]->_nodes.back() ); - tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() ); - dumpMove( tgtNode ); + _eToSmooth[1].first = _eToSmooth[1].second = _eos.size(); - SMDS_FacePosition* pos = static_cast( tgtNode->GetPosition() ); - pos->SetUParameter( newUV.X() ); - pos->SetVParameter( newUV.Y() ); - } - } - return true; + for ( int i = _eos.size() - 1; i >= _eToSmooth[0].second; --i ) + { + if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) + { + if ( needSmoothing( _leOnV[1]._cosin, + _eos[i]->_len * leOnV[1]->_lenFactor, _curveLen * ( 1.-_leParams[i] )) || + isToSmooth( i )) + _eos[i]->Set( _LayerEdge::TO_SMOOTH ); + else + break; + } + _eToSmooth[1].first = i; + } +} + +//================================================================================ +/*! + * \brief Check if iE-th _LayerEdge needs smoothing + */ +//================================================================================ + +bool _Smoother1D::isToSmooth( int iE ) +{ + SMESH_NodeXYZ pi( _eos[iE]->_nodes[0] ); + SMESH_NodeXYZ p0( _eos[iE]->_2neibors->srcNode(0) ); + SMESH_NodeXYZ p1( _eos[iE]->_2neibors->srcNode(1) ); + gp_XYZ seg0 = pi - p0; + gp_XYZ seg1 = p1 - pi; + gp_XYZ tangent = seg0 + seg1; + double tangentLen = tangent.Modulus(); + double segMinLen = Min( seg0.Modulus(), seg1.Modulus() ); + if ( tangentLen < std::numeric_limits::min() ) + return false; + tangent /= tangentLen; + + for ( size_t i = 0; i < _eos[iE]->_neibors.size(); ++i ) + { + _LayerEdge* ne = _eos[iE]->_neibors[i]; + if ( !ne->Is( _LayerEdge::TO_SMOOTH ) || + ne->_nodes.size() < 2 || + ne->_nodes[0]->GetPosition()->GetDim() != 2 ) + continue; + gp_XYZ edgeVec = SMESH_NodeXYZ( ne->_nodes.back() ) - SMESH_NodeXYZ( ne->_nodes[0] ); + double proj = edgeVec * tangent; + if ( needSmoothing( 1., proj, segMinLen )) + return true; + } + return false; +} + +//================================================================================ +/*! + * \brief smooth _LayerEdge's on a staight EDGE or circular EDGE + */ +//================================================================================ + +bool _Smoother1D::smoothAnalyticEdge( _SolidData& data, + Handle(ShapeAnalysis_Surface)& surface, + const TopoDS_Face& F, + SMESH_MesherHelper& helper) +{ + if ( !isAnalytic() ) return false; + + size_t iFrom = 0, iTo = _eos._edges.size(); + + if ( _anaCurve->IsKind( STANDARD_TYPE( Geom_Line ))) + { + if ( F.IsNull() ) // 3D + { + SMESH_TNodeXYZ pSrc0( _eos._edges[iFrom]->_2neibors->srcNode(0) ); + SMESH_TNodeXYZ pSrc1( _eos._edges[iTo-1]->_2neibors->srcNode(1) ); + //const gp_XYZ lineDir = pSrc1 - pSrc0; + //_LayerEdge* vLE0 = getLEdgeOnV( 0 ); + //_LayerEdge* vLE1 = getLEdgeOnV( 1 ); + // bool shiftOnly = ( vLE0->Is( _LayerEdge::NORMAL_UPDATED ) || + // vLE0->Is( _LayerEdge::BLOCKED ) || + // vLE1->Is( _LayerEdge::NORMAL_UPDATED ) || + // vLE1->Is( _LayerEdge::BLOCKED )); + for ( int iEnd = 0; iEnd < 2; ++iEnd ) + { + iFrom = _eToSmooth[ iEnd ].first, iTo = _eToSmooth[ iEnd ].second; + if ( iFrom >= iTo ) continue; + SMESH_TNodeXYZ p0( _eos[iFrom]->_2neibors->tgtNode(0) ); + SMESH_TNodeXYZ p1( _eos[iTo-1]->_2neibors->tgtNode(1) ); + double param0 = ( iFrom == 0 ) ? 0. : _leParams[ iFrom-1 ]; + double param1 = _leParams[ iTo ]; + for ( size_t i = iFrom; i < iTo; ++i ) + { + _LayerEdge* edge = _eos[i]; + SMDS_MeshNode* tgtNode = const_cast( edge->_nodes.back() ); + double param = ( _leParams[i] - param0 ) / ( param1 - param0 ); + gp_XYZ newPos = p0 * ( 1. - param ) + p1 * param; + + // if ( shiftOnly || edge->Is( _LayerEdge::NORMAL_UPDATED )) + // { + // gp_XYZ curPos = SMESH_TNodeXYZ ( tgtNode ); + // double shift = ( lineDir * ( newPos - pSrc0 ) - + // lineDir * ( curPos - pSrc0 )); + // newPos = curPos + lineDir * shift / lineDir.SquareModulus(); + // } + if ( edge->Is( _LayerEdge::BLOCKED )) + { + SMESH_TNodeXYZ pSrc( edge->_nodes[0] ); + double curThick = pSrc.SquareDistance( tgtNode ); + double newThink = ( pSrc - newPos ).SquareModulus(); + if ( newThink > curThick ) + continue; + } + edge->_pos.back() = newPos; + tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() ); + dumpMove( tgtNode ); + } + } + } + else // 2D + { + _LayerEdge* eV0 = getLEdgeOnV( 0 ); + _LayerEdge* eV1 = getLEdgeOnV( 1 ); + gp_XY uvV0 = eV0->LastUV( F, *data.GetShapeEdges( eV0 )); + gp_XY uvV1 = eV1->LastUV( F, *data.GetShapeEdges( eV1 )); + if ( eV0->_nodes.back() == eV1->_nodes.back() ) // closed edge + { + int iPeriodic = helper.GetPeriodicIndex(); + if ( iPeriodic == 1 || iPeriodic == 2 ) + { + uvV1.SetCoord( iPeriodic, helper.GetOtherParam( uvV1.Coord( iPeriodic ))); + if ( uvV0.Coord( iPeriodic ) > uvV1.Coord( iPeriodic )) + std::swap( uvV0, uvV1 ); + } + } + for ( int iEnd = 0; iEnd < 2; ++iEnd ) + { + iFrom = _eToSmooth[ iEnd ].first, iTo = _eToSmooth[ iEnd ].second; + if ( iFrom >= iTo ) continue; + _LayerEdge* e0 = _eos[iFrom]->_2neibors->_edges[0]; + _LayerEdge* e1 = _eos[iTo-1]->_2neibors->_edges[1]; + gp_XY uv0 = ( e0 == eV0 ) ? uvV0 : e0->LastUV( F, _eos ); + gp_XY uv1 = ( e1 == eV1 ) ? uvV1 : e1->LastUV( F, _eos ); + double param0 = ( iFrom == 0 ) ? 0. : _leParams[ iFrom-1 ]; + double param1 = _leParams[ iTo ]; + gp_XY rangeUV = uv1 - uv0; + for ( size_t i = iFrom; i < iTo; ++i ) + { + if ( _eos[i]->Is( _LayerEdge::BLOCKED )) continue; + double param = ( _leParams[i] - param0 ) / ( param1 - param0 ); + gp_XY newUV = uv0 + param * rangeUV; + + gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() ); + SMDS_MeshNode* tgtNode = const_cast( _eos[i]->_nodes.back() ); + tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() ); + dumpMove( tgtNode ); + + if ( SMDS_FacePositionPtr pos = tgtNode->GetPosition() ) // NULL if F is noShrink + { + pos->SetUParameter( newUV.X() ); + pos->SetVParameter( newUV.Y() ); + } + + gp_XYZ newUV0( newUV.X(), newUV.Y(), 0 ); + + if ( !_eos[i]->Is( _LayerEdge::SMOOTHED )) + { + _eos[i]->Set( _LayerEdge::SMOOTHED ); // to check in refine() (IPAL54237) + if ( _eos[i]->_pos.size() > 2 ) + { + // modify previous positions to make _LayerEdge less sharply bent + vector& uvVec = _eos[i]->_pos; + const gp_XYZ uvShift = newUV0 - uvVec.back(); + const double len2 = ( uvVec.back() - uvVec[ 0 ] ).SquareModulus(); + int iPrev = uvVec.size() - 2; + while ( iPrev > 0 ) + { + double r = ( uvVec[ iPrev ] - uvVec[0] ).SquareModulus() / len2; + uvVec[ iPrev ] += uvShift * r; + --iPrev; + } + } + } + _eos[i]->_pos.back() = newUV0; + } + } + } + return true; } if ( _anaCurve->IsKind( STANDARD_TYPE( Geom_Circle ))) @@ -5536,9 +6260,10 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData& data, if ( uLast < 0 ) uLast += 2 * M_PI; - for ( size_t i = iFrom; i < iTo; ++i ) + for ( size_t i = 0; i < _eos.size(); ++i ) { - if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue; + if ( _eos[i]->Is( _LayerEdge::BLOCKED )) continue; + //if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) continue; double u = uLast * _leParams[i]; gp_Pnt p = ElCLib::Value( u, newCirc ); _eos._edges[i]->_pos.back() = p.XYZ(); @@ -5570,9 +6295,10 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData& data, gp_Ax2d axis( center, vec0 ); gp_Circ2d circ( axis, radius ); - for ( size_t i = iFrom; i < iTo; ++i ) + for ( size_t i = 0; i < _eos.size(); ++i ) { - if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue; + if ( _eos[i]->Is( _LayerEdge::BLOCKED )) continue; + //if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) continue; double newU = uLast * _leParams[i]; gp_Pnt2d newUV = ElCLib::Value( newU, circ ); _eos._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 ); @@ -5582,9 +6308,12 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData& data, tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() ); dumpMove( tgtNode ); - SMDS_FacePosition* pos = static_cast( tgtNode->GetPosition() ); - pos->SetUParameter( newUV.X() ); - pos->SetVParameter( newUV.Y() ); + if ( SMDS_FacePositionPtr pos = tgtNode->GetPosition() ) // NULL if F is noShrink + { + pos->SetUParameter( newUV.X() ); + pos->SetVParameter( newUV.Y() ); + } + _eos[i]->Set( _LayerEdge::SMOOTHED ); // to check in refine() (IPAL54237) } } return true; @@ -5599,15 +6328,17 @@ bool _Smoother1D::smoothAnalyticEdge( _SolidData& data, */ //================================================================================ -bool _Smoother1D::smoothComplexEdge( _SolidData& data, +bool _Smoother1D::smoothComplexEdge( _SolidData& /*data*/, Handle(ShapeAnalysis_Surface)& surface, const TopoDS_Face& F, - SMESH_MesherHelper& helper) + SMESH_MesherHelper& /*helper*/) { if ( _offPoints.empty() ) return false; + // ---------------------------------------------- // move _offPoints along normals of _LayerEdge's + // ---------------------------------------------- _LayerEdge* e[2] = { getLEdgeOnV(0), getLEdgeOnV(1) }; if ( e[0]->Is( _LayerEdge::NORMAL_UPDATED )) @@ -5651,18 +6382,23 @@ bool _Smoother1D::smoothComplexEdge( _SolidData& data, } } + // ----------------------------------------------------------------- // project tgt nodes of extreme _LayerEdge's to the offset segments + // ----------------------------------------------------------------- - if ( e[0]->Is( _LayerEdge::NORMAL_UPDATED )) _iSeg[0] = 0; - if ( e[1]->Is( _LayerEdge::NORMAL_UPDATED )) _iSeg[1] = _offPoints.size()-2; + const int updatedOrBlocked = _LayerEdge::NORMAL_UPDATED | _LayerEdge::BLOCKED; + if ( e[0]->Is( updatedOrBlocked )) _iSeg[0] = 0; + if ( e[1]->Is( updatedOrBlocked )) _iSeg[1] = _offPoints.size()-2; gp_Pnt pExtreme[2], pProj[2]; + bool isProjected[2]; for ( int is2nd = 0; is2nd < 2; ++is2nd ) { pExtreme[ is2nd ] = SMESH_TNodeXYZ( e[is2nd]->_nodes.back() ); int i = _iSeg[ is2nd ]; int di = is2nd ? -1 : +1; - bool projected = false; + bool & projected = isProjected[ is2nd ]; + projected = false; double uOnSeg, distMin = Precision::Infinite(), dist, distPrev = 0; int nbWorse = 0; do { @@ -5709,42 +6445,52 @@ bool _Smoother1D::smoothComplexEdge( _SolidData& data, gp_Vec vDiv0( pExtreme[0], pProj[0] ); gp_Vec vDiv1( pExtreme[1], pProj[1] ); double d0 = vDiv0.Magnitude(); - double d1 = vDiv1.Magnitude(); - if ( e[0]->_normal * vDiv0.XYZ() < 0 ) e[0]->_len += d0; - else e[0]->_len -= d0; - if ( e[1]->_normal * vDiv1.XYZ() < 0 ) e[1]->_len += d1; - else e[1]->_len -= d1; + double d1 = isProjected[1] ? vDiv1.Magnitude() : 0; + if ( e[0]->Is( _LayerEdge::BLOCKED )) { + if ( e[0]->_normal * vDiv0.XYZ() < 0 ) e[0]->_len += d0; + else e[0]->_len -= d0; + } + if ( e[1]->Is( _LayerEdge::BLOCKED )) { + if ( e[1]->_normal * vDiv1.XYZ() < 0 ) e[1]->_len += d1; + else e[1]->_len -= d1; + } + // --------------------------------------------------------------------------------- // compute normalized length of the offset segments located between the projections + // --------------------------------------------------------------------------------- + + // temporary replace extreme _offPoints by pExtreme + gp_XYZ opXYZ[2] = { _offPoints[ _iSeg[0] ]._xyz, + _offPoints[ _iSeg[1]+1 ]._xyz }; + _offPoints[ _iSeg[0] ]._xyz = pExtreme[0].XYZ(); + _offPoints[ _iSeg[1]+ 1]._xyz = pExtreme[1].XYZ(); size_t iSeg = 0, nbSeg = _iSeg[1] - _iSeg[0] + 1; vector< double > len( nbSeg + 1 ); len[ iSeg++ ] = 0; - len[ iSeg++ ] = pProj[ 0 ].Distance( _offPoints[ _iSeg[0]+1 ]._xyz )/* * e[0]->_lenFactor*/; + len[ iSeg++ ] = pProj[ 0 ].Distance( _offPoints[ _iSeg[0]+1 ]._xyz ); for ( size_t i = _iSeg[0]+1; i <= _iSeg[1]; ++i, ++iSeg ) { len[ iSeg ] = len[ iSeg-1 ] + _offPoints[i].Distance( _offPoints[i+1] ); } - len[ nbSeg ] -= pProj[ 1 ].Distance( _offPoints[ _iSeg[1]+1 ]._xyz )/* * e[1]->_lenFactor*/; + // if ( isProjected[ 1 ]) + // len[ nbSeg ] -= pProj[ 1 ].Distance( _offPoints[ _iSeg[1]+1 ]._xyz ); + // else + // len[ nbSeg ] += pExtreme[ 1 ].Distance( _offPoints[ _iSeg[1]+1 ]._xyz ); - // d0 *= e[0]->_lenFactor; - // d1 *= e[1]->_lenFactor; double fullLen = len.back() - d0 - d1; for ( iSeg = 0; iSeg < len.size(); ++iSeg ) len[iSeg] = ( len[iSeg] - d0 ) / fullLen; - // temporary replace extreme _offPoints by pExtreme - gp_XYZ op[2] = { _offPoints[ _iSeg[0] ]._xyz, - _offPoints[ _iSeg[1]+1 ]._xyz }; - _offPoints[ _iSeg[0] ]._xyz = pExtreme[0].XYZ(); - _offPoints[ _iSeg[1]+ 1]._xyz = pExtreme[1].XYZ(); - + // ------------------------------------------------------------- // distribute tgt nodes of _LayerEdge's between the projections + // ------------------------------------------------------------- iSeg = 0; - for ( size_t i = 0; i < _eos._edges.size(); ++i ) + for ( size_t i = 0; i < _eos.size(); ++i ) { - if ( _eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue; + if ( _eos[i]->Is( _LayerEdge::BLOCKED )) continue; + //if ( !_eos[i]->Is( _LayerEdge::TO_SMOOTH )) continue; while ( iSeg+2 < len.size() && _leParams[i] > len[ iSeg+1 ] ) iSeg++; double r = ( _leParams[i] - len[ iSeg ]) / ( len[ iSeg+1 ] - len[ iSeg ]); @@ -5753,23 +6499,23 @@ bool _Smoother1D::smoothComplexEdge( _SolidData& data, if ( surface.IsNull() ) { - _eos._edges[i]->_pos.back() = p; + _eos[i]->_pos.back() = p; } else // project a new node position to a FACE { - gp_Pnt2d uv ( _eos._edges[i]->_pos.back().X(), _eos._edges[i]->_pos.back().Y() ); + gp_Pnt2d uv ( _eos[i]->_pos.back().X(), _eos[i]->_pos.back().Y() ); gp_Pnt2d uv2( surface->NextValueOfUV( uv, p, fTol )); p = surface->Value( uv2 ).XYZ(); - _eos._edges[i]->_pos.back().SetCoord( uv2.X(), uv2.Y(), 0 ); + _eos[i]->_pos.back().SetCoord( uv2.X(), uv2.Y(), 0 ); } - SMDS_MeshNode* tgtNode = const_cast( _eos._edges[i]->_nodes.back() ); + SMDS_MeshNode* tgtNode = const_cast( _eos[i]->_nodes.back() ); tgtNode->setXYZ( p.X(), p.Y(), p.Z() ); dumpMove( tgtNode ); } - _offPoints[ _iSeg[0] ]._xyz = op[0]; - _offPoints[ _iSeg[1]+1 ]._xyz = op[1]; + _offPoints[ _iSeg[0] ]._xyz = opXYZ[0]; + _offPoints[ _iSeg[1]+1 ]._xyz = opXYZ[1]; return true; } @@ -5804,15 +6550,27 @@ void _Smoother1D::prepare(_SolidData& data) double fullLen = _leParams.back() + pPrev.Distance( SMESH_TNodeXYZ( getLEdgeOnV(1)->_nodes[0])); for ( size_t i = 0; i < _leParams.size()-1; ++i ) _leParams[i] = _leParams[i+1] / fullLen; + _leParams.back() = 1.; } + _LayerEdge* leOnV[2] = { getLEdgeOnV(0), getLEdgeOnV(1) }; + + // get cosin to use in findEdgesToSmooth() + _edgeDir[0] = getEdgeDir( E, leOnV[0]->_nodes[0], data.GetHelper() ); + _edgeDir[1] = getEdgeDir( E, leOnV[1]->_nodes[0], data.GetHelper() ); + _leOnV[0]._cosin = Abs( leOnV[0]->_cosin ); + _leOnV[1]._cosin = Abs( leOnV[1]->_cosin ); + if ( _eos._sWOL.IsNull() ) // 3D + for ( int iEnd = 0; iEnd < 2; ++iEnd ) + _leOnV[iEnd]._cosin = Abs( _edgeDir[iEnd].Normalized() * leOnV[iEnd]->_normal ); + if ( isAnalytic() ) return; // divide E to have offset segments with low deflection BRepAdaptor_Curve c3dAdaptor( E ); - const double curDeflect = 0.1; //0.3; // 0.01; // Curvature deflection - const double angDeflect = 0.1; //0.2; // 0.09; // Angular deflection + const double curDeflect = 0.1; //0.01; // Curvature deflection == |p1p2|*sin(p1p2,p1pM) + const double angDeflect = 0.1; //0.09; // Angular deflection == sin(p1pM,pMp2) GCPnts_TangentialDeflection discret(c3dAdaptor, angDeflect, curDeflect); if ( discret.NbPoints() <= 2 ) { @@ -5822,17 +6580,40 @@ void _Smoother1D::prepare(_SolidData& data) const double u0 = c3dAdaptor.FirstParameter(); gp_Pnt p; gp_Vec tangent; - _offPoints.resize( discret.NbPoints() ); - for ( size_t i = 0; i < _offPoints.size(); i++ ) + if ( discret.NbPoints() >= (int) _eos.size() + 2 ) { - double u = discret.Parameter( i+1 ); - c3dAdaptor.D1( u, p, tangent ); - _offPoints[i]._xyz = p.XYZ(); - _offPoints[i]._edgeDir = tangent.XYZ(); - _offPoints[i]._param = GCPnts_AbscissaPoint::Length( c3dAdaptor, u0, u ) / _curveLen; + _offPoints.resize( discret.NbPoints() ); + for ( size_t i = 0; i < _offPoints.size(); i++ ) + { + double u = discret.Parameter( i+1 ); + c3dAdaptor.D1( u, p, tangent ); + _offPoints[i]._xyz = p.XYZ(); + _offPoints[i]._edgeDir = tangent.XYZ(); + _offPoints[i]._param = GCPnts_AbscissaPoint::Length( c3dAdaptor, u0, u ) / _curveLen; + } } + else + { + std::vector< double > params( _eos.size() + 2 ); - _LayerEdge* leOnV[2] = { getLEdgeOnV(0), getLEdgeOnV(1) }; + params[0] = data.GetHelper().GetNodeU( E, leOnV[0]->_nodes[0] ); + params.back() = data.GetHelper().GetNodeU( E, leOnV[1]->_nodes[0] ); + for ( size_t i = 0; i < _eos.size(); i++ ) + params[i+1] = data.GetHelper().GetNodeU( E, _eos[i]->_nodes[0] ); + + if ( params[1] > params[ _eos.size() ] ) + std::reverse( params.begin() + 1, params.end() - 1 ); + + _offPoints.resize( _eos.size() + 2 ); + for ( size_t i = 0; i < _offPoints.size(); i++ ) + { + const double u = params[i]; + c3dAdaptor.D1( u, p, tangent ); + _offPoints[i]._xyz = p.XYZ(); + _offPoints[i]._edgeDir = tangent.XYZ(); + _offPoints[i]._param = GCPnts_AbscissaPoint::Length( c3dAdaptor, u0, u ) / _curveLen; + } + } // set _2edges _offPoints [0]._2edges.set( &_leOnV[0], &_leOnV[0], 0.5, 0.5 ); @@ -5871,11 +6652,14 @@ void _Smoother1D::prepare(_SolidData& data) int iLBO = _offPoints.size() - 2; // last but one - _edgeDir[0] = getEdgeDir( E, leOnV[0]->_nodes[0], data.GetHelper() ); - _edgeDir[1] = getEdgeDir( E, leOnV[1]->_nodes[0], data.GetHelper() ); - - _leOnV[ 0 ]._normal = getNormalNormal( leOnV[0]->_normal, _edgeDir[0] ); - _leOnV[ 1 ]._normal = getNormalNormal( leOnV[1]->_normal, _edgeDir[1] ); + if ( leOnV[ 0 ]->Is( _LayerEdge::MULTI_NORMAL )) + _leOnV[ 0 ]._normal = getNormalNormal( _eos._edges[1]->_normal, _edgeDir[0] ); + else + _leOnV[ 0 ]._normal = getNormalNormal( leOnV[0]->_normal, _edgeDir[0] ); + if ( leOnV[ 1 ]->Is( _LayerEdge::MULTI_NORMAL )) + _leOnV[ 1 ]._normal = getNormalNormal( _eos._edges.back()->_normal, _edgeDir[1] ); + else + _leOnV[ 1 ]._normal = getNormalNormal( leOnV[1]->_normal, _edgeDir[1] ); _leOnV[ 0 ]._len = 0; _leOnV[ 1 ]._len = 0; _leOnV[ 0 ]._lenFactor = _offPoints[1 ]._2edges._edges[1]->_lenFactor; @@ -5909,7 +6693,7 @@ void _Smoother1D::prepare(_SolidData& data) //================================================================================ /*! - * \brief set _normal of _leOnV[is2nd] to be normal to the EDGE + * \brief return _normal of _leOnV[is2nd] normal to the EDGE */ //================================================================================ @@ -5920,9 +6704,46 @@ gp_XYZ _Smoother1D::getNormalNormal( const gp_XYZ & normal, gp_XYZ norm = edgeDir ^ cross; double size = norm.Modulus(); + // if ( size == 0 ) // MULTI_NORMAL _LayerEdge + // return gp_XYZ( 1e-100, 1e-100, 1e-100 ); + + if ( size < 1e-5 ) // normal || edgeDir (almost) at inflation along EDGE (bos #20643) + { + const _LayerEdge* le = _eos._edges[ _eos._edges.size() / 2 ]; + const gp_XYZ& leNorm = le->_normal; + + cross = leNorm ^ edgeDir; + norm = edgeDir ^ cross; + size = norm.Modulus(); + } + return norm / size; } +//================================================================================ +/*! + * \brief Writes a script creating a mesh composed of _offPoints + */ +//================================================================================ + +void _Smoother1D::offPointsToPython() const +{ + const char* fname = "/tmp/offPoints.py"; + cout << "exec(open('"<_edges.size(); ++i ) - { - eos->_edges[i]->SetSmooLen( Precision::Infinite() ); - } - SMESH_subMeshIteratorPtr smIt = eos->_subMesh->getDependsOnIterator(/*includeSelf=*/false); - while ( smIt->more() ) // loop on sub-shapes of the FACE - { - _EdgesOnShape* eoe = GetShapeEdges( smIt->next()->GetId() ); - if ( !eoe ) continue; - - vector<_LayerEdge*>& eE = eoe->_edges; - for ( size_t iE = 0; iE < eE.size(); ++iE ) // loop on _LayerEdge's on EDGE or VERTEX - { - if ( eE[iE]->_cosin <= theMinSmoothCosin ) - continue; + // for ( size_t i = 0; i < eos->_edges.size(); ++i ) + // { + // eos->_edges[i]->SetSmooLen( Precision::Infinite() ); + // } + // SMESH_subMeshIteratorPtr smIt = eos->_subMesh->getDependsOnIterator(/*includeSelf=*/false); + // while ( smIt->more() ) // loop on sub-shapes of the FACE + // { + // _EdgesOnShape* eoe = GetShapeEdges( smIt->next()->GetId() ); + // if ( !eoe ) continue; - SMDS_ElemIteratorPtr segIt = eE[iE]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Edge); - while ( segIt->more() ) - { - const SMDS_MeshElement* seg = segIt->next(); - if ( !eos->_subMesh->DependsOn( seg->getshapeId() )) - continue; - if ( seg->GetNode(0) != eE[iE]->_nodes[0] ) - continue; // not to check a seg twice - for ( size_t iN = 0; iN < eE[iE]->_neibors.size(); ++iN ) - { - _LayerEdge* eN = eE[iE]->_neibors[iN]; - if ( eN->_nodes[0]->getshapeId() != eos->_shapeID ) - continue; - double dist = SMESH_MeshAlgos::GetDistance( seg, SMESH_TNodeXYZ( eN->_nodes[0] )); - double smooLen = getSmoothingThickness( eE[iE]->_cosin, dist ); - eN->SetSmooLen( Min( smooLen, eN->GetSmooLen() )); - eN->Set( _LayerEdge::NEAR_BOUNDARY ); - } - } - } - } + // vector<_LayerEdge*>& eE = eoe->_edges; + // for ( size_t iE = 0; iE < eE.size(); ++iE ) // loop on _LayerEdge's on EDGE or VERTEX + // { + // if ( eE[iE]->_cosin <= theMinSmoothCosin ) + // continue; + + // SMDS_ElemIteratorPtr segIt = eE[iE]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Edge); + // while ( segIt->more() ) + // { + // const SMDS_MeshElement* seg = segIt->next(); + // if ( !eos->_subMesh->DependsOn( seg->getshapeId() )) + // continue; + // if ( seg->GetNode(0) != eE[iE]->_nodes[0] ) + // continue; // not to check a seg twice + // for ( size_t iN = 0; iN < eE[iE]->_neibors.size(); ++iN ) + // { + // _LayerEdge* eN = eE[iE]->_neibors[iN]; + // if ( eN->_nodes[0]->getshapeId() != eos->_shapeID ) + // continue; + // double dist = SMESH_MeshAlgos::GetDistance( seg, SMESH_TNodeXYZ( eN->_nodes[0] )); + // double smooLen = getSmoothingThickness( eE[iE]->_cosin, dist ); + // eN->SetSmooLen( Min( smooLen, eN->GetSmooLen() )); + // eN->Set( _LayerEdge::NEAR_BOUNDARY ); + // } + // } + // } + // } } // if ( eos->ShapeType() == TopAbs_FACE ) for ( size_t i = 0; i < eos->_edges.size(); ++i ) @@ -6104,11 +6925,12 @@ void _SolidData::PrepareEdgesToSmoothOnFace( _EdgesOnShape* eos, bool substitute avgLen /= edge->_simplices.size(); if (( edge->_curvature = _Curvature::New( avgNormProj, avgLen ))) { + edge->Set( _LayerEdge::SMOOTHED_C1 ); isCurved = true; - SMDS_FacePosition* fPos = dynamic_cast( edge->_nodes[0]->GetPosition() ); + SMDS_FacePositionPtr fPos = edge->_nodes[0]->GetPosition(); if ( !fPos ) for ( size_t iS = 0; iS < edge->_simplices.size() && !fPos; ++iS ) - fPos = dynamic_cast( edge->_simplices[iS]._nPrev->GetPosition() ); + fPos = edge->_simplices[iS]._nPrev->GetPosition(); if ( fPos ) edge->_curvature->_uv.SetCoord( fPos->GetUParameter(), fPos->GetVParameter() ); } @@ -6122,20 +6944,98 @@ void _SolidData::PrepareEdgesToSmoothOnFace( _EdgesOnShape* eos, bool substitute eos->_edgeForOffset = 0; double maxCosin = -1; + //bool hasNoShrink = false; for ( TopExp_Explorer eExp( eos->_shape, TopAbs_EDGE ); eExp.More(); eExp.Next() ) { _EdgesOnShape* eoe = GetShapeEdges( eExp.Current() ); if ( !eoe || eoe->_edges.empty() ) continue; + // if ( eos->GetData()._noShrinkShapes.count( eoe->_shapeID )) + // hasNoShrink = true; + vector<_LayerEdge*>& eE = eoe->_edges; _LayerEdge* e = eE[ eE.size() / 2 ]; - if ( e->_cosin > maxCosin ) + if ( !e->Is( _LayerEdge::RISKY_SWOL ) && e->_cosin > maxCosin ) { eos->_edgeForOffset = e; maxCosin = e->_cosin; } + + if ( !eoe->_sWOL.IsNull() ) + for ( _LayerEdge* le : eoe->_edges ) + if ( le->Is( _LayerEdge::RISKY_SWOL ) && e->_cosin > 0 ) + { + // make _neibors on FACE be smoothed after le->Is( BLOCKED ) + for ( _LayerEdge* neibor : le->_neibors ) + { + int shapeDim = neibor->BaseShapeDim(); + if ( shapeDim == 2 ) + neibor->Set( _LayerEdge::NEAR_BOUNDARY ); // on FACE + else if ( shapeDim == 0 ) + neibor->Set( _LayerEdge::RISKY_SWOL ); // on VERTEX + + if ( !neibor->_curvature ) + { + gp_XY uv = helper.GetNodeUV( F, neibor->_nodes[0] ); + neibor->_curvature = _Factory::NewCurvature(); + neibor->_curvature->_r = 0; + neibor->_curvature->_k = 0; + neibor->_curvature->_h2lenRatio = 0; + neibor->_curvature->_uv = uv; + } + } + } + } // loop on EDGEs + + // Try to initialize _Mapper2D + + // if ( hasNoShrink ) + // return; + + SMDS_ElemIteratorPtr fIt = eos->_subMesh->GetSubMeshDS()->GetElements(); + if ( !fIt->more() || fIt->next()->NbCornerNodes() != 4 ) + return; + + // get EDGEs of quadrangle bottom + std::list< TopoDS_Edge > edges; + std::list< int > nbEdgesInWire; + int nbWire = SMESH_Block::GetOrderedEdges( F, edges, nbEdgesInWire ); + if ( nbWire != 1 || nbEdgesInWire.front() < 4 ) + return; + const SMDS_MeshNode* node; + while ( true ) // make edges start at a corner VERTEX + { + node = SMESH_Algo::VertexNode( helper.IthVertex( 0, edges.front() ), helper.GetMeshDS() ); + if ( node && helper.IsCornerOfStructure( node, eos->_subMesh->GetSubMeshDS(), helper )) + break; + edges.pop_front(); + if ( edges.empty() ) + return; + } + std::list< TopoDS_Edge >::iterator edgeIt = edges.begin(); + while ( true ) // make edges finish at a corner VERTEX + { + node = SMESH_Algo::VertexNode( helper.IthVertex( 1, *edgeIt ), helper.GetMeshDS() ); + ++edgeIt; + if ( node && helper.IsCornerOfStructure( node, eos->_subMesh->GetSubMeshDS(), helper )) + { + edges.erase( edgeIt, edges.end() ); + break; + } + if ( edgeIt == edges.end() ) + return; } - } + + // get structure of nodes + TParam2ColumnMap param2ColumnMap; + if ( !helper.LoadNodeColumns( param2ColumnMap, F, edges, helper.GetMeshDS() )) + return; + + eos->_mapper2D = new _Mapper2D( param2ColumnMap, eos->GetData()._n2eMap ); + + } // if eos is of curved FACE + + return; } //================================================================================ @@ -6176,7 +7076,7 @@ void _SolidData::AddShapesToSmooth( const set< _EdgesOnShape* >& eosToSmooth, */ //================================================================================ -void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelper& helper ) +void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelper& /*helper*/ ) { // find intersection of neighbor _LayerEdge's to limit _maxLen // according to local curvature (IPAL52648) @@ -6199,7 +7099,7 @@ void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelp if ( eI->_nodes[0]->GetID() < eN->_nodes[0]->GetID() ) // treat this pair once { _EdgesOnShape* eosN = data.GetShapeEdges( eN ); - limitMaxLenByCurvature( eI, eN, eosI, *eosN, helper ); + limitMaxLenByCurvature( eI, eN, eosI, *eosN, eosI._hyp.ToSmooth() ); } } } @@ -6213,7 +7113,7 @@ void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelp for ( size_t i = 1; i < eosI._edges.size(); ++i ) { _LayerEdge* eI = eosI._edges[i]; - limitMaxLenByCurvature( eI, e0, eosI, eosI, helper ); + limitMaxLenByCurvature( eI, e0, eosI, eosI, eosI._hyp.ToSmooth() ); e0 = eI; } } @@ -6226,12 +7126,17 @@ void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelp */ //================================================================================ -void _ViscousBuilder::limitMaxLenByCurvature( _LayerEdge* e1, - _LayerEdge* e2, - _EdgesOnShape& eos1, - _EdgesOnShape& eos2, - SMESH_MesherHelper& helper ) +void _ViscousBuilder::limitMaxLenByCurvature( _LayerEdge* e1, + _LayerEdge* e2, + _EdgesOnShape& /*eos1*/, + _EdgesOnShape& /*eos2*/, + const bool /*isSmoothable*/ ) { + if (( e1->_nodes[0]->GetPosition()->GetDim() != + e2->_nodes[0]->GetPosition()->GetDim() ) && + ( e1->_cosin < 0.75 )) + return; // angle > 90 deg at e1 + gp_XYZ plnNorm = e1->_normal ^ e2->_normal; double norSize = plnNorm.SquareModulus(); if ( norSize < std::numeric_limits::min() ) @@ -6251,9 +7156,10 @@ void _ViscousBuilder::limitMaxLenByCurvature( _LayerEdge* e1, double ovl = ( u1 * e1->_normal * dir12 - u2 * e2->_normal * dir12 ) / dir12.SquareModulus(); if ( ovl > theSmoothThickToElemSizeRatio ) - { - e1->_maxLen = Min( e1->_maxLen, 0.75 * u1 / e1->_lenFactor ); - e2->_maxLen = Min( e2->_maxLen, 0.75 * u2 / e2->_lenFactor ); + { + const double coef = 0.75; + e1->SetMaxLen( Min( e1->_maxLen, coef * u1 / e1->_lenFactor )); + e2->SetMaxLen( Min( e2->_maxLen, coef * u2 / e2->_lenFactor )); } } } @@ -6274,6 +7180,7 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper& _EdgesOnShape& eos = data._edgesOnShape[iS]; if ( eos._edges.empty() ) continue; if ( eos.ShapeType() != TopAbs_EDGE && eos.ShapeType() != TopAbs_VERTEX ) continue; + if ( !eos._sWOL.IsNull() ) continue; // PAL23566 for ( size_t i = 0; i < eos._edges.size(); ++i ) { @@ -6320,7 +7227,7 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper& src2->GetID() < edge->_nodes[0]->GetID() ) continue; // avoid using same segment twice - // a _LayerEdge containg tgt2 + // a _LayerEdge containing tgt2 _LayerEdge* neiborEdge = edge->_2neibors->_edges[j]; _TmpMeshFaceOnEdge* f = new _TmpMeshFaceOnEdge( edge, neiborEdge, --_tmpFaceID ); @@ -6399,9 +7306,9 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper& ( f->_le2->IsOnEdge() && f->_le2->_2neibors->include( edge ))) continue; } dist1 = dist2 = Precision::Infinite(); - if ( !edge->SegTriaInter( lastSegment, f->_nn[0], f->_nn[1], f->_nn[2], dist1, eps )) + if ( !edge->SegTriaInter( lastSegment, f->n(0), f->n(1), f->n(2), dist1, eps )) dist1 = Precision::Infinite(); - if ( !edge->SegTriaInter( lastSegment, f->_nn[3], f->_nn[2], f->_nn[0], dist2, eps )) + if ( !edge->SegTriaInter( lastSegment, f->n(3), f->n(2), f->n(0), dist2, eps )) dist2 = Precision::Infinite(); if (( dist1 > segLen ) && ( dist2 > segLen )) continue; @@ -6409,7 +7316,7 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper& if ( edge->IsOnEdge() ) { // skip perpendicular EDGEs - gp_Vec fSegDir = SMESH_TNodeXYZ( f->_nn[0] ) - SMESH_TNodeXYZ( f->_nn[3] ); + gp_Vec fSegDir = SMESH_TNodeXYZ( f->n(0) ) - SMESH_TNodeXYZ( f->n(3) ); bool isParallel = ( isLessAngle( eSegDir0, fSegDir, angle45 ) || isLessAngle( eSegDir1, fSegDir, angle45 ) || isLessAngle( eSegDir0, fSegDir.Reversed(), angle45 ) || @@ -6428,7 +7335,7 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper& // else // { // double shortLen = 0.75 * ( Min( dist1, dist2 ) / edge->_lenFactor ); - // edge->_maxLen = Min( shortLen, edge->_maxLen ); + // edge->SetMaxLen( Min( shortLen, edge->_maxLen )); // } } @@ -6458,6 +7365,80 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper& } } +//================================================================================ +/*! + * \brief Find _LayerEdge's located on boundary of a convex FACE whose normal + * will be updated at each inflation step + */ +//================================================================================ + +void _ViscousBuilder::findEdgesToUpdateNormalNearConvexFace( _ConvexFace & convFace, + _SolidData& data, + SMESH_MesherHelper& helper ) +{ + const TGeomID convFaceID = getMeshDS()->ShapeToIndex( convFace._face ); + const double preci = BRep_Tool::Tolerance( convFace._face ); + Handle(ShapeAnalysis_Surface) surface = helper.GetSurface( convFace._face ); + + bool edgesToUpdateFound = false; + + map< TGeomID, _EdgesOnShape* >::iterator id2eos = convFace._subIdToEOS.begin(); + for ( ; id2eos != convFace._subIdToEOS.end(); ++id2eos ) + { + _EdgesOnShape& eos = * id2eos->second; + if ( !eos._sWOL.IsNull() ) continue; + if ( !eos._hyp.ToSmooth() ) continue; + for ( size_t i = 0; i < eos._edges.size(); ++i ) + { + _LayerEdge* ledge = eos._edges[ i ]; + if ( ledge->Is( _LayerEdge::UPD_NORMAL_CONV )) continue; // already checked + if ( ledge->Is( _LayerEdge::MULTI_NORMAL )) continue; // not inflatable + + gp_XYZ tgtPos = ( SMESH_NodeXYZ( ledge->_nodes[0] ) + + ledge->_normal * ledge->_lenFactor * ledge->_maxLen ); + + // the normal must be updated if distance from tgtPos to surface is less than + // target thickness + + // find an initial UV for search of a projection of tgtPos to surface + const SMDS_MeshNode* nodeInFace = 0; + SMDS_ElemIteratorPtr fIt = ledge->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face); + while ( fIt->more() && !nodeInFace ) + { + const SMDS_MeshElement* f = fIt->next(); + if ( convFaceID != f->getshapeId() ) continue; + + SMDS_ElemIteratorPtr nIt = f->nodesIterator(); + while ( nIt->more() && !nodeInFace ) + { + const SMDS_MeshElement* n = nIt->next(); + if ( n->getshapeId() == convFaceID ) + nodeInFace = static_cast< const SMDS_MeshNode* >( n ); + } + } + if ( !nodeInFace ) + continue; + gp_XY uv = helper.GetNodeUV( convFace._face, nodeInFace ); + + // projection + surface->NextValueOfUV( uv, tgtPos, preci ); + double dist = surface->Gap(); + if ( dist < 0.95 * ledge->_maxLen ) + { + ledge->Set( _LayerEdge::UPD_NORMAL_CONV ); + if ( !ledge->_curvature ) ledge->_curvature = _Factory::NewCurvature(); + ledge->_curvature->_uv.SetCoord( uv.X(), uv.Y() ); + edgesToUpdateFound = true; + } + } + } + + if ( !convFace._isTooCurved && edgesToUpdateFound ) + { + data._convexFaces.insert( make_pair( convFaceID, convFace )).first->second; + } +} + //================================================================================ /*! * \brief Modify normals of _LayerEdge's on EDGE's to avoid intersection with @@ -6468,7 +7449,7 @@ void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper& bool _ViscousBuilder::updateNormals( _SolidData& data, SMESH_MesherHelper& helper, int stepNb, - double stepSize) + double /*stepSize*/) { updateNormalsOfC1Vertices( data ); @@ -6565,8 +7546,8 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, // compute new _normals for ( size_t i = 0; i < intEdgesDist.size(); ++i ) { - _LayerEdge* edge2 = intEdgesDist[i].first; - double distWgt = edge1->_len / intEdgesDist[i].second; + _LayerEdge* edge2 = intEdgesDist[i].first; + double distWgt = edge1->_len / intEdgesDist[i].second; // if ( edge1->Is( _LayerEdge::BLOCKED ) && // edge2->Is( _LayerEdge::BLOCKED )) continue; if ( edge2->Is( _LayerEdge::MARKED )) continue; @@ -6604,12 +7585,17 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, e2neIt = edge2newEdge.insert( make_pair( edge1, zeroEdge )).first; e2neIt->second._normal += distWgt * newNormal; e2neIt->second._cosin = newCos; - e2neIt->second._maxLen = 0.7 * minIntDist / edge1->_lenFactor; + e2neIt->second.SetMaxLen( 0.7 * minIntDist / edge1->_lenFactor ); if ( iter > 0 && sgn1 * sgn2 < 0 && edge1->_cosin < 0 ) e2neIt->second._normal += dir2; + e2neIt = edge2newEdge.insert( make_pair( edge2, zeroEdge )).first; e2neIt->second._normal += distWgt * newNormal; - e2neIt->second._cosin = edge2->_cosin; + if ( Precision::IsInfinite( zeroEdge._maxLen )) + { + e2neIt->second._cosin = edge2->_cosin; + e2neIt->second.SetMaxLen( 1.3 * minIntDist / edge1->_lenFactor ); + } if ( iter > 0 && sgn1 * sgn2 < 0 && edge2->_cosin < 0 ) e2neIt->second._normal += dir1; } @@ -6626,9 +7612,10 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, for ( e2neIt = edge2newEdge.begin(); e2neIt != edge2newEdge.end(); ++e2neIt ) { _LayerEdge* edge = e2neIt->first; - if ( edge->Is( _LayerEdge::BLOCKED )) continue; _LayerEdge& newEdge = e2neIt->second; _EdgesOnShape* eos = data.GetShapeEdges( edge ); + if ( edge->Is( _LayerEdge::BLOCKED ) && newEdge._maxLen > edge->_len ) + continue; // Check if a new _normal is OK: newEdge._normal.Normalize(); @@ -6637,7 +7624,7 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, if ( newEdge._maxLen < edge->_len && iter > 0 ) // limit _maxLen { edge->InvalidateStep( stepNb + 1, *eos, /*restoreLength=*/true ); - edge->_maxLen = newEdge._maxLen; + edge->SetMaxLen( newEdge._maxLen ); edge->SetNewLength( newEdge._maxLen, *eos, helper ); } continue; // the new _normal is bad @@ -6657,12 +7644,11 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, else // edge inflates along a FACE { TopoDS_Shape V = helper.GetSubShapeByNode( edge->_nodes[0], getMeshDS() ); - PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE ); + PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE, &eos->_sWOL ); while ( const TopoDS_Shape* E = eIt->next() ) { - if ( !helper.IsSubShape( *E, /*FACE=*/eos->_sWOL )) - continue; - gp_Vec edgeDir = getEdgeDir( TopoDS::Edge( *E ), TopoDS::Vertex( V )); + gp_Vec edgeDir = getEdgeDir( TopoDS::Edge( *E ), TopoDS::Vertex( V ), + eos->_hyp.Get1stLayerThickness() ); double angle = edgeDir.Angle( newEdge._normal ); // [0,PI] if ( angle < M_PI / 2 ) shapesToSmooth.insert( data.GetShapeEdges( *E )); @@ -6804,7 +7790,7 @@ bool _ViscousBuilder::isNewNormalOk( _SolidData& data, //================================================================================ bool _ViscousBuilder::updateNormalsOfSmoothed( _SolidData& data, - SMESH_MesherHelper& helper, + SMESH_MesherHelper& /*helper*/, const int nbSteps, const double stepSize ) { @@ -6912,6 +7898,8 @@ bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData& data, for ( ; id2face != data._convexFaces.end(); ++id2face ) { _ConvexFace & convFace = (*id2face).second; + convFace._normalsFixedOnBorders = false; // to update at each inflation step + if ( convFace._normalsFixed ) continue; // already fixed if ( convFace.CheckPrisms() ) @@ -7274,6 +8262,59 @@ bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData& data, return true; } +//================================================================================ +/*! + * \brief Return max curvature of a FACE + */ +//================================================================================ + +double _ConvexFace::GetMaxCurvature( _SolidData& data, + _EdgesOnShape& eof, + BRepLProp_SLProps& surfProp, + SMESH_MesherHelper& helper) +{ + double maxCurvature = 0; + + TopoDS_Face F = TopoDS::Face( eof._shape ); + + const int nbTestPnt = 5; + const double oriFactor = ( F.Orientation() == TopAbs_REVERSED ? +1. : -1. ); + SMESH_subMeshIteratorPtr smIt = eof._subMesh->getDependsOnIterator(/*includeSelf=*/true); + while ( smIt->more() ) + { + SMESH_subMesh* sm = smIt->next(); + const TGeomID subID = sm->GetId(); + + // find _LayerEdge's of a sub-shape + _EdgesOnShape* eos; + if (( eos = data.GetShapeEdges( subID ))) + this->_subIdToEOS.insert( make_pair( subID, eos )); + else + continue; + + // check concavity and curvature and limit data._stepSize + const double minCurvature = + 1. / ( eos->_hyp.GetTotalThickness() * ( 1 + theThickToIntersection )); + size_t iStep = Max( 1, eos->_edges.size() / nbTestPnt ); + for ( size_t i = 0; i < eos->_edges.size(); i += iStep ) + { + gp_XY uv = helper.GetNodeUV( F, eos->_edges[ i ]->_nodes[0] ); + surfProp.SetParameters( uv.X(), uv.Y() ); + if ( surfProp.IsCurvatureDefined() ) + { + double curvature = Max( surfProp.MaxCurvature() * oriFactor, + surfProp.MinCurvature() * oriFactor ); + maxCurvature = Max( maxCurvature, curvature ); + + if ( curvature > minCurvature ) + this->_isTooCurved = true; + } + } + } // loop on sub-shapes of the FACE + + return maxCurvature; +} + //================================================================================ /*! * \brief Finds a center of curvature of a surface at a _LayerEdge @@ -7482,7 +8523,7 @@ bool _LayerEdge::FindIntersection( SMESH_ElementSearcher& searcher, gp_XYZ _LayerEdge::PrevCheckPos( _EdgesOnShape* eos ) const { - size_t i = Is( NORMAL_UPDATED ) ? _pos.size()-2 : 0; + size_t i = Is( NORMAL_UPDATED ) && IsOnFace() ? _pos.size()-2 : 0; if ( !eos || eos->_sWOL.IsNull() ) return _pos[ i ]; @@ -7556,13 +8597,14 @@ gp_Ax1 _LayerEdge::LastSegment(double& segLen, _EdgesOnShape& eos) const //================================================================================ /*! - * \brief Return the last position of the target node on a FACE. + * \brief Return the last (or \a which) position of the target node on a FACE. * \param [in] F - the FACE this _LayerEdge is inflated along + * \param [in] which - index of position * \return gp_XY - result UV */ //================================================================================ -gp_XY _LayerEdge::LastUV( const TopoDS_Face& F, _EdgesOnShape& eos ) const +gp_XY _LayerEdge::LastUV( const TopoDS_Face& F, _EdgesOnShape& eos, int which ) const { if ( F.IsSame( eos._sWOL )) // F is my FACE return gp_XY( _pos.back().X(), _pos.back().Y() ); @@ -7571,7 +8613,7 @@ gp_XY _LayerEdge::LastUV( const TopoDS_Face& F, _EdgesOnShape& eos ) const return gp_XY( 1e100, 1e100 ); // _sWOL is EDGE of F; _pos.back().X() is the last U on the EDGE - double f, l, u = _pos.back().X(); + double f, l, u = _pos[ which < 0 ? _pos.size()-1 : which ].X(); Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( TopoDS::Edge(eos._sWOL), F, f,l); if ( !C2d.IsNull() && f <= u && u <= l ) return C2d->Value( u ).XY(); @@ -7649,7 +8691,7 @@ bool _LayerEdge::SegTriaInter( const gp_Ax1& lastSegment, * \param [in] eov - EOS of the VERTEX * \param [in] eos - EOS of the FACE * \param [in] step - inflation step - * \param [in,out] badSmooEdges - not untangled _LayerEdge's + * \param [in,out] badSmooEdges - tangled _LayerEdge's */ //================================================================================ @@ -7711,7 +8753,7 @@ void _LayerEdge::MoveNearConcaVer( const _EdgesOnShape* eov, prevPosV = surface.Value( prevPosV.X(), prevPosV.Y() ).XYZ(); } - SMDS_FacePosition* fPos; + SMDS_FacePositionPtr fPos; //double r = 1. - Min( 0.9, step / 10. ); for ( set< _LayerEdge* >::iterator e = edges.begin(); e != edges.end(); ++e ) { @@ -7725,9 +8767,9 @@ void _LayerEdge::MoveNearConcaVer( const _EdgesOnShape* eov, // set _curvature to make edgeF updated by putOnOffsetSurface() if ( !edgeF->_curvature ) - if (( fPos = dynamic_cast( edgeF->_nodes[0]->GetPosition() ))) + if (( fPos = edgeF->_nodes[0]->GetPosition() )) { - edgeF->_curvature = new _Curvature; + edgeF->_curvature = _Factory::NewCurvature(); edgeF->_curvature->_r = 0; edgeF->_curvature->_k = 0; edgeF->_curvature->_h2lenRatio = 0; @@ -8200,13 +9242,13 @@ int _LayerEdge::Smooth(const int step, const bool isConcaveFace, bool findBest ) //================================================================================ /*! - * \brief Chooses a smoothing technic giving a position most close to an initial one. + * \brief Chooses a smoothing technique giving a position most close to an initial one. * For a correct result, _simplices must contain nodes lying on geometry. */ //================================================================================ void _LayerEdge::ChooseSmooFunction( const set< TGeomID >& concaveVertices, - const TNode2Edge& n2eMap) + const TNode2Edge& /*n2eMap*/) { if ( _smooFunction ) return; @@ -8235,7 +9277,7 @@ void _LayerEdge::ChooseSmooFunction( const set< TGeomID >& concaveVertices, } } - // // this coice is done only if ( !concaveVertices.empty() ) for Grids/smesh/bugs_19/X1 + // // this choice is done only if ( !concaveVertices.empty() ) for Grids/smesh/bugs_19/X1 // // where the nodes are smoothed too far along a sphere thus creating // // inverted _simplices // double dist[theNbSmooFuns]; @@ -8399,7 +9441,7 @@ gp_XYZ _LayerEdge::smoothAngular() else norm += cross; } - catch (Standard_Failure) { // if |cross| == 0. + catch (Standard_Failure&) { // if |cross| == 0. } } gp_XYZ vec = newPos - pN; @@ -8411,7 +9453,7 @@ gp_XYZ _LayerEdge::smoothAngular() //================================================================================ /*! - * \brief Computes a new node position using weigthed node positions + * \brief Computes a new node position using weighted node positions */ //================================================================================ @@ -8485,7 +9527,7 @@ gp_XYZ _LayerEdge::smoothNefPolygon() { gp_XYZ newPos(0,0,0); - // get a plane to seach a solution on + // get a plane to search a solution on vector< gp_XYZ > vecs( _simplices.size() + 1 ); size_t i; @@ -8676,7 +9718,7 @@ gp_XYZ _LayerEdge::smoothNefPolygon() { ////////////////////////////////// NEW gp_XYZ newPos(0,0,0); - // get a plane to seach a solution on + // get a plane to search a solution on size_t i; gp_XYZ center(0,0,0); @@ -8892,6 +9934,8 @@ void _LayerEdge::SetNewLength( double len, _EdgesOnShape& eos, SMESH_MesherHelpe Block( eos.GetData() ); } const double lenDelta = len - _len; + // if ( lenDelta < 0 ) + // return; if ( lenDelta < len * 1e-3 ) { Block( eos.GetData() ); @@ -8935,46 +9979,13 @@ void _LayerEdge::SetNewLength( double len, _EdgesOnShape& eos, SMESH_MesherHelpe _pos.push_back( newXYZ ); if ( !eos._sWOL.IsNull() ) - { - double distXYZ[4]; - bool uvOK = false; - if ( eos.SWOLType() == TopAbs_EDGE ) - { - double u = Precision::Infinite(); // to force projection w/o distance check - uvOK = helper.CheckNodeU( TopoDS::Edge( eos._sWOL ), n, u, - /*tol=*/2*lenDelta, /*force=*/true, distXYZ ); - _pos.back().SetCoord( u, 0, 0 ); - if ( _nodes.size() > 1 && uvOK ) - { - SMDS_EdgePosition* pos = static_cast( n->GetPosition() ); - pos->SetUParameter( u ); - } - } - else // TopAbs_FACE - { - gp_XY uv( Precision::Infinite(), 0 ); - uvOK = helper.CheckNodeUV( TopoDS::Face( eos._sWOL ), n, uv, - /*tol=*/2*lenDelta, /*force=*/true, distXYZ ); - _pos.back().SetCoord( uv.X(), uv.Y(), 0 ); - if ( _nodes.size() > 1 && uvOK ) - { - SMDS_FacePosition* pos = static_cast( n->GetPosition() ); - pos->SetUParameter( uv.X() ); - pos->SetVParameter( uv.Y() ); - } - } - if ( uvOK ) - { - n->setXYZ( distXYZ[1], distXYZ[2], distXYZ[3]); - } - else + if ( !UpdatePositionOnSWOL( n, 2*lenDelta, eos, helper )) { n->setXYZ( oldXYZ.X(), oldXYZ.Y(), oldXYZ.Z() ); _pos.pop_back(); Block( eos.GetData() ); return; } - } _len = len; @@ -8983,13 +9994,57 @@ void _LayerEdge::SetNewLength( double len, _EdgesOnShape& eos, SMESH_MesherHelpe { for ( size_t i = 0; i < _neibors.size(); ++i ) //if ( _len > _neibors[i]->GetSmooLen() ) - _neibors[i]->Set( MOVED ); + _neibors[i]->Set( MOVED ); Set( MOVED ); } dumpMove( n ); //debug } + +//================================================================================ +/*! + * \brief Update last position on SWOL by projecting node on SWOL +*/ +//================================================================================ + +bool _LayerEdge::UpdatePositionOnSWOL( SMDS_MeshNode* n, + double tol, + _EdgesOnShape& eos, + SMESH_MesherHelper& helper ) +{ + double distXYZ[4]; + bool uvOK = false; + if ( eos.SWOLType() == TopAbs_EDGE ) + { + double u = Precision::Infinite(); // to force projection w/o distance check + uvOK = helper.CheckNodeU( TopoDS::Edge( eos._sWOL ), n, u, tol, /*force=*/true, distXYZ ); + _pos.back().SetCoord( u, 0, 0 ); + if ( _nodes.size() > 1 && uvOK ) + { + SMDS_EdgePositionPtr pos = n->GetPosition(); + pos->SetUParameter( u ); + } + } + else // TopAbs_FACE + { + gp_XY uv( Precision::Infinite(), 0 ); + uvOK = helper.CheckNodeUV( TopoDS::Face( eos._sWOL ), n, uv, tol, /*force=*/true, distXYZ ); + _pos.back().SetCoord( uv.X(), uv.Y(), 0 ); + if ( _nodes.size() > 1 && uvOK ) + { + SMDS_FacePositionPtr pos = n->GetPosition(); + pos->SetUParameter( uv.X() ); + pos->SetVParameter( uv.Y() ); + } + } + if ( uvOK ) + { + n->setXYZ( distXYZ[1], distXYZ[2], distXYZ[3]); + } + return uvOK; +} + //================================================================================ /*! * \brief Set BLOCKED flag and propagate limited _maxLen to _neibors @@ -9001,7 +10056,12 @@ void _LayerEdge::Block( _SolidData& data ) //if ( Is( BLOCKED )) return; Set( BLOCKED ); - _maxLen = _len; + SMESH_Comment msg( "#BLOCK shape="); + msg << data.GetShapeEdges( this )->_shapeID + << ", nodes " << _nodes[0]->GetID() << ", " << _nodes.back()->GetID(); + dumpCmd( msg + " -- BEGIN"); + + SetMaxLen( _len ); std::queue<_LayerEdge*> queue; queue.push( this ); @@ -9023,18 +10083,21 @@ void _LayerEdge::Block( _SolidData& data ) minDist = Min( pSrc.SquareDistance( pTgtN ), minDist ); minDist = Min( pTgt.SquareDistance( pSrcN ), minDist ); double newMaxLen = edge->_maxLen + 0.5 * Sqrt( minDist ); - if ( edge->_nodes[0]->getshapeId() == neibor->_nodes[0]->getshapeId() ) + //if ( edge->_nodes[0]->getshapeId() == neibor->_nodes[0]->getshapeId() ) viscous_layers_00/A3 { - newMaxLen *= edge->_lenFactor / neibor->_lenFactor; + //newMaxLen *= edge->_lenFactor / neibor->_lenFactor; + // newMaxLen *= Min( edge->_lenFactor / neibor->_lenFactor, + // neibor->_lenFactor / edge->_lenFactor ); } if ( neibor->_maxLen > newMaxLen ) { - neibor->_maxLen = newMaxLen; + neibor->SetMaxLen( newMaxLen ); if ( neibor->_maxLen < neibor->_len ) { _EdgesOnShape* eos = data.GetShapeEdges( neibor ); + int lastStep = neibor->Is( BLOCKED ) ? 1 : 0; while ( neibor->_len > neibor->_maxLen && - neibor->NbSteps() > 1 ) + neibor->NbSteps() > lastStep ) neibor->InvalidateStep( neibor->NbSteps(), *eos, /*restoreLength=*/true ); neibor->SetNewLength( neibor->_maxLen, *eos, data.GetHelper() ); //neibor->Block( data ); @@ -9043,6 +10106,7 @@ void _LayerEdge::Block( _SolidData& data ) } } } + dumpCmd( msg + " -- END"); } //================================================================================ @@ -9065,7 +10129,7 @@ void _LayerEdge::InvalidateStep( size_t curStep, const _EdgesOnShape& eos, bool TopLoc_Location loc; if ( eos.SWOLType() == TopAbs_EDGE ) { - SMDS_EdgePosition* pos = static_cast( n->GetPosition() ); + SMDS_EdgePositionPtr pos = n->GetPosition(); pos->SetUParameter( nXYZ.X() ); double f,l; Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( eos._sWOL ), loc, f,l); @@ -9073,7 +10137,7 @@ void _LayerEdge::InvalidateStep( size_t curStep, const _EdgesOnShape& eos, bool } else { - SMDS_FacePosition* pos = static_cast( n->GetPosition() ); + SMDS_FacePositionPtr pos = n->GetPosition(); pos->SetUParameter( nXYZ.X() ); pos->SetVParameter( nXYZ.Y() ); Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(eos._sWOL), loc ); @@ -9085,9 +10149,15 @@ void _LayerEdge::InvalidateStep( size_t curStep, const _EdgesOnShape& eos, bool if ( restoreLength ) { - _len -= ( nXYZ.XYZ() - curXYZ ).Modulus() / _lenFactor; + if ( NbSteps() == 0 ) + _len = 0.; + else if ( IsOnFace() && Is( MOVED )) + _len = ( nXYZ.XYZ() - SMESH_NodeXYZ( _nodes[0] )) * _normal; + else + _len -= ( nXYZ.XYZ() - curXYZ ).Modulus() / _lenFactor; } } + return; } //================================================================================ @@ -9123,10 +10193,23 @@ void _LayerEdge::SmoothPos( const vector< double >& segLen, const double tol ) int iSmoothed = GetSmoothedPos( tol ); if ( !iSmoothed ) return; - //if ( 1 || Is( DISTORTED )) + gp_XYZ normal = _normal; + if ( Is( NORMAL_UPDATED )) { - gp_XYZ normal = _normal; - if ( Is( NORMAL_UPDATED )) + double minDot = 1; + for ( size_t i = 0; i < _neibors.size(); ++i ) + { + if ( _neibors[i]->IsOnFace() ) + { + double dot = _normal * _neibors[i]->_normal; + if ( dot < minDot ) + { + normal = _neibors[i]->_normal; + minDot = dot; + } + } + } + if ( minDot == 1. ) for ( size_t i = 1; i < _pos.size(); ++i ) { normal = _pos[i] - _pos[0]; @@ -9137,44 +10220,70 @@ void _LayerEdge::SmoothPos( const vector< double >& segLen, const double tol ) break; } } - const double r = 0.2; - for ( int iter = 0; iter < 50; ++iter ) + } + const double r = 0.2; + for ( int iter = 0; iter < 50; ++iter ) + { + double minDot = 1; + for ( size_t i = Max( 1, iSmoothed-1-iter ); i < _pos.size()-1; ++i ) { - double minDot = 1; - for ( size_t i = Max( 1, iSmoothed-1-iter ); i < _pos.size()-1; ++i ) - { - gp_XYZ midPos = 0.5 * ( _pos[i-1] + _pos[i+1] ); - gp_XYZ newPos = ( 1-r ) * midPos + r * _pos[i]; - _pos[i] = newPos; - double midLen = 0.5 * ( segLen[i-1] + segLen[i+1] ); - double newLen = ( 1-r ) * midLen + r * segLen[i]; - const_cast< double& >( segLen[i] ) = newLen; - // check angle between normal and (_pos[i+1], _pos[i] ) - gp_XYZ posDir = _pos[i+1] - _pos[i]; - double size = posDir.SquareModulus(); - if ( size > RealSmall() ) - minDot = Min( minDot, ( normal * posDir ) * ( normal * posDir ) / size ); - } - if ( minDot > 0.5 * 0.5 ) - break; + gp_XYZ midPos = 0.5 * ( _pos[i-1] + _pos[i+1] ); + gp_XYZ newPos = ( 1-r ) * midPos + r * _pos[i]; + _pos[i] = newPos; + double midLen = 0.5 * ( segLen[i-1] + segLen[i+1] ); + double newLen = ( 1-r ) * midLen + r * segLen[i]; + const_cast< double& >( segLen[i] ) = newLen; + // check angle between normal and (_pos[i+1], _pos[i] ) + gp_XYZ posDir = _pos[i+1] - _pos[i]; + double size = posDir.SquareModulus(); + if ( size > RealSmall() ) + minDot = Min( minDot, ( normal * posDir ) * ( normal * posDir ) / size ); } + if ( minDot > 0.5 * 0.5 ) + break; } - // else - // { - // for ( size_t i = 1; i < _pos.size()-1; ++i ) - // { - // if ((int) i < iSmoothed && ( segLen[i] / segLen.back() < 0.5 )) - // continue; - - // double wgt = segLen[i] / segLen.back(); - // gp_XYZ normPos = _pos[0] + _normal * wgt * _len; - // gp_XYZ tgtPos = ( 1 - wgt ) * _pos[0] + wgt * _pos.back(); - // gp_XYZ newPos = ( 1 - wgt ) * normPos + wgt * tgtPos; - // _pos[i] = newPos; - // } - // } + return; +} + +//================================================================================ +/*! + * \brief Print flags + */ +//================================================================================ + +std::string _LayerEdge::DumpFlags() const +{ + SMESH_Comment dump; + for ( int flag = 1; flag < 0x1000000; flag *= 2 ) + if ( _flags & flag ) + { + EFlags f = (EFlags) flag; + switch ( f ) { + case TO_SMOOTH: dump << "TO_SMOOTH"; break; + case MOVED: dump << "MOVED"; break; + case SMOOTHED: dump << "SMOOTHED"; break; + case DIFFICULT: dump << "DIFFICULT"; break; + case ON_CONCAVE_FACE: dump << "ON_CONCAVE_FACE"; break; + case BLOCKED: dump << "BLOCKED"; break; + case INTERSECTED: dump << "INTERSECTED"; break; + case NORMAL_UPDATED: dump << "NORMAL_UPDATED"; break; + case UPD_NORMAL_CONV: dump << "UPD_NORMAL_CONV"; break; + case MARKED: dump << "MARKED"; break; + case MULTI_NORMAL: dump << "MULTI_NORMAL"; break; + case NEAR_BOUNDARY: dump << "NEAR_BOUNDARY"; break; + case SMOOTHED_C1: dump << "SMOOTHED_C1"; break; + case DISTORTED: dump << "DISTORTED"; break; + case RISKY_SWOL: dump << "RISKY_SWOL"; break; + case SHRUNK: dump << "SHRUNK"; break; + case UNUSED_FLAG: dump << "UNUSED_FLAG"; break; + } + dump << " "; + } + cout << dump << endl; + return dump; } + //================================================================================ /*! * \brief Create layers of prisms @@ -9194,7 +10303,7 @@ bool _ViscousBuilder::refine(_SolidData& data) double f,l, u = 0; gp_XY uv; vector< gp_XYZ > pos3D; - bool isOnEdge; + bool isOnEdge, isTooConvexFace = false; TGeomID prevBaseId = -1; TNode2Edge* n2eMap = 0; TNode2Edge::iterator n2e; @@ -9233,11 +10342,11 @@ bool _ViscousBuilder::refine(_SolidData& data) surface = helper.GetSurface( geomFace ); // propagate _toSmooth back to _eosC1, which was unset in findShapesToSmooth() for ( size_t i = 0; i < eos._eosC1.size(); ++i ) - { eos._eosC1[ i ]->_toSmooth = true; - for ( size_t j = 0; j < eos._eosC1[i]->_edges.size(); ++j ) - eos._eosC1[i]->_edges[j]->Set( _LayerEdge::SMOOTHED_C1 ); - } + + isTooConvexFace = false; + if ( _ConvexFace* cf = data.GetConvexFace( eos._shapeID )) + isTooConvexFace = cf->_isTooCurved; } vector< double > segLen; @@ -9253,8 +10362,8 @@ bool _ViscousBuilder::refine(_SolidData& data) if ( eos._sWOL.IsNull() ) { bool useNormal = true; - bool usePos = false; - bool smoothed = false; + bool usePos = false; + bool smoothed = false; double preci = 0.1 * edge._len; if ( eos._toSmooth && edge._pos.size() > 2 ) { @@ -9262,8 +10371,7 @@ bool _ViscousBuilder::refine(_SolidData& data) } if ( smoothed ) { - if ( !surface.IsNull() && - !data._convexFaces.count( eos._shapeID )) // edge smoothed on FACE + if ( !surface.IsNull() && !isTooConvexFace ) // edge smoothed on FACE { useNormal = usePos = false; gp_Pnt2d uv = helper.GetNodeUV( geomFace, edge._nodes[0] ); @@ -9338,14 +10446,31 @@ bool _ViscousBuilder::refine(_SolidData& data) } else if ( eos._isRegularSWOL ) // usual SWOL { - for ( size_t j = 1; j < edge._pos.size(); ++j ) - segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus(); + if ( edge.Is( _LayerEdge::SMOOTHED )) + { + SMESH_NodeXYZ p0( edge._nodes[0] ); + for ( size_t j = 1; j < edge._pos.size(); ++j ) + { + gp_XYZ pj = surface->Value( edge._pos[j].X(), edge._pos[j].Y() ).XYZ(); + segLen[j] = ( pj - p0 ) * edge._normal; + } + } + else + { + for ( size_t j = 1; j < edge._pos.size(); ++j ) + segLen[j] = segLen[j-1] + (edge._pos[j-1] - edge._pos[j] ).Modulus(); + } } - else if ( !surface.IsNull() ) // SWOL surface with singularities + else // SWOL is surface with singularities or irregularly parametrized curve { pos3D.resize( edge._pos.size() ); - for ( size_t j = 0; j < edge._pos.size(); ++j ) - pos3D[j] = surface->Value( edge._pos[j].X(), edge._pos[j].Y() ).XYZ(); + + if ( !surface.IsNull() ) + for ( size_t j = 0; j < edge._pos.size(); ++j ) + pos3D[j] = surface->Value( edge._pos[j].X(), edge._pos[j].Y() ).XYZ(); + else if ( !curve.IsNull() ) + for ( size_t j = 0; j < edge._pos.size(); ++j ) + pos3D[j] = curve->Value( edge._pos[j].X() ).XYZ(); for ( size_t j = 1; j < edge._pos.size(); ++j ) segLen[j] = segLen[j-1] + ( pos3D[j-1] - pos3D[j] ).Modulus(); @@ -9377,41 +10502,32 @@ bool _ViscousBuilder::refine(_SolidData& data) if ( n2eMap && (( n2e = n2eMap->find( edge._nodes[0] )) != n2eMap->end() )) { edgeOnSameNode = n2e->second; - useExistingPos = ( edgeOnSameNode->_len < edge._len ); + useExistingPos = ( edgeOnSameNode->_len < edge._len || + segLen[0] == segLen.back() ); // too short inflation step (bos #20643) const gp_XYZ& otherTgtPos = edgeOnSameNode->_pos.back(); SMDS_PositionPtr lastPos = tgtNode->GetPosition(); if ( isOnEdge ) { - SMDS_EdgePosition* epos = static_cast( lastPos ); + SMDS_EdgePositionPtr epos = lastPos; epos->SetUParameter( otherTgtPos.X() ); } else { - SMDS_FacePosition* fpos = static_cast( lastPos ); + SMDS_FacePositionPtr fpos = lastPos; fpos->SetUParameter( otherTgtPos.X() ); fpos->SetVParameter( otherTgtPos.Y() ); } } - // calculate height of the first layer - double h0; - const double T = segLen.back(); //data._hyp.GetTotalThickness(); - const double f = eos._hyp.GetStretchFactor(); - const int N = eos._hyp.GetNumberLayers(); - const double fPowN = pow( f, N ); - if ( fPowN - 1 <= numeric_limits::min() ) - h0 = T / N; - else - h0 = T * ( f - 1 )/( fPowN - 1 ); - - const double zeroLen = std::numeric_limits::min(); // create intermediate nodes - double hSum = 0, hi = h0/f; + const double h0 = eos._hyp.Get1stLayerThickness( segLen.back() ); + const double zeroLen = std::numeric_limits::min(); + double hSum = 0, hi = h0/eos._hyp.GetStretchFactor(); size_t iSeg = 1; for ( size_t iStep = 1; iStep < edge._nodes.size(); ++iStep ) { // compute an intermediate position - hi *= f; + hi *= eos._hyp.GetStretchFactor(); hSum += hi; while ( hSum > segLen[iSeg] && iSeg < segLen.size()-1 ) ++iSeg; @@ -9476,7 +10592,7 @@ bool _ViscousBuilder::refine(_SolidData& data) u = helper.GetNodeU( geomEdge, node ); pos = curve->Value( u ).Transformed(loc); - SMDS_EdgePosition* epos = static_cast( node->GetPosition() ); + SMDS_EdgePositionPtr epos = node->GetPosition(); epos->SetUParameter( u ); } else @@ -9486,7 +10602,7 @@ bool _ViscousBuilder::refine(_SolidData& data) uv = helper.GetNodeUV( geomFace, node ); pos = surface->Value( uv ); - SMDS_FacePosition* fpos = static_cast( node->GetPosition() ); + SMDS_FacePositionPtr fpos = node->GetPosition(); fpos->SetUParameter( uv.X() ); fpos->SetVParameter( uv.Y() ); } @@ -9527,7 +10643,6 @@ bool _ViscousBuilder::refine(_SolidData& data) set< vector* > nnSet; set< int > degenEdgeInd; vector degenVols; - vector isRiskySWOL; TopExp_Explorer exp( data._solid, TopAbs_FACE ); for ( ; exp.More(); exp.Next() ) @@ -9535,6 +10650,11 @@ bool _ViscousBuilder::refine(_SolidData& data) const TGeomID faceID = getMeshDS()->ShapeToIndex( exp.Current() ); if ( data._ignoreFaceIds.count( faceID )) continue; + _EdgesOnShape* eos = data.GetShapeEdges( faceID ); + SMDS_MeshGroup* group = StdMeshers_ViscousLayers::CreateGroup( eos->_hyp.GetGroupName(), + *helper.GetMesh(), + SMDSAbs_Volume ); + std::vector< const SMDS_MeshElement* > vols; const bool isReversedFace = data._reversedFaceIds.count( faceID ); SMESHDS_SubMesh* fSubM = getMeshDS()->MeshElements( exp.Current() ); SMDS_ElemIteratorPtr fIt = fSubM->GetElements(); @@ -9545,7 +10665,6 @@ bool _ViscousBuilder::refine(_SolidData& data) nnVec.resize( nbNodes ); nnSet.clear(); degenEdgeInd.clear(); - isRiskySWOL.resize( nbNodes ); size_t maxZ = 0, minZ = std::numeric_limits::max(); SMDS_NodeIteratorPtr nIt = face->nodeIterator(); for ( int iN = 0; iN < nbNodes; ++iN ) @@ -9556,7 +10675,6 @@ bool _ViscousBuilder::refine(_SolidData& data) nnVec[ i ] = & edge->_nodes; maxZ = std::max( maxZ, nnVec[ i ]->size() ); minZ = std::min( minZ, nnVec[ i ]->size() ); - //isRiskySWOL[ i ] = edge->Is( _LayerEdge::RISKY_SWOL ); if ( helper.HasDegeneratedEdges() ) nnSet.insert( nnVec[ i ]); @@ -9567,14 +10685,20 @@ bool _ViscousBuilder::refine(_SolidData& data) if ( 0 < nnSet.size() && nnSet.size() < 3 ) continue; + vols.clear(); + const SMDS_MeshElement* vol; + switch ( nbNodes ) { case 3: // TRIA { // PENTA for ( size_t iZ = 1; iZ < minZ; ++iZ ) - helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1], - (*nnVec[0])[iZ], (*nnVec[1])[iZ], (*nnVec[2])[iZ]); + { + vol = helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], (*nnVec[2])[iZ-1], + (*nnVec[0])[iZ], (*nnVec[1])[iZ], (*nnVec[2])[iZ]); + vols.push_back( vol ); + } for ( size_t iZ = minZ; iZ < maxZ; ++iZ ) { @@ -9587,16 +10711,18 @@ bool _ViscousBuilder::refine(_SolidData& data) int i2 = *degenEdgeInd.begin(); int i0 = helper.WrapIndex( i2 - 1, nbNodes ); int i1 = helper.WrapIndex( i2 + 1, nbNodes ); - helper.AddVolume( (*nnVec[i0])[iZ-1], (*nnVec[i1])[iZ-1], - (*nnVec[i1])[iZ ], (*nnVec[i0])[iZ ], (*nnVec[i2]).back()); + vol = helper.AddVolume( (*nnVec[i0])[iZ-1], (*nnVec[i1])[iZ-1], + (*nnVec[i1])[iZ ], (*nnVec[i0])[iZ ], (*nnVec[i2]).back()); + vols.push_back( vol ); } else // TETRA { int i3 = !degenEdgeInd.count(0) ? 0 : !degenEdgeInd.count(1) ? 1 : 2; - helper.AddVolume( (*nnVec[ 0 ])[ i3 == 0 ? iZ-1 : nnVec[0]->size()-1 ], - (*nnVec[ 1 ])[ i3 == 1 ? iZ-1 : nnVec[1]->size()-1 ], - (*nnVec[ 2 ])[ i3 == 2 ? iZ-1 : nnVec[2]->size()-1 ], - (*nnVec[ i3 ])[ iZ ]); + vol = helper.AddVolume( (*nnVec[ 0 ])[ i3 == 0 ? iZ-1 : nnVec[0]->size()-1 ], + (*nnVec[ 1 ])[ i3 == 1 ? iZ-1 : nnVec[1]->size()-1 ], + (*nnVec[ 2 ])[ i3 == 2 ? iZ-1 : nnVec[2]->size()-1 ], + (*nnVec[ i3 ])[ iZ ]); + vols.push_back( vol ); } } break; // TRIA @@ -9605,10 +10731,13 @@ bool _ViscousBuilder::refine(_SolidData& data) { // HEX for ( size_t iZ = 1; iZ < minZ; ++iZ ) - helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], - (*nnVec[2])[iZ-1], (*nnVec[3])[iZ-1], - (*nnVec[0])[iZ], (*nnVec[1])[iZ], - (*nnVec[2])[iZ], (*nnVec[3])[iZ]); + { + vol = helper.AddVolume( (*nnVec[0])[iZ-1], (*nnVec[1])[iZ-1], + (*nnVec[2])[iZ-1], (*nnVec[3])[iZ-1], + (*nnVec[0])[iZ], (*nnVec[1])[iZ], + (*nnVec[2])[iZ], (*nnVec[3])[iZ]); + vols.push_back( vol ); + } for ( size_t iZ = minZ; iZ < maxZ; ++iZ ) { @@ -9627,9 +10756,9 @@ bool _ViscousBuilder::refine(_SolidData& data) int i0 = helper.WrapIndex( i3 + 1, nbNodes ); int i1 = helper.WrapIndex( i0 + 1, nbNodes ); - const SMDS_MeshElement* vol = - helper.AddVolume( nnVec[i3]->back(), (*nnVec[i0])[iZ], (*nnVec[i0])[iZ-1], - nnVec[i2]->back(), (*nnVec[i1])[iZ], (*nnVec[i1])[iZ-1]); + vol = helper.AddVolume( nnVec[i3]->back(), (*nnVec[i0])[iZ], (*nnVec[i0])[iZ-1], + nnVec[i2]->back(), (*nnVec[i1])[iZ], (*nnVec[i1])[iZ-1]); + vols.push_back( vol ); if ( !ok && vol ) degenVols.push_back( vol ); } @@ -9637,15 +10766,15 @@ bool _ViscousBuilder::refine(_SolidData& data) default: // degen HEX { - const SMDS_MeshElement* vol = - helper.AddVolume( nnVec[0]->size() > iZ-1 ? (*nnVec[0])[iZ-1] : nnVec[0]->back(), - nnVec[1]->size() > iZ-1 ? (*nnVec[1])[iZ-1] : nnVec[1]->back(), - nnVec[2]->size() > iZ-1 ? (*nnVec[2])[iZ-1] : nnVec[2]->back(), - nnVec[3]->size() > iZ-1 ? (*nnVec[3])[iZ-1] : nnVec[3]->back(), - nnVec[0]->size() > iZ ? (*nnVec[0])[iZ] : nnVec[0]->back(), - nnVec[1]->size() > iZ ? (*nnVec[1])[iZ] : nnVec[1]->back(), - nnVec[2]->size() > iZ ? (*nnVec[2])[iZ] : nnVec[2]->back(), - nnVec[3]->size() > iZ ? (*nnVec[3])[iZ] : nnVec[3]->back()); + vol = helper.AddVolume( nnVec[0]->size() > iZ-1 ? (*nnVec[0])[iZ-1] : nnVec[0]->back(), + nnVec[1]->size() > iZ-1 ? (*nnVec[1])[iZ-1] : nnVec[1]->back(), + nnVec[2]->size() > iZ-1 ? (*nnVec[2])[iZ-1] : nnVec[2]->back(), + nnVec[3]->size() > iZ-1 ? (*nnVec[3])[iZ-1] : nnVec[3]->back(), + nnVec[0]->size() > iZ ? (*nnVec[0])[iZ] : nnVec[0]->back(), + nnVec[1]->size() > iZ ? (*nnVec[1])[iZ] : nnVec[1]->back(), + nnVec[2]->size() > iZ ? (*nnVec[2])[iZ] : nnVec[2]->back(), + nnVec[3]->size() > iZ ? (*nnVec[3])[iZ] : nnVec[3]->back()); + vols.push_back( vol ); degenVols.push_back( vol ); } } @@ -9656,6 +10785,11 @@ bool _ViscousBuilder::refine(_SolidData& data) return error("Not supported type of element", data._index); } // switch ( nbNodes ) + + if ( group ) + for ( size_t i = 0; i < vols.size(); ++i ) + group->Add( vols[ i ]); + } // while ( fIt->more() ) } // loop on FACEs @@ -9664,16 +10798,491 @@ bool _ViscousBuilder::refine(_SolidData& data) SMESH_ComputeErrorPtr& err = _mesh->GetSubMesh( data._solid )->GetComputeError(); if ( !err || err->IsOK() ) { - err.reset( new SMESH_ComputeError( COMPERR_WARNING, - "Degenerated volumes created" )); - err->myBadElements.insert( err->myBadElements.end(), - degenVols.begin(),degenVols.end() ); + SMESH_BadInputElements* badElems = + new SMESH_BadInputElements( getMeshDS(), COMPERR_WARNING, "Bad quality volumes created" ); + badElems->myBadElements.insert( badElems->myBadElements.end(), + degenVols.begin(),degenVols.end() ); + err.reset( badElems ); } } return true; } +namespace VISCOUS_3D +{ + struct ShrinkFace; + //-------------------------------------------------------------------------------- + /*! + * \brief Pair of periodic FACEs + */ + struct PeriodicFaces + { + typedef StdMeshers_ProjectionUtils::TrsfFinder3D Trsf; + + ShrinkFace* _shriFace[2]; + TNodeNodeMap _nnMap; + Trsf _trsf; + + PeriodicFaces( ShrinkFace* sf1, ShrinkFace* sf2 ): _shriFace{ sf1, sf2 } {} + bool IncludeShrunk( const TopoDS_Face& face, const TopTools_MapOfShape& shrunkFaces ) const; + bool MoveNodes( const TopoDS_Face& tgtFace ); + void Clear() { _nnMap.clear(); } + bool IsEmpty() const { return _nnMap.empty(); } + }; + + //-------------------------------------------------------------------------------- + /*! + * \brief Shrink FACE data used to find periodic FACEs + */ + struct ShrinkFace + { + // ................................................................................ + struct BndPart //!< part of FACE boundary, either shrink or no-shrink + { + bool _isShrink, _isReverse; + int _nbSegments; + AverageHyp* _hyp; + std::vector< SMESH_NodeXYZ > _nodes; + TopAbs_ShapeEnum _vertSWOLType[2]; // shrink part includes VERTEXes + AverageHyp* _vertHyp[2]; + double _edgeWOLLen[2]; // length of wol EDGE + double _tol; // to compare _edgeWOLLen's + + BndPart(): + _isShrink(0), _isReverse(0), _nbSegments(0), _hyp(0), + _vertSWOLType{ TopAbs_WIRE, TopAbs_WIRE }, _vertHyp{ 0, 0 }, _edgeWOLLen{ 0., 0.} + {} + + bool IsEqualLengthEWOL( const BndPart& other ) const + { + return ( std::abs( _edgeWOLLen[0] - other._edgeWOLLen[0] ) < _tol && + std::abs( _edgeWOLLen[1] - other._edgeWOLLen[1] ) < _tol ); + } + + bool operator==( const BndPart& other ) const + { + return ( _isShrink == other._isShrink && + _nbSegments == other._nbSegments && + _nodes.size() == other._nodes.size() && + vertSWOLType1() == other.vertSWOLType1() && + vertSWOLType2() == other.vertSWOLType2() && + (( !_isShrink ) || + ( *_hyp == *other._hyp && + vertHyp1() == other.vertHyp1() && + vertHyp2() == other.vertHyp2() && + IsEqualLengthEWOL( other ))) + ); + } + bool CanAppend( const BndPart& other ) + { + return ( _isShrink == other._isShrink && + (( !_isShrink ) || + ( *_hyp == *other._hyp && + *_hyp == vertHyp2() && + vertHyp2() == other.vertHyp1() )) + ); + } + void Append( const BndPart& other ) + { + _nbSegments += other._nbSegments; + bool hasCommonNode = ( _nodes.back()->GetID() == other._nodes.front()->GetID() ); + _nodes.insert( _nodes.end(), other._nodes.begin() + hasCommonNode, other._nodes.end() ); + _vertSWOLType[1] = other._vertSWOLType[1]; + if ( _isShrink ) { + _vertHyp[1] = other._vertHyp[1]; + _edgeWOLLen[1] = other._edgeWOLLen[1]; + } + } + const SMDS_MeshNode* Node(size_t i) const + { + return _nodes[ _isReverse ? ( _nodes.size() - 1 - i ) : i ]._node; + } + void Reverse() { _isReverse = !_isReverse; } + const TopAbs_ShapeEnum& vertSWOLType1() const { return _vertSWOLType[ _isReverse ]; } + const TopAbs_ShapeEnum& vertSWOLType2() const { return _vertSWOLType[ !_isReverse ]; } + const AverageHyp& vertHyp1() const { return *(_vertHyp[ _isReverse ]); } + const AverageHyp& vertHyp2() const { return *(_vertHyp[ !_isReverse ]); } + }; + // ................................................................................ + + SMESH_subMesh* _subMesh; + _SolidData* _data1; + _SolidData* _data2; + + std::list< BndPart > _boundary; + int _boundarySize, _nbBoundaryParts; + + void Init( SMESH_subMesh* sm, _SolidData* sd1, _SolidData* sd2 ) + { + _subMesh = sm; _data1 = sd1; _data2 = sd2; + } + bool IsSame( const TopoDS_Face& face ) const + { + return _subMesh->GetSubShape().IsSame( face ); + } + bool IsShrunk( const TopTools_MapOfShape& shrunkFaces ) const + { + return shrunkFaces.Contains( _subMesh->GetSubShape() ); + } + + //================================================================================ + /*! + * Check if meshes on two FACEs are equal + */ + bool IsPeriodic( ShrinkFace& other, PeriodicFaces& periodic ) + { + if ( !IsSameNbElements( other )) + return false; + + this->SetBoundary(); + other.SetBoundary(); + if ( this->_boundarySize != other._boundarySize || + this->_nbBoundaryParts != other._nbBoundaryParts ) + return false; + + for ( int isReverse = 0; isReverse < 2; ++isReverse ) + { + if ( isReverse ) + Reverse( _boundary ); + + // check boundaries + bool equalBoundary = false; + for ( int iP = 0; iP < _nbBoundaryParts && !equalBoundary; ++iP ) + { + if ( ! ( equalBoundary = ( this->_boundary == other._boundary ))) + // set first part at end + _boundary.splice( _boundary.end(), _boundary, _boundary.begin() ); + } + if ( !equalBoundary ) + continue; + + // check connectivity + std::set elemsThis, elemsOther; + this->GetElements( elemsThis ); + other.GetElements( elemsOther ); + SMESH_MeshEditor::Sew_Error err = + SMESH_MeshEditor::FindMatchingNodes( elemsThis, elemsOther, + this->_boundary.front().Node(0), + other._boundary.front().Node(0), + this->_boundary.front().Node(1), + other._boundary.front().Node(1), + periodic._nnMap ); + if ( err != SMESH_MeshEditor::SEW_OK ) + continue; + + // check node positions + std::vector< gp_XYZ > srcPnts, tgtPnts; + this->GetBoundaryPoints( srcPnts ); + other.GetBoundaryPoints( tgtPnts ); + if ( !periodic._trsf.Solve( srcPnts, tgtPnts )) { + continue; + } + double tol = std::numeric_limits::max(); // tolerance by segment size + for ( size_t i = 1; i < srcPnts.size(); ++i ) { + tol = Min( tol, ( srcPnts[i-1] - srcPnts[i] ).SquareModulus() ); + } + tol = 0.01 * Sqrt( tol ); + for ( BndPart& boundary : _boundary ) { // tolerance by VL thickness + if ( boundary._isShrink ) + tol = Min( tol, boundary._hyp->Get1stLayerThickness() / 50. ); + } + bool nodeCoincide = true; + TNodeNodeMap::iterator n2n = periodic._nnMap.begin(); + for ( ; n2n != periodic._nnMap.end() && nodeCoincide; ++n2n ) + { + SMESH_NodeXYZ nSrc = n2n->first; + SMESH_NodeXYZ nTgt = n2n->second; + gp_XYZ pTgt = periodic._trsf.Transform( nSrc ); + nodeCoincide = (( pTgt - nTgt ).SquareModulus() < tol * tol ); + } + if ( nodeCoincide ) + return true; + } + return false; + } + + bool IsSameNbElements( ShrinkFace& other ) // check number of mesh faces + { + SMESHDS_SubMesh* sm1 = this->_subMesh->GetSubMeshDS(); + SMESHDS_SubMesh* sm2 = other._subMesh->GetSubMeshDS(); + return ( sm1->NbElements() == sm2->NbElements() && + sm1->NbNodes() == sm2->NbNodes() ); + } + + void Reverse( std::list< BndPart >& boundary ) + { + boundary.reverse(); + for ( std::list< BndPart >::iterator part = boundary.begin(); part != boundary.end(); ++part ) + part->Reverse(); + } + + void SetBoundary() + { + if ( !_boundary.empty() ) + return; + + TopoDS_Face F = TopoDS::Face( _subMesh->GetSubShape() ); + if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD ); + std::list< TopoDS_Edge > edges; + std::list< int > nbEdgesInWire; + /*int nbWires =*/ SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire); + + // std::list< TopoDS_Edge >::iterator edgesEnd = edges.end(); + // if ( nbWires > 1 ) { + // edgesEnd = edges.begin(); + // std::advance( edgesEnd, nbEdgesInWire.front() ); + // } + StdMeshers_FaceSide fSide( F, edges, _subMesh->GetFather(), + /*fwd=*/true, /*skipMedium=*/true ); + _boundarySize = fSide.NbSegments(); + + //TopoDS_Vertex vv[2]; + //std::list< TopoDS_Edge >::iterator edgeIt = edges.begin(); + for ( int iE = 0; iE < nbEdgesInWire.front(); ++iE ) + { + BndPart bndPart; + + std::vector nodes = fSide.GetOrderedNodes( iE ); + bndPart._nodes.assign( nodes.begin(), nodes.end() ); + bndPart._nbSegments = bndPart._nodes.size() - 1; + + _EdgesOnShape* eos = _data1->GetShapeEdges( fSide.EdgeID( iE )); + + bndPart._isShrink = ( eos->SWOLType() == TopAbs_FACE ); + if ( bndPart._isShrink ) + if (( _data1->_noShrinkShapes.count( eos->_shapeID )) || + ( _data2 && _data2->_noShrinkShapes.count( eos->_shapeID ))) + bndPart._isShrink = false; + + if ( bndPart._isShrink ) + { + bndPart._hyp = & eos->_hyp; + _EdgesOnShape* eov[2] = { _data1->GetShapeEdges( fSide.FirstVertex( iE )), + _data1->GetShapeEdges( fSide.LastVertex ( iE )) }; + for ( int iV = 0; iV < 2; ++iV ) + { + bndPart._vertHyp [iV] = & eov[iV]->_hyp; + bndPart._vertSWOLType[iV] = eov[iV]->SWOLType(); + if ( _data1->_noShrinkShapes.count( eov[iV]->_shapeID )) + bndPart._vertSWOLType[iV] = TopAbs_SHAPE; + if ( _data2 && bndPart._vertSWOLType[iV] != TopAbs_SHAPE ) + { + eov[iV] = _data2->GetShapeEdges( iV ? fSide.LastVertex(iE) : fSide.FirstVertex(iE )); + if ( _data2->_noShrinkShapes.count( eov[iV]->_shapeID )) + bndPart._vertSWOLType[iV] = TopAbs_SHAPE; + else if ( eov[iV]->SWOLType() > bndPart._vertSWOLType[iV] ) + bndPart._vertSWOLType[iV] = eov[iV]->SWOLType(); + } + } + bndPart._edgeWOLLen[0] = fSide.EdgeLength( iE - 1 ); + bndPart._edgeWOLLen[1] = fSide.EdgeLength( iE + 1 ); + + bndPart._tol = std::numeric_limits::max(); // tolerance by segment size + for ( size_t i = 1; i < bndPart._nodes.size(); ++i ) + bndPart._tol = Min( bndPart._tol, + ( bndPart._nodes[i-1] - bndPart._nodes[i] ).SquareModulus() ); + } + + if ( _boundary.empty() || ! _boundary.back().CanAppend( bndPart )) + _boundary.push_back( bndPart ); + else + _boundary.back().Append( bndPart ); + } + + _nbBoundaryParts = _boundary.size(); + if ( _nbBoundaryParts > 1 && _boundary.front()._isShrink == _boundary.back()._isShrink ) + { + _boundary.back().Append( _boundary.front() ); + _boundary.pop_front(); + --_nbBoundaryParts; + } + } + + void GetElements( std::set& theElems) + { + if ( SMESHDS_SubMesh* sm = _subMesh->GetSubMeshDS() ) + for ( SMDS_ElemIteratorPtr fIt = sm->GetElements(); fIt->more(); ) + theElems.insert( theElems.end(), fIt->next() ); + + return ; + } + + void GetBoundaryPoints( std::vector< gp_XYZ >& points ) + { + points.reserve( _boundarySize ); + size_t nb = _boundary.rbegin()->_nodes.size(); + smIdType lastID = _boundary.rbegin()->Node( nb - 1 )->GetID(); + std::list< BndPart >::const_iterator part = _boundary.begin(); + for ( ; part != _boundary.end(); ++part ) + { + size_t nb = part->_nodes.size(); + size_t iF = 0; + size_t iR = nb - 1; + size_t* i = part->_isReverse ? &iR : &iF; + if ( part->_nodes[ *i ]->GetID() == lastID ) + ++iF, --iR; + for ( ; iF < nb; ++iF, --iR ) + points.push_back( part->_nodes[ *i ]); + --iF, ++iR; + lastID = part->_nodes[ *i ]->GetID(); + } + } + }; // struct ShrinkFace + + //-------------------------------------------------------------------------------- + /*! + * \brief Periodic FACEs + */ + struct Periodicity + { + std::vector< ShrinkFace > _shrinkFaces; + std::vector< PeriodicFaces > _periodicFaces; + + PeriodicFaces* GetPeriodic( const TopoDS_Face& face, const TopTools_MapOfShape& shrunkFaces ) + { + for ( size_t i = 0; i < _periodicFaces.size(); ++i ) + if ( _periodicFaces[ i ].IncludeShrunk( face, shrunkFaces )) + return & _periodicFaces[ i ]; + return 0; + } + void ClearPeriodic( const TopoDS_Face& face ) + { + for ( size_t i = 0; i < _periodicFaces.size(); ++i ) + if ( _periodicFaces[ i ]._shriFace[0]->IsSame( face ) || + _periodicFaces[ i ]._shriFace[1]->IsSame( face )) + _periodicFaces[ i ].Clear(); + } + }; + + //================================================================================ + /*! + * Check if a pair includes the given FACE and the other FACE is already shrunk + */ + bool PeriodicFaces::IncludeShrunk( const TopoDS_Face& face, + const TopTools_MapOfShape& shrunkFaces ) const + { + if ( IsEmpty() ) return false; + return (( _shriFace[0]->IsSame( face ) && _shriFace[1]->IsShrunk( shrunkFaces )) || + ( _shriFace[1]->IsSame( face ) && _shriFace[0]->IsShrunk( shrunkFaces ))); + } + + //================================================================================ + /*! + * Make equal meshes on periodic faces by moving corresponding nodes + */ + bool PeriodicFaces::MoveNodes( const TopoDS_Face& tgtFace ) + { + int iTgt = _shriFace[1]->IsSame( tgtFace ); + int iSrc = 1 - iTgt; + + _SolidData* dataSrc = _shriFace[iSrc]->_data1; + _SolidData* dataTgt = _shriFace[iTgt]->_data1; + + Trsf * trsf = & _trsf, trsfInverse; + if ( iSrc != 0 ) + { + trsfInverse = _trsf; + if ( !trsfInverse.Invert()) + return false; + trsf = &trsfInverse; + } + SMESHDS_Mesh* meshDS = dataSrc->GetHelper().GetMeshDS(); + + dumpFunction(SMESH_Comment("periodicMoveNodes_F") + << _shriFace[iSrc]->_subMesh->GetId() << "_F" + << _shriFace[iTgt]->_subMesh->GetId() ); + TNode2Edge::iterator n2e; + TNodeNodeMap::iterator n2n = _nnMap.begin(); + for ( ; n2n != _nnMap.end(); ++n2n ) + { + const SMDS_MeshNode* const* nn = & n2n->first; + const SMDS_MeshNode* nSrc = nn[ iSrc ]; + const SMDS_MeshNode* nTgt = nn[ iTgt ]; + + if (( nSrc->GetPosition()->GetDim() == 2 ) || + (( n2e = dataSrc->_n2eMap.find( nSrc )) == dataSrc->_n2eMap.end() )) + { + SMESH_NodeXYZ pSrc = nSrc; + gp_XYZ pTgt = trsf->Transform( pSrc ); + meshDS->MoveNode( nTgt, pTgt.X(), pTgt.Y(), pTgt.Z() ); + } + else + { + _LayerEdge* leSrc = n2e->second; + n2e = dataTgt->_n2eMap.find( nTgt ); + if ( n2e == dataTgt->_n2eMap.end() ) + break; + _LayerEdge* leTgt = n2e->second; + if ( leSrc->_nodes.size() != leTgt->_nodes.size() ) + break; + for ( size_t iN = 1; iN < leSrc->_nodes.size(); ++iN ) + { + SMESH_NodeXYZ pSrc = leSrc->_nodes[ iN ]; + gp_XYZ pTgt = trsf->Transform( pSrc ); + meshDS->MoveNode( leTgt->_nodes[ iN ], pTgt.X(), pTgt.Y(), pTgt.Z() ); + + dumpMove( leTgt->_nodes[ iN ]); + } + } + } + bool done = ( n2n == _nnMap.end() ); + debugMsg( "PeriodicFaces::MoveNodes " + << _shriFace[iSrc]->_subMesh->GetId() << " -> " + << _shriFace[iTgt]->_subMesh->GetId() << " -- " + << ( done ? "DONE" : "FAIL")); + dumpFunctionEnd(); + + return done; + } +} // namespace VISCOUS_3D; Periodicity part + + +//================================================================================ +/*! + * \brief Find FACEs to shrink, that are equally meshed before shrink (i.e. periodic) + * and should remain equal after shrink + */ +//================================================================================ + +void _ViscousBuilder::findPeriodicFaces() +{ + // make map of (ids of FACEs to shrink mesh on) to (list of _SolidData containing + // _LayerEdge's inflated along FACE or EDGE) + std::map< TGeomID, std::list< _SolidData* > > id2sdMap; + for ( size_t i = 0 ; i < _sdVec.size(); ++i ) + { + _SolidData& data = _sdVec[i]; + std::map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin(); + for (; s2s != data._shrinkShape2Shape.end(); ++s2s ) + if ( s2s->second.ShapeType() == TopAbs_FACE ) + id2sdMap[ getMeshDS()->ShapeToIndex( s2s->second )].push_back( &data ); + } + + _periodicity.reset( new Periodicity ); + _periodicity->_shrinkFaces.resize( id2sdMap.size() ); + + std::map< TGeomID, std::list< _SolidData* > >::iterator id2sdIt = id2sdMap.begin(); + for ( size_t i = 0; i < id2sdMap.size(); ++i, ++id2sdIt ) + { + _SolidData* sd1 = id2sdIt->second.front(); + _SolidData* sd2 = id2sdIt->second.back(); + _periodicity->_shrinkFaces[ i ].Init( _mesh->GetSubMeshContaining( id2sdIt->first ), sd1, sd2 ); + } + + for ( size_t i1 = 0; i1 < _periodicity->_shrinkFaces.size(); ++i1 ) + for ( size_t i2 = i1 + 1; i2 < _periodicity->_shrinkFaces.size(); ++i2 ) + { + PeriodicFaces pf( & _periodicity->_shrinkFaces[ i1 ], + & _periodicity->_shrinkFaces[ i2 ]); + if ( pf._shriFace[0]->IsPeriodic( *pf._shriFace[1], pf )) + { + _periodicity->_periodicFaces.push_back( pf ); + } + } + return; +} + //================================================================================ /*! * \brief Shrink 2D mesh on faces to let space for inflated layers @@ -9690,11 +11299,11 @@ bool _ViscousBuilder::shrink(_SolidData& theData) _SolidData& data = _sdVec[i]; map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin(); for (; s2s != data._shrinkShape2Shape.end(); ++s2s ) - if ( s2s->second.ShapeType() == TopAbs_FACE && !_shrinkedFaces.Contains( s2s->second )) + if ( s2s->second.ShapeType() == TopAbs_FACE && !_shrunkFaces.Contains( s2s->second )) { f2sdMap[ getMeshDS()->ShapeToIndex( s2s->second )].push_back( &data ); - // Put mesh faces on the shrinked FACE to the proxy sub-mesh to avoid + // Put mesh faces on the shrunk FACE to the proxy sub-mesh to avoid // usage of mesh faces made in addBoundaryElements() by the 3D algo or // by StdMeshers_QuadToTriaAdaptor if ( SMESHDS_SubMesh* smDS = getMeshDS()->MeshElements( s2s->second )) @@ -9733,7 +11342,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData) vector< _EdgesOnShape* > subEOS; vector< _LayerEdge* > lEdges; - // loop on FACEs to srink mesh on + // loop on FACEs to shrink mesh on map< TGeomID, list< _SolidData* > >::iterator f2sd = f2sdMap.begin(); for ( ; f2sd != f2sdMap.end(); ++f2sd ) { @@ -9746,21 +11355,31 @@ bool _ViscousBuilder::shrink(_SolidData& theData) continue; _SolidData& data = *dataList.front(); + _SolidData* data2 = dataList.size() > 1 ? dataList.back() : 0; const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first )); SMESH_subMesh* sm = _mesh->GetSubMesh( F ); SMESHDS_SubMesh* smDS = sm->GetSubMeshDS(); Handle(Geom_Surface) surface = BRep_Tool::Surface( F ); - _shrinkedFaces.Add( F ); + _shrunkFaces.Add( F ); helper.SetSubShape( F ); + // ============================== + // Use periodicity to move nodes + // ============================== + + PeriodicFaces* periodic = _periodicity->GetPeriodic( F, _shrunkFaces ); + bool movedByPeriod = ( periodic && periodic->MoveNodes( F )); + // =========================== // Prepare data for shrinking // =========================== - // Collect nodes to smooth, they are marked at the beginning of this method + // Collect nodes to smooth (they are marked at the beginning of this method) vector < const SMDS_MeshNode* > smoothNodes; + + if ( !movedByPeriod ) { SMDS_NodeIteratorPtr nIt = smDS->GetNodes(); while ( nIt->more() ) @@ -9797,15 +11416,21 @@ bool _ViscousBuilder::shrink(_SolidData& theData) if ( data._noShrinkShapes.count( subID )) continue; _EdgesOnShape* eos = data.GetShapeEdges( subID ); - if ( !eos || eos->_sWOL.IsNull() ) continue; - + if ( !eos || eos->_sWOL.IsNull() ) + if ( data2 ) // check in adjacent SOLID + { + eos = data2->GetShapeEdges( subID ); + if ( !eos || eos->_sWOL.IsNull() ) + continue; + } subEOS.push_back( eos ); - for ( size_t i = 0; i < eos->_edges.size(); ++i ) - { - lEdges.push_back( eos->_edges[ i ] ); - prepareEdgeToShrink( *eos->_edges[ i ], *eos, helper, smDS ); - } + if ( !movedByPeriod ) + for ( size_t i = 0; i < eos->_edges.size(); ++i ) + { + lEdges.push_back( eos->_edges[ i ] ); + prepareEdgeToShrink( *eos->_edges[ i ], *eos, helper, smDS ); + } } } @@ -9858,7 +11483,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData) { const SMDS_MeshNode* n = smoothNodes[i]; nodesToSmooth[ i ]._node = n; - // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices + // src nodes must be already replaced by tgt nodes to have tgt nodes in _simplices _Simplex::GetSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, 0, sortSimplices); // fix up incorrect uv of nodes on the FACE helper.GetNodeUV( F, n, 0, &isOkUV); @@ -9877,13 +11502,16 @@ bool _ViscousBuilder::shrink(_SolidData& theData) if ( eos.SWOLType() == TopAbs_EDGE ) { SMESH_subMesh* edgeSM = _mesh->GetSubMesh( eos._sWOL ); - _Shrinker1D& srinker = e2shrMap[ edgeSM->GetId() ]; - eShri1D.insert( & srinker ); - srinker.AddEdge( eos._edges[0], eos, helper ); VISCOUS_3D::ToClearSubWithMain( edgeSM, data._solid ); - // restore params of nodes on EGDE if the EDGE has been already - // srinked while srinking other FACE - srinker.RestoreParams(); + if ( !movedByPeriod ) + { + _Shrinker1D& shrinker = e2shrMap[ edgeSM->GetId() ]; + eShri1D.insert( & shrinker ); + shrinker.AddEdge( eos._edges[0], eos, helper ); + // restore params of nodes on EDGE if the EDGE has been already + // shrunk while shrinking other FACE + shrinker.RestoreParams(); + } } for ( size_t i = 0; i < eos._edges.size(); ++i ) { @@ -9898,7 +11526,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData) } bool toFixTria = false; // to improve quality of trias by diagonal swap - if ( isConcaveFace ) + if ( isConcaveFace && !movedByPeriod ) { const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles(); if ( hasTria != hasQuad ) { @@ -9917,24 +11545,24 @@ bool _ViscousBuilder::shrink(_SolidData& theData) // Perform shrinking // ================== - bool shrinked = true; + bool shrunk = !movedByPeriod; int nbBad, shriStep=0, smooStep=0; _SmoothNode::SmoothType smoothType = isConcaveFace ? _SmoothNode::ANGULAR : _SmoothNode::LAPLACIAN; SMESH_Comment errMsg; - while ( shrinked ) + while ( shrunk ) { shriStep++; // Move boundary nodes (actually just set new UV) // ----------------------------------------------- dumpFunction(SMESH_Comment("moveBoundaryOnF")<first<<"_st"<SetNewLength2d( surface, F, eos, helper ); + shrunk |= eos._edges[i]->SetNewLength2d( surface, F, eos, helper ); } } dumpFunctionEnd(); @@ -10025,7 +11653,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData) // } // } - } // while ( shrinked ) + } // while ( shrunk ) if ( !errMsg.empty() ) // Try to re-compute the shrink FACE { @@ -10063,6 +11691,8 @@ bool _ViscousBuilder::shrink(_SolidData& theData) getMeshDS()->RemoveFreeNode( n, smDS, /*fromGroups=*/false ); } } + _periodicity->ClearPeriodic( F ); + // restore position and UV of target nodes gp_Pnt p; for ( size_t iS = 0; iS < subEOS.size(); ++iS ) @@ -10070,17 +11700,18 @@ bool _ViscousBuilder::shrink(_SolidData& theData) { _LayerEdge* edge = subEOS[iS]->_edges[i]; SMDS_MeshNode* tgtNode = const_cast< SMDS_MeshNode*& >( edge->_nodes.back() ); - if ( edge->_pos.empty() ) continue; + if ( edge->_pos.empty() || + edge->Is( _LayerEdge::SHRUNK )) continue; if ( subEOS[iS]->SWOLType() == TopAbs_FACE ) { - SMDS_FacePosition* pos = static_cast( tgtNode->GetPosition() ); + SMDS_FacePositionPtr pos = tgtNode->GetPosition(); pos->SetUParameter( edge->_pos[0].X() ); pos->SetVParameter( edge->_pos[0].Y() ); p = surface->Value( edge->_pos[0].X(), edge->_pos[0].Y() ); } else { - SMDS_EdgePosition* pos = static_cast( tgtNode->GetPosition() ); + SMDS_EdgePositionPtr pos = tgtNode->GetPosition(); pos->SetUParameter( edge->_pos[0].Coord( U_TGT )); p = BRepAdaptor_Curve( TopoDS::Edge( subEOS[iS]->_sWOL )).Value( pos->GetUParameter() ); } @@ -10135,6 +11766,8 @@ bool _ViscousBuilder::shrink(_SolidData& theData) { _EdgesOnShape& eos = * subEOS[ iS ]; if ( eos.ShapeType() != TopAbs_EDGE ) continue; + if ( eos.size() == 0 ) + continue; const TopoDS_Edge& E = TopoDS::Edge( eos._shape ); data.SortOnEdge( E, eos._edges ); @@ -10157,6 +11790,8 @@ bool _ViscousBuilder::shrink(_SolidData& theData) uvPtVec[ i ].param = helper.GetNodeU( E, edges[i]->_nodes[0] ); uvPtVec[ i ].SetUV( helper.GetNodeUV( F, edges[i]->_nodes.back() )); } + // if ( edges.empty() ) + // continue; BRep_Tool::Range( E, uvPtVec[0].param, uvPtVec.back().param ); StdMeshers_FaceSide fSide( uvPtVec, F, E, _mesh ); StdMeshers_ViscousLayers2D::SetProxyMeshOfEdge( fSide ); @@ -10188,7 +11823,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData) return error( errMsg ); } // end of re-meshing in case of failed smoothing - else + else if ( !movedByPeriod ) { // No wrongly shaped faces remain; final smooth. Set node XYZ. bool isStructuredFixed = false; @@ -10226,11 +11861,13 @@ bool _ViscousBuilder::shrink(_SolidData& theData) // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh VISCOUS_3D::ToClearSubWithMain( sm, data._solid ); + if ( data2 ) + VISCOUS_3D::ToClearSubWithMain( sm, data2->_solid ); - } // loop on FACES to srink mesh on + } // loop on FACES to shrink mesh on - // Replace source nodes by target nodes in shrinked mesh edges + // Replace source nodes by target nodes in shrunk mesh edges map< int, _Shrinker1D >::iterator e2shr = e2shrMap.begin(); for ( ; e2shr != e2shrMap.end(); ++e2shr ) @@ -10248,7 +11885,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData) bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, _EdgesOnShape& eos, SMESH_MesherHelper& helper, - const SMESHDS_SubMesh* faceSubMesh) + const SMESHDS_SubMesh* /*faceSubMesh*/) { const SMDS_MeshNode* srcNode = edge._nodes[0]; const SMDS_MeshNode* tgtNode = edge._nodes.back(); @@ -10273,7 +11910,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, edge._pos[0].SetCoord( tgtUV.X(), tgtUV.Y(), 0 ); // set UV of source node to target node - SMDS_FacePosition* pos = static_cast( tgtNode->GetPosition() ); + SMDS_FacePositionPtr pos = tgtNode->GetPosition(); pos->SetUParameter( srcUV.X() ); pos->SetVParameter( srcUV.Y() ); } @@ -10302,9 +11939,17 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, if ( !n2 ) return error(SMESH_Comment("Wrongly meshed EDGE ") << getMeshDS()->ShapeToIndex( E )); + if ( n2 == tgtNode || // for 3D_mesh_GHS3D_01/B1 + n2 == edge._nodes[1] ) // bos #20643 + { + // shrunk by other SOLID + edge.Set( _LayerEdge::SHRUNK ); // ??? + return true; + } + double uSrc = helper.GetNodeU( E, srcNode, n2 ); double uTgt = helper.GetNodeU( E, tgtNode, srcNode ); - double u2 = helper.GetNodeU( E, n2, srcNode ); + double u2 = helper.GetNodeU( E, n2, srcNode ); //edge._pos.clear(); @@ -10323,7 +11968,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, edge._simplices[0]._nPrev = n2; // set U of source node to the target node - SMDS_EdgePosition* pos = static_cast( tgtNode->GetPosition() ); + SMDS_EdgePositionPtr pos = tgtNode->GetPosition(); pos->SetUParameter( uSrc ); } return true; @@ -10356,7 +12001,7 @@ void _ViscousBuilder::restoreNoShrink( _LayerEdge& edge ) const TopLoc_Location loc; Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( S ), loc, f, l ); if ( curve.IsNull() ) return; - SMDS_EdgePosition* ePos = static_cast( srcNode->GetPosition() ); + SMDS_EdgePositionPtr ePos = srcNode->GetPosition(); p = curve->Value( ePos->GetUParameter() ); break; } @@ -10374,7 +12019,7 @@ void _ViscousBuilder::restoreNoShrink( _LayerEdge& edge ) const //================================================================================ /*! - * \brief Try to fix triangles with high aspect ratio by swaping diagonals + * \brief Try to fix triangles with high aspect ratio by swapping diagonals */ //================================================================================ @@ -10557,8 +12202,8 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface, continue; // simplex of quadrangle created by addBoundaryElements() // find intersection of 2 lines: curUV-tgtUV and that connecting simplex nodes - gp_XY uvN1 = helper.GetNodeUV( F, _simplices[i]._nPrev ); - gp_XY uvN2 = helper.GetNodeUV( F, _simplices[i]._nNext ); + gp_XY uvN1 = helper.GetNodeUV( F, _simplices[i]._nPrev, tgtNode ); + gp_XY uvN2 = helper.GetNodeUV( F, _simplices[i]._nNext, tgtNode ); gp_XY dirN = uvN2 - uvN1; double det = uvDir.Crossed( dirN ); if ( Abs( det ) < std::numeric_limits::min() ) continue; @@ -10582,7 +12227,7 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface, { return true; } - SMDS_FacePosition* pos = static_cast( tgtNode->GetPosition() ); + SMDS_FacePositionPtr pos = tgtNode->GetPosition(); pos->SetUParameter( newUV.X() ); pos->SetVParameter( newUV.Y() ); @@ -10590,13 +12235,15 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface, gp_Pnt p = surface->Value( newUV.X(), newUV.Y() ); tgtNode->setXYZ( p.X(), p.Y(), p.Z() ); dumpMove( tgtNode ); +#else + if ( surface.IsNull() ) {} #endif } else // _sWOL is TopAbs_EDGE { const TopoDS_Edge& E = TopoDS::Edge( eos._sWOL ); const SMDS_MeshNode* n2 = _simplices[0]._nPrev; - SMDS_EdgePosition* tgtPos = static_cast( tgtNode->GetPosition() ); + SMDS_EdgePositionPtr tgtPos = tgtNode->GetPosition(); const double u2 = helper.GetNodeU( E, n2, tgtNode ); const double uSrc = _pos[0].Coord( U_SRC ); @@ -10716,7 +12363,7 @@ bool _SmoothNode::Smooth(int& nbBad, return false; } - SMDS_FacePosition* pos = static_cast( _node->GetPosition() ); + SMDS_FacePositionPtr pos = _node->GetPosition(); pos->SetUParameter( newPos.X() ); pos->SetVParameter( newPos.Y() ); @@ -10736,7 +12383,7 @@ bool _SmoothNode::Smooth(int& nbBad, //================================================================================ /*! - * \brief Computes new UV using angle based smoothing technic + * \brief Computes new UV using angle based smoothing technique */ //================================================================================ @@ -10794,34 +12441,6 @@ gp_XY _SmoothNode::computeAngularPos(vector& uv, return newPos; } -//================================================================================ -/*! - * \brief Delete _SolidData - */ -//================================================================================ - -_SolidData::~_SolidData() -{ - TNode2Edge::iterator n2e = _n2eMap.begin(); - for ( ; n2e != _n2eMap.end(); ++n2e ) - { - _LayerEdge* & e = n2e->second; - if ( e ) - { - delete e->_curvature; - if ( e->_2neibors ) - delete e->_2neibors->_plnNorm; - delete e->_2neibors; - } - delete e; - e = 0; - } - _n2eMap.clear(); - - delete _helper; - _helper = 0; -} - //================================================================================ /*! * \brief Keep a _LayerEdge inflated along the EDGE @@ -10839,7 +12458,7 @@ void _Shrinker1D::AddEdge( const _LayerEdge* e, _done = false; } // check _LayerEdge - if ( e == _edges[0] || e == _edges[1] ) + if ( e == _edges[0] || e == _edges[1] || e->_nodes.size() < 2 ) return; if ( eos.SWOLType() != TopAbs_EDGE ) throw SALOME_Exception(LOCALIZED("Wrong _LayerEdge is added")); @@ -10853,11 +12472,18 @@ void _Shrinker1D::AddEdge( const _LayerEdge* e, double u = helper.GetNodeU( _geomEdge, e->_nodes[0], e->_nodes.back()); _edges[ u < 0.5*(f+l) ? 0 : 1 ] = e; - // Update _nodes + // Check if the nodes are already shrunk by another SOLID const SMDS_MeshNode* tgtNode0 = TgtNode( 0 ); const SMDS_MeshNode* tgtNode1 = TgtNode( 1 ); + _done = (( tgtNode0 && tgtNode0->NbInverseElements( SMDSAbs_Edge ) == 2 ) || + ( tgtNode1 && tgtNode1->NbInverseElements( SMDSAbs_Edge ) == 2 )); + if ( _done ) + _nodes.resize( 1, nullptr ); + + // Update _nodes + if ( _nodes.empty() ) { SMESHDS_SubMesh * eSubMesh = helper.GetMeshDS()->MeshElements( _geomEdge ); @@ -10868,7 +12494,7 @@ void _Shrinker1D::AddEdge( const _LayerEdge* e, GeomAdaptor_Curve aCurve(C, f,l); const double totLen = GCPnts_AbscissaPoint::Length(aCurve, f, l); - int nbExpectNodes = eSubMesh->NbNodes(); + smIdType nbExpectNodes = eSubMesh->NbNodes(); _initU .reserve( nbExpectNodes ); _normPar.reserve( nbExpectNodes ); _nodes .reserve( nbExpectNodes ); @@ -10876,9 +12502,18 @@ void _Shrinker1D::AddEdge( const _LayerEdge* e, while ( nIt->more() ) { const SMDS_MeshNode* node = nIt->next(); + + // skip refinement nodes if ( node->NbInverseElements(SMDSAbs_Edge) == 0 || node == tgtNode0 || node == tgtNode1 ) - continue; // refinement nodes + continue; + bool hasMarkedFace = false; + SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face); + while ( fIt->more() && !hasMarkedFace ) + hasMarkedFace = fIt->next()->isMarked(); + if ( !hasMarkedFace ) + continue; + _nodes.push_back( node ); _initU.push_back( helper.GetNodeU( _geomEdge, node )); double len = GCPnts_AbscissaPoint::Length(aCurve, f, _initU.back()); @@ -10917,6 +12552,8 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper) double f,l; if ( set3D || _done ) { + dumpFunction(SMESH_Comment("shrink1D_E") << helper.GetMeshDS()->ShapeToIndex( _geomEdge )<< + "_F" << helper.GetSubShapeID() ); Handle(Geom_Curve) C = BRep_Tool::Curve(_geomEdge, f,l); GeomAdaptor_Curve aCurve(C, f,l); @@ -10929,17 +12566,18 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper) for ( size_t i = 0; i < _nodes.size(); ++i ) { if ( !_nodes[i] ) continue; - if ( _initU[i] < f || l < _initU[i] ) continue; double len = totLen * _normPar[i]; GCPnts_AbscissaPoint discret( aCurve, len, f ); if ( !discret.IsDone() ) return throw SALOME_Exception(LOCALIZED("GCPnts_AbscissaPoint failed")); double u = discret.Parameter(); - SMDS_EdgePosition* pos = static_cast( _nodes[i]->GetPosition() ); + SMDS_EdgePositionPtr pos = _nodes[i]->GetPosition(); pos->SetUParameter( u ); gp_Pnt p = C->Value( u ); const_cast< SMDS_MeshNode*>( _nodes[i] )->setXYZ( p.X(), p.Y(), p.Z() ); + dumpMove( _nodes[i] ); } + dumpFunctionEnd(); } else { @@ -10952,9 +12590,8 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper) for ( size_t i = 0; i < _nodes.size(); ++i ) { if ( !_nodes[i] ) continue; - if ( _initU[i] < f || l < _initU[i] ) continue; double u = f * ( 1-_normPar[i] ) + l * _normPar[i]; - SMDS_EdgePosition* pos = static_cast( _nodes[i]->GetPosition() ); + SMDS_EdgePositionPtr pos = _nodes[i]->GetPosition(); pos->SetUParameter( u ); } } @@ -10972,7 +12609,7 @@ void _Shrinker1D::RestoreParams() for ( size_t i = 0; i < _nodes.size(); ++i ) { if ( !_nodes[i] ) continue; - SMDS_EdgePosition* pos = static_cast( _nodes[i]->GetPosition() ); + SMDS_EdgePositionPtr pos = _nodes[i]->GetPosition(); pos->SetUParameter( _initU[i] ); } _done = false; @@ -10980,7 +12617,7 @@ void _Shrinker1D::RestoreParams() //================================================================================ /*! - * \brief Replace source nodes by target nodes in shrinked mesh edges + * \brief Replace source nodes by target nodes in shrunk mesh edges */ //================================================================================ @@ -11013,6 +12650,140 @@ void _Shrinker1D::SwapSrcTgtNodes( SMESHDS_Mesh* mesh ) } } +//================================================================================ +/*! + * \brief Setup quadPoints + */ +//================================================================================ + +_Mapper2D::_Mapper2D( const TParam2ColumnMap & param2ColumnMap, const TNode2Edge& n2eMap ) +{ + size_t i, iSize = _quadPoints.iSize = param2ColumnMap.size(); + size_t j, jSize = _quadPoints.jSize = param2ColumnMap.begin()->second.size(); + if ( _quadPoints.iSize < 3 || + _quadPoints.jSize < 3 ) + return; + _quadPoints.uv_grid.resize( iSize * jSize ); + + // set nodes + i = 0; + for ( auto & u_columnNodes : param2ColumnMap ) + { + for ( j = 0; j < u_columnNodes.second.size(); ++j ) + _quadPoints.UVPt( i, j ).node = u_columnNodes.second[ j ]; + ++i; + } + + // compute x parameter on borders + uvPnt( 0, 0 ).x = 0; + uvPnt( 0, jSize-1 ).x = 0; + gp_Pnt p0, pPrev0 = SMESH_NodeXYZ( uvPnt( 0, 0 ).node ); + gp_Pnt p1, pPrev1 = SMESH_NodeXYZ( uvPnt( 0, jSize-1 ).node ); + for ( i = 1; i < iSize; ++i ) + { + p0 = SMESH_NodeXYZ( uvPnt( i, 0 ).node ); + p1 = SMESH_NodeXYZ( uvPnt( i, jSize-1 ).node ); + uvPnt( i, 0 ).x = uvPnt( i-1, 0 ).x + p0.Distance( pPrev0 ); + uvPnt( i, jSize-1 ).x = uvPnt( i-1, jSize-1 ).x + p1.Distance( pPrev1 ); + pPrev0 = p0; + pPrev1 = p1; + } + for ( i = 1; i < iSize-1; ++i ) + { + uvPnt( i, 0 ).x /= uvPnt( iSize-1, 0 ).x; + uvPnt( i, jSize-1 ).x /= uvPnt( iSize-1, jSize-1 ).x; + uvPnt( i, 0 ).y = 0; + uvPnt( i, jSize-1 ).y = 1; + } + + // compute y parameter on borders + uvPnt( 0, 0 ).y = 0; + uvPnt( iSize-1, 0 ).y = 0; + pPrev0 = SMESH_NodeXYZ( uvPnt( 0, 0 ).node ); + pPrev1 = SMESH_NodeXYZ( uvPnt( iSize-1, 0 ).node ); + for ( j = 1; j < jSize; ++j ) + { + p0 = SMESH_NodeXYZ( uvPnt( 0, j ).node ); + p1 = SMESH_NodeXYZ( uvPnt( iSize-1, j ).node ); + uvPnt( 0, j ).y = uvPnt( 0, j-1 ).y + p0.Distance( pPrev0 ); + uvPnt( iSize-1, j ).y = uvPnt( iSize-1, j-1 ).y + p1.Distance( pPrev1 ); + pPrev0 = p0; + pPrev1 = p1; + } + for ( j = 1; j < jSize-1; ++j ) + { + uvPnt( 0, j ).y /= uvPnt( 0, jSize-1 ).y; + uvPnt( iSize-1, j ).y /= uvPnt( iSize-1, jSize-1 ).y; + uvPnt( 0, j ).x = 0; + uvPnt( iSize-1, j ).x = 1; + } + + // compute xy of internal nodes + for ( i = 1; i < iSize-1; ++i ) + { + const double x0 = uvPnt( i, 0 ).x; + const double x1 = uvPnt( i, jSize-1 ).x; + for ( j = 1; j < jSize-1; ++j ) + { + const double y0 = uvPnt( 0, j ).y; + const double y1 = uvPnt( iSize-1, j ).y; + double x = (x0 + y0 * (x1 - x0)) / (1 - (y1 - y0) * (x1 - x0)); + double y = y0 + x * (y1 - y0); + uvPnt( i, j ).x = x; + uvPnt( i, j ).y = y; + } + } + + // replace base nodes with target ones + for ( i = 0; i < iSize; ++i ) + for ( j = 0; j < jSize; ++j ) + { + auto n2e = n2eMap.find( uvPnt( i, j ).node ); + uvPnt( i, j ).node = n2e->second->_nodes.back(); + } + + return; +} + +//================================================================================ +/*! + * \brief Compute positions of nodes of 2D structured mesh using TFI + */ +//================================================================================ + +bool _Mapper2D::ComputeNodePositions() +{ + if ( _quadPoints.uv_grid.empty() ) + return true; + + size_t i, iSize = _quadPoints.iSize; + size_t j, jSize = _quadPoints.jSize; + + SMESH_NodeXYZ a0 ( uvPnt( 0, 0 ).node ); + SMESH_NodeXYZ a1 ( uvPnt( iSize-1, 0 ).node ); + SMESH_NodeXYZ a2 ( uvPnt( iSize-1, jSize-1 ).node ); + SMESH_NodeXYZ a3 ( uvPnt( 0, jSize-1 ).node ); + + for ( i = 1; i < iSize-1; ++i ) + { + SMESH_NodeXYZ p0 ( uvPnt( i, 0 ).node ); + SMESH_NodeXYZ p2 ( uvPnt( i, jSize-1 ).node ); + for ( j = 1; j < jSize-1; ++j ) + { + SMESH_NodeXYZ p1 ( uvPnt( iSize-1, j ).node ); + SMESH_NodeXYZ p3 ( uvPnt( 0, j ).node ); + double x = uvPnt( i, j ).x; + double y = uvPnt( i, j ).y; + + gp_XYZ p = SMESH_MesherHelper::calcTFI( x, y, a0,a1,a2,a3, p0,p1,p2,p3 ); + const_cast< SMDS_MeshNode* >( uvPnt( i, j ).node )->setXYZ( p.X(), p.Y(), p.Z() ); + + dumpMove( uvPnt( i, j ).node ); + } + } + return true; +} + //================================================================================ /*! * \brief Creates 2D and 1D elements on boundaries of new prisms @@ -11033,7 +12804,8 @@ bool _ViscousBuilder::addBoundaryElements(_SolidData& data) for ( int iE = 1; iE <= geomEdges.Extent(); ++iE ) { const TopoDS_Edge& E = TopoDS::Edge( geomEdges(iE)); - if ( data._noShrinkShapes.count( getMeshDS()->ShapeToIndex( E ))) + const TGeomID edgeID = getMeshDS()->ShapeToIndex( E ); + if ( data._noShrinkShapes.count( edgeID )) continue; // Get _LayerEdge's based on E @@ -11058,11 +12830,11 @@ bool _ViscousBuilder::addBoundaryElements(_SolidData& data) const SMDS_MeshNode* tgtN0 = ledges[0]->_nodes.back(); const SMDS_MeshNode* tgtN1 = ledges[1]->_nodes.back(); int nbSharedPyram = 0; - SMDS_ElemIteratorPtr vIt = tgtN0->GetInverseElementIterator(SMDSAbs_Volume); + SMDS_ElemIteratorPtr vIt = tgtN1->GetInverseElementIterator(SMDSAbs_Volume); while ( vIt->more() ) { const SMDS_MeshElement* v = vIt->next(); - nbSharedPyram += int( v->GetNodeIndex( tgtN1 ) >= 0 ); + nbSharedPyram += int( v->GetNodeIndex( tgtN0 ) >= 0 ); } if ( nbSharedPyram > 1 ) continue; // not free border of the pyramid @@ -11082,10 +12854,9 @@ bool _ViscousBuilder::addBoundaryElements(_SolidData& data) // Find out orientation and type of face to create bool reverse = false, isOnFace; - - map< TGeomID, TopoDS_Shape >::iterator e2f = - data._shrinkShape2Shape.find( getMeshDS()->ShapeToIndex( E )); TopoDS_Shape F; + + map< TGeomID, TopoDS_Shape >::iterator e2f = data._shrinkShape2Shape.find( edgeID ); if (( isOnFace = ( e2f != data._shrinkShape2Shape.end() ))) { F = e2f->second.Oriented( TopAbs_FORWARD ); @@ -11095,17 +12866,12 @@ bool _ViscousBuilder::addBoundaryElements(_SolidData& data) if ( helper.IsReversedSubMesh( TopoDS::Face(F) )) reverse = !reverse; } - else + else if ( !data._ignoreFaceIds.count( e2f->first )) { // find FACE with layers sharing E - PShapeIteratorPtr fIt = helper.GetAncestors( E, *_mesh, TopAbs_FACE ); - while ( fIt->more() && F.IsNull() ) - { - const TopoDS_Shape* pF = fIt->next(); - if ( helper.IsSubShape( *pF, data._solid) && - !data._ignoreFaceIds.count( e2f->first )) - F = *pF; - } + PShapeIteratorPtr fIt = helper.GetAncestors( E, *_mesh, TopAbs_FACE, &data._solid ); + if ( fIt->more() ) + F = *( fIt->next() ); } // Find the sub-mesh to add new faces SMESHDS_SubMesh* sm = 0; @@ -11116,18 +12882,44 @@ bool _ViscousBuilder::addBoundaryElements(_SolidData& data) if ( !sm ) return error("error in addBoundaryElements()", data._index); + // Find a proxy sub-mesh of the FACE of an adjacent SOLID, which will use the new boundary + // faces for 3D meshing (PAL23414) + SMESHDS_SubMesh* adjSM = 0; + if ( isOnFace ) + { + const TGeomID faceID = sm->GetID(); + PShapeIteratorPtr soIt = helper.GetAncestors( F, *_mesh, TopAbs_SOLID ); + while ( const TopoDS_Shape* solid = soIt->next() ) + if ( !solid->IsSame( data._solid )) + { + size_t iData = _solids.FindIndex( *solid ) - 1; + if ( iData < _sdVec.size() && + _sdVec[ iData ]._ignoreFaceIds.count( faceID ) && + _sdVec[ iData ]._shrinkShape2Shape.count( edgeID ) == 0 ) + { + SMESH_ProxyMesh::SubMesh* proxySub = + _sdVec[ iData ]._proxyMesh->getFaceSubM( TopoDS::Face( F ), /*create=*/false); + if ( proxySub && proxySub->NbElements() > 0 ) + adjSM = proxySub; + } + } + } + // Make faces const int dj1 = reverse ? 0 : 1; const int dj2 = reverse ? 1 : 0; + vector< const SMDS_MeshElement*> ff; // new faces row + SMESHDS_Mesh* m = getMeshDS(); for ( size_t j = 1; j < ledges.size(); ++j ) { vector< const SMDS_MeshNode*>& nn1 = ledges[j-dj1]->_nodes; vector< const SMDS_MeshNode*>& nn2 = ledges[j-dj2]->_nodes; + ff.resize( std::max( nn1.size(), nn2.size() ), NULL ); if ( nn1.size() == nn2.size() ) { if ( isOnFace ) for ( size_t z = 1; z < nn1.size(); ++z ) - sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] )); + sm->AddElement( ff[z-1] = m->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] )); else for ( size_t z = 1; z < nn1.size(); ++z ) sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z] )); @@ -11136,7 +12928,7 @@ bool _ViscousBuilder::addBoundaryElements(_SolidData& data) { if ( isOnFace ) for ( size_t z = 1; z < nn2.size(); ++z ) - sm->AddElement( getMeshDS()->AddFace( nn1[0], nn2[z-1], nn2[z] )); + sm->AddElement( ff[z-1] = m->AddFace( nn1[0], nn2[z-1], nn2[z] )); else for ( size_t z = 1; z < nn2.size(); ++z ) sm->AddElement( new SMDS_FaceOfNodes( nn1[0], nn2[z-1], nn2[z] )); @@ -11145,11 +12937,19 @@ bool _ViscousBuilder::addBoundaryElements(_SolidData& data) { if ( isOnFace ) for ( size_t z = 1; z < nn1.size(); ++z ) - sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[0], nn1[z] )); + sm->AddElement( ff[z-1] = m->AddFace( nn1[z-1], nn2[0], nn1[z] )); else for ( size_t z = 1; z < nn1.size(); ++z ) sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[0], nn2[z] )); } + + if ( adjSM ) // add faces to a proxy SM of the adjacent SOLID + { + for ( size_t z = 0; z < ff.size(); ++z ) + if ( ff[ z ]) + adjSM->AddElement( ff[ z ]); + ff.clear(); + } } // Make edges @@ -11160,7 +12960,7 @@ bool _ViscousBuilder::addBoundaryElements(_SolidData& data) if ( eos && eos->SWOLType() == TopAbs_EDGE ) { vector< const SMDS_MeshNode*>& nn = edge->_nodes; - if ( nn.size() < 2 || nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() ) + if ( nn.size() < 2 || nn[1]->NbInverseElements( SMDSAbs_Edge ) >= 2 ) continue; helper.SetSubShape( eos->_sWOL ); helper.SetElementsOnShape( true );