X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ViscousLayers.cxx;h=124b61e1e66e93fdabe4468d1de23cc8247b1e2c;hp=73536f8eb82acd8d75eaa3fe8f0689f663d98659;hb=4cd2499bddcd3da3ec8900fe825bc98669b789b5;hpb=e201ae8a2a7d37a84c195803ee97237be92f2f9b diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index 73536f8eb..124b61e1e 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -1,20 +1,20 @@ -// Copyright (C) 2007-2011 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2012 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 -// License as published by the Free Software Foundation; either -// version 2.1 of the License. +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License. // -// This library is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// Lesser General Public License for more details. +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. // -// You should have received a copy of the GNU Lesser General Public -// License along with this library; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // -// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com // // File : StdMeshers_ViscousLayers.cxx @@ -32,21 +32,34 @@ #include "SMESHDS_Hypothesis.hxx" #include "SMESH_Algo.hxx" #include "SMESH_ComputeError.hxx" +#include "SMESH_ControlsDef.hxx" #include "SMESH_Gen.hxx" #include "SMESH_Group.hxx" #include "SMESH_Mesh.hxx" #include "SMESH_MesherHelper.hxx" +#include "SMESH_ProxyMesh.hxx" #include "SMESH_subMesh.hxx" #include "SMESH_subMeshEventListener.hxx" -#include "SMESH_ProxyMesh.hxx" #include "utilities.h" +#include #include +#include +#include +#include #include +#include +#include +#include #include +#include #include +#include +#include #include +#include +#include #include #include #include @@ -58,18 +71,18 @@ #include #include #include -#include #include #include -#include +#include +#include //#define __myDEBUG using namespace std; //================================================================================ -namespace VISCOUS +namespace VISCOUS_3D { typedef int TGeomID; @@ -106,13 +119,15 @@ namespace VISCOUS //-------------------------------------------------------------------------------- /*! * \brief Listener of events of 3D sub-meshes computed with viscous layers. - * It is used to clear an inferior dim sub-mesh modified by viscous layers + * It is used to clear an inferior dim sub-meshes modified by viscous layers */ - class _SrinkShapeListener : SMESH_subMeshEventListener + class _ShrinkShapeListener : SMESH_subMeshEventListener { - _SrinkShapeListener(): SMESH_subMeshEventListener(/*isDeletable=*/false) {} - static SMESH_subMeshEventListener* Get() { static _SrinkShapeListener l; return &l; } + _ShrinkShapeListener() + : SMESH_subMeshEventListener(/*isDeletable=*/false, + "StdMeshers_ViscousLayers::_ShrinkShapeListener") {} public: + static SMESH_subMeshEventListener* Get() { static _ShrinkShapeListener l; return &l; } virtual void ProcessEvent(const int event, const int eventType, SMESH_subMesh* solidSM, @@ -124,23 +139,6 @@ namespace VISCOUS SMESH_subMeshEventListener::ProcessEvent(event,eventType,solidSM,data,hyp); } } - static void ToClearSubMeshWithSolid( SMESH_subMesh* sm, - const TopoDS_Shape& solid) - { - SMESH_subMesh* solidSM = sm->GetFather()->GetSubMesh( solid ); - SMESH_subMeshEventListenerData* data = solidSM->GetEventListenerData( Get()); - if ( data ) - { - if ( find( data->mySubMeshes.begin(), data->mySubMeshes.end(), sm ) == - data->mySubMeshes.end()) - data->mySubMeshes.push_back( sm ); - } - else - { - data = SMESH_subMeshEventListenerData::MakeData( /*dependent=*/sm ); - sm->SetEventListener( Get(), data, /*whereToListenTo=*/solidSM ); - } - } }; //-------------------------------------------------------------------------------- /*! @@ -150,7 +148,9 @@ namespace VISCOUS */ class _ViscousListener : SMESH_subMeshEventListener { - _ViscousListener(): SMESH_subMeshEventListener(/*isDeletable=*/false) {} + _ViscousListener(): + SMESH_subMeshEventListener(/*isDeletable=*/false, + "StdMeshers_ViscousLayers::_ViscousListener") {} static SMESH_subMeshEventListener* Get() { static _ViscousListener l; return &l; } public: virtual void ProcessEvent(const int event, @@ -188,6 +188,32 @@ namespace VISCOUS } }; + //================================================================================ + /*! + * \brief sets a sub-mesh event listener to clear sub-meshes of sub-shapes of + * the main shape when sub-mesh of the main shape is cleared, + * for example to clear sub-meshes of FACEs when sub-mesh of a SOLID + * is cleared + */ + //================================================================================ + + void ToClearSubWithMain( SMESH_subMesh* sub, const TopoDS_Shape& main) + { + SMESH_subMesh* mainSM = sub->GetFather()->GetSubMesh( main ); + SMESH_subMeshEventListenerData* data = + mainSM->GetEventListenerData( _ShrinkShapeListener::Get()); + if ( data ) + { + if ( find( data->mySubMeshes.begin(), data->mySubMeshes.end(), sub ) == + data->mySubMeshes.end()) + data->mySubMeshes.push_back( sub ); + } + else + { + data = SMESH_subMeshEventListenerData::MakeData( /*dependent=*/sub ); + sub->SetEventListener( _ShrinkShapeListener::Get(), data, /*whereToListenTo=*/mainSM ); + } + } //-------------------------------------------------------------------------------- /*! * \brief Simplex (triangle or tetrahedron) based on 1 (tria) or 2 (tet) nodes of @@ -214,17 +240,22 @@ namespace VISCOUS - M[0][2]*M[1][1]*M[2][0]); return determinant > 1e-100; } - bool IsForward(const gp_XY& tgtUV, - const TopoDS_Face& face, - SMESH_MesherHelper& helper, - const double refSign) const + bool IsForward(const gp_XY& tgtUV, + const SMDS_MeshNode* smoothedNode, + const TopoDS_Face& face, + SMESH_MesherHelper& helper, + const double refSign) const { - gp_XY prevUV = helper.GetNodeUV( face, _nPrev ); - gp_XY nextUV = helper.GetNodeUV( face, _nNext ); + gp_XY prevUV = helper.GetNodeUV( face, _nPrev, smoothedNode ); + gp_XY nextUV = helper.GetNodeUV( face, _nNext, smoothedNode ); gp_Vec2d v1( tgtUV, prevUV ), v2( tgtUV, nextUV ); double d = v1 ^ v2; return d*refSign > 1e-100; } + bool IsNeighbour(const _Simplex& other) const + { + return _nPrev == other._nNext || _nNext == other._nPrev; + } }; //-------------------------------------------------------------------------------- /*! @@ -234,6 +265,7 @@ namespace VISCOUS { double _r; // radius double _k; // factor to correct node smoothed position + double _h2lenRatio; // avgNormProj / (2*avgDist) public: static _Curvature* New( double avgNormProj, double avgDist ) { @@ -244,10 +276,12 @@ namespace VISCOUS c->_r = avgDist * avgDist / avgNormProj; c->_k = avgDist * avgDist / c->_r / c->_r; c->_k *= ( c->_r < 0 ? 1/1.1 : 1.1 ); // not to be too restrictive + c->_h2lenRatio = avgNormProj / ( avgDist + avgDist ); } return c; } double lenDelta(double len) const { return _k * ( _r + len ); } + double lenDeltaByDist(double dist) const { return dist * _h2lenRatio; } }; struct _LayerEdge; //-------------------------------------------------------------------------------- @@ -267,6 +301,10 @@ namespace VISCOUS gp_XYZ* _plnNorm; _2NearEdges() { _nodes[0]=_nodes[1]=0; _plnNorm = 0; } + void reverse() { + std::swap( _nodes[0], _nodes[1] ); + std::swap( _wgt[0], _wgt[1] ); + } }; //-------------------------------------------------------------------------------- /*! @@ -321,6 +359,14 @@ namespace VISCOUS void Copy( _LayerEdge& other, SMESH_MesherHelper& helper ); void SetCosin( double cosin ); }; + struct _LayerEdgeCmp + { + bool operator () (const _LayerEdge* e1, const _LayerEdge* e2) const + { + const bool cmpNodes = ( e1 && e2 && e1->_nodes.size() && e2->_nodes.size() ); + return cmpNodes ? ( e1->_nodes[0]->GetID() < e2->_nodes[0]->GetID()) : ( e1 < e2 ); + } + }; //-------------------------------------------------------------------------------- typedef map< const SMDS_MeshNode*, _LayerEdge*, TIDCompare > TNode2Edge; @@ -353,9 +399,8 @@ namespace VISCOUS // FACE's WOL, srink on which is forbiden due to algo on the adjacent SOLID set< TGeomID > _noShrinkFaces; - // end index in _edges of _LayerEdge's based on EDGE (map key) to - // FACE (maybe NULL) they are inflated along - //map< int, TopoDS_Face > _endEdge2Face; + // to + map< TGeomID,Handle(Geom_Curve)> _edge2curve; // end indices in _edges of _LayerEdge on one shape to smooth vector< int > _endEdgeToSmooth; @@ -368,6 +413,13 @@ namespace VISCOUS const StdMeshers_ViscousLayers* h=0, _MeshOfSolid* m=0) :_solid(s), _hyp(h), _proxyMesh(m) {} ~_SolidData(); + + Handle(Geom_Curve) CurveForSmooth( const TopoDS_Edge& E, + const int iFrom, + const int iTo, + Handle(Geom_Surface)& surface, + const TopoDS_Face& F, + SMESH_MesherHelper& helper); }; //-------------------------------------------------------------------------------- /*! @@ -383,6 +435,7 @@ namespace VISCOUS Handle(Geom_Surface)& surface, SMESH_MesherHelper& helper, const double refSign, + bool isCentroidal, bool set3D); }; //-------------------------------------------------------------------------------- @@ -416,7 +469,8 @@ namespace VISCOUS _SolidData& data); void getSimplices( const SMDS_MeshNode* node, vector<_Simplex>& simplices, const set& ingnoreShapes, - const _SolidData* dataToCheckOri = 0); + const _SolidData* dataToCheckOri = 0, + const bool toSort = false); bool sortEdges( _SolidData& data, vector< vector<_LayerEdge*> >& edgesByGeom); void limitStepSize( _SolidData& data, @@ -424,13 +478,20 @@ namespace VISCOUS const double cosin); void limitStepSize( _SolidData& data, const double minSize); bool inflate(_SolidData& data); - bool smoothAndCheck(_SolidData& data, int nbSteps, double & distToIntersection); + bool smoothAndCheck(_SolidData& data, const int nbSteps, double & distToIntersection); + bool smoothAnalyticEdge( _SolidData& data, + const int iFrom, + const int iTo, + Handle(Geom_Surface)& surface, + const TopoDS_Face& F, + SMESH_MesherHelper& helper); bool updateNormals( _SolidData& data, SMESH_MesherHelper& helper ); bool refine(_SolidData& data); bool shrink(); bool prepareEdgeToShrink( _LayerEdge& edge, const TopoDS_Face& F, SMESH_MesherHelper& helper, const SMESHDS_SubMesh* faceSubMesh ); + void fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& helper); bool addBoundaryElements(); bool error( const string& text, int solidID=-1 ); @@ -478,7 +539,8 @@ namespace VISCOUS virtual SMDSAbs_ElementType GetType() const { return SMDSAbs_Face; } virtual vtkIdType GetVtkType() const { return -1; } virtual SMDSAbs_EntityType GetEntityType() const { return SMDSEntity_Last; } - virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const + virtual SMDSAbs_GeometryType GetGeomType() const { return SMDSGeom_TRIANGLE; } +virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const { return SMDS_ElemIteratorPtr( new SMDS_NodeVectorElemIterator( _nn.begin(), _nn.end()));} }; //-------------------------------------------------------------------------------- @@ -497,7 +559,7 @@ namespace VISCOUS _nn[3]=_le2->_nodes[0]; } }; -} // namespace VISCOUS +} // namespace VISCOUS_3D //================================================================================ // StdMeshers_ViscousLayers hypothesis @@ -509,10 +571,15 @@ StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH _name = StdMeshers_ViscousLayers::GetHypType(); _param_algo_dim = -3; // auxiliary hyp used by 3D algos } // -------------------------------------------------------------------------------- -void StdMeshers_ViscousLayers::SetIgnoreFaces(const std::vector& faceIds) +void StdMeshers_ViscousLayers::SetBndShapesToIgnore(const std::vector& faceIds) { - if ( faceIds != _ignoreFaceIds ) - _ignoreFaceIds = faceIds, NotifySubMeshesHypothesisModification(); + if ( faceIds != _ignoreBndShapeIds ) + _ignoreBndShapeIds = faceIds, NotifySubMeshesHypothesisModification(); +} // -------------------------------------------------------------------------------- +bool StdMeshers_ViscousLayers::IsIgnoredShape(const int shapeID) const +{ + return ( find( _ignoreBndShapeIds.begin(), _ignoreBndShapeIds.end(), shapeID ) + != _ignoreBndShapeIds.end() ); } // -------------------------------------------------------------------------------- void StdMeshers_ViscousLayers::SetTotalThickness(double thickness) { @@ -534,7 +601,7 @@ StdMeshers_ViscousLayers::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theShape, const bool toMakeN2NMap) const { - using namespace VISCOUS; + using namespace VISCOUS_3D; _ViscousBuilder bulder; SMESH_ComputeErrorPtr err = bulder.Compute( theMesh, theShape ); if ( err && !err->IsOK() ) @@ -570,17 +637,17 @@ std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save) save << " " << _nbLayers << " " << _thickness << " " << _stretchFactor - << " " << _ignoreFaceIds.size(); - for ( unsigned i = 0; i < _ignoreFaceIds.size(); ++i ) - save << " " << _ignoreFaceIds[i]; + << " " << _ignoreBndShapeIds.size(); + for ( unsigned i = 0; i < _ignoreBndShapeIds.size(); ++i ) + save << " " << _ignoreBndShapeIds[i]; return save; } // -------------------------------------------------------------------------------- std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load) { int nbFaces, faceID; load >> _nbLayers >> _thickness >> _stretchFactor >> nbFaces; - while ( _ignoreFaceIds.size() < nbFaces && load >> faceID ) - _ignoreFaceIds.push_back( faceID ); + while ( _ignoreBndShapeIds.size() < nbFaces && load >> faceID ) + _ignoreBndShapeIds.push_back( faceID ); return load; } // -------------------------------------------------------------------------------- bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh* theMesh, @@ -660,38 +727,102 @@ namespace while ( eIt->more()) { const TopoDS_Edge* e = static_cast( eIt->next() ); - if ( helper.IsSubShape( *e, F ) && BRep_Tool::Curve( *e, loc,f,l)) + if ( helper.IsSubShape( *e, F ) && !BRep_Tool::Curve( *e, loc,f,l).IsNull() ) edges.push_back( *e ); } gp_XYZ dir(0,0,0); if ( !( ok = ( edges.size() > 0 ))) return dir; // get average dir of edges going fromV - gp_Vec edgeDir; - for ( unsigned i = 0; i < edges.size(); ++i ) - { - edgeDir = getEdgeDir( edges[i], fromV ); - double size2 = edgeDir.SquareMagnitude(); - if ( size2 > numeric_limits::min() ) - edgeDir /= sqrt( size2 ); - else - ok = false; - dir += edgeDir.XYZ(); - } + gp_XYZ edgeDir; + //if ( edges.size() > 1 ) + for ( unsigned i = 0; i < edges.size(); ++i ) + { + edgeDir = getEdgeDir( edges[i], fromV ); + double size2 = edgeDir.SquareModulus(); + if ( size2 > numeric_limits::min() ) + edgeDir /= sqrt( size2 ); + else + ok = false; + dir += edgeDir; + } gp_XYZ fromEdgeDir = getFaceDir( F, edges[0], node, helper, ok ); - if ( edges.size() == 1 || dir.SquareModulus() < 1e-10) + if ( edges.size() == 1 ) dir = fromEdgeDir; + else if ( dir.SquareModulus() < 0.1 ) // ~< 20 degrees + dir = fromEdgeDir + getFaceDir( F, edges[1], node, helper, ok ); else if ( dir * fromEdgeDir < 0 ) dir *= -1; if ( ok ) { //dir /= edges.size(); if ( cosin ) { - double angle = edgeDir.Angle( dir ); + double angle = gp_Vec( edgeDir ).Angle( dir ); *cosin = cos( angle ); } } return dir; } + //================================================================================ + /*! + * \brief Returns true if a FACE is bound by a concave EDGE + */ + //================================================================================ + + bool isConcave( const TopoDS_Face& F, SMESH_MesherHelper& helper ) + { + gp_Vec2d drv1, drv2; + gp_Pnt2d p; + TopExp_Explorer eExp( F.Oriented( TopAbs_FORWARD ), TopAbs_EDGE ); + for ( ; eExp.More(); eExp.Next() ) + { + const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() ); + if ( BRep_Tool::Degenerated( E )) continue; + // check if 2D curve is concave + BRepAdaptor_Curve2d curve( E, F ); + const int nbIntervals = curve.NbIntervals( GeomAbs_C2 ); + TColStd_Array1OfReal intervals(1, nbIntervals + 1 ); + curve.Intervals( intervals, GeomAbs_C2 ); + bool isConvex = true; + for ( int i = 1; i <= nbIntervals && isConvex; ++i ) + { + double u1 = intervals( i ); + double u2 = intervals( i+1 ); + curve.D2( 0.5*( u1+u2 ), p, drv1, drv2 ); + double cross = drv2 ^ drv1; + if ( E.Orientation() == TopAbs_REVERSED ) + cross = -cross; + isConvex = ( cross < 1e-9 ); + } + // check if concavity is strong enough to care about it + //const double maxAngle = 5 * Standard_PI180; + if ( !isConvex ) + { + //cout << "Concave FACE " << helper.GetMeshDS()->ShapeToIndex( F ) << endl; + return true; + // map< double, const SMDS_MeshNode* > u2nodes; + // if ( !SMESH_Algo::GetSortedNodesOnEdge( helper.GetMeshDS(), E, + // /*ignoreMedium=*/true, u2nodes)) + // continue; + // map< double, const SMDS_MeshNode* >::iterator u2n = u2nodes.begin(); + // gp_Pnt2d uvPrev = helper.GetNodeUV( F, u2n->second ); + // double uPrev = u2n->first; + // for ( ++u2n; u2n != u2nodes.end(); ++u2n ) + // { + // gp_Pnt2d uv = helper.GetNodeUV( F, u2n->second ); + // gp_Vec2d segmentDir( uvPrev, uv ); + // curve.D1( uPrev, p, drv1 ); + // try { + // if ( fabs( segmentDir.Angle( drv1 )) > maxAngle ) + // return true; + // } + // catch ( ... ) {} + // uvPrev = uv; + // uPrev = u2n->first; + // } + } + } + return false; + } //-------------------------------------------------------------------------------- // DEBUG. Dump intermediate node positions into a python script #ifdef __myDEBUG @@ -703,35 +834,42 @@ namespace py = new ofstream(fname); *py << "from smesh import *" << endl << "meshSO = GetCurrentStudy().FindObjectID('0:1:2:3')" << endl - << "mesh = Mesh( meshSO.GetObject()._narrow( SMESH.SMESH_Mesh ))"<GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE || curN->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) { - double dist = SMESH_MeshEditor::TNodeXYZ( face->GetNode(i)).Distance( nextN ); + double dist = SMESH_TNodeXYZ( face->GetNode(i)).Distance( nextN ); if ( dist < minSize ) minSize = dist, iN = i; } @@ -1308,7 +1473,7 @@ void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize) if ( data._stepSizeNodes[0] ) { double dist = - SMESH_MeshEditor::TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]); + SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]); data._stepSizeCoeff = data._stepSize / dist; } } @@ -1329,7 +1494,7 @@ bool _ViscousBuilder::sortEdges( _SolidData& data, list< TGeomID > shapesToSmooth; SMESH_MesherHelper helper( *_mesh ); - bool ok; + bool ok = true; for ( unsigned iS = 0; iS < edgesByGeom.size(); ++iS ) { @@ -1550,8 +1715,8 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, } case SMDS_TOP_VERTEX: { TopoDS_Vertex V = TopoDS::Vertex( helper.GetSubShapeByNode( node, getMeshDS())); - gp_Vec inFaceDir = getFaceDir( F, V, node, helper, normOK); - double angle = inFaceDir.Angle( edge._normal ); // [0,PI] + gp_XYZ inFaceDir = getFaceDir( F, V, node, helper, normOK); + double angle = gp_Vec( inFaceDir).Angle( edge._normal ); // [0,PI] edge._cosin = cos( angle ); //cout << "Cosin on VERTEX " << edge._cosin << " node " << node->GetID() << endl; break; @@ -1595,7 +1760,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, } else { - edge._pos.push_back( SMESH_MeshEditor::TNodeXYZ( node )); + edge._pos.push_back( SMESH_TNodeXYZ( node )); if ( posType == SMDS_TOP_FACE ) { @@ -1603,7 +1768,7 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, double avgNormProj = 0, avgLen = 0; for ( unsigned i = 0; i < edge._simplices.size(); ++i ) { - gp_XYZ vec = edge._pos.back() - SMESH_MeshEditor::TNodeXYZ( edge._simplices[i]._nPrev ); + gp_XYZ vec = edge._pos.back() - SMESH_TNodeXYZ( edge._simplices[i]._nPrev ); avgNormProj += edge._normal * vec; avgLen += vec.Modulus(); } @@ -1693,9 +1858,9 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1, if ( _nodes[0]->GetPosition()->GetTypeOfPosition() != SMDS_TOP_EDGE ) return; - gp_XYZ pos = SMESH_MeshEditor::TNodeXYZ( _nodes[0] ); - gp_XYZ vec1 = pos - SMESH_MeshEditor::TNodeXYZ( n1 ); - gp_XYZ vec2 = pos - SMESH_MeshEditor::TNodeXYZ( n2 ); + gp_XYZ pos = SMESH_TNodeXYZ( _nodes[0] ); + gp_XYZ vec1 = pos - SMESH_TNodeXYZ( n1 ); + gp_XYZ vec2 = pos - SMESH_TNodeXYZ( n2 ); // Set _curvature @@ -1783,14 +1948,14 @@ void _LayerEdge::SetCosin( double cosin ) void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node, vector<_Simplex>& simplices, const set& ingnoreShapes, - const _SolidData* dataToCheckOri) + const _SolidData* dataToCheckOri, + const bool toSort) { - SMESH_MeshEditor editor( _mesh ); SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face); while ( fIt->more() ) { const SMDS_MeshElement* f = fIt->next(); - const TGeomID shapeInd = editor.FindShape( f ); + const TGeomID shapeInd = f->getshapeId(); if ( ingnoreShapes.count( shapeInd )) continue; const int nbNodes = f->NbCornerNodes(); int srcInd = f->GetNodeIndex( node ); @@ -1800,7 +1965,25 @@ void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node, std::swap( nPrev, nNext ); simplices.push_back( _Simplex( nPrev, nNext )); } - simplices.resize( simplices.size() ); + + if ( toSort ) + { + vector<_Simplex> sortedSimplices( simplices.size() ); + sortedSimplices[0] = simplices[0]; + int nbFound = 0; + for ( size_t i = 1; i < simplices.size(); ++i ) + { + for ( size_t j = 1; j < simplices.size(); ++j ) + if ( sortedSimplices[i-1]._nNext == simplices[j]._nPrev ) + { + sortedSimplices[i] = simplices[j]; + nbFound++; + break; + } + } + if ( nbFound == simplices.size() - 1 ) + simplices.swap( sortedSimplices ); + } } //================================================================================ @@ -1836,7 +2019,7 @@ void _ViscousBuilder::makeGroupOfLE() for ( unsigned j = 0 ; j < _sdVec[i]._edges.size(); ++j ) { _LayerEdge& edge = *_sdVec[i]._edges[j]; - SMESH_MeshEditor::TNodeXYZ nXYZ( edge._nodes[0] ); + SMESH_TNodeXYZ nXYZ( edge._nodes[0] ); nXYZ += edge._normal * _sdVec[i]._stepSize; dumpCmd(SMESH_Comment("mesh.AddEdge([ ") <GetID() << ", mesh.AddNode( " << nXYZ.X()<<","<< nXYZ.Y()<<","<< nXYZ.Z()<<")])"); @@ -1891,7 +2074,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) { if ( data._edges[i]->IsOnEdge() ) continue; data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon ); - if ( geomSize > intersecDist ) + if ( geomSize > intersecDist && intersecDist > 0 ) geomSize = intersecDist; } if ( data._stepSize > 0.3 * geomSize ) @@ -1969,7 +2152,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) limitStepSize( data, 0.25 * distToIntersection ); if ( data._stepSizeNodes[0] ) data._stepSize = data._stepSizeCoeff * - SMESH_MeshEditor::TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]); + SMESH_TNodeXYZ(data._stepSizeNodes[0]).Distance(data._stepSizeNodes[1]); } if (nbSteps == 0 ) @@ -1985,7 +2168,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) //================================================================================ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, - int nbSteps, + const int nbSteps, double & distToIntersection) { if ( data._endEdgeToSmooth.empty() ) @@ -2019,20 +2202,25 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, TGeomID sInd = data._edges[ iBeg ]->_nodes[0]->getshapeId(); if ( data._edges[ iBeg ]->IsOnEdge() ) - { + { dumpFunction(SMESH_Comment("smooth")<SmoothOnEdge(surface, F, helper); + + // try a simple solution on an analytic EDGE + if ( !smoothAnalyticEdge( data, iBeg, iEnd, surface, F, helper )) + { + // smooth on EDGE's + int step = 0; + do { + moved = false; + for ( int i = iBeg; i < iEnd; ++i ) + { + moved |= data._edges[i]->SmoothOnEdge(surface, F, helper); + } + dumpCmd( SMESH_Comment("# end step ")<_nodes.back() ); + SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() ); for ( unsigned j = 0; j < edge->_simplices.size(); ++j ) if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ )) { @@ -2083,14 +2271,23 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, distToIntersection = Precision::Infinite(); double dist; - const SMDS_MeshElement* intFace = 0, *closestFace = 0; + const SMDS_MeshElement* intFace = 0; +#ifdef __myDEBUG + const SMDS_MeshElement* closestFace = 0; int iLE = 0; +#endif for ( unsigned i = 0; i < data._edges.size(); ++i ) { if ( data._edges[i]->FindIntersection( *searcher, dist, data._epsilon, &intFace )) return false; if ( distToIntersection > dist ) - distToIntersection = dist, closestFace = intFace, iLE = i; + { + distToIntersection = dist; +#ifdef __myDEBUG + iLE = i; + closestFace = intFace; +#endif + } } #ifdef __myDEBUG if ( closestFace ) @@ -2106,6 +2303,267 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, return true; } +//================================================================================ +/*! + * \brief Return a curve of the EDGE to be used for smoothing and arrange + * _LayerEdge's to be in a consequent order + */ +//================================================================================ + +Handle(Geom_Curve) _SolidData::CurveForSmooth( const TopoDS_Edge& E, + const int iFrom, + const int iTo, + Handle(Geom_Surface)& surface, + const TopoDS_Face& F, + SMESH_MesherHelper& helper) +{ + TGeomID eIndex = helper.GetMeshDS()->ShapeToIndex( E ); + + map< TGeomID, Handle(Geom_Curve)>::iterator i2curve = _edge2curve.find( eIndex ); + + if ( i2curve == _edge2curve.end() ) + { + // sort _LayerEdge's by position on the EDGE + { + map< double, _LayerEdge* > u2edge; + for ( int i = iFrom; i < iTo; ++i ) + u2edge.insert( make_pair( helper.GetNodeU( E, _edges[i]->_nodes[0] ), _edges[i] )); + + ASSERT( u2edge.size() == iTo - iFrom ); + map< double, _LayerEdge* >::iterator u2e = u2edge.begin(); + for ( int i = iFrom; i < iTo; ++i, ++u2e ) + _edges[i] = u2e->second; + + // set _2neibors according to the new order + for ( int i = iFrom; i < iTo-1; ++i ) + if ( _edges[i]->_2neibors->_nodes[1] != _edges[i+1]->_nodes.back() ) + _edges[i]->_2neibors->reverse(); + if ( u2edge.size() > 1 && + _edges[iTo-1]->_2neibors->_nodes[0] != _edges[iTo-2]->_nodes.back() ) + _edges[iTo-1]->_2neibors->reverse(); + } + + SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( eIndex ); + + TopLoc_Location loc; double f,l; + + Handle(Geom_Line) line; + Handle(Geom_Circle) circle; + bool isLine, isCirc; + if ( F.IsNull() ) // 3D case + { + // check if the EDGE is a line + Handle(Geom_Curve) curve = BRep_Tool::Curve( E, loc, f, l); + if ( curve->IsKind( STANDARD_TYPE( Geom_TrimmedCurve ))) + curve = Handle(Geom_TrimmedCurve)::DownCast( curve )->BasisCurve(); + + line = Handle(Geom_Line)::DownCast( curve ); + circle = Handle(Geom_Circle)::DownCast( curve ); + isLine = (!line.IsNull()); + isCirc = (!circle.IsNull()); + + if ( !isLine && !isCirc ) // Check if the EDGE is close to a line + { + Bnd_B3d bndBox; + SMDS_NodeIteratorPtr nIt = smDS->GetNodes(); + while ( nIt->more() ) + bndBox.Add( SMESH_TNodeXYZ( nIt->next() )); + gp_XYZ size = bndBox.CornerMax() - bndBox.CornerMin(); + + SMESH_TNodeXYZ p0( _edges[iFrom]->_2neibors->_nodes[0] ); + SMESH_TNodeXYZ p1( _edges[iFrom]->_2neibors->_nodes[1] ); + const double lineTol = 1e-2 * ( p0 - p1 ).Modulus(); + for ( int i = 0; i < 3 && !isLine; ++i ) + isLine = ( size.Coord( i+1 ) <= lineTol ); + } + if ( !isLine && !isCirc && iTo-iFrom > 2) // Check if the EDGE is close to a circle + { + // TODO + } + } + else // 2D case + { + // check if the EDGE is a line + Handle(Geom2d_Curve) curve = BRep_Tool::CurveOnSurface( E, F, f, l); + if ( curve->IsKind( STANDARD_TYPE( Geom2d_TrimmedCurve ))) + curve = Handle(Geom2d_TrimmedCurve)::DownCast( curve )->BasisCurve(); + + Handle(Geom2d_Line) line2d = Handle(Geom2d_Line)::DownCast( curve ); + Handle(Geom2d_Circle) circle2d = Handle(Geom2d_Circle)::DownCast( curve ); + isLine = (!line2d.IsNull()); + isCirc = (!circle2d.IsNull()); + + if ( !isLine && !isCirc) // Check if the EDGE is close to a line + { + Bnd_B2d bndBox; + SMDS_NodeIteratorPtr nIt = smDS->GetNodes(); + while ( nIt->more() ) + bndBox.Add( helper.GetNodeUV( F, nIt->next() )); + gp_XY size = bndBox.CornerMax() - bndBox.CornerMin(); + + const double lineTol = 1e-2 * sqrt( bndBox.SquareExtent() ); + for ( int i = 0; i < 2 && !isLine; ++i ) + isLine = ( size.Coord( i+1 ) <= lineTol ); + } + if ( !isLine && !isCirc && iTo-iFrom > 2) // Check if the EDGE is close to a circle + { + // TODO + } + if ( isLine ) + { + line = new Geom_Line( gp::OX() ); // only type does matter + } + else if ( isCirc ) + { + gp_Pnt2d p = circle2d->Location(); + gp_Ax2 ax( gp_Pnt( p.X(), p.Y(), 0), gp::DX()); + circle = new Geom_Circle( ax, 1.); // only center position does matter + } + } + + Handle(Geom_Curve)& res = _edge2curve[ eIndex ]; + if ( isLine ) + res = line; + else if ( isCirc ) + res = circle; + + return res; + } + return i2curve->second; +} + +//================================================================================ +/*! + * \brief smooth _LayerEdge's on a staight EDGE or circular EDGE + */ +//================================================================================ + +bool _ViscousBuilder::smoothAnalyticEdge( _SolidData& data, + const int iFrom, + const int iTo, + Handle(Geom_Surface)& surface, + const TopoDS_Face& F, + SMESH_MesherHelper& helper) +{ + TopoDS_Shape S = helper.GetSubShapeByNode( data._edges[ iFrom ]->_nodes[0], + helper.GetMeshDS()); + TopoDS_Edge E = TopoDS::Edge( S ); + + Handle(Geom_Curve) curve = data.CurveForSmooth( E, iFrom, iTo, surface, F, helper ); + if ( curve.IsNull() ) return false; + + // compute a relative length of segments + vector< double > len( iTo-iFrom+1 ); + { + double curLen, prevLen = len[0] = 1.0; + for ( int i = iFrom; i < iTo; ++i ) + { + curLen = prevLen * data._edges[i]->_2neibors->_wgt[0] / data._edges[i]->_2neibors->_wgt[1]; + len[i-iFrom+1] = len[i-iFrom] + curLen; + prevLen = curLen; + } + } + + if ( curve->IsKind( STANDARD_TYPE( Geom_Line ))) + { + if ( F.IsNull() ) // 3D + { + SMESH_TNodeXYZ p0( data._edges[iFrom]->_2neibors->_nodes[0]); + SMESH_TNodeXYZ p1( data._edges[iTo-1]->_2neibors->_nodes[1]); + for ( int i = iFrom; i < iTo; ++i ) + { + double r = len[i-iFrom] / len.back(); + gp_XYZ newPos = p0 * ( 1. - r ) + p1 * r; + data._edges[i]->_pos.back() = newPos; + SMDS_MeshNode* tgtNode = const_cast( data._edges[i]->_nodes.back() ); + tgtNode->setXYZ( newPos.X(), newPos.Y(), newPos.Z() ); + dumpMove( tgtNode ); + } + } + else + { + gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]); + gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]); + if ( data._edges[iFrom]->_2neibors->_nodes[0] == + data._edges[iTo-1]->_2neibors->_nodes[1] ) // closed edge + { + int iPeriodic = helper.GetPeriodicIndex(); + if ( iPeriodic == 1 || iPeriodic == 2 ) + { + uv1.SetCoord( iPeriodic, helper.GetOtherParam( uv1.Coord( iPeriodic ))); + if ( uv0.Coord( iPeriodic ) > uv1.Coord( iPeriodic )) + std::swap( uv0, uv1 ); + } + } + const gp_XY rangeUV = uv1 - uv0; + for ( int i = iFrom; i < iTo; ++i ) + { + double r = len[i-iFrom] / len.back(); + gp_XY newUV = uv0 + r * rangeUV; + data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 ); + + gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() ); + SMDS_MeshNode* tgtNode = const_cast( data._edges[i]->_nodes.back() ); + 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() ); + } + } + return true; + } + + if ( curve->IsKind( STANDARD_TYPE( Geom_Circle ))) + { + Handle(Geom_Circle) circle = Handle(Geom_Circle)::DownCast( curve ); + gp_Pnt center3D = circle->Location(); + + if ( F.IsNull() ) // 3D + { + return false; // TODO ??? + } + else // 2D + { + const gp_XY center( center3D.X(), center3D.Y() ); + + gp_XY uv0 = helper.GetNodeUV( F, data._edges[iFrom]->_2neibors->_nodes[0]); + gp_XY uvM = helper.GetNodeUV( F, data._edges[iFrom]->_nodes.back()); + gp_XY uv1 = helper.GetNodeUV( F, data._edges[iTo-1]->_2neibors->_nodes[1]); + gp_Vec2d vec0( center, uv0 ); + gp_Vec2d vecM( center, uvM ); + gp_Vec2d vec1( center, uv1 ); + double uLast = vec0.Angle( vec1 ); // -PI - +PI + double uMidl = vec0.Angle( vecM ); + if ( uLast * uMidl < 0. ) + uLast += ( uMidl > 0 ? +2. : -2. ) * M_PI; + const double radius = 0.5 * ( vec0.Magnitude() + vec1.Magnitude() ); + + gp_Ax2d axis( center, vec0 ); + gp_Circ2d circ( axis, radius ); + for ( int i = iFrom; i < iTo; ++i ) + { + double newU = uLast * len[i-iFrom] / len.back(); + gp_Pnt2d newUV = ElCLib::Value( newU, circ ); + data._edges[i]->_pos.back().SetCoord( newUV.X(), newUV.Y(), 0 ); + + gp_Pnt newPos = surface->Value( newUV.X(), newUV.Y() ); + SMDS_MeshNode* tgtNode = const_cast( data._edges[i]->_nodes.back() ); + 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() ); + } + } + return true; + } + + return false; +} + //================================================================================ /*! * \brief Modify normals of _LayerEdge's on EDGE's to avoid intersection with @@ -2177,7 +2635,8 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, // 1) Find intersections double dist; const SMDS_MeshElement* face; - map< _LayerEdge*, set< _LayerEdge* > > edge2CloseEdge; + typedef map< _LayerEdge*, set< _LayerEdge*, _LayerEdgeCmp >, _LayerEdgeCmp > TLEdge2LEdgeSet; + TLEdge2LEdgeSet edge2CloseEdge; const double eps = data._epsilon * data._epsilon; for ( unsigned i = 0; i < data._edges.size(); ++i ) @@ -2187,7 +2646,7 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, if ( edge->FindIntersection( *searcher, dist, eps, &face )) { const TmpMeshFaceOnEdge* f = (const TmpMeshFaceOnEdge*) face; - set< _LayerEdge* > & ee = edge2CloseEdge[ edge ]; + set< _LayerEdge*, _LayerEdgeCmp > & ee = edge2CloseEdge[ edge ]; ee.insert( f->_le1 ); ee.insert( f->_le2 ); if ( f->_le1->IsOnEdge() && f->_le1->_sWOL.IsNull() ) @@ -2203,12 +2662,12 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, { dumpFunction(SMESH_Comment("updateNormals")< >::iterator e2ee = edge2CloseEdge.begin(); + TLEdge2LEdgeSet::iterator e2ee = edge2CloseEdge.begin(); for ( ; e2ee != edge2CloseEdge.end(); ++e2ee ) { _LayerEdge* edge1 = e2ee->first; _LayerEdge* edge2 = 0; - set< _LayerEdge* >& ee = e2ee->second; + set< _LayerEdge*, _LayerEdgeCmp >& ee = e2ee->second; // find EDGEs the edges reside TopoDS_Edge E1, E2; @@ -2216,7 +2675,7 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, if ( S.ShapeType() != TopAbs_EDGE ) continue; // TODO: find EDGE by VERTEX E1 = TopoDS::Edge( S ); - set< _LayerEdge* >::iterator eIt = ee.begin(); + set< _LayerEdge*, _LayerEdgeCmp >::iterator eIt = ee.begin(); while ( E2.IsNull() && eIt != ee.end()) { _LayerEdge* e2 = *eIt++; @@ -2405,7 +2864,7 @@ bool _LayerEdge::FindIntersection( SMESH_ElementSearcher& searcher, } if ( intFound ) { - if ( dist < segLen*(1.01)) + if ( dist < segLen*(1.01) && dist > -(_len-segLen) ) segmentIntersected = true; if ( distance > dist ) distance = dist, iFace = j; @@ -2463,7 +2922,7 @@ gp_Ax1 _LayerEdge::LastSegment(double& segLen) const gp_Ax1 segDir; if ( iPrev < 0 ) { - segDir.SetLocation( SMESH_MeshEditor::TNodeXYZ( _nodes[0] )); + segDir.SetLocation( SMESH_TNodeXYZ( _nodes[0] )); segDir.SetDirection( _normal ); segLen = 0; } @@ -2484,7 +2943,7 @@ gp_Ax1 _LayerEdge::LastSegment(double& segLen) const Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face(_sWOL), loc ); pPrev = surface->Value( pPrev.X(), pPrev.Y() ).Transformed( loc ); } - dir = SMESH_MeshEditor::TNodeXYZ( _nodes.back() ) - pPrev.XYZ(); + dir = SMESH_TNodeXYZ( _nodes.back() ) - pPrev.XYZ(); } segDir.SetLocation( pPrev ); segDir.SetDirection( dir ); @@ -2514,9 +2973,9 @@ bool _LayerEdge::SegTriaInter( const gp_Ax1& lastSegment, gp_XYZ orig = lastSegment.Location().XYZ(); gp_XYZ dir = lastSegment.Direction().XYZ(); - SMESH_MeshEditor::TNodeXYZ vert0( n0 ); - SMESH_MeshEditor::TNodeXYZ vert1( n1 ); - SMESH_MeshEditor::TNodeXYZ vert2( n2 ); + SMESH_TNodeXYZ vert0( n0 ); + SMESH_TNodeXYZ vert1( n1 ); + SMESH_TNodeXYZ vert2( n2 ); /* calculate distance from vert0 to ray origin */ gp_XYZ tvec = orig - vert0; @@ -2597,18 +3056,19 @@ bool _LayerEdge::SmoothOnEdge(Handle(Geom_Surface)& surface, ASSERT( IsOnEdge() ); SMDS_MeshNode* tgtNode = const_cast( _nodes.back() ); - SMESH_MeshEditor::TNodeXYZ oldPos( tgtNode ); + SMESH_TNodeXYZ oldPos( tgtNode ); double dist01, distNewOld; - SMESH_MeshEditor::TNodeXYZ p0( _2neibors->_nodes[0]); - SMESH_MeshEditor::TNodeXYZ p1( _2neibors->_nodes[1]); + SMESH_TNodeXYZ p0( _2neibors->_nodes[0]); + SMESH_TNodeXYZ p1( _2neibors->_nodes[1]); dist01 = p0.Distance( _2neibors->_nodes[1] ); gp_Pnt newPos = p0 * _2neibors->_wgt[0] + p1 * _2neibors->_wgt[1]; double lenDelta = 0; if ( _curvature ) { - lenDelta = _curvature->lenDelta( _len ); + //lenDelta = _curvature->lenDelta( _len ); + lenDelta = _curvature->lenDeltaByDist( dist01 ); newPos.ChangeCoord() += _normal * lenDelta; } @@ -2619,7 +3079,7 @@ bool _LayerEdge::SmoothOnEdge(Handle(Geom_Surface)& surface, if ( _2neibors->_plnNorm ) { // put newPos on the plane defined by source node and _plnNorm - gp_XYZ new2src = SMESH_MeshEditor::TNodeXYZ( _nodes[0] ) - newPos.XYZ(); + gp_XYZ new2src = SMESH_TNodeXYZ( _nodes[0] ) - newPos.XYZ(); double new2srcProj = (*_2neibors->_plnNorm) * new2src; newPos.ChangeCoord() += (*_2neibors->_plnNorm) * new2srcProj; } @@ -2665,7 +3125,7 @@ bool _LayerEdge::Smooth(int& badNb) // compute new position for the last _pos gp_XYZ newPos (0,0,0); for ( unsigned i = 0; i < _simplices.size(); ++i ) - newPos += SMESH_MeshEditor::TNodeXYZ( _simplices[i]._nPrev ); + newPos += SMESH_TNodeXYZ( _simplices[i]._nPrev ); newPos /= _simplices.size(); if ( _curvature ) @@ -2686,7 +3146,7 @@ bool _LayerEdge::Smooth(int& badNb) // } // count quality metrics (orientation) of tetras around _tgtNode int nbOkBefore = 0; - SMESH_MeshEditor::TNodeXYZ tgtXYZ( _nodes.back() ); + SMESH_TNodeXYZ tgtXYZ( _nodes.back() ); for ( unsigned i = 0; i < _simplices.size(); ++i ) nbOkBefore += _simplices[i].IsForward( _nodes[0], &tgtXYZ ); @@ -2699,7 +3159,7 @@ bool _LayerEdge::Smooth(int& badNb) SMDS_MeshNode* n = const_cast< SMDS_MeshNode* >( _nodes.back() ); - _len -= prevPos.Distance(SMESH_MeshEditor::TNodeXYZ( n )); + _len -= prevPos.Distance(SMESH_TNodeXYZ( n )); _len += prevPos.Distance(newPos); n->setXYZ( newPos.X(), newPos.Y(), newPos.Z()); @@ -2727,7 +3187,7 @@ void _LayerEdge::SetNewLength( double len, SMESH_MesherHelper& helper ) } SMDS_MeshNode* n = const_cast< SMDS_MeshNode*>( _nodes.back() ); - SMESH_MeshEditor::TNodeXYZ oldXYZ( n ); + SMESH_TNodeXYZ oldXYZ( n ); gp_XYZ nXYZ = oldXYZ + _normal * ( len - _len ) * _lenFactor; n->setXYZ( nXYZ.X(), nXYZ.Y(), nXYZ.Z() ); @@ -2938,6 +3398,15 @@ bool _ViscousBuilder::refine(_SolidData& data) } } + if ( !getMeshDS()->IsEmbeddedMode() ) + // Log node movement + for ( unsigned i = 0; i < data._edges.size(); ++i ) + { + _LayerEdge& edge = *data._edges[i]; + SMESH_TNodeXYZ p ( edge._nodes.back() ); + getMeshDS()->MoveNode( p._node, p.X(), p.Y(), p.Z() ); + } + // TODO: make quadratic prisms and polyhedrons(?) helper.SetElementsOnShape(true); @@ -3017,14 +3486,16 @@ bool _ViscousBuilder::shrink() SMDS_ElemIteratorPtr fIt = smDS->GetElements(); while ( fIt->more() ) proxySub->AddElement( fIt->next() ); + // as a result 3D algo will use elements from proxySub and not from smDS } } } SMESH_MesherHelper helper( *_mesh ); + helper.ToFixNodeParameters( true ); // EDGE's to shrink - map< int, _Shrinker1D > e2shrMap; + map< TGeomID, _Shrinker1D > e2shrMap; // loop on FACES to srink mesh on map< TGeomID, _SolidData* >::iterator f2sd = f2sdMap.begin(); @@ -3033,7 +3504,6 @@ bool _ViscousBuilder::shrink() _SolidData& data = *f2sd->second; TNode2Edge& n2eMap = data._n2eMap; const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( f2sd->first )); - const bool reverse = ( data._reversedFaceIds.count( f2sd->first )); Handle(Geom_Surface) surface = BRep_Tool::Surface(F); @@ -3046,8 +3516,8 @@ bool _ViscousBuilder::shrink() // Prepare data for shrinking // =========================== - // Collect nodes to smooth as src nodes are not yet replaced by tgt ones - // and thus all nodes on FACE connected to 2d elements are to be smoothed + // Collect nodes to smooth, as src nodes are not yet replaced by tgt ones + // and thus all nodes on a FACE connected to 2d elements are to be smoothed vector < const SMDS_MeshNode* > smoothNodes; { SMDS_NodeIteratorPtr nIt = smDS->GetNodes(); @@ -3061,12 +3531,15 @@ bool _ViscousBuilder::shrink() // Find out face orientation double refSign = 1; const set ignoreShapes; + bool isOkUV; if ( !smoothNodes.empty() ) { - gp_XY uv = helper.GetNodeUV( F, smoothNodes[0] ); vector<_Simplex> simplices; getSimplices( smoothNodes[0], simplices, ignoreShapes ); - if ( simplices[0].IsForward(uv, F, helper,refSign) != (!reverse)) + helper.GetNodeUV( F, simplices[0]._nPrev, 0, &isOkUV ); // fix UV of silpmex nodes + helper.GetNodeUV( F, simplices[0]._nNext, 0, &isOkUV ); + gp_XY uv = helper.GetNodeUV( F, smoothNodes[0], 0, &isOkUV ); + if ( !simplices[0].IsForward(uv, smoothNodes[0], F, helper,refSign) ) refSign = -1; } @@ -3114,6 +3587,9 @@ bool _ViscousBuilder::shrink() } } + // find out if a FACE is concave + const bool isConcaveFace = isConcave( F, helper ); + // Create _SmoothNode's on face F vector< _SmoothNode > nodesToSmooth( smoothNodes.size() ); { @@ -3123,7 +3599,9 @@ bool _ViscousBuilder::shrink() const SMDS_MeshNode* n = smoothNodes[i]; nodesToSmooth[ i ]._node = n; // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices - getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes ); + getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, NULL, isConcaveFace ); + // fix up incorrect uv of nodes on the FACE + helper.GetNodeUV( F, n, 0, &isOkUV); dumpMove( n ); } dumpFunctionEnd(); @@ -3142,6 +3620,7 @@ bool _ViscousBuilder::shrink() _Shrinker1D& srinker = e2shrMap[ edgeIndex ]; eShri1D.insert( & srinker ); srinker.AddEdge( edge, helper ); + VISCOUS_3D::ToClearSubWithMain( _mesh->GetSubMesh( edge->_sWOL ), data._solid ); // restore params of nodes on EGDE if the EDGE has been already // srinked while srinking another FACE srinker.RestoreParams(); @@ -3166,8 +3645,6 @@ bool _ViscousBuilder::shrink() shrinked |= lEdges[i]->SetNewLength2d( surface,F,helper ); } dumpFunctionEnd(); - if ( !shrinked ) - break; // Move nodes on EDGE's set< _Shrinker1D* >::iterator shr = eShri1D.begin(); @@ -3188,7 +3665,9 @@ bool _ViscousBuilder::shrink() moved = false; for ( unsigned i = 0; i < nodesToSmooth.size(); ++i ) { - moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,/*set3D=*/false ); + moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign, + /*isCentroidal=*/isConcaveFace, + /*set3D=*/isConcaveFace); } if ( badNb < oldBadNb ) nbNoImpSteps = 0; @@ -3200,18 +3679,44 @@ bool _ViscousBuilder::shrink() if ( badNb > 0 ) return error(SMESH_Comment("Can't shrink 2D mesh on face ") << f2sd->first ); } - // No wrongly shaped faces remain; final smooth. Set node XYZ - for ( int st = 3; st; --st ) + // No wrongly shaped faces remain; final smooth. Set node XYZ. + // First, find out a needed quality of smoothing (high for quadrangles only) + bool highQuality; + { + const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles(); + if ( hasTria != hasQuad ) { + highQuality = hasQuad; + } + else { + set nbNodesSet; + SMDS_ElemIteratorPtr fIt = smDS->GetElements(); + while ( fIt->more() && nbNodesSet.size() < 2 ) + nbNodesSet.insert( fIt->next()->NbCornerNodes() ); + highQuality = ( *nbNodesSet.begin() == 4 ); + } + } + if ( !highQuality && isConcaveFace ) + fixBadFaces( F, helper ); // fix narrow faces by swaping diagonals + for ( int st = highQuality ? 10 : 3; st; --st ) { dumpFunction(SMESH_Comment("shrinkFace")<first<<"_st"<<++smooStep); // debug for ( unsigned i = 0; i < nodesToSmooth.size(); ++i ) - nodesToSmooth[i].Smooth( badNb,surface,helper,refSign,/*set3D=*/st==1 ); + nodesToSmooth[i].Smooth( badNb,surface,helper,refSign, + /*isCentroidal=*/isConcaveFace,/*set3D=*/st==1 ); dumpFunctionEnd(); } - // Set event listener to clear FACE sub-mesh together with SOLID sub-mesh - _SrinkShapeListener::ToClearSubMeshWithSolid( sm, data._solid ); + // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh + VISCOUS_3D::ToClearSubWithMain( sm, data._solid ); + + if ( !getMeshDS()->IsEmbeddedMode() ) + // Log node movement + for ( unsigned i = 0; i < nodesToSmooth.size(); ++i ) + { + SMESH_TNodeXYZ p ( nodesToSmooth[i]._node ); + getMeshDS()->MoveNode( nodesToSmooth[i]._node, p.X(), p.Y(), p.Z() ); + } - }// loop on FACES to srink mesh on + } // loop on FACES to srink mesh on // Replace source nodes by target nodes in shrinked mesh edges @@ -3397,6 +3902,124 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, // return true; } +//================================================================================ +/*! + * \brief Try to fix triangles with high aspect ratio by swaping diagonals + */ +//================================================================================ + +void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& helper) +{ + SMESH::Controls::AspectRatio qualifier; + SMESH::Controls::TSequenceOfXYZ points(3), points1(3), points2(3); + const double maxAspectRatio = 4.; + + // find bad triangles + + vector< const SMDS_MeshElement* > badTrias; + vector< double > badAspects; + SMESHDS_SubMesh* sm = helper.GetMeshDS()->MeshElements( F ); + SMDS_ElemIteratorPtr fIt = sm->GetElements(); + while ( fIt->more() ) + { + const SMDS_MeshElement * f = fIt->next(); + if ( f->NbCornerNodes() != 3 ) continue; + for ( int iP = 0; iP < 3; ++iP ) points(iP+1) = SMESH_TNodeXYZ( f->GetNode(iP)); + double aspect = qualifier.GetValue( points ); + if ( aspect > maxAspectRatio ) + { + badTrias.push_back( f ); + badAspects.push_back( aspect ); + } + } + if ( badTrias.empty() ) + return; + + // find couples of faces to swap diagonal + + typedef pair < const SMDS_MeshElement* , const SMDS_MeshElement* > T2Trias; + vector< T2Trias > triaCouples; + + TIDSortedElemSet involvedFaces, emptySet; + for ( size_t iTia = 0; iTia < badTrias.size(); ++iTia ) + { + T2Trias trias [3]; + double aspRatio [3]; + int i1, i2, i3; + + involvedFaces.insert( badTrias[iTia] ); + for ( int iP = 0; iP < 3; ++iP ) + points(iP+1) = SMESH_TNodeXYZ( badTrias[iTia]->GetNode(iP)); + + // find triangles adjacent to badTrias[iTia] with better aspect ratio after diag-swaping + int bestCouple = -1; + for ( int iSide = 0; iSide < 3; ++iSide ) + { + const SMDS_MeshNode* n1 = badTrias[iTia]->GetNode( iSide ); + const SMDS_MeshNode* n2 = badTrias[iTia]->GetNode(( iSide+1 ) % 3 ); + trias [iSide].first = badTrias[iTia]; + trias [iSide].second = SMESH_MeshEditor::FindFaceInSet( n1, n2, emptySet, involvedFaces, + & i1, & i2 ); + if ( ! trias[iSide].second || trias[iSide].second->NbCornerNodes() != 3 ) + continue; + + // aspect ratio of an adjacent tria + for ( int iP = 0; iP < 3; ++iP ) + points2(iP+1) = SMESH_TNodeXYZ( trias[iSide].second->GetNode(iP)); + double aspectInit = qualifier.GetValue( points2 ); + + // arrange nodes as after diag-swaping + if ( helper.WrapIndex( i1+1, 3 ) == i2 ) + i3 = helper.WrapIndex( i1-1, 3 ); + else + i3 = helper.WrapIndex( i1+1, 3 ); + points1 = points; + points1( 1+ iSide ) = points2( 1+ i3 ); + points2( 1+ i2 ) = points1( 1+ ( iSide+2 ) % 3 ); + + // aspect ratio after diag-swaping + aspRatio[ iSide ] = qualifier.GetValue( points1 ) + qualifier.GetValue( points2 ); + if ( aspRatio[ iSide ] > aspectInit + badAspects[ iTia ] ) + continue; + + if ( bestCouple < 0 || aspRatio[ bestCouple ] > aspRatio[ iSide ] ) + bestCouple = iSide; + } + + if ( bestCouple >= 0 ) + { + triaCouples.push_back( trias[bestCouple] ); + involvedFaces.insert ( trias[bestCouple].second ); + } + else + { + involvedFaces.erase( badTrias[iTia] ); + } + } + if ( triaCouples.empty() ) + return; + + // swap diagonals + + SMESH_MeshEditor editor( helper.GetMesh() ); + dumpFunction(SMESH_Comment("beforeSwapDiagonals_F")< uv( _simplices.size() ); + for ( size_t i = 0; i < _simplices.size(); ++i ) + uv[i] = helper.GetNodeUV( face, _simplices[i]._nPrev, _node ); + // compute new UV for the node gp_XY newPos (0,0); - for ( unsigned i = 0; i < _simplices.size(); ++i ) - newPos += helper.GetNodeUV( face, _simplices[i]._nPrev ); - newPos /= _simplices.size(); + if ( isCentroidal && _simplices.size() > 3 ) + { + // average centers of diagonals wieghted with their reciprocal lengths + if ( _simplices.size() == 4 ) + { + double w1 = 1. / ( uv[2]-uv[0] ).SquareModulus(); + double w2 = 1. / ( uv[3]-uv[1] ).SquareModulus(); + newPos = ( w1 * ( uv[2]+uv[0] ) + w2 * ( uv[3]+uv[1] )) / ( w1+w2 ) / 2; + } + else + { + double sumWeight = 0; + int nb = _simplices.size() == 4 ? 2 : _simplices.size(); + for ( int i = 0; i < nb; ++i ) + { + int iFrom = i + 2; + int iTo = i + _simplices.size() - 1; + for ( int j = iFrom; j < iTo; ++j ) + { + int i2 = SMESH_MesherHelper::WrapIndex( j, _simplices.size() ); + double w = 1. / ( uv[i]-uv[i2] ).SquareModulus(); + sumWeight += w; + newPos += w * ( uv[i]+uv[i2] ); + } + } + newPos /= 2 * sumWeight; + } + } + else + { + // Laplacian smooth + isCentroidal = false; + for ( size_t i = 0; i < _simplices.size(); ++i ) + newPos += uv[i]; + newPos /= _simplices.size(); + } // count quality metrics (orientation) of triangles around the node int nbOkBefore = 0; gp_XY tgtUV = helper.GetNodeUV( face, _node ); for ( unsigned i = 0; i < _simplices.size(); ++i ) - nbOkBefore += _simplices[i].IsForward( tgtUV, face, helper, refSign ); + nbOkBefore += _simplices[i].IsForward( tgtUV, _node, face, helper, refSign ); int nbOkAfter = 0; for ( unsigned i = 0; i < _simplices.size(); ++i ) - nbOkAfter += _simplices[i].IsForward( newPos, face, helper, refSign ); + nbOkAfter += _simplices[i].IsForward( newPos, _node, face, helper, refSign ); if ( nbOkAfter < nbOkBefore ) + { + // if ( isCentroidal ) + // return Smooth( badNb, surface, helper, refSign, !isCentroidal, set3D ); + badNb += _simplices.size() - nbOkBefore; return false; + } SMDS_FacePosition* pos = static_cast( _node->GetPosition() ); pos->SetUParameter( newPos.X() ); @@ -3597,10 +4264,10 @@ void _Shrinker1D::AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper ) return; TopLoc_Location loc; Handle(Geom_Curve) C = BRep_Tool::Curve(E, loc, f,l); - GeomAdaptor_Curve aCurve(C); + GeomAdaptor_Curve aCurve(C, f,l); const double totLen = GCPnts_AbscissaPoint::Length(aCurve, f, l); - int nbExpectNodes = eSubMesh->NbNodes() - e->_nodes.size(); + int nbExpectNodes = eSubMesh->NbNodes(); _initU .reserve( nbExpectNodes ); _normPar.reserve( nbExpectNodes ); _nodes .reserve( nbExpectNodes ); @@ -3651,7 +4318,7 @@ void _Shrinker1D::Compute(bool set3D, SMESH_MesherHelper& helper) if ( set3D || _done ) { Handle(Geom_Curve) C = BRep_Tool::Curve(E, f,l); - GeomAdaptor_Curve aCurve(C); + GeomAdaptor_Curve aCurve(C, f,l); if ( _edges[0] ) f = helper.GetNodeU( E, _edges[0]->_nodes.back(), _nodes[0] ); @@ -3708,6 +4375,7 @@ void _Shrinker1D::RestoreParams() } _done = false; } + //================================================================================ /*! * \brief Replace source nodes by target nodes in shrinked mesh edges @@ -3801,16 +4469,18 @@ bool _ViscousBuilder::addBoundaryElements() // Find out orientation and type of face to create - bool reverse = false, tria = false, isOnFace; + bool reverse = false, isOnFace; map< TGeomID, TopoDS_Shape >::iterator e2f = data._shrinkShape2Shape.find( getMeshDS()->ShapeToIndex( E )); TopoDS_Shape F; - if ( isOnFace = ( e2f != data._shrinkShape2Shape.end() )) + if (( isOnFace = ( e2f != data._shrinkShape2Shape.end() ))) { F = e2f->second.Oriented( TopAbs_FORWARD ); reverse = ( helper.GetSubShapeOri( F, E ) == TopAbs_REVERSED ); if ( helper.GetSubShapeOri( data._solid, F ) == TopAbs_REVERSED ) + reverse = !reverse, F.Reverse(); + if ( helper.IsReversedSubMesh( TopoDS::Face(F) )) reverse = !reverse; } else @@ -3824,8 +4494,6 @@ bool _ViscousBuilder::addBoundaryElements() !_ignoreShapeIds.count( e2f->first )) F = *pF; } - - tria = true; } // Find the sub-mesh to add new faces SMESHDS_SubMesh* sm = 0; @@ -3839,19 +4507,33 @@ bool _ViscousBuilder::addBoundaryElements() // Make faces const int dj1 = reverse ? 0 : 1; const int dj2 = reverse ? 1 : 0; - vector newFaces; - newFaces.reserve(( ledges.size() - 1 ) * (ledges[0]->_nodes.size() - 1 )); for ( unsigned j = 1; j < ledges.size(); ++j ) { vector< const SMDS_MeshNode*>& nn1 = ledges[j-dj1]->_nodes; vector< const SMDS_MeshNode*>& nn2 = ledges[j-dj2]->_nodes; if ( isOnFace ) - for ( unsigned z = 1; z < nn1.size(); ++z ) + for ( size_t z = 1; z < nn1.size(); ++z ) sm->AddElement( getMeshDS()->AddFace( nn1[z-1], nn2[z-1], nn2[z], nn1[z] )); else - for ( unsigned z = 1; z < nn1.size(); ++z ) + for ( size_t z = 1; z < nn1.size(); ++z ) sm->AddElement( new SMDS_FaceOfNodes( nn1[z-1], nn2[z-1], nn2[z], nn1[z])); } + + // Make edges + for ( int isFirst = 0; isFirst < 2; ++isFirst ) + { + _LayerEdge* edge = isFirst ? ledges.front() : ledges.back(); + if ( !edge->_sWOL.IsNull() && edge->_sWOL.ShapeType() == TopAbs_EDGE ) + { + vector< const SMDS_MeshNode*>& nn = edge->_nodes; + if ( nn[1]->GetInverseElementIterator( SMDSAbs_Edge )->more() ) + continue; + helper.SetSubShape( edge->_sWOL ); + helper.SetElementsOnShape( true ); + for ( size_t z = 1; z < nn.size(); ++z ) + helper.AddEdge( nn[z-1], nn[z] ); + } + } } }