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=4326b4f7aa2c63f46addb3800a7c256bffc34fec;hb=4cd2499bddcd3da3ec8900fe825bc98669b789b5;hpb=5acab3d3271d9516ca568577b1465693f5db96dd diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index 4326b4f7a..124b61e1e 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -1,4 +1,4 @@ -// 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 @@ -15,6 +15,7 @@ // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// // File : StdMeshers_ViscousLayers.cxx // Created : Wed Dec 1 15:15:34 2010 @@ -81,7 +82,7 @@ using namespace std; //================================================================================ -namespace VISCOUS +namespace VISCOUS_3D { typedef int TGeomID; @@ -120,11 +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(): 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, @@ -136,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 ); - } - } }; //-------------------------------------------------------------------------------- /*! @@ -162,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, @@ -200,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 @@ -251,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 ) { @@ -261,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; //-------------------------------------------------------------------------------- @@ -522,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()));} }; //-------------------------------------------------------------------------------- @@ -541,7 +559,7 @@ namespace VISCOUS _nn[3]=_le2->_nodes[0]; } }; -} // namespace VISCOUS +} // namespace VISCOUS_3D //================================================================================ // StdMeshers_ViscousLayers hypothesis @@ -553,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 != _ignoreBndShapeIds ) + _ignoreBndShapeIds = faceIds, NotifySubMeshesHypothesisModification(); +} // -------------------------------------------------------------------------------- +bool StdMeshers_ViscousLayers::IsIgnoredShape(const int shapeID) const { - if ( faceIds != _ignoreFaceIds ) - _ignoreFaceIds = faceIds, NotifySubMeshesHypothesisModification(); + return ( find( _ignoreBndShapeIds.begin(), _ignoreBndShapeIds.end(), shapeID ) + != _ignoreBndShapeIds.end() ); } // -------------------------------------------------------------------------------- void StdMeshers_ViscousLayers::SetTotalThickness(double thickness) { @@ -578,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() ) @@ -614,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, @@ -710,27 +733,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 ); } } @@ -843,7 +869,7 @@ namespace #endif } -using namespace VISCOUS; +using namespace VISCOUS_3D; //================================================================================ /*! @@ -970,6 +996,9 @@ SMESH_ComputeErrorPtr _ViscousBuilder::Compute(SMESH_Mesh& theMesh, { if ( ! makeLayer(_sdVec[i]) ) return _error; + + if ( _sdVec[i]._edges.size() == 0 ) + continue; if ( ! inflate(_sdVec[i]) ) return _error; @@ -1042,7 +1071,7 @@ bool _ViscousBuilder::findFacesWithLayers() vector ignoreFaces; for ( unsigned i = 0; i < _sdVec.size(); ++i ) { - vector ids = _sdVec[i]._hyp->GetIgnoreFaces(); + vector ids = _sdVec[i]._hyp->GetBndShapesToIgnore(); for ( unsigned i = 0; i < ids.size(); ++i ) { const TopoDS_Shape& s = getMeshDS()->IndexToShape( ids[i] ); @@ -1067,7 +1096,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 ); } } @@ -1362,7 +1391,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; @@ -1686,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; @@ -2045,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 ) @@ -2498,22 +2527,21 @@ bool _ViscousBuilder::smoothAnalyticEdge( _SolidData& data, 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 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.*M_PI; // 0.0 - 2*PI - if ( uMidl < 0 ) uMidl += 2.*M_PI; - const bool sense = ( uMidl < uLast ); + 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, sense ); + 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(); @@ -2836,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; @@ -3039,7 +3067,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; } @@ -3369,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); @@ -3454,9 +3492,10 @@ bool _ViscousBuilder::shrink() } 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(); @@ -3581,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(); @@ -3626,7 +3666,8 @@ bool _ViscousBuilder::shrink() for ( unsigned i = 0; i < nodesToSmooth.size(); ++i ) { moved |= nodesToSmooth[i].Smooth( badNb,surface,helper,refSign, - /*isCentroidal=*/isConcaveFace,/*set3D=*/false ); + /*isCentroidal=*/isConcaveFace, + /*set3D=*/isConcaveFace); } if ( badNb < oldBadNb ) nbNoImpSteps = 0; @@ -3643,12 +3684,10 @@ bool _ViscousBuilder::shrink() bool highQuality; { const bool hasTria = _mesh->NbTriangles(), hasQuad = _mesh->NbQuadrangles(); - if ( hasTria != hasQuad ) - { + if ( hasTria != hasQuad ) { highQuality = hasQuad; } - else - { + else { set nbNodesSet; SMDS_ElemIteratorPtr fIt = smDS->GetElements(); while ( fIt->more() && nbNodesSet.size() < 2 ) @@ -3667,7 +3706,15 @@ bool _ViscousBuilder::shrink() 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 + 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 @@ -4220,7 +4267,7 @@ void _Shrinker1D::AddEdge( const _LayerEdge* e, SMESH_MesherHelper& helper ) 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 ); @@ -4432,6 +4479,8 @@ bool _ViscousBuilder::addBoundaryElements() 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 @@ -4463,12 +4512,28 @@ bool _ViscousBuilder::addBoundaryElements() 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] ); + } + } } }