X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ViscousLayers.cxx;h=b2f6855b01c065ec0245b23b6a550c94e1552f25;hp=d4f0f678f3e188a36dd7dbf1c4ece2aaad1dd120;hb=050aa87698efb9b123985f0bd69aab20ad42c728;hpb=0de979b4949cfda1ed932273abb9a92ecc71a4ec diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index d4f0f678f..b2f6855b0 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -45,6 +45,8 @@ #include "StdMeshers_FaceSide.hxx" #include +#include +#include #include #include #include @@ -54,12 +56,14 @@ #include #include #include +#include #include #include #include #include #include #include +#include #include #include #include @@ -322,7 +326,7 @@ namespace VISCOUS_3D gp_XYZ _normal; // to solid surface vector _pos; // points computed during inflation - double _len; // length achived with the last step + double _len; // length achived with the last inflation step double _cosin; // of angle (_normal ^ surface) double _lenFactor; // to compute _len taking _cosin into account @@ -361,7 +365,7 @@ namespace VISCOUS_3D const double& epsilon) const; gp_Ax1 LastSegment(double& segLen) const; bool IsOnEdge() const { return _2neibors; } - void Copy( _LayerEdge& other, SMESH_MesherHelper& helper ); + gp_XYZ Copy( _LayerEdge& other, SMESH_MesherHelper& helper ); void SetCosin( double cosin ); }; struct _LayerEdgeCmp @@ -393,9 +397,13 @@ namespace VISCOUS_3D const SMDS_MeshNode* _stepSizeNodes[2]; TNode2Edge _n2eMap; + // map to find _n2eMap of another _SolidData by a shrink shape shared by two _SolidData's + map< TGeomID, TNode2Edge* > _s2neMap; // edges of _n2eMap. We keep same data in two containers because // iteration over the map is 5 time longer than over the vector vector< _LayerEdge* > _edges; + // edges on EDGE's with null _sWOL, whose _simplices are used to stop inflation + vector< _LayerEdge* > _simplexTestEdges; // key: an id of shape (EDGE or VERTEX) shared by a FACE with // layers and a FACE w/o layers @@ -478,6 +486,9 @@ namespace VISCOUS_3D bool makeLayer(_SolidData& data); bool setEdgeData(_LayerEdge& edge, const set& subIds, SMESH_MesherHelper& helper, _SolidData& data); + gp_XYZ getWeigthedNormal( const SMDS_MeshNode* n, + std::pair< TGeomID, gp_XYZ > fId2Normal[], + const int nbFaces ); bool findNeiborsOnEdge(const _LayerEdge* edge, const SMDS_MeshNode*& n1, const SMDS_MeshNode*& n2, @@ -486,6 +497,8 @@ namespace VISCOUS_3D const set& ingnoreShapes, const _SolidData* dataToCheckOri = 0, const bool toSort = false); + void findSimplexTestEdges( _SolidData& data, + vector< vector<_LayerEdge*> >& edgesByGeom); bool sortEdges( _SolidData& data, vector< vector<_LayerEdge*> >& edgesByGeom); void limitStepSize( _SolidData& data, @@ -793,23 +806,29 @@ namespace // get average dir of edges going fromV gp_XYZ edgeDir; //if ( edges.size() > 1 ) - for ( size_t 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; - } + for ( size_t 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 = 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; + try { + if ( edges.size() == 1 ) + dir = fromEdgeDir; + else if ( dir.SquareModulus() < 0.1 ) // ~< 20 degrees + dir = fromEdgeDir.Normalized() + getFaceDir( F, edges[1], node, helper, ok ).Normalized(); + else if ( dir * fromEdgeDir < 0 ) + dir *= -1; + } + catch ( Standard_Failure ) + { + ok = false; + } if ( ok ) { //dir /= edges.size(); @@ -1388,6 +1407,19 @@ bool _ViscousBuilder::findFacesWithLayers() } } + // add FACEs of other SOLIDs to _ignoreFaceIds + for ( size_t i = 0; i < _sdVec.size(); ++i ) + { + shapes.Clear(); + TopExp::MapShapes(_sdVec[i]._solid, TopAbs_FACE, shapes); + + for ( exp.Init( _mesh->GetShapeToMesh(), TopAbs_FACE ); exp.More(); exp.Next() ) + { + if ( !shapes.Contains( exp.Current() )) + _sdVec[i]._ignoreFaceIds.insert( getMeshDS()->ShapeToIndex( exp.Current() )); + } + } + return true; } @@ -1415,7 +1447,6 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) } // make a map to find new nodes on sub-shapes shared with other SOLID - map< TGeomID, TNode2Edge* > s2neMap; map< TGeomID, TNode2Edge* >::iterator s2ne; map< TGeomID, TopoDS_Shape >::iterator s2s = data._shrinkShape2Shape.begin(); for (; s2s != data._shrinkShape2Shape.end(); ++s2s ) @@ -1428,7 +1459,7 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) if ( s2s2 != _sdVec[i]._shrinkShape2Shape.end() && *s2s == *s2s2 && !_sdVec[i]._n2eMap.empty() ) { - s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap )); + data._s2neMap.insert( make_pair( shapeInd, &_sdVec[i]._n2eMap )); break; } } @@ -1480,21 +1511,23 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) const int shapeID = n->getshapeId(); edgesByGeom[ shapeID ].push_back( edge ); + SMESH_TNodeXYZ xyz( n ); + // set edge data or find already refined _LayerEdge and get data from it if ( n->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE && - ( s2ne = s2neMap.find( shapeID )) != s2neMap.end() && + ( s2ne = data._s2neMap.find( shapeID )) != data._s2neMap.end() && ( n2e2 = (*s2ne).second->find( n )) != s2ne->second->end()) { _LayerEdge* foundEdge = (*n2e2).second; - edge->Copy( *foundEdge, helper ); - // location of the last node is modified but we can restore - // it by node position on _sWOL stored by the node + gp_XYZ lastPos = edge->Copy( *foundEdge, helper ); + foundEdge->_pos.push_back( lastPos ); + // location of the last node is modified and we restore it by foundEdge->_pos.back() const_cast< SMDS_MeshNode* > - ( edge->_nodes.back() )->setXYZ( n->X(), n->Y(), n->Z() ); + ( edge->_nodes.back() )->setXYZ( xyz.X(), xyz.Y(), xyz.Z() ); } else { - edge->_nodes.push_back( helper.AddNode( n->X(), n->Y(), n->Z() )); + edge->_nodes.push_back( helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() )); if ( !setEdgeData( *edge, subIds, helper, data )) return false; } @@ -1521,12 +1554,14 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) if ( data._stepSize < 1. ) data._epsilon *= data._stepSize; - // Put _LayerEdge's into the vector data._edges + // fill data._simplexTestEdges + findSimplexTestEdges( data, edgesByGeom ); + // Put _LayerEdge's into the vector data._edges if ( !sortEdges( data, edgesByGeom )) return false; - // Set target nodes into _Simplex and _2NearEdges + // Set target nodes into _Simplex and _2NearEdges of _LayerEdge's TNode2Edge::iterator n2e; for ( size_t i = 0; i < data._edges.size(); ++i ) { @@ -1541,13 +1576,13 @@ bool _ViscousBuilder::makeLayer(_SolidData& data) n = (*n2e).second->_nodes.back(); data._edges[i]->_2neibors->_edges[j] = n2e->second; } - else - for ( size_t j = 0; j < data._edges[i]->_simplices.size(); ++j ) - { - _Simplex& s = data._edges[i]->_simplices[j]; - s._nNext = data._n2eMap[ s._nNext ]->_nodes.back(); - s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back(); - } + //else + for ( size_t j = 0; j < data._edges[i]->_simplices.size(); ++j ) + { + _Simplex& s = data._edges[i]->_simplices[j]; + s._nNext = data._n2eMap[ s._nNext ]->_nodes.back(); + s._nPrev = data._n2eMap[ s._nPrev ]->_nodes.back(); + } } dumpFunctionEnd(); @@ -1595,7 +1630,7 @@ void _ViscousBuilder::limitStepSize( _SolidData& data, */ //================================================================================ -void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize) +void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize ) { if ( minSize < data._stepSize ) { @@ -1609,6 +1644,96 @@ void _ViscousBuilder::limitStepSize( _SolidData& data, const double minSize) } } +//================================================================================ +/*! + * Fill data._simplexTestEdges. These _LayerEdge's are used to stop inflation + * in the case where there are no _LayerEdge's on a curved convex FACE, + * as e.g. on a fillet surface with no internal nodes - issue 22580, + * so that collision of viscous internal faces is not detected by check of + * intersection of _LayerEdge's with the viscous internal faces. + */ +//================================================================================ + +void _ViscousBuilder::findSimplexTestEdges( _SolidData& data, + vector< vector<_LayerEdge*> >& edgesByGeom) +{ + data._simplexTestEdges.clear(); + + SMESH_MesherHelper helper( *_mesh ); + + vector< vector<_LayerEdge*> * > ledgesOnEdges; + set< const SMDS_MeshNode* > usedNodes; + + const double minCurvature = 1. / data._hyp->GetTotalThickness(); + + for ( size_t iS = 0; iS < edgesByGeom.size(); ++iS ) + { + // look for a FACE with layers and w/o _LayerEdge's + const vector<_LayerEdge*>& eS = edgesByGeom[iS]; + if ( !eS.empty() ) continue; + const TopoDS_Shape& S = getMeshDS()->IndexToShape( iS ); + if ( S.IsNull() || S.ShapeType() != TopAbs_FACE ) continue; + if ( data._ignoreFaceIds.count( iS )) continue; + + const TopoDS_Face& F = TopoDS::Face( S ); + + // look for _LayerEdge's on EDGEs with null _sWOL + ledgesOnEdges.clear(); + TopExp_Explorer eExp( F, TopAbs_EDGE ); + for ( ; eExp.More(); eExp.Next() ) + { + TGeomID iE = getMeshDS()->ShapeToIndex( eExp.Current() ); + vector<_LayerEdge*>& eE = edgesByGeom[iE]; + if ( !eE.empty() && eE[0]->_sWOL.IsNull() ) + ledgesOnEdges.push_back( & eE ); + } + if ( ledgesOnEdges.empty() ) continue; + + // check FACE convexity + const _LayerEdge* le = ledgesOnEdges[0]->back(); + gp_XY uv = helper.GetNodeUV( F, le->_nodes[0] ); + BRepAdaptor_Surface surf( F ); + BRepLProp_SLProps surfProp( surf, uv.X(), uv.Y(), 2, 1e-6 ); + if ( !surfProp.IsCurvatureDefined() ) + continue; + double surfCurvature = Max( Abs( surfProp.MaxCurvature() ), + Abs( surfProp.MinCurvature() )); + if ( surfCurvature < minCurvature ) + continue; + gp_Dir minDir, maxDir; + surfProp.CurvatureDirections( maxDir, minDir ); + if ( F.Orientation() == TopAbs_REVERSED ) { + maxDir.Reverse(); minDir.Reverse(); + } + const gp_XYZ& inDir = le->_normal; + if ( inDir * maxDir.XYZ() < 0 && + inDir * minDir.XYZ() < 0 ) + continue; + + limitStepSize( data, 0.9 / surfCurvature ); + + // add _simplices to the _LayerEdge's + for ( size_t iE = 0; iE < ledgesOnEdges.size(); ++iE ) + { + const vector<_LayerEdge*>& ledges = *ledgesOnEdges[iE]; + for ( size_t iLE = 0; iLE < ledges.size(); ++iLE ) + { + _LayerEdge* ledge = ledges[iLE]; + const SMDS_MeshNode* srcNode = ledge->_nodes[0]; + if ( !usedNodes.insert( srcNode ).second ) continue; + + getSimplices( srcNode, ledge->_simplices, data._ignoreFaceIds, &data ); + for ( size_t i = 0; i < ledge->_simplices.size(); ++i ) + { + usedNodes.insert( ledge->_simplices[i]._nPrev ); + usedNodes.insert( ledge->_simplices[i]._nNext ); + } + data._simplexTestEdges.push_back( ledge ); + } + } + } +} + //================================================================================ /*! * \brief Separate shapes (and _LayerEdge's on them) to smooth from the rest ones @@ -1630,7 +1755,7 @@ bool _ViscousBuilder::sortEdges( _SolidData& data, { vector<_LayerEdge*>& eS = edgesByGeom[iS]; if ( eS.empty() ) continue; - TopoDS_Shape S = getMeshDS()->IndexToShape( iS ); + const TopoDS_Shape& S = getMeshDS()->IndexToShape( iS ); bool needSmooth = false; switch ( S.ShapeType() ) { @@ -1803,61 +1928,62 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, set::iterator id = faceIds.begin(); TopoDS_Face F; + std::pair< TGeomID, gp_XYZ > id2Norm[20]; for ( ; id != faceIds.end(); ++id ) { const TopoDS_Shape& s = getMeshDS()->IndexToShape( *id ); if ( s.IsNull() || s.ShapeType() != TopAbs_FACE || !subIds.count( *id )) continue; - totalNbFaces++; F = TopoDS::Face( s ); - // IDEA: if there is a problem with finding a normal, - // we can compute an area-weighted sum of normals of all faces sharing the node gp_XY uv = helper.GetNodeUV( F, node, 0, &normOK ); Handle(Geom_Surface) surface = BRep_Tool::Surface( F ); - surface->D1( uv.X(), uv.Y(), p, du,dv ); - geomNorm = du ^ dv; - double size2 = geomNorm.SquareMagnitude(); - if ( size2 < 1e-10 ) // singularity { - SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face); - while ( fIt->more() ) + gp_Dir normal; + if ( GeomLib::NormEstim( surface, uv, 1e-10, normal ) < 3 ) { - const SMDS_MeshElement* f = fIt->next(); - if ( editor.FindShape( f ) == *id ) + geomNorm = normal; + normOK = true; + } + else // hard singularity + { + SMDS_ElemIteratorPtr fIt = node->GetInverseElementIterator(SMDSAbs_Face); + while ( fIt->more() ) { - SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) geomNorm.XYZ(), /*normalized=*/false ); - size2 = geomNorm.SquareMagnitude(); - break; + const SMDS_MeshElement* f = fIt->next(); + if ( editor.FindShape( f ) == *id ) + { + SMESH_MeshAlgos::FaceNormal( f, (gp_XYZ&) geomNorm.XYZ(), /*normalized=*/false ); + if ( helper.IsReversedSubMesh( F )) + geomNorm.Reverse(); + break; + } } + double size2 = geomNorm.SquareMagnitude(); + if ( size2 > numeric_limits::min() ) + geomNorm /= sqrt( size2 ); + else + normOK = false; } - // double ddu = 0, ddv = 0; - // if ( du.SquareMagnitude() > dv.SquareMagnitude() ) - // ddu = 1e-3; - // else - // ddv = 1e-3; - // surface->D1( uv.X()+ddu, uv.Y()+ddv, p, du,dv ); - // geomNorm = du ^ dv; - // size2 = geomNorm.SquareMagnitude(); - // if ( size2 < 1e-10 ) - // { - // surface->D1( uv.X()-ddu, uv.Y()-ddv, p, du,dv ); - // geomNorm = du ^ dv; - // size2 = geomNorm.SquareMagnitude(); - // } } - if ( size2 > numeric_limits::min() ) - geomNorm /= sqrt( size2 ); - else - normOK = false; if ( helper.GetSubShapeOri( data._solid, F ) != TopAbs_REVERSED ) geomNorm.Reverse(); + id2Norm[ totalNbFaces ].first = *id; + id2Norm[ totalNbFaces ].second = geomNorm.XYZ(); + totalNbFaces++; edge._normal += geomNorm.XYZ(); } if ( totalNbFaces == 0 ) return error(SMESH_Comment("Can't get normal to node ") << node->GetID(), data._index); - edge._normal /= totalNbFaces; + if ( totalNbFaces < 3 ) + { + //edge._normal /= totalNbFaces; + } + else + { + edge._normal = getWeigthedNormal( node, id2Norm, totalNbFaces ); + } switch ( posType ) { @@ -1959,6 +2085,87 @@ bool _ViscousBuilder::setEdgeData(_LayerEdge& edge, return true; } +//================================================================================ +/*! + * \brief Return normal at a node weighted with angles taken by FACEs + * \param [in] n - the node + * \param [in] fId2Normal - FACE ids and normals + * \param [in] nbFaces - nb of FACEs meeting at the node + * \return gp_XYZ - computed normal + */ +//================================================================================ + +gp_XYZ _ViscousBuilder::getWeigthedNormal( const SMDS_MeshNode* n, + std::pair< TGeomID, gp_XYZ > fId2Normal[], + const int nbFaces ) +{ + gp_XYZ resNorm(0,0,0); + TopoDS_Shape V = SMESH_MesherHelper::GetSubShapeByNode( n, getMeshDS() ); + if ( V.ShapeType() != TopAbs_VERTEX ) + { + for ( int i = 0; i < nbFaces; ++i ) + resNorm += fId2Normal[i].second / nbFaces ; + return resNorm; + } + + double angles[30]; + for ( int i = 0; i < nbFaces; ++i ) + { + const TopoDS_Face& F = TopoDS::Face( getMeshDS()->IndexToShape( fId2Normal[i].first )); + + // look for two EDGEs shared by F and other FACEs within fId2Normal + TopoDS_Edge ee[2]; + int nbE = 0; + PShapeIteratorPtr eIt = SMESH_MesherHelper::GetAncestors( V, *_mesh, TopAbs_EDGE ); + while ( const TopoDS_Shape* E = eIt->next() ) + { + if ( !SMESH_MesherHelper::IsSubShape( *E, F )) + continue; + bool isSharedEdge = false; + for ( int j = 0; j < nbFaces && !isSharedEdge; ++j ) + { + if ( i == j ) continue; + const TopoDS_Shape& otherF = getMeshDS()->IndexToShape( fId2Normal[j].first ); + isSharedEdge = SMESH_MesherHelper::IsSubShape( *E, otherF ); + } + if ( !isSharedEdge ) + continue; + ee[ nbE ] = TopoDS::Edge( *E ); + ee[ nbE ].Orientation( SMESH_MesherHelper::GetSubShapeOri( F, *E )); + if ( ++nbE == 2 ) + break; + } + + // get an angle between the two EDGEs + angles[i] = 0; + if ( nbE < 1 ) continue; + if ( nbE == 1 ) + { + ee[ 1 ] == ee[ 0 ]; + } + else + { + TopoDS_Vertex v10 = SMESH_MesherHelper::IthVertex( 1, ee[ 0 ]); + TopoDS_Vertex v01 = SMESH_MesherHelper::IthVertex( 0, ee[ 1 ]); + if ( !v10.IsSame( v01 )) + std::swap( ee[0], ee[1] ); + } + angles[i] = SMESH_MesherHelper::GetAngle( ee[0], ee[1], F ); + } + + // compute a weighted normal + double sumAngle = 0; + for ( int i = 0; i < nbFaces; ++i ) + { + angles[i] = ( angles[i] > 2*M_PI ) ? 0 : M_PI - angles[i]; + sumAngle += angles[i]; + } + for ( int i = 0; i < nbFaces; ++i ) + resNorm += angles[i] / sumAngle * fId2Normal[i].second; + + return resNorm; +} + //================================================================================ /*! * \brief Find 2 neigbor nodes of a node on EDGE @@ -2062,7 +2269,7 @@ void _LayerEdge::SetDataByNeighbors( const SMDS_MeshNode* n1, */ //================================================================================ -void _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper ) +gp_XYZ _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper ) { _nodes = other._nodes; _normal = other._normal; @@ -2074,16 +2281,25 @@ void _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper ) _curvature = 0; std::swap( _curvature, other._curvature ); _2neibors = 0; std::swap( _2neibors, other._2neibors ); + gp_XYZ lastPos( 0,0,0 ); if ( _sWOL.ShapeType() == TopAbs_EDGE ) { double u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes[0] ); _pos.push_back( gp_XYZ( u, 0, 0)); + + u = helper.GetNodeU( TopoDS::Edge( _sWOL ), _nodes.back() ); + lastPos.SetX( u ); } else // TopAbs_FACE { gp_XY uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes[0]); _pos.push_back( gp_XYZ( uv.X(), uv.Y(), 0)); + + uv = helper.GetNodeUV( TopoDS::Face( _sWOL ), _nodes.back() ); + lastPos.SetX( uv.X() ); + lastPos.SetY( uv.Y() ); } + return lastPos; } //================================================================================ @@ -2095,7 +2311,7 @@ void _LayerEdge::Copy( _LayerEdge& other, SMESH_MesherHelper& helper ) void _LayerEdge::SetCosin( double cosin ) { _cosin = cosin; - _lenFactor = ( _cosin > 0.1 ) ? 1./sqrt(1-_cosin*_cosin) : 1.0; + _lenFactor = ( Abs( _cosin ) > 0.1 ) ? 1./sqrt(1-_cosin*_cosin) : 1.0; } //================================================================================ @@ -2387,8 +2603,8 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, else { // smooth on FACE's - int step = 0, badNb = 0; moved = true; - while (( ++step <= 5 && moved ) || improved ) + int step = 0, stepLimit = 5, badNb = 0; moved = true; + while (( ++step <= stepLimit && moved ) || improved ) { dumpFunction(SMESH_Comment("smooth")<Smooth(badNb); improved = ( badNb < oldBadNb ); + // issue 22576. no bad faces but still there are intersections to fix + if ( improved && badNb == 0 ) + stepLimit = step + 3; + dumpFunctionEnd(); } if ( badNb > 0 ) @@ -2423,6 +2643,24 @@ bool _ViscousBuilder::smoothAndCheck(_SolidData& data, } } // loop on shapes to smooth + // Check orientation of simplices of _simplexTestEdges + for ( size_t i = 0; i < data._simplexTestEdges.size(); ++i ) + { + const _LayerEdge* edge = data._simplexTestEdges[i]; + SMESH_TNodeXYZ tgtXYZ( edge->_nodes.back() ); + for ( size_t j = 0; j < edge->_simplices.size(); ++j ) + if ( !edge->_simplices[j].IsForward( edge->_nodes[0], &tgtXYZ )) + { +#ifdef __myDEBUG + cout << "Bad simplex of _simplexTestEdges (" + << " "<< edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID() + << " "<< edge->_simplices[j]._nPrev->GetID() + << " "<< edge->_simplices[j]._nNext->GetID() << " )" << endl; +#endif + return false; + } + } + // Check if the last segments of _LayerEdge intersects 2D elements; // checked elements are either temporary faces or faces on surfaces w/o the layers @@ -2875,19 +3113,19 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, if ( FF1[0].IsNull() || FF2[0].IsNull() ) continue; -// // get a new normal for edge1 + // get a new normal for edge1 bool ok; gp_Vec dir1 = edge1->_normal, dir2 = edge2->_normal; if ( edge1->_cosin < 0 ) dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ).Normalized(); if ( edge2->_cosin < 0 ) dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok ).Normalized(); - // gp_Vec dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ); -// gp_Vec dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok2 ); -// double wgt1 = ( edge1->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 ); -// double wgt2 = ( edge2->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 ); -// gp_Vec newNorm = wgt1 * dir1 + wgt2 * dir2; -// newNorm.Normalize(); + // gp_Vec dir1 = getFaceDir( FF1[0], E1, edge1->_nodes[0], helper, ok ); + // gp_Vec dir2 = getFaceDir( FF2[0], E2, edge2->_nodes[0], helper, ok2 ); + // double wgt1 = ( edge1->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 ); + // double wgt2 = ( edge2->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 ); + // gp_Vec newNorm = wgt1 * dir1 + wgt2 * dir2; + // newNorm.Normalize(); double wgt1 = ( edge1->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 ); double wgt2 = ( edge2->_cosin + 1 ) / ( edge1->_cosin + edge2->_cosin + 2 ); @@ -2901,7 +3139,7 @@ bool _ViscousBuilder::updateNormals( _SolidData& data, n1 = edge1->_2neibors->_edges[0]->_nodes[0]; n2 = edge1->_2neibors->_edges[1]->_nodes[0]; //if ( !findNeiborsOnEdge( edge1, n1, n2, data )) - //continue; + // continue; edge1->SetDataByNeighbors( n1, n2, helper ); gp_Vec dirInFace; if ( edge1->_cosin < 0 ) @@ -3437,10 +3675,14 @@ bool _ViscousBuilder::refine(_SolidData& data) Handle(Geom_Surface) surface; TopoDS_Edge geomEdge; TopoDS_Face geomFace; + TopoDS_Shape prevSWOL; TopLoc_Location loc; - double f,l, u/*, distXYZ[4]*/; + double f,l, u; gp_XY uv; bool isOnEdge; + TGeomID prevBaseId = -1; + TNode2Edge* n2eMap = 0; + TNode2Edge::iterator n2e; for ( size_t i = 0; i < data._edges.size(); ++i ) { @@ -3460,27 +3702,46 @@ bool _ViscousBuilder::refine(_SolidData& data) edge._nodes[1] = 0; edge._nodes.back() = tgtNode; } - if ( !edge._sWOL.IsNull() ) + // get data of a shrink shape + if ( !edge._sWOL.IsNull() && edge._sWOL != prevSWOL ) { isOnEdge = ( edge._sWOL.ShapeType() == TopAbs_EDGE ); - // restore position of the last node -// gp_Pnt p; if ( isOnEdge ) { geomEdge = TopoDS::Edge( edge._sWOL ); - curve = BRep_Tool::Curve( geomEdge, loc, f,l); -// double u = helper.GetNodeU( tgtNode ); -// p = curve->Value( u ); + curve = BRep_Tool::Curve( geomEdge, loc, f,l); } else { geomFace = TopoDS::Face( edge._sWOL ); - surface = BRep_Tool::Surface( geomFace, loc ); -// gp_XY uv = helper.GetNodeUV( tgtNode ); -// p = surface->Value( uv.X(), uv.Y() ); + surface = BRep_Tool::Surface( geomFace, loc ); + } + prevSWOL = edge._sWOL; + } + // restore shapePos of the last node by already treated _LayerEdge of another _SolidData + const TGeomID baseShapeId = edge._nodes[0]->getshapeId(); + if ( baseShapeId != prevBaseId ) + { + map< TGeomID, TNode2Edge* >::iterator s2ne = data._s2neMap.find( baseShapeId ); + n2eMap = ( s2ne == data._s2neMap.end() ) ? 0 : n2eMap = s2ne->second; + prevBaseId = baseShapeId; + } + if ( n2eMap && (( n2e = n2eMap->find( edge._nodes[0] )) != n2eMap->end() )) + { + _LayerEdge* foundEdge = n2e->second; + const gp_XYZ& foundPos = foundEdge->_pos.back(); + SMDS_PositionPtr lastPos = tgtNode->GetPosition(); + if ( isOnEdge ) + { + SMDS_EdgePosition* epos = static_cast( lastPos ); + epos->SetUParameter( foundPos.X() ); + } + else + { + SMDS_FacePosition* fpos = static_cast( lastPos ); + fpos->SetUParameter( foundPos.X() ); + fpos->SetVParameter( foundPos.Y() ); } -// p.Transform( loc ); -// const_cast< SMDS_MeshNode* >( tgtNode )->setXYZ( p.X(), p.Y(), p.Z() ); } // calculate height of the first layer double h0; @@ -3518,12 +3779,14 @@ bool _ViscousBuilder::refine(_SolidData& data) if ( isOnEdge ) { u = pos.X(); - pos = curve->Value( u ).Transformed(loc); + if ( !node ) + pos = curve->Value( u ).Transformed(loc); } else { uv.SetCoord( pos.X(), pos.Y() ); - pos = surface->Value( pos.X(), pos.Y() ).Transformed(loc); + if ( !node ) + pos = surface->Value( pos.X(), pos.Y() ).Transformed(loc); } } // create or update the node @@ -3551,11 +3814,18 @@ bool _ViscousBuilder::refine(_SolidData& data) { u = 0.5 * ( u + helper.GetNodeU( geomEdge, node )); pos = curve->Value( u ).Transformed(loc); + + SMDS_EdgePosition* epos = static_cast( node->GetPosition() ); + epos->SetUParameter( u ); } else { uv = 0.5 * ( uv + helper.GetNodeUV( geomFace, node )); pos = surface->Value( uv.X(), uv.Y()).Transformed(loc); + + SMDS_FacePosition* fpos = static_cast( node->GetPosition() ); + fpos->SetUParameter( uv.X() ); + fpos->SetVParameter( uv.Y() ); } } node->setXYZ( pos.X(), pos.Y(), pos.Z() );