X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ViscousLayers.cxx;h=a2a44a23829f391426f09303ed72a60a62c9baef;hb=de78fd2b92839c0b70fca18939f243d9b201ce1a;hp=b5a076c5e3edc5ca459963484f9c044197a9d31d;hpb=bcfa36bbd041d724ff56560fb606261235d5f9e1;p=modules%2Fsmesh.git diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index b5a076c5e..a2a44a238 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -95,6 +95,12 @@ namespace VISCOUS_3D enum UIndex { U_TGT = 1, U_SRC, LEN_TGT }; const double theMinSmoothCosin = 0.1; + const double theSmoothThickToElemSizeRatio = 0.3; + + bool needSmoothing( double cosin, double tgtThick, double elemSize ) + { + return cosin * tgtThick > theSmoothThickToElemSizeRatio * elemSize; + } /*! * \brief SMESH_ProxyMesh computed by _ViscousBuilder for a SOLID. @@ -425,7 +431,7 @@ namespace VISCOUS_3D set _reversedFaceIds; set _ignoreFaceIds; // WOL FACEs and FACEs of other SOLIDS - double _stepSize, _stepSizeCoeff; + double _stepSize, _stepSizeCoeff, _geomSize; const SMDS_MeshNode* _stepSizeNodes[2]; TNode2Edge _n2eMap; // nodes and _LayerEdge's based on them @@ -589,6 +595,7 @@ namespace VISCOUS_3D const bool toSort = false); void findSimplexTestEdges( _SolidData& data, vector< vector<_LayerEdge*> >& edgesByGeom); + void computeGeomSize( _SolidData& data ); bool sortEdges( _SolidData& data, vector< vector<_LayerEdge*> >& edgesByGeom); void limitStepSizeByCurvature( _SolidData& data ); @@ -724,6 +731,7 @@ namespace VISCOUS_3D return _surface->Value( uv.X(), uv.Y() ).XYZ(); } }; + } // namespace VISCOUS_3D @@ -849,22 +857,32 @@ namespace gp_Vec dir; double f,l; gp_Pnt p; Handle(Geom_Curve) c = BRep_Tool::Curve( E, f, l ); + if ( c.IsNull() ) return gp_XYZ( 1e100, 1e100, 1e100 ); double u = helper.GetNodeU( E, atNode ); c->D1( u, p, dir ); return dir.XYZ(); } //-------------------------------------------------------------------------------- + gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV, + const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok, + double* cosin=0); + //-------------------------------------------------------------------------------- gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Edge& fromE, const SMDS_MeshNode* node, SMESH_MesherHelper& helper, bool& ok) { + double f,l; + Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l ); + if ( c.IsNull() ) + { + TopoDS_Vertex v = helper.IthVertex( 0, fromE ); + return getFaceDir( F, v, node, helper, ok ); + } gp_XY uv = helper.GetNodeUV( F, node, 0, &ok ); Handle(Geom_Surface) surface = BRep_Tool::Surface( F ); gp_Pnt p; gp_Vec du, dv, norm; surface->D1( uv.X(),uv.Y(), p, du,dv ); norm = du ^ dv; - double f,l; - Handle(Geom_Curve) c = BRep_Tool::Curve( fromE, f, l ); double u = helper.GetNodeU( fromE, node, 0, &ok ); c->D1( u, p, du ); TopAbs_Orientation o = helper.GetSubShapeOri( F.Oriented(TopAbs_FORWARD), fromE); @@ -888,7 +906,7 @@ namespace //-------------------------------------------------------------------------------- gp_XYZ getFaceDir( const TopoDS_Face& F, const TopoDS_Vertex& fromV, const SMDS_MeshNode* node, SMESH_MesherHelper& helper, - bool& ok, double* cosin=0) + bool& ok, double* cosin) { TopoDS_Face faceFrw = F; faceFrw.Orientation( TopAbs_FORWARD ); @@ -1019,6 +1037,59 @@ namespace } return false; } + //================================================================================ + /*! + * \brief Computes mimimal distance of face in-FACE nodes from an EDGE + * \param [in] face - the mesh face to treat + * \param [in] nodeOnEdge - a node on the EDGE + * \param [out] faceSize - the computed distance + * \return bool - true if faceSize computed + */ + //================================================================================ + + bool getDistFromEdge( const SMDS_MeshElement* face, + const SMDS_MeshNode* nodeOnEdge, + double & faceSize ) + { + faceSize = Precision::Infinite(); + bool done = false; + + int nbN = face->NbCornerNodes(); + int iOnE = face->GetNodeIndex( nodeOnEdge ); + int iNext[2] = { SMESH_MesherHelper::WrapIndex( iOnE+1, nbN ), + SMESH_MesherHelper::WrapIndex( iOnE-1, nbN ) }; + const SMDS_MeshNode* nNext[2] = { face->GetNode( iNext[0] ), + face->GetNode( iNext[1] ) }; + gp_XYZ segVec, segEnd = SMESH_TNodeXYZ( nodeOnEdge ); // segment on EDGE + double segLen = -1.; + // look for two neighbor not in-FACE nodes of face + for ( int i = 0; i < 2; ++i ) + { + if ( nNext[i]->GetPosition()->GetDim() != 2 && + nNext[i]->GetID() < nodeOnEdge->GetID() ) + { + // look for an in-FACE node + for ( int iN = 0; iN < nbN; ++iN ) + { + if ( iN == iOnE || iN == iNext[i] ) + continue; + SMESH_TNodeXYZ pInFace = face->GetNode( iN ); + gp_XYZ v = pInFace - segEnd; + if ( segLen < 0 ) + { + segVec = SMESH_TNodeXYZ( nNext[i] ) - segEnd; + segLen = segVec.Modulus(); + } + double distToSeg = v.Crossed( segVec ).Modulus() / segLen; + faceSize = Min( faceSize, distToSeg ); + done = true; + } + segLen = -1; + } + } + return done; + } + //-------------------------------------------------------------------------------- // DEBUG. Dump intermediate node positions into a python script // HOWTO use: run python commands written in a console to see @@ -1643,22 +1714,43 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) while ( eIt->more() ) { const SMDS_MeshElement* face = eIt->next(); + double faceMaxCosin = -1; + _LayerEdge* maxCosinEdge = 0; + int nbDegenNodes = 0; + newNodes.resize( face->NbCornerNodes() ); - double faceMaxCosin = -1; - _LayerEdge* maxCosinEdge = 0; - for ( int i = 0 ; i < face->NbCornerNodes(); ++i ) + for ( size_t i = 0 ; i < newNodes.size(); ++i ) { - const SMDS_MeshNode* n = face->GetNode(i); + const SMDS_MeshNode* n = face->GetNode( i ); + const int shapeID = n->getshapeId(); + const bool onDegenShap = helper.IsDegenShape( shapeID ); + const bool onDegenEdge = ( onDegenShap && n->GetPosition()->GetDim() == 1 ); + if ( onDegenShap ) + { + if ( onDegenEdge ) + { + // substitute n on a degenerated EDGE with a node on a corresponding VERTEX + const TopoDS_Shape& E = getMeshDS()->IndexToShape( shapeID ); + TopoDS_Vertex V = helper.IthVertex( 0, TopoDS::Edge( E )); + if ( const SMDS_MeshNode* vN = SMESH_Algo::VertexNode( V, getMeshDS() )) { + n = vN; + nbDegenNodes++; + } + } + else + { + nbDegenNodes++; + } + } TNode2Edge::iterator n2e = data._n2eMap.insert( make_pair( n, (_LayerEdge*)0 )).first; if ( !(*n2e).second ) { // add a _LayerEdge _LayerEdge* edge = new _LayerEdge(); - n2e->second = edge; edge->_nodes.push_back( n ); - const int shapeID = n->getshapeId(); - const bool noShrink = data._noShrinkShapes.count( shapeID ); + n2e->second = edge; edgesByGeom[ shapeID ].push_back( edge ); + const bool noShrink = data._noShrinkShapes.count( shapeID ); SMESH_TNodeXYZ xyz( n ); @@ -1693,7 +1785,13 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) } } newNodes[ i ] = n2e->second->_nodes.back(); + + if ( onDegenEdge ) + data._n2eMap.insert( make_pair( face->GetNode( i ), n2e->second )); } + if ( newNodes.size() - nbDegenNodes < 2 ) + continue; + // create a temporary face const SMDS_MeshElement* newFace = new _TmpMeshFace( newNodes, --_tmpFaceID, face->getshapeId() ); @@ -1702,6 +1800,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) // compute inflation step size by min size of element on a convex surface if ( faceMaxCosin > theMinSmoothCosin ) limitStepSize( data, face, maxCosinEdge ); + } // loop on 2D elements on a FACE } // loop on FACEs of a SOLID @@ -1721,16 +1820,17 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) const SMDS_MeshNode* nn[2]; for ( size_t i = 0; i < data._edges.size(); ++i ) { - if ( data._edges[i]->IsOnEdge() ) + _LayerEdge* edge = data._edges[i]; + if ( edge->IsOnEdge() ) { // get neighbor nodes - bool hasData = ( data._edges[i]->_2neibors->_edges[0] ); + bool hasData = ( edge->_2neibors->_edges[0] ); if ( hasData ) // _LayerEdge is a copy of another one { - nn[0] = data._edges[i]->_2neibors->srcNode(0); - nn[1] = data._edges[i]->_2neibors->srcNode(1); + nn[0] = edge->_2neibors->srcNode(0); + nn[1] = edge->_2neibors->srcNode(1); } - else if ( !findNeiborsOnEdge( data._edges[i], nn[0],nn[1], data )) + else if ( !findNeiborsOnEdge( edge, nn[0],nn[1], data )) { return false; } @@ -1739,18 +1839,37 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) { if (( n2e = data._n2eMap.find( nn[j] )) == data._n2eMap.end() ) return error("_LayerEdge not found by src node", data._index); - data._edges[i]->_2neibors->_edges[j] = n2e->second; + edge->_2neibors->_edges[j] = n2e->second; } if ( !hasData ) - data._edges[i]->SetDataByNeighbors( nn[0], nn[1], helper); + edge->SetDataByNeighbors( nn[0], nn[1], helper); } - for ( size_t j = 0; j < data._edges[i]->_simplices.size(); ++j ) + for ( size_t j = 0; j < edge->_simplices.size(); ++j ) { - _Simplex& s = data._edges[i]->_simplices[j]; + _Simplex& s = edge->_simplices[j]; s._nNext = data._n2eMap[ s._nNext ]->_nodes.back(); s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back(); } + + // For an _LayerEdge on a degenerated EDGE, copy some data from + // a corresponding _LayerEdge on a VERTEX + // (issue 52453, pb on a downloaded SampleCase2-Tet-netgen-mephisto.hdf) + if ( helper.IsDegenShape( edge->_nodes[0]->getshapeId() )) + { + // Generally we should not get here + const TopoDS_Shape& E = getMeshDS()->IndexToShape( edge->_nodes[0]->getshapeId() ); + if ( E.ShapeType() != TopAbs_EDGE ) + continue; + TopoDS_Vertex V = helper.IthVertex( 0, TopoDS::Edge( E )); + const SMDS_MeshNode* vN = SMESH_Algo::VertexNode( V, getMeshDS() ); + if (( n2e = data._n2eMap.find( vN )) == data._n2eMap.end() ) + continue; + const _LayerEdge* vEdge = n2e->second; + edge->_normal = vEdge->_normal; + edge->_lenFactor = vEdge->_lenFactor; + edge->_cosin = vEdge->_cosin; + } } dumpFunctionEnd(); @@ -1948,8 +2067,12 @@ void _ViscousBuilder::limitStepSizeByCurvature( _SolidData& data ) bool _ViscousBuilder::sortEdges( _SolidData& data, vector< vector<_LayerEdge*> >& edgesByGeom) { + // define allowed thickness + computeGeomSize( data ); // compute data._geomSize + const double tgtThick = Min( 0.5 * data._geomSize, data._hyp->GetTotalThickness() ); + // Find shapes needing smoothing; such a shape has _LayerEdge._normal on it's - // boundry inclined at a sharp angle to the shape + // boundry inclined to the shape at a sharp angle list< TGeomID > shapesToSmooth; @@ -1966,31 +2089,33 @@ bool _ViscousBuilder::sortEdges( _SolidData& data, { case TopAbs_EDGE: { - bool isShrinkEdge = !eS[0]->_sWOL.IsNull(); + if ( SMESH_Algo::isDegenerated( TopoDS::Edge( S ))) + break; + //bool isShrinkEdge = !eS[0]->_sWOL.IsNull(); for ( TopoDS_Iterator vIt( S ); vIt.More() && !needSmooth; vIt.Next() ) { TGeomID iV = getMeshDS()->ShapeToIndex( vIt.Value() ); vector<_LayerEdge*>& eV = edgesByGeom[ iV ]; if ( eV.empty() ) continue; - // double cosin = eV[0]->_cosin; - // bool badCosin = - // ( !eV[0]->_sWOL.IsNull() && ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE || !isShrinkEdge)); - // if ( badCosin ) - // { - // gp_Vec dir1, dir2; - // if ( eV[0]->_sWOL.ShapeType() == TopAbs_EDGE ) - // dir1 = getEdgeDir( TopoDS::Edge( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() )); - // else - // dir1 = getFaceDir( TopoDS::Face( eV[0]->_sWOL ), TopoDS::Vertex( vIt.Value() ), - // eV[0]->_nodes[0], helper, ok); - // dir2 = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() )); - // double angle = dir1.Angle( dir2 ); - // cosin = cos( angle ); - // } gp_Vec eDir = getEdgeDir( TopoDS::Edge( S ), TopoDS::Vertex( vIt.Value() )); double angle = eDir.Angle( eV[0]->_normal ); double cosin = Cos( angle ); - needSmooth = ( cosin > theMinSmoothCosin ); + if ( cosin > theMinSmoothCosin ) + { + // compare tgtThick with the length of an end segment + SMDS_ElemIteratorPtr eIt = eV[0]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Edge); + while ( eIt->more() ) + { + const SMDS_MeshElement* endSeg = eIt->next(); + if ( endSeg->getshapeId() == iS ) + { + double segLen = + SMESH_TNodeXYZ( endSeg->GetNode(0) ).Distance( endSeg->GetNode(1 )); + needSmooth = needSmoothing( cosin, tgtThick, segLen ); + break; + } + } + } } break; } @@ -2001,25 +2126,38 @@ bool _ViscousBuilder::sortEdges( _SolidData& data, TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() ); vector<_LayerEdge*>& eE = edgesByGeom[ iE ]; if ( eE.empty() ) continue; - if ( eE[0]->_sWOL.IsNull() ) + // TopLoc_Location loc; + // Handle(Geom_Surface) surface = BRep_Tool::Surface( TopoDS::Face( S ), loc ); + // bool isPlane = GeomLib_IsPlanarSurface( surface ).IsPlanar(); + //if ( eE[0]->_sWOL.IsNull() ) { + double faceSize; for ( size_t i = 0; i < eE.size() && !needSmooth; ++i ) - needSmooth = ( eE[i]->_cosin > theMinSmoothCosin ); - } - else - { - const TopoDS_Face& F1 = TopoDS::Face( S ); - const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL ); - const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() ); - for ( size_t i = 0; i < eE.size() && !needSmooth; ++i ) - { - gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok ); - gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok ); - double angle = dir1.Angle( dir2 ); - double cosin = cos( angle ); - needSmooth = ( cosin > theMinSmoothCosin ); - } + if ( eE[i]->_cosin > theMinSmoothCosin ) + { + SMDS_ElemIteratorPtr fIt = eE[i]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face); + while ( fIt->more() && !needSmooth ) + { + const SMDS_MeshElement* face = fIt->next(); + if ( getDistFromEdge( face, eE[i]->_nodes[0], faceSize )) + needSmooth = needSmoothing( eE[i]->_cosin, tgtThick, faceSize ); + } + } } + // else + // { + // const TopoDS_Face& F1 = TopoDS::Face( S ); + // const TopoDS_Face& F2 = TopoDS::Face( eE[0]->_sWOL ); + // const TopoDS_Edge& E = TopoDS::Edge( eExp.Current() ); + // for ( size_t i = 0; i < eE.size() && !needSmooth; ++i ) + // { + // gp_Vec dir1 = getFaceDir( F1, E, eE[i]->_nodes[0], helper, ok ); + // gp_Vec dir2 = getFaceDir( F2, E, eE[i]->_nodes[0], helper, ok ); + // double angle = dir1.Angle( ); + // double cosin = cos( angle ); + // needSmooth = ( cosin > theMinSmoothCosin ); + // } + // } } break; } @@ -2027,6 +2165,7 @@ bool _ViscousBuilder::sortEdges( _SolidData& data, continue; default:; } + if ( needSmooth ) { if ( S.ShapeType() == TopAbs_EDGE ) shapesToSmooth.push_front( iS ); @@ -2354,9 +2493,19 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node, Handle(Geom_Surface) surface = BRep_Tool::Surface( face ); int pointKind = GeomLib::NormEstim( surface, uv, 1e-5, normal ); enum { REGULAR = 0, QUASYSINGULAR, CONICAL, IMPOSSIBLE }; + + if ( pointKind == IMPOSSIBLE && + node->GetPosition()->GetDim() == 2 ) // node inside the FACE + { + // probably NormEstim() failed due to a too high tolerance + pointKind = GeomLib::NormEstim( surface, uv, 1e-20, normal ); + isOK = ( pointKind < IMPOSSIBLE ); + } if ( pointKind < IMPOSSIBLE ) { - if ( pointKind != REGULAR && !shiftInside ) + if ( pointKind != REGULAR && + !shiftInside && + node->GetPosition()->GetDim() < 2 ) // FACE boundary { gp_XYZ normShift = getFaceNormal( node, face, helper, isOK, /*shiftInside=*/true ); if ( normShift * normal.XYZ() < 0. ) @@ -2364,7 +2513,8 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node, } isOK = true; } - else // hard singularity, to call with shiftInside=true ? + + if ( !isOK ) // hard singularity, to call with shiftInside=true ? { const TGeomID faceID = helper.GetMeshDS()->ShapeToIndex( face ); @@ -2377,7 +2527,9 @@ gp_XYZ _ViscousBuilder::getFaceNormal(const SMDS_MeshNode* node, isOK = SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) normal.XYZ(), /*normalized=*/true ); if ( isOK ) { - if ( helper.IsReversedSubMesh( face )) + TopoDS_Face ff = face; + ff.Orientation( TopAbs_FORWARD ); + if ( helper.IsReversedSubMesh( ff )) normal.Reverse(); break; } @@ -2529,11 +2681,11 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1, // Set _curvature - double sumLen = vec1.Modulus() + vec2.Modulus(); + double sumLen = vec1.Modulus() + vec2.Modulus(); _2neibors->_wgt[0] = 1 - vec1.Modulus() / sumLen; _2neibors->_wgt[1] = 1 - vec2.Modulus() / sumLen; double avgNormProj = 0.5 * ( _normal * vec1 + _normal * vec2 ); - double avgLen = 0.5 * ( vec1.Modulus() + vec2.Modulus() ); + double avgLen = 0.5 * ( vec1.Modulus() + vec2.Modulus() ); if ( _curvature ) delete _curvature; _curvature = _Curvature::New( avgNormProj, avgLen ); // if ( _curvature ) @@ -2547,10 +2699,13 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1, if ( _sWOL.IsNull() ) { TopoDS_Shape S = helper.GetSubShapeByNode( _nodes[0], helper.GetMeshDS() ); - gp_XYZ dirE = getEdgeDir( TopoDS::Edge( S ), _nodes[0], helper ); + TopoDS_Edge E = TopoDS::Edge( S ); + // if ( SMESH_Algo::isDegenerated( E )) + // return; + gp_XYZ dirE = getEdgeDir( E, _nodes[0], helper ); gp_XYZ plnNorm = dirE ^ _normal; - double proj0 = plnNorm * vec1; - double proj1 = plnNorm * vec2; + double proj0 = plnNorm * vec1; + double proj1 = plnNorm * vec2; if ( fabs( proj0 ) > 1e-10 || fabs( proj1 ) > 1e-10 ) { if ( _2neibors->_plnNorm ) delete _2neibors->_plnNorm; @@ -2717,6 +2872,31 @@ void _ViscousBuilder::makeGroupOfLE() #endif } +//================================================================================ +/*! + * \brief Find maximal _LayerEdge length (layer thickness) limited by geometry + */ +//================================================================================ + +void _ViscousBuilder::computeGeomSize( _SolidData& data ) +{ + data._geomSize = Precision::Infinite(); + double intersecDist; + auto_ptr searcher + ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), + data._proxyMesh->GetFaces( data._solid )) ); + + TNode2Edge::iterator n2e = data._n2eMap.begin(), n2eEnd = data._n2eMap.end(); + for ( ; n2e != n2eEnd; ++n2e ) + { + _LayerEdge* edge = n2e->second; + if ( edge->IsOnEdge() ) continue; + edge->FindIntersection( *searcher, intersecDist, data._epsilon ); + if ( data._geomSize > intersecDist && intersecDist > 0 ) + data._geomSize = intersecDist; + } +} + //================================================================================ /*! * \brief Increase length of _LayerEdge's to reach the required thickness of layers @@ -2729,19 +2909,8 @@ 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; - auto_ptr searcher - ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), - data._proxyMesh->GetFaces( data._solid )) ); - for ( size_t i = 0; i < data._edges.size(); ++i ) - { - if ( data._edges[i]->IsOnEdge() ) continue; - data._edges[i]->FindIntersection( *searcher, intersecDist, data._epsilon ); - if ( geomSize > intersecDist && intersecDist > 0 ) - geomSize = intersecDist; - } - if ( data._stepSize > 0.3 * geomSize ) - limitStepSize( data, 0.3 * geomSize ); + if ( data._stepSize > 0.3 * data._geomSize ) + limitStepSize( data, 0.3 * data._geomSize ); const double tgtThick = data._hyp->GetTotalThickness(); if ( data._stepSize > tgtThick ) @@ -2750,7 +2919,7 @@ bool _ViscousBuilder::inflate(_SolidData& data) if ( data._stepSize < 1. ) data._epsilon = data._stepSize * 1e-7; - debugMsg( "-- geomSize = " << geomSize << ", stepSize = " << data._stepSize ); + debugMsg( "-- geomSize = " << data._geomSize << ", stepSize = " << data._stepSize ); double avgThick = 0, curThick = 0, distToIntersection = Precision::Infinite(); int nbSteps = 0, nbRepeats = 0; @@ -4819,6 +4988,7 @@ bool _ViscousBuilder::refine(_SolidData& data) helper.SetElementsOnShape(true); vector< vector* > nnVec; + set< vector* > nnSet; set< int > degenEdgeInd; vector degenVols; @@ -4834,20 +5004,26 @@ bool _ViscousBuilder::refine(_SolidData& data) const SMDS_MeshElement* face = fIt->next(); const int nbNodes = face->NbCornerNodes(); nnVec.resize( nbNodes ); + nnSet.clear(); degenEdgeInd.clear(); int nbZ = 0; - SMDS_ElemIteratorPtr nIt = face->nodesIterator(); + SMDS_NodeIteratorPtr nIt = face->nodeIterator(); for ( int iN = 0; iN < nbNodes; ++iN ) { - const SMDS_MeshNode* n = static_cast( nIt->next() ); + const SMDS_MeshNode* n = nIt->next(); nnVec[ iN ] = & data._n2eMap[ n ]->_nodes; if ( nnVec[ iN ]->size() < 2 ) degenEdgeInd.insert( iN ); else nbZ = nnVec[ iN ]->size(); + + if ( helper.HasDegeneratedEdges() ) + nnSet.insert( nnVec[ iN ]); } if ( nbZ == 0 ) continue; + if ( 0 < nnSet.size() && nnSet.size() < 3 ) + continue; switch ( nbNodes ) { @@ -6057,6 +6233,8 @@ bool _ViscousBuilder::addBoundaryElements() { SMESH_MesherHelper helper( *_mesh ); + vector< const SMDS_MeshNode* > faceNodes; + for ( size_t i = 0; i < _sdVec.size(); ++i ) { _SolidData& data = _sdVec[i]; @@ -6099,8 +6277,13 @@ bool _ViscousBuilder::addBoundaryElements() if ( nbSharedPyram > 1 ) continue; // not free border of the pyramid - if ( getMeshDS()->FindFace( ledges[0]->_nodes[0], ledges[0]->_nodes[1], - ledges[1]->_nodes[0], ledges[1]->_nodes[1])) + faceNodes.clear(); + faceNodes.push_back( ledges[0]->_nodes[0] ); + faceNodes.push_back( ledges[1]->_nodes[0] ); + if ( ledges[0]->_nodes.size() > 1 ) faceNodes.push_back( ledges[0]->_nodes[1] ); + if ( ledges[1]->_nodes.size() > 1 ) faceNodes.push_back( ledges[1]->_nodes[1] ); + + if ( getMeshDS()->FindElement( faceNodes, SMDSAbs_Face, /*noMedium=*/true)) continue; // faces already created } for ( ++u2n; u2n != u2nodes.end(); ++u2n )