X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ViscousLayers.cxx;h=680708d7ea81b1578e264bf8870c206a848ef67a;hp=e6d99c74e66faa45c40ff2f019a26c62c11dd2d3;hb=d9f4b53e489dd5857db264ede6acded7b076c9f1;hpb=04f997252152407f9180e03f0af428ab2ca6f4be diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index e6d99c74e..680708d7e 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-2022 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,9 +98,12 @@ #include #include #include +#include #ifdef _DEBUG_ +#ifndef WIN32 #define __myDEBUG +#endif //#define __NOT_INVALIDATE_BAD_SMOOTH //#define __NODES_AT_POS #endif @@ -200,8 +208,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 +377,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 +387,7 @@ namespace VISCOUS_3D struct _LayerEdge; struct _EdgesOnShape; struct _Smoother1D; + struct _Mapper2D; typedef map< const SMDS_MeshNode*, _LayerEdge*, TIDCompare > TNode2Edge; //-------------------------------------------------------------------------------- @@ -453,6 +447,10 @@ namespace VISCOUS_3D 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, @@ -503,7 +501,7 @@ namespace VISCOUS_3D 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 @@ -580,6 +578,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; } @@ -620,12 +619,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; } @@ -634,9 +651,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; }; //-------------------------------------------------------------------------------- @@ -660,11 +687,14 @@ 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 @@ -678,8 +708,12 @@ 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(); }; //-------------------------------------------------------------------------------- @@ -739,6 +773,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 @@ -755,7 +790,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 @@ -764,7 +799,7 @@ 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; @@ -782,8 +817,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 ); @@ -884,6 +919,7 @@ namespace VISCOUS_3D const double refSign ); }; struct PyDump; + struct Periodicity; //-------------------------------------------------------------------------------- /*! * \brief Builder of viscous layers @@ -907,13 +943,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, @@ -996,15 +1034,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; - PyDump* _pyDump; + int _tmpFaceID; + PyDump* _pyDump; }; //-------------------------------------------------------------------------------- /*! @@ -1078,7 +1117,7 @@ namespace VISCOUS_3D 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); @@ -1092,30 +1131,45 @@ namespace VISCOUS_3D 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 { @@ -1123,17 +1177,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; @@ -1144,10 +1202,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(); } }; @@ -1205,6 +1263,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 @@ -1212,10 +1291,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 @@ -1247,6 +1327,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, @@ -1301,6 +1390,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) @@ -1313,14 +1405,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; @@ -1346,22 +1445,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(); } //-------------------------------------------------------------------------------- @@ -1378,8 +1532,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) @@ -1420,7 +1574,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 ); @@ -1452,7 +1606,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 ); @@ -1472,15 +1626,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 { @@ -1584,7 +1738,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 @@ -1674,8 +1828,8 @@ 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 ostream* py; int theNbPyFunc; @@ -1684,12 +1838,12 @@ namespace VISCOUS_3D PyDump(SMESH_Mesh& m) { int tag = 3 + m.GetId(); const char* fname = "/tmp/viscous.py"; - cout << "execfile('"<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; @@ -2139,7 +2320,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 )) @@ -2175,7 +2356,7 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith) 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 ]; + TopoDS_Shape fWOL = FF[ ignore[0] ? 0 : 1 ]; // FACE w/o layers // add EDGE to maps if ( !fWOL.IsNull()) @@ -2255,7 +2436,7 @@ bool _ViscousBuilder::findFacesWithLayers(const bool onlyWith) } } - // Add to _noShrinkShapes sub-shapes of FACE's that can't be shrinked since + // 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 ) @@ -2480,17 +2661,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(); @@ -2512,8 +2682,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; @@ -2524,34 +2697,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() ) @@ -2585,11 +2743,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 ); @@ -2619,13 +2777,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; @@ -2641,7 +2798,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 @@ -2675,6 +2832,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; @@ -2845,6 +3038,7 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) cnvFace._face = F; cnvFace._normalsFixed = false; cnvFace._isTooCurved = false; + cnvFace._normalsFixedOnBorders = false; double maxCurvature = cnvFace.GetMaxCurvature( data, eof, surfProp, helper ); if ( maxCurvature > 0 ) @@ -3028,13 +3222,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 ); @@ -3082,7 +3276,7 @@ bool _ViscousBuilder::findShapesToSmooth( _SolidData& data ) } - // 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; @@ -3212,7 +3406,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() )); } @@ -3227,12 +3421,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 @@ -3242,6 +3446,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 @@ -3265,6 +3516,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 = @@ -3313,19 +3565,18 @@ void _ViscousBuilder::setShapeData( _EdgesOnShape& eos, { SMESHDS_SubMesh* smDS = sm->GetSubMeshDS(); if ( !smDS ) return; - eos._faceNormals.resize( smDS->NbElements() ); + 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 { @@ -3351,7 +3602,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 ) { @@ -3365,9 +3616,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 ) { @@ -3377,6 +3628,17 @@ bool _EdgesOnShape::GetNormal( const SMDS_MeshElement* face, gp_Vec& norm ) return ok; } +//================================================================================ +/*! + * \brief EdgesOnShape destructor + */ +//================================================================================ + +_EdgesOnShape::~_EdgesOnShape() +{ + delete _edgeSmoother; + delete _mapper2D; +} //================================================================================ /*! @@ -3391,12 +3653,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 @@ -3456,13 +3719,14 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, 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 { @@ -3551,30 +3815,37 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, break; } case TopAbs_VERTEX: { + 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 ) - { - getFaceDir( TopoDS::Face( eos._sWOL ), TopoDS::Vertex( eos._shape ), - node, helper, normOK, &edge._cosin ); - } - else 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 = 1; iF < totalNbFaces; ++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; + 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: @@ -3599,13 +3870,26 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, // 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 ) @@ -3623,7 +3907,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 @@ -3638,10 +3922,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 ); @@ -3669,7 +3957,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 } @@ -3808,7 +4096,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 ); @@ -4008,7 +4296,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 { @@ -4085,7 +4373,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 */ //================================================================================ @@ -4130,7 +4418,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 */ //================================================================================ @@ -4155,7 +4471,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() @@ -4199,8 +4514,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 ) @@ -4229,12 +4547,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; } //================================================================================ @@ -4415,7 +4744,7 @@ 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; + 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 ) @@ -4606,7 +4935,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) if ( eos._edges[i]->_nodes.size() > 1 ) avgThick += Min( 1., eos._edges[i]->_len / shapeTgtThick ); else - avgThick += shapeTgtThick; + avgThick += 1; nbActiveEdges += ( ! eos._edges[i]->Is( _LayerEdge::BLOCKED )); } } @@ -4787,24 +5116,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 ) { @@ -4861,7 +5212,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() ) @@ -5062,7 +5413,7 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, continue; // ignore intersection with intFace of an adjacent FACE - if ( dist > 0.1 * eos._edges[i]->_len ) + if ( dist > 0.01 * eos._edges[i]->_len ) { bool toIgnore = false; if ( eos._toSmooth ) @@ -5098,8 +5449,11 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, // intersection not ignored - if ( toBlockInfaltion && - dist < ( eos._edges[i]->_len * theThickToIntersection )) + double minDist = 0; + if ( eos._edges[i]->_maxLen < 0.99 * eos._hyp.GetTotalThickness() ) // limited length + minDist = eos._edges[i]->_len * theThickToIntersection; + + if ( toBlockInfaltion && dist < minDist ) { 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() ); @@ -5138,7 +5491,9 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, } // loop on data._edgesOnShape if ( !is1stBlocked ) + { dumpFunctionEnd(); + } if ( closestFace && le ) { @@ -5332,14 +5687,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; - offsetMaker.PerformByJoin( eos._shape, -offset, Precision::Confusion() ); + offsetMaker.PerformByJoin( eos._shape, -eos._offsetValue, Precision::Confusion() ); if ( !offsetMaker.IsDone() ) return; TopExp_Explorer fExp( offsetMaker.Shape(), TopAbs_FACE ); @@ -5351,7 +5706,7 @@ void _ViscousBuilder::makeOffsetSurface( _EdgesOnShape& eos, SMESH_MesherHelper& eos._offsetSurf = new ShapeAnalysis_Surface( surf ); } - catch ( Standard_Failure ) + catch ( Standard_Failure& ) { } } @@ -5391,6 +5746,7 @@ 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]; @@ -5402,12 +5758,23 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, 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; @@ -5436,6 +5803,35 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, { 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()" ); + } + } } } @@ -5449,12 +5845,13 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, break; if ( i < eos._edges.size() ) { - dumpFunction(SMESH_Comment("putOnOffsetSurface_S") << 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(); } @@ -5470,7 +5867,7 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, SMESH_subMeshIteratorPtr smIt = eos._subMesh->getDependsOnIterator(/*includeSelf=*/false); while ( smIt->more() ) { - SMESH_subMesh* sm = smIt->next(); + 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; @@ -5479,6 +5876,28 @@ void _ViscousBuilder::putOnOffsetSurface( _EdgesOnShape& eos, } 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 ); + } + } } //================================================================================ @@ -5781,9 +6200,11 @@ 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() ); + } gp_XYZ newUV0( newUV.X(), newUV.Y(), 0 ); @@ -5893,10 +6314,11 @@ 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) } } @@ -5912,10 +6334,10 @@ 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; @@ -6144,6 +6566,7 @@ void _Smoother1D::prepare(_SolidData& data) _edgeDir[1] = getEdgeDir( E, leOnV[1]->_nodes[0], data.GetHelper() ); _leOnV[0]._cosin = Abs( leOnV[0]->_cosin ); _leOnV[1]._cosin = Abs( leOnV[1]->_cosin ); + _leOnV[0]._flags = _leOnV[1]._flags = 0; if ( _eos._sWOL.IsNull() ) // 3D for ( int iEnd = 0; iEnd < 2; ++iEnd ) _leOnV[iEnd]._cosin = Abs( _edgeDir[iEnd].Normalized() * leOnV[iEnd]->_normal ); @@ -6153,7 +6576,7 @@ void _Smoother1D::prepare(_SolidData& data) // divide E to have offset segments with low deflection BRepAdaptor_Curve c3dAdaptor( E ); - const double curDeflect = 0.1; //0.01; // Curvature deflection == |p1p2]*sin(p1p2,p1pM) + 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 ) @@ -6291,8 +6714,18 @@ gp_XYZ _Smoother1D::getNormalNormal( const gp_XYZ & normal, // if ( size == 0 ) // MULTI_NORMAL _LayerEdge // return gp_XYZ( 1e-100, 1e-100, 1e-100 ); - return norm / size; -} + 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; +} //================================================================================ /*! @@ -6303,7 +6736,7 @@ gp_XYZ _Smoother1D::getNormalNormal( const gp_XYZ & normal, void _Smoother1D::offPointsToPython() const { const char* fname = "/tmp/offPoints.py"; - cout << "execfile('"<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() ); } @@ -6518,20 +6951,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; } //================================================================================ @@ -6572,7 +7083,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) @@ -6624,9 +7135,9 @@ void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelp void _ViscousBuilder::limitMaxLenByCurvature( _LayerEdge* e1, _LayerEdge* e2, - _EdgesOnShape& eos1, - _EdgesOnShape& eos2, - const bool isSmoothable ) + _EdgesOnShape& /*eos1*/, + _EdgesOnShape& /*eos2*/, + const bool /*isSmoothable*/ ) { if (( e1->_nodes[0]->GetPosition()->GetDim() != e2->_nodes[0]->GetPosition()->GetDim() ) && @@ -6676,6 +7187,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 ) { @@ -6801,9 +7313,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; @@ -6811,7 +7323,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 ) || @@ -6921,7 +7433,7 @@ void _ViscousBuilder::findEdgesToUpdateNormalNearConvexFace( _ConvexFace & if ( dist < 0.95 * ledge->_maxLen ) { ledge->Set( _LayerEdge::UPD_NORMAL_CONV ); - if ( !ledge->_curvature ) ledge->_curvature = new _Curvature; + if ( !ledge->_curvature ) ledge->_curvature = _Factory::NewCurvature(); ledge->_curvature->_uv.SetCoord( uv.X(), uv.Y() ); edgesToUpdateFound = true; } @@ -6944,7 +7456,7 @@ void _ViscousBuilder::findEdgesToUpdateNormalNearConvexFace( _ConvexFace & bool _ViscousBuilder::updateNormals( _SolidData& data, SMESH_MesherHelper& helper, int stepNb, - double stepSize) + double /*stepSize*/) { updateNormalsOfC1Vertices( data ); @@ -7109,7 +7621,7 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, _LayerEdge* edge = e2neIt->first; _LayerEdge& newEdge = e2neIt->second; _EdgesOnShape* eos = data.GetShapeEdges( edge ); - if ( edge->Is( _LayerEdge::BLOCKED && newEdge._maxLen > edge->_len )) + if ( edge->Is( _LayerEdge::BLOCKED ) && newEdge._maxLen > edge->_len ) continue; // Check if a new _normal is OK: @@ -7142,7 +7654,8 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE, &eos->_sWOL ); while ( const TopoDS_Shape* E = eIt->next() ) { - 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 )); @@ -7284,7 +7797,7 @@ bool _ViscousBuilder::isNewNormalOk( _SolidData& data, //================================================================================ bool _ViscousBuilder::updateNormalsOfSmoothed( _SolidData& data, - SMESH_MesherHelper& helper, + SMESH_MesherHelper& /*helper*/, const int nbSteps, const double stepSize ) { @@ -8247,7 +8760,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 ) { @@ -8261,9 +8774,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; @@ -8736,13 +9249,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; @@ -8771,7 +9284,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]; @@ -8935,7 +9448,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; @@ -8947,7 +9460,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 */ //================================================================================ @@ -9428,6 +9941,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() ); @@ -9471,46 +9986,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; @@ -9519,13 +10001,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 @@ -9610,7 +10136,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); @@ -9618,7 +10144,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 ); @@ -9942,11 +10468,16 @@ bool _ViscousBuilder::refine(_SolidData& data) 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(); @@ -9978,41 +10509,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; @@ -10077,7 +10599,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 @@ -10087,7 +10609,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() ); } @@ -10135,6 +10657,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(); @@ -10165,14 +10692,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 ) { @@ -10185,16 +10718,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 @@ -10203,10 +10738,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 ) { @@ -10225,9 +10763,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 ); } @@ -10235,15 +10773,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 ); } } @@ -10254,6 +10792,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 @@ -10262,16 +10805,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, - "Bad quality 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 @@ -10288,11 +11306,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 )) @@ -10331,7 +11349,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 ) { @@ -10351,15 +11369,24 @@ bool _ViscousBuilder::shrink(_SolidData& theData) 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) vector < const SMDS_MeshNode* > smoothNodes; + + if ( !movedByPeriod ) { SMDS_NodeIteratorPtr nIt = smDS->GetNodes(); while ( nIt->more() ) @@ -10405,11 +11432,12 @@ bool _ViscousBuilder::shrink(_SolidData& theData) } 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 ); + } } } @@ -10481,13 +11509,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 ) { @@ -10502,7 +11533,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 ) { @@ -10521,24 +11552,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(); @@ -10629,7 +11660,7 @@ bool _ViscousBuilder::shrink(_SolidData& theData) // } // } - } // while ( shrinked ) + } // while ( shrunk ) if ( !errMsg.empty() ) // Try to re-compute the shrink FACE { @@ -10667,6 +11698,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 ) @@ -10678,14 +11711,14 @@ bool _ViscousBuilder::shrink(_SolidData& theData) 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() ); } @@ -10740,6 +11773,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 ); @@ -10762,6 +11797,19 @@ 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 ( uvPtVec[ 0 ].node == uvPtVec.back().node && // closed + helper.IsSeamShape( uvPtVec[ 0 ].node->GetShapeID() )) + { + uvPtVec[ 0 ].SetUV( helper.GetNodeUV( F, + edges[0]->_nodes.back(), + edges[1]->_nodes.back() )); + size_t i = edges.size() - 1; + uvPtVec[ i ].SetUV( helper.GetNodeUV( F, + edges[i ]->_nodes.back(), + edges[i-1]->_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 ); @@ -10776,8 +11824,12 @@ bool _ViscousBuilder::shrink(_SolidData& theData) smDS->Clear(); // compute the mesh on the FACE + TopTools_IndexedMapOfShape allowed(1); + allowed.Add( sm->GetSubShape() ); + sm->SetAllowedSubShapes( &allowed ); sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); sm->ComputeStateEngine( SMESH_subMesh::COMPUTE_SUBMESH ); + sm->SetAllowedSubShapes( nullptr ); // re-fill proxy sub-meshes of the FACE for ( size_t i = 0 ; i < _sdVec.size(); ++i ) @@ -10793,7 +11845,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; @@ -10834,10 +11886,10 @@ bool _ViscousBuilder::shrink(_SolidData& theData) 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 ) @@ -10855,7 +11907,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(); @@ -10880,7 +11932,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() ); } @@ -10909,6 +11961,14 @@ 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 ); @@ -10930,7 +11990,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; @@ -10963,7 +12023,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; } @@ -11164,8 +12224,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; @@ -11189,7 +12249,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() ); @@ -11197,13 +12257,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 ); @@ -11323,7 +12385,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() ); @@ -11343,7 +12405,7 @@ bool _SmoothNode::Smooth(int& nbBad, //================================================================================ /*! - * \brief Computes new UV using angle based smoothing technic + * \brief Computes new UV using angle based smoothing technique */ //================================================================================ @@ -11368,17 +12430,16 @@ gp_XY _SmoothNode::computeAngularPos(vector& uv, edgeSize.back() = edgeSize.front(); gp_XY newPos(0,0); - //int nbEdges = 0; - double sumSize = 0; + double sumWgt = 0; for ( size_t i = 1; i < edgeDir.size(); ++i ) { - if ( edgeDir[i-1].X() > 1. ) continue; - int i1 = i-1; + const int i1 = i-1; + if ( edgeDir[i1].X() > 1. ) continue; while ( edgeDir[i].X() > 1. && ++i < edgeDir.size() ); if ( i == edgeDir.size() ) break; gp_XY p = uv[i]; gp_XY norm1( -edgeDir[i1].Y(), edgeDir[i1].X() ); - gp_XY norm2( -edgeDir[i].Y(), edgeDir[i].X() ); + gp_XY norm2( -edgeDir[i ].Y(), edgeDir[i ].X() ); gp_XY bisec = norm1 + norm2; double bisecSize = bisec.Modulus(); if ( bisecSize < numeric_limits::min() ) @@ -11388,47 +12449,19 @@ gp_XY _SmoothNode::computeAngularPos(vector& uv, } bisec /= bisecSize; - gp_XY dirToN = uvToFix - p; - double distToN = dirToN.Modulus(); + gp_XY dirToN = uvToFix - p; + double distToN = bisec * dirToN; if ( bisec * dirToN < 0 ) distToN = -distToN; - newPos += ( p + bisec * distToN ) * ( edgeSize[i1] + edgeSize[i] ); - //++nbEdges; - sumSize += edgeSize[i1] + edgeSize[i]; + double wgt = edgeSize[i1] + edgeSize[i]; + newPos += ( p + bisec * distToN ) * wgt; + sumWgt += wgt; } - newPos /= /*nbEdges * */sumSize; + newPos /= sumWgt; 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 @@ -11460,11 +12493,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 ); @@ -11475,7 +12515,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 ); @@ -11533,6 +12573,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); @@ -11550,11 +12592,13 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper) 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 { @@ -11568,7 +12612,7 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper) { if ( !_nodes[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 ); } } @@ -11586,7 +12630,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; @@ -11594,7 +12638,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 */ //================================================================================ @@ -11627,6 +12671,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