X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ViscousLayers.cxx;h=e164282c6c0dd7b27d77feea80eb52c062eccdb8;hb=1f44408d18cb7ea84d857b4c3f0868e08c778cb6;hp=2ea12f46f9f493261488698571c710d05db25b63;hpb=e3409934a4fb5496d8d687026fc7f1ee2d40c735;p=modules%2Fsmesh.git diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index 2ea12f46f..e164282c6 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-2011 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 @@ -43,9 +43,18 @@ #include "utilities.h" #include +#include +#include +#include #include +#include +#include +#include #include +#include #include +#include +#include #include #include #include @@ -62,7 +71,7 @@ #include #include -#include +#include #include //#define __myDEBUG @@ -107,7 +116,7 @@ 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 { @@ -268,6 +277,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] ); + } }; //-------------------------------------------------------------------------------- /*! @@ -362,9 +375,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; @@ -377,6 +389,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); }; //-------------------------------------------------------------------------------- /*! @@ -433,7 +452,13 @@ 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(); @@ -878,7 +903,7 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh& theMesh, return _error; addBoundaryElements(); - + makeGroupOfLE(); // debug return _error; @@ -886,7 +911,7 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh& theMesh, //================================================================================ /*! - * \brief Finds SOLIDs to compute using viscous layers. Fill _sdVec + * \brief Finds SOLIDs to compute using viscous layers. Fills _sdVec */ //================================================================================ @@ -1338,7 +1363,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 ) { @@ -1994,7 +2019,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) //================================================================================ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, - int nbSteps, + const int nbSteps, double & distToIntersection) { if ( data._endEdgeToSmooth.empty() ) @@ -2028,21 +2053,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 )) + { + dumpFunction(SMESH_Comment("smooth")<SmoothOnEdge(surface, F, helper); + } + dumpCmd( SMESH_Comment("# end step ")<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 ) @@ -2115,6 +2153,253 @@ 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() ); + } + } + 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]); + for ( int i = iFrom; i < iTo; ++i ) + { + double r = len[i-iFrom] / len.back(); + gp_XY newUV = uv0 * ( 1. - r ) + uv1 * r; + 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() ); + + 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 < 0 ) uLast += 2*PI; // 0.0 - 2*PI + if ( uMidl < 0 ) uMidl += 2*PI; + const bool sense = ( uMidl < uLast ); + const double radius = 0.5 * ( vec0.Magnitude() + vec1.Magnitude() ); + + gp_Ax2d axis( center, vec0 ); + gp_Circ2d circ ( axis, radius, sense ); + 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() ); + + 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 @@ -2226,7 +2511,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++; @@ -3027,6 +3312,7 @@ 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 } } } @@ -3056,8 +3342,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(); @@ -3125,6 +3411,7 @@ bool _ViscousBuilder::shrink() } // Create _SmoothNode's on face F + bool isOkUV; vector< _SmoothNode > nodesToSmooth( smoothNodes.size() ); { dumpFunction(SMESH_Comment("beforeShrinkFace")<first); // debug @@ -3134,6 +3421,8 @@ bool _ViscousBuilder::shrink() nodesToSmooth[ i ]._node = n; // src nodes must be replaced by tgt nodes to have tgt nodes in _simplices getSimplices( n, nodesToSmooth[ i ]._simplices, ignoreShapes ); + // fix up incorrect uv of nodes on the FACE + helper.GetNodeUV( F, n, 0, &isOkUV); dumpMove( n ); } dumpFunctionEnd(); @@ -3210,15 +3499,32 @@ 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 ); + } + } + for ( int st = highQuality ? 8 : 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 ); dumpFunctionEnd(); } - // Set event listener to clear FACE sub-mesh together with SOLID sub-mesh + // Set an event listener to clear FACE sub-mesh together with SOLID sub-mesh _SrinkShapeListener::ToClearSubMeshWithSolid( sm, data._solid ); }// loop on FACES to srink mesh on @@ -3607,7 +3913,7 @@ 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(); @@ -3661,7 +3967,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] ); @@ -3718,6 +4024,7 @@ void _Shrinker1D::RestoreParams() } _done = false; } + //================================================================================ /*! * \brief Replace source nodes by target nodes in shrinked mesh edges @@ -3811,12 +4118,12 @@ 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 ); @@ -3834,8 +4141,6 @@ bool _ViscousBuilder::addBoundaryElements() !_ignoreShapeIds.count( e2f->first )) F = *pF; } - - tria = true; } // Find the sub-mesh to add new faces SMESHDS_SubMesh* sm = 0; @@ -3849,8 +4154,6 @@ 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;