X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ViscousLayers.cxx;h=5ec3c8fa19de78041c42cc09802e152f5a0cc28e;hb=4031dc980843986775a48d3f9eae8642be249887;hp=a855c24a35a5a3182d9c9bcdc3f96b054ef9a581;hpb=1067ffa6e7e5c394e3a1b17219d8b355a57607cd;p=modules%2Fsmesh.git diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index a855c24a3..5ec3c8fa1 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,13 +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 #include #include @@ -224,8 +223,11 @@ namespace VISCOUS_3D 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] = @@ -431,12 +433,18 @@ namespace VISCOUS_3D //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 ); }; //-------------------------------------------------------------------------------- /*! @@ -566,20 +574,17 @@ virtual SMDS_ElemIteratorPtr elementsIterator(SMDSAbs_ElementType type) const // 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::SetBndShapesToIgnore(const std::vector& faceIds) +void StdMeshers_ViscousLayers::SetBndShapes(const std::vector& faceIds, bool toIgnore) { - if ( faceIds != _ignoreBndShapeIds ) - _ignoreBndShapeIds = faceIds, NotifySubMeshesHypothesisModification(); -} // -------------------------------------------------------------------------------- -bool StdMeshers_ViscousLayers::IsIgnoredShape(const int shapeID) const -{ - return ( find( _ignoreBndShapeIds.begin(), _ignoreBndShapeIds.end(), shapeID ) - != _ignoreBndShapeIds.end() ); + if ( faceIds != _shapeIds ) + _shapeIds = faceIds, NotifySubMeshesHypothesisModification(); + if ( _isToIgnoreShapes != toIgnore ) + _isToIgnoreShapes = toIgnore, NotifySubMeshesHypothesisModification(); } // -------------------------------------------------------------------------------- void StdMeshers_ViscousLayers::SetTotalThickness(double thickness) { @@ -637,17 +642,22 @@ std::ostream & StdMeshers_ViscousLayers::SaveTo(std::ostream & save) save << " " << _nbLayers << " " << _thickness << " " << _stretchFactor - << " " << _ignoreBndShapeIds.size(); - for ( unsigned i = 0; i < _ignoreBndShapeIds.size(); ++i ) - save << " " << _ignoreBndShapeIds[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 ( _ignoreBndShapeIds.size() < nbFaces && load >> faceID ) - _ignoreBndShapeIds.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, @@ -1071,7 +1081,7 @@ bool _ViscousBuilder::findFacesWithLayers() vector ignoreFaces; for ( unsigned i = 0; i < _sdVec.size(); ++i ) { - vector ids = _sdVec[i]._hyp->GetBndShapesToIgnore(); + vector ids = _sdVec[i]._hyp->GetBndShapes(); for ( unsigned i = 0; i < ids.size(); ++i ) { const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[i] ); @@ -1096,7 +1106,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 ); } } @@ -1961,9 +1971,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 ) @@ -2067,14 +2078,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 ) @@ -2265,9 +2276,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; @@ -2627,10 +2638,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; @@ -3550,7 +3561,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; @@ -3594,12 +3605,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 ); @@ -3634,6 +3646,8 @@ bool _ViscousBuilder::shrink() bool shrinked = true; int badNb, shriStep=0, smooStep=0; + _SmoothNode::SmoothType smoothType + = isConcaveFace ? _SmoothNode::CENTROIDAL : _SmoothNode::LAPLACIAN; while ( shrinked ) { // Move boundary nodes (actually just set new UV) @@ -3647,6 +3661,7 @@ bool _ViscousBuilder::shrink() dumpFunctionEnd(); // Move nodes on EDGE's + // (XYZ is set as soon as a needed length reached in SetNewLength2d()) set< _Shrinker1D* >::iterator shr = eShri1D.begin(); for ( ; shr != eShri1D.end(); ++shr ) (*shr)->Compute( /*set3D=*/false, helper ); @@ -3666,8 +3681,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; @@ -3679,32 +3693,26 @@ 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. - // 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 VISCOUS_3D::ToClearSubWithMain( sm, data._solid ); @@ -3752,6 +3760,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; @@ -3781,7 +3790,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 ) { @@ -3958,8 +3967,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; @@ -4057,11 +4066,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(); @@ -4085,22 +4096,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() ); @@ -4122,7 +4133,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() ); @@ -4134,7 +4145,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 ) @@ -4159,13 +4193,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(); @@ -4207,6 +4241,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 @@ -4480,7 +4574,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