X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ViscousLayers.cxx;h=6a6d871c43d65d25128c00a3f466aa903662e998;hp=c875ca104a4f622be907664d77f9c5f2329f6597;hb=e2a638a02647937174d628352feeef9640c14ec0;hpb=f5016d85b7b4b88623723027a1585c6414c4dc66 diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index c875ca104..6a6d871c4 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2012 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2013 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 @@ -36,12 +36,12 @@ #include "SMESH_Gen.hxx" #include "SMESH_Group.hxx" #include "SMESH_Mesh.hxx" +#include "SMESH_MeshAlgos.hxx" #include "SMESH_MesherHelper.hxx" #include "SMESH_ProxyMesh.hxx" #include "SMESH_subMesh.hxx" #include "SMESH_subMeshEventListener.hxx" - -#include "utilities.h" +#include "StdMeshers_FaceSide.hxx" #include #include @@ -82,7 +82,7 @@ using namespace std; //================================================================================ -namespace VISCOUS +namespace VISCOUS_3D { typedef int TGeomID; @@ -121,13 +121,13 @@ namespace VISCOUS * \brief Listener of events of 3D sub-meshes computed with 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() + _ShrinkShapeListener() : SMESH_subMeshEventListener(/*isDeletable=*/false, - "StdMeshers_ViscousLayers::_SrinkShapeListener") {} - static SMESH_subMeshEventListener* Get() { static _SrinkShapeListener l; return &l; } + "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, @@ -139,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 ); - } - } }; //-------------------------------------------------------------------------------- /*! @@ -205,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 @@ -215,8 +224,11 @@ namespace VISCOUS struct _Simplex { const SMDS_MeshNode *_nPrev, *_nNext; // nodes on a smoothed mesh surface - _Simplex(const SMDS_MeshNode* nPrev=0, const SMDS_MeshNode* nNext=0) - : _nPrev(nPrev), _nNext(nNext) {} + const SMDS_MeshNode *_nOpp; // in 2D case, a node opposite to a smoothed node in QUAD + _Simplex(const SMDS_MeshNode* nPrev=0, + const SMDS_MeshNode* nNext=0, + const SMDS_MeshNode* nOpp=0) + : _nPrev(nPrev), _nNext(nNext), _nOpp(nOpp) {} bool IsForward(const SMDS_MeshNode* nSrc, const gp_XYZ* pntTgt) const { const double M[3][3] = @@ -256,6 +268,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 ) { @@ -266,10 +279,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; //-------------------------------------------------------------------------------- @@ -419,12 +434,18 @@ namespace VISCOUS //vector _nodesAround; vector<_Simplex> _simplices; // for quality check + enum SmoothType { LAPLACIAN, CENTROIDAL, ANGULAR }; + bool Smooth(int& badNb, Handle(Geom_Surface)& surface, SMESH_MesherHelper& helper, const double refSign, - bool isCentroidal, + SmoothType how, bool set3D); + + gp_XY computeAngularPos(vector& uv, + const gp_XY& uvToFix, + const double refSign ); }; //-------------------------------------------------------------------------------- /*! @@ -547,22 +568,24 @@ virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const _nn[3]=_le2->_nodes[0]; } }; -} // namespace VISCOUS +} // namespace VISCOUS_3D //================================================================================ // StdMeshers_ViscousLayers hypothesis // StdMeshers_ViscousLayers::StdMeshers_ViscousLayers(int hypId, int studyId, SMESH_Gen* gen) :SMESH_Hypothesis(hypId, studyId, gen), - _nbLayers(1), _thickness(1), _stretchFactor(1) + _isToIgnoreShapes(18), _nbLayers(1), _thickness(1), _stretchFactor(1) { _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::SetBndShapes(const std::vector& faceIds, bool toIgnore) { - if ( faceIds != _ignoreFaceIds ) - _ignoreFaceIds = faceIds, NotifySubMeshesHypothesisModification(); + if ( faceIds != _shapeIds ) + _shapeIds = faceIds, NotifySubMeshesHypothesisModification(); + if ( _isToIgnoreShapes != toIgnore ) + _isToIgnoreShapes = toIgnore, NotifySubMeshesHypothesisModification(); } // -------------------------------------------------------------------------------- void StdMeshers_ViscousLayers::SetTotalThickness(double thickness) { @@ -584,7 +607,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() ) @@ -620,17 +643,22 @@ 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]; + << " " << _shapeIds.size(); + for ( unsigned i = 0; i < _shapeIds.size(); ++i ) + save << " " << _shapeIds[i]; + save << " " << !_isToIgnoreShapes; // negate to keep the behavior in old studies. return save; } // -------------------------------------------------------------------------------- std::istream & StdMeshers_ViscousLayers::LoadFrom(std::istream & load) { - int nbFaces, faceID; + int nbFaces, faceID, shapeToTreat; load >> _nbLayers >> _thickness >> _stretchFactor >> nbFaces; - while ( _ignoreFaceIds.size() < nbFaces && load >> faceID ) - _ignoreFaceIds.push_back( faceID ); + while ( _shapeIds.size() < nbFaces && load >> faceID ) + _shapeIds.push_back( faceID ); + if ( load >> shapeToTreat ) + _isToIgnoreShapes = !shapeToTreat; + else + _isToIgnoreShapes = true; // old behavior return load; } // -------------------------------------------------------------------------------- bool StdMeshers_ViscousLayers::SetParametersByMesh(const SMESH_Mesh* theMesh, @@ -716,27 +744,30 @@ namespace 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 ); } } @@ -750,13 +781,15 @@ namespace bool isConcave( const TopoDS_Face& F, SMESH_MesherHelper& helper ) { + // if ( helper.Count( F, TopAbs_WIRE, /*useMap=*/false) > 1 ) + // return true; 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; + if ( SMESH_Algo::isDegenerated( E )) continue; // check if 2D curve is concave BRepAdaptor_Curve2d curve( E, F ); const int nbIntervals = curve.NbIntervals( GeomAbs_C2 ); @@ -768,10 +801,10 @@ namespace double u1 = intervals( i ); double u2 = intervals( i+1 ); curve.D2( 0.5*( u1+u2 ), p, drv1, drv2 ); - double cross = drv2 ^ drv1; + double cross = drv2 * drv1; //drv2 ^ drv1; if ( E.Orientation() == TopAbs_REVERSED ) cross = -cross; - isConvex = ( cross < 1e-9 ); + isConvex = ( cross > -1e-9 ); } // check if concavity is strong enough to care about it //const double maxAngle = 5 * Standard_PI180; @@ -801,6 +834,26 @@ namespace // } } } + // check angles at VERTEXes + TError error; + TSideVector wires = StdMeshers_FaceSide::GetFaceWires( F, *helper.GetMesh(), 0, error ); + for ( size_t iW = 0; iW < wires.size(); ++iW ) + { + const int nbEdges = wires[iW]->NbEdges(); + if ( nbEdges < 2 && SMESH_Algo::isDegenerated( wires[iW]->Edge(0))) + continue; + for ( int iE1 = 0; iE1 < nbEdges; ++iE1 ) + { + if ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE1 ))) continue; + int iE2 = ( iE1 + 1 ) % nbEdges; + while ( SMESH_Algo::isDegenerated( wires[iW]->Edge( iE2 ))) + iE2 = ( iE2 + 1 ) % nbEdges; + double angle = helper.GetAngle( wires[iW]->Edge( iE1 ), + wires[iW]->Edge( iE2 ), F ); + if ( angle < -5. * M_PI / 180. ) + return true; + } + } return false; } //-------------------------------------------------------------------------------- @@ -812,9 +865,11 @@ namespace const char* fname = "/tmp/viscous.py"; cout << "execfile('"< ignoreFaces; for ( unsigned i = 0; i < _sdVec.size(); ++i ) { - vector ids = _sdVec[i]._hyp->GetIgnoreFaces(); + vector ids = _sdVec[i]._hyp->GetBndShapes(); for ( unsigned i = 0; i < ids.size(); ++i ) { const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[i] ); @@ -1076,7 +1131,7 @@ bool _ViscousBuilder::findFacesWithLayers() { _ignoreShapeIds.insert( faceInd ); ignoreFaces.push_back( exp.Current() ); - if ( SMESH_Algo::IsReversedSubMesh( TopoDS::Face( exp.Current() ), getMeshDS())) + if ( helper.IsReversedSubMesh( TopoDS::Face( exp.Current() ))) _sdVec[i]._reversedFaceIds.insert( faceInd ); } } @@ -1371,7 +1426,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) if ( data._stepSize < 1. ) data._epsilon *= data._stepSize; - // Put _LayerEdge's into a vector + // Put _LayerEdge's into the vector data._edges if ( !sortEdges( data, edgesByGeom )) return false; @@ -1695,8 +1750,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; @@ -1941,9 +1996,10 @@ void _ViscousBuilder::getSimplices( const SMDS_MeshNode* node, int srcInd = f->GetNodeIndex( node ); const SMDS_MeshNode* nPrev = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd-1, nbNodes )); const SMDS_MeshNode* nNext = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+1, nbNodes )); + const SMDS_MeshNode* nOpp = f->GetNode( SMESH_MesherHelper::WrapIndex( srcInd+2, nbNodes )); if ( dataToCheckOri && dataToCheckOri->_reversedFaceIds.count( shapeInd )) std::swap( nPrev, nNext ); - simplices.push_back( _Simplex( nPrev, nNext )); + simplices.push_back( _Simplex( nPrev, nNext, nOpp )); } if ( toSort ) @@ -2047,14 +2103,14 @@ bool _ViscousBuilder::inflate(_SolidData& data) // Limit inflation step size by geometry size found by itersecting // normals of _LayerEdge's with mesh faces double geomSize = Precision::Infinite(), intersecDist; - SMESH_MeshEditor editor( _mesh ); auto_ptr searcher - ( editor.GetElementSearcher( data._proxyMesh->GetFaces( data._solid )) ); + ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), + data._proxyMesh->GetFaces( data._solid )) ); for ( unsigned i = 0; i < data._edges.size(); ++i ) { 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 ) @@ -2245,9 +2301,9 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, // Check if the last segments of _LayerEdge intersects 2D elements; // checked elements are either temporary faces or faces on surfaces w/o the layers - SMESH_MeshEditor editor( _mesh ); auto_ptr searcher - ( editor.GetElementSearcher( data._proxyMesh->GetFaces( data._solid )) ); + ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), + data._proxyMesh->GetFaces( data._solid )) ); distToIntersection = Precision::Infinite(); double dist; @@ -2502,6 +2558,10 @@ bool _ViscousBuilder::smoothAnalyticEdge( _SolidData& data, if ( F.IsNull() ) // 3D { + if ( data._edges[iFrom]->_2neibors->_nodes[0] == + data._edges[iTo-1]->_2neibors->_nodes[1] ) + return true; // closed EDGE - nothing to do + return false; // TODO ??? } else // 2D @@ -2607,10 +2667,10 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, // 1) to find and fix intersection // 2) to check that no new intersection appears as result of 1) - SMESH_MeshEditor editor( _mesh ); SMDS_ElemIteratorPtr fIt( new SMDS_ElementVectorIterator( tmpFaces.begin(), tmpFaces.end())); - auto_ptr searcher ( editor.GetElementSearcher( fIt )); + auto_ptr searcher + ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), fIt )); // 1) Find intersections double dist; @@ -2844,7 +2904,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; @@ -3047,7 +3107,8 @@ bool _LayerEdge::SmoothOnEdge(Handle(Geom_Surface)& surface, double lenDelta = 0; if ( _curvature ) { - lenDelta = _curvature->lenDelta( _len ); + //lenDelta = _curvature->lenDelta( _len ); + lenDelta = _curvature->lenDeltaByDist( dist01 ); newPos.ChangeCoord() += _normal * lenDelta; } @@ -3529,7 +3590,7 @@ bool _ViscousBuilder::shrink() sm->getDependsOnIterator(/*includeSelf=*/false, /*complexShapeFirst=*/false); while ( subIt->more() ) { - SMESH_subMesh* sub = subIt->next(); + SMESH_subMesh* sub = subIt->next(); SMESHDS_SubMesh* subDS = sub->GetSubMeshDS(); if ( subDS->NbNodes() == 0 || !n2eMap.count( subDS->GetNodes()->next() )) continue; @@ -3573,12 +3634,13 @@ bool _ViscousBuilder::shrink() vector< _SmoothNode > nodesToSmooth( smoothNodes.size() ); { dumpFunction(SMESH_Comment("beforeShrinkFace")<first); // debug + const bool sortSimplices = isConcaveFace; for ( unsigned i = 0; i < smoothNodes.size(); ++i ) { 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, NULL, isConcaveFace ); + getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes, NULL, sortSimplices ); // fix up incorrect uv of nodes on the FACE helper.GetNodeUV( F, n, 0, &isOkUV); dumpMove( n ); @@ -3599,6 +3661,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(); @@ -3612,11 +3675,14 @@ bool _ViscousBuilder::shrink() bool shrinked = true; int badNb, shriStep=0, smooStep=0; + _SmoothNode::SmoothType smoothType + = isConcaveFace ? _SmoothNode::CENTROIDAL : _SmoothNode::LAPLACIAN; while ( shrinked ) { + shriStep++; // Move boundary nodes (actually just set new UV) // ----------------------------------------------- - dumpFunction(SMESH_Comment("moveBoundaryOnF")<first<<"_st"<first<<"_st"<::iterator shr = eShri1D.begin(); for ( ; shr != eShri1D.end(); ++shr ) (*shr)->Compute( /*set3D=*/false, helper ); @@ -3644,8 +3711,7 @@ bool _ViscousBuilder::shrink() for ( unsigned i = 0; i < nodesToSmooth.size(); ++i ) { moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign, - /*isCentroidal=*/isConcaveFace, - /*set3D=*/isConcaveFace); + smoothType, /*set3D=*/isConcaveFace); } if ( badNb < oldBadNb ) nbNoImpSteps = 0; @@ -3656,35 +3722,31 @@ bool _ViscousBuilder::shrink() } if ( badNb > 0 ) return error(SMESH_Comment("Can't shrink 2D mesh on face ") << f2sd->first ); + if ( shriStep > 200 ) + return error(SMESH_Comment("Infinite loop at shrinking 2D mesh on face ") << f2sd->first ); } + // No wrongly shaped faces remain; final smooth. Set node XYZ. - // First, find out a needed quality of smoothing (high for quadrangles only) - bool highQuality; + bool isStructuredFixed = false; + if ( SMESH_2D_Algo* algo = dynamic_cast( sm->GetAlgo() )) + isStructuredFixed = algo->FixInternalNodes( *data._proxyMesh, F ); + if ( !isStructuredFixed ) { - 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 ( isConcaveFace ) + fixBadFaces( F, helper ); // fix narrow faces by swapping 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, + smoothType,/*set3D=*/st==1 ); + } + dumpFunctionEnd(); } } - 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, - /*isCentroidal=*/isConcaveFace,/*set3D=*/st==1 ); - dumpFunctionEnd(); - } // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh - _SrinkShapeListener::ToClearSubMeshWithSolid( sm, data._solid ); + VISCOUS_3D::ToClearSubWithMain( sm, data._solid ); if ( !getMeshDS()->IsEmbeddedMode() ) // Log node movement @@ -3730,6 +3792,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, double uvLen = uvDir.Magnitude(); uvDir /= uvLen; edge._normal.SetCoord( uvDir.X(),uvDir.Y(), 0); + edge._len = uvLen; // IMPORTANT to have src nodes NOT yet REPLACED by tgt nodes in shrinked faces vector faces; @@ -3759,7 +3822,7 @@ bool _ViscousBuilder::prepareEdgeToShrink( _LayerEdge& edge, } multimap< double, const SMDS_MeshNode* >::iterator p2n = proj2node.begin(), p2nEnd; - const double minProj = p2n->first; + const double minProj = p2n->first; const double projThreshold = 1.1 * uvLen; if ( minProj > projThreshold ) { @@ -3936,8 +3999,8 @@ void _ViscousBuilder::fixBadFaces(const TopoDS_Face& F, SMESH_MesherHelper& help 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 ); + trias [iSide].second = SMESH_MeshAlgos::FindFaceInSet( n1, n2, emptySet, involvedFaces, + & i1, & i2 ); if ( ! trias[iSide].second || trias[iSide].second->NbCornerNodes() != 3 ) continue; @@ -4035,11 +4098,13 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface, double proj = uvDirN * uvDir * kSafe; if ( proj < stepSize && proj > minStepSize ) stepSize = proj; + else if ( proj < minStepSize ) + stepSize = minStepSize; } } gp_Pnt2d newUV; - if ( stepSize == uvLen ) + if ( uvLen - stepSize < _len / 20. ) { newUV = tgtUV; _pos.clear(); @@ -4063,22 +4128,22 @@ bool _LayerEdge::SetNewLength2d( Handle(Geom_Surface)& surface, { TopoDS_Edge E = TopoDS::Edge( _sWOL ); const SMDS_MeshNode* n2 = _simplices[0]._nPrev; + SMDS_EdgePosition* tgtPos = static_cast( tgtNode->GetPosition() ); const double u2 = helper.GetNodeU( E, n2, tgtNode ); const double uSrc = _pos[0].Coord( U_SRC ); const double lenTgt = _pos[0].Coord( LEN_TGT ); double newU = _pos[0].Coord( U_TGT ); - if ( lenTgt < 0.99 * fabs( uSrc-u2 )) + if ( lenTgt < 0.99 * fabs( uSrc-u2 )) // n2 got out of src-tgt range { _pos.clear(); } else { - newU = 0.1 * uSrc + 0.9 * u2; + newU = 0.1 * tgtPos->GetUParameter() + 0.9 * u2; } - SMDS_EdgePosition* pos = static_cast( tgtNode->GetPosition() ); - pos->SetUParameter( newU ); + tgtPos->SetUParameter( newU ); #ifdef __myDEBUG gp_XY newUV = helper.GetNodeUV( F, tgtNode, _nodes[0]); gp_Pnt p = surface->Value( newUV.X(), newUV.Y() ); @@ -4100,7 +4165,7 @@ bool _SmoothNode::Smooth(int& badNb, Handle(Geom_Surface)& surface, SMESH_MesherHelper& helper, const double refSign, - bool isCentroidal, + SmoothType how, bool set3D) { const TopoDS_Face& face = TopoDS::Face( helper.GetSubShape() ); @@ -4112,7 +4177,30 @@ bool _SmoothNode::Smooth(int& badNb, // compute new UV for the node gp_XY newPos (0,0); - if ( isCentroidal && _simplices.size() > 3 ) +/* if ( how == ANGULAR && _simplices.size() == 4 ) + { + vector corners; corners.reserve(4); + for ( size_t i = 0; i < _simplices.size(); ++i ) + if ( _simplices[i]._nOpp ) + corners.push_back( helper.GetNodeUV( face, _simplices[i]._nOpp, _node )); + if ( corners.size() == 4 ) + { + newPos = helper.calcTFI + ( 0.5, 0.5, + corners[0], corners[1], corners[2], corners[3], + uv[1], uv[2], uv[3], uv[0] ); + } + // vector p( _simplices.size() * 2 + 1 ); + // p.clear(); + // for ( size_t i = 0; i < _simplices.size(); ++i ) + // { + // p.push_back( uv[i] ); + // if ( _simplices[i]._nOpp ) + // p.push_back( helper.GetNodeUV( face, _simplices[i]._nOpp, _node )); + // } + // newPos = computeAngularPos( p, helper.GetNodeUV( face, _node ), refSign ); + } + else*/ if ( how == CENTROIDAL && _simplices.size() > 3 ) { // average centers of diagonals wieghted with their reciprocal lengths if ( _simplices.size() == 4 ) @@ -4137,13 +4225,13 @@ bool _SmoothNode::Smooth(int& badNb, newPos += w * ( uv[i]+uv[i2] ); } } - newPos /= 2 * sumWeight; + newPos /= 2 * sumWeight; // 2 is to get a middle between uv's } } else { // Laplacian smooth - isCentroidal = false; + //isCentroidal = false; for ( size_t i = 0; i < _simplices.size(); ++i ) newPos += uv[i]; newPos /= _simplices.size(); @@ -4185,6 +4273,66 @@ bool _SmoothNode::Smooth(int& badNb, return ( (tgtUV-newPos).SquareModulus() > 1e-10 ); } +//================================================================================ +/*! + * \brief Computes new UV using angle based smoothing technic + */ +//================================================================================ + +gp_XY _SmoothNode::computeAngularPos(vector& uv, + const gp_XY& uvToFix, + const double refSign) +{ + uv.push_back( uv.front() ); + + vector< gp_XY > edgeDir( uv.size() ); + vector< double > edgeSize( uv.size() ); + for ( size_t i = 1; i < edgeDir.size(); ++i ) + { + edgeDir[i-1] = uv[i] - uv[i-1]; + edgeSize[i-1] = edgeDir[i-1].Modulus(); + if ( edgeSize[i-1] < numeric_limits::min() ) + edgeDir[i-1].SetX( 100 ); + else + edgeDir[i-1] /= edgeSize[i-1] * refSign; + } + edgeDir.back() = edgeDir.front(); + edgeSize.back() = edgeSize.front(); + + gp_XY newPos(0,0); + int nbEdges = 0; + double sumSize = 0; + for ( size_t i = 1; i < edgeDir.size(); ++i ) + { + if ( edgeDir[i-1].X() > 1. ) continue; + int i1 = i-1; + 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 bisec = norm1 + norm2; + double bisecSize = bisec.Modulus(); + if ( bisecSize < numeric_limits::min() ) + { + bisec = -edgeDir[i1] + edgeDir[i]; + bisecSize = bisec.Modulus(); + } + bisec /= bisecSize; + + gp_XY dirToN = uvToFix - p; + double distToN = dirToN.Modulus(); + if ( bisec * dirToN < 0 ) + distToN = -distToN; + + newPos += ( p + bisec * distToN ) * ( edgeSize[i1] + edgeSize[i] ); + ++nbEdges; + sumSize += edgeSize[i1] + edgeSize[i]; + } + newPos /= /*nbEdges * */sumSize; + return newPos; +} + //================================================================================ /*! * \brief Delete _SolidData @@ -4458,7 +4606,7 @@ bool _ViscousBuilder::addBoundaryElements() reverse = ( helper.GetSubShapeOri( F, E ) == TopAbs_REVERSED ); if ( helper.GetSubShapeOri( data._solid, F ) == TopAbs_REVERSED ) reverse = !reverse, F.Reverse(); - if ( SMESH_Algo::IsReversedSubMesh( TopoDS::Face(F), getMeshDS() )) + if ( helper.IsReversedSubMesh( TopoDS::Face(F) )) reverse = !reverse; } else