+ const double u0 = c3dAdaptor.FirstParameter();
+ gp_Pnt p; gp_Vec tangent;
+ _offPoints.resize( discret.NbPoints() );
+ for ( size_t i = 0; i < _offPoints.size(); i++ )
+ {
+ double u = discret.Parameter( i+1 );
+ c3dAdaptor.D1( u, p, tangent );
+ _offPoints[i]._xyz = p.XYZ();
+ _offPoints[i]._edgeDir = tangent.XYZ();
+ _offPoints[i]._param = GCPnts_AbscissaPoint::Length( c3dAdaptor, u0, u ) / _curveLen;
+ }
+
+ _LayerEdge* leOnV[2] = { getLEdgeOnV(0), getLEdgeOnV(1) };
+
+ // set _2edges
+ _offPoints [0]._2edges.set( &_leOnV[0], &_leOnV[0], 0.5, 0.5 );
+ _offPoints.back()._2edges.set( &_leOnV[1], &_leOnV[1], 0.5, 0.5 );
+ _2NearEdges tmp2edges;
+ tmp2edges._edges[1] = _eos._edges[0];
+ _leOnV[0]._2neibors = & tmp2edges;
+ _leOnV[0]._nodes = leOnV[0]->_nodes;
+ _leOnV[1]._nodes = leOnV[1]->_nodes;
+ _LayerEdge* eNext, *ePrev = & _leOnV[0];
+ for ( size_t iLE = 0, i = 1; i < _offPoints.size()-1; i++ )
+ {
+ // find _LayerEdge's located before and after an offset point
+ // (_eos._edges[ iLE ] is next after ePrev)
+ while ( iLE < _eos._edges.size() && _offPoints[i]._param > _leParams[ iLE ] )
+ ePrev = _eos._edges[ iLE++ ];
+ eNext = ePrev->_2neibors->_edges[1];
+
+ gp_Pnt p0 = SMESH_TNodeXYZ( ePrev->_nodes[0] );
+ gp_Pnt p1 = SMESH_TNodeXYZ( eNext->_nodes[0] );
+ double r = p0.Distance( _offPoints[i]._xyz ) / p0.Distance( p1 );
+ _offPoints[i]._2edges.set( ePrev, eNext, 1-r, r );
+ }
+
+ // replace _LayerEdge's on VERTEX by _leOnV in _offPoints._2edges
+ for ( size_t i = 0; i < _offPoints.size(); i++ )
+ if ( _offPoints[i]._2edges._edges[0] == leOnV[0] )
+ _offPoints[i]._2edges._edges[0] = & _leOnV[0];
+ else break;
+ for ( size_t i = _offPoints.size()-1; i > 0; i-- )
+ if ( _offPoints[i]._2edges._edges[1] == leOnV[1] )
+ _offPoints[i]._2edges._edges[1] = & _leOnV[1];
+ else break;
+
+ // set _normal of _leOnV[0] and _leOnV[1] to be normal to the EDGE
+
+ int iLBO = _offPoints.size() - 2; // last but one
+
+ _edgeDir[0] = getEdgeDir( E, leOnV[0]->_nodes[0], data.GetHelper() );
+ _edgeDir[1] = getEdgeDir( E, leOnV[1]->_nodes[0], data.GetHelper() );
+
+ _leOnV[ 0 ]._normal = getNormalNormal( leOnV[0]->_normal, _edgeDir[0] );
+ _leOnV[ 1 ]._normal = getNormalNormal( leOnV[1]->_normal, _edgeDir[1] );
+ _leOnV[ 0 ]._len = 0;
+ _leOnV[ 1 ]._len = 0;
+ _leOnV[ 0 ]._lenFactor = _offPoints[1 ]._2edges._edges[1]->_lenFactor;
+ _leOnV[ 1 ]._lenFactor = _offPoints[iLBO]._2edges._edges[0]->_lenFactor;
+
+ _iSeg[0] = 0;
+ _iSeg[1] = _offPoints.size()-2;
+
+ // initialize OffPnt::_len
+ for ( size_t i = 0; i < _offPoints.size(); ++i )
+ _offPoints[i]._len = 0;
+
+ if ( _eos._edges[0]->NbSteps() > 1 ) // already inflated several times, init _xyz
+ {
+ _leOnV[0]._len = leOnV[0]->_len;
+ _leOnV[1]._len = leOnV[1]->_len;
+ for ( size_t i = 0; i < _offPoints.size(); i++ )
+ {
+ _LayerEdge* e0 = _offPoints[i]._2edges._edges[0];
+ _LayerEdge* e1 = _offPoints[i]._2edges._edges[1];
+ const double w0 = _offPoints[i]._2edges._wgt[0];
+ const double w1 = _offPoints[i]._2edges._wgt[1];
+ double avgLen = ( e0->_len * w0 + e1->_len * w1 );
+ gp_XYZ avgXYZ = ( SMESH_TNodeXYZ( e0->_nodes.back() ) * w0 +
+ SMESH_TNodeXYZ( e1->_nodes.back() ) * w1 );
+ _offPoints[i]._xyz = avgXYZ;
+ _offPoints[i]._len = avgLen;
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief set _normal of _leOnV[is2nd] to be normal to the EDGE
+ */
+//================================================================================
+
+gp_XYZ _Smoother1D::getNormalNormal( const gp_XYZ & normal,
+ const gp_XYZ& edgeDir)
+{
+ gp_XYZ cross = normal ^ edgeDir;
+ gp_XYZ norm = edgeDir ^ cross;
+ double size = norm.Modulus();
+
+ return norm / size;
+}
+
+//================================================================================
+/*!
+ * \brief Sort _LayerEdge's by a parameter on a given EDGE
+ */
+//================================================================================
+
+void _SolidData::SortOnEdge( const TopoDS_Edge& E,
+ vector< _LayerEdge* >& edges)
+{
+ map< double, _LayerEdge* > u2edge;
+ for ( size_t i = 0; i < edges.size(); ++i )
+ u2edge.insert( u2edge.end(),
+ make_pair( _helper->GetNodeU( E, edges[i]->_nodes[0] ), edges[i] ));
+
+ ASSERT( u2edge.size() == edges.size() );
+ map< double, _LayerEdge* >::iterator u2e = u2edge.begin();
+ for ( size_t i = 0; i < edges.size(); ++i, ++u2e )
+ edges[i] = u2e->second;
+
+ Sort2NeiborsOnEdge( edges );
+}
+
+//================================================================================
+/*!
+ * \brief Set _2neibors according to the order of _LayerEdge on EDGE
+ */
+//================================================================================
+
+void _SolidData::Sort2NeiborsOnEdge( vector< _LayerEdge* >& edges )
+{
+ if ( edges.size() < 2 || !edges[0]->_2neibors ) return;
+
+ for ( size_t i = 0; i < edges.size()-1; ++i )
+ if ( edges[i]->_2neibors->tgtNode(1) != edges[i+1]->_nodes.back() )
+ edges[i]->_2neibors->reverse();
+
+ const size_t iLast = edges.size() - 1;
+ if ( edges.size() > 1 &&
+ edges[iLast]->_2neibors->tgtNode(0) != edges[iLast-1]->_nodes.back() )
+ edges[iLast]->_2neibors->reverse();
+}
+
+//================================================================================
+/*!
+ * \brief Return _EdgesOnShape* corresponding to the shape
+ */
+//================================================================================
+
+_EdgesOnShape* _SolidData::GetShapeEdges(const TGeomID shapeID )
+{
+ if ( shapeID < (int)_edgesOnShape.size() &&
+ _edgesOnShape[ shapeID ]._shapeID == shapeID )
+ return _edgesOnShape[ shapeID ]._subMesh ? & _edgesOnShape[ shapeID ] : 0;
+
+ for ( size_t i = 0; i < _edgesOnShape.size(); ++i )
+ if ( _edgesOnShape[i]._shapeID == shapeID )
+ return _edgesOnShape[i]._subMesh ? & _edgesOnShape[i] : 0;
+
+ return 0;
+}
+
+//================================================================================
+/*!
+ * \brief Return _EdgesOnShape* corresponding to the shape
+ */
+//================================================================================
+
+_EdgesOnShape* _SolidData::GetShapeEdges(const TopoDS_Shape& shape )
+{
+ SMESHDS_Mesh* meshDS = _proxyMesh->GetMesh()->GetMeshDS();
+ return GetShapeEdges( meshDS->ShapeToIndex( shape ));
+}
+
+//================================================================================
+/*!
+ * \brief Prepare data of the _LayerEdge for smoothing on FACE
+ */
+//================================================================================
+
+void _SolidData::PrepareEdgesToSmoothOnFace( _EdgesOnShape* eos, bool substituteSrcNodes )
+{
+ SMESH_MesherHelper helper( *_proxyMesh->GetMesh() );
+
+ set< TGeomID > vertices;
+ TopoDS_Face F;
+ if ( eos->ShapeType() == TopAbs_FACE )
+ {
+ // check FACE concavity and get concave VERTEXes
+ F = TopoDS::Face( eos->_shape );
+ if ( isConcave( F, helper, &vertices ))
+ _concaveFaces.insert( eos->_shapeID );
+
+ // set eos._eosConcaVer
+ eos->_eosConcaVer.clear();
+ eos->_eosConcaVer.reserve( vertices.size() );
+ for ( set< TGeomID >::iterator v = vertices.begin(); v != vertices.end(); ++v )
+ {
+ _EdgesOnShape* eov = GetShapeEdges( *v );
+ if ( eov && eov->_edges.size() == 1 )
+ {
+ eos->_eosConcaVer.push_back( eov );
+ for ( size_t i = 0; i < eov->_edges[0]->_neibors.size(); ++i )
+ eov->_edges[0]->_neibors[i]->Set( _LayerEdge::DIFFICULT );
+ }
+ }
+
+ // SetSmooLen() to _LayerEdge's on FACE
+ for ( size_t i = 0; i < eos->_edges.size(); ++i )
+ {
+ eos->_edges[i]->SetSmooLen( Precision::Infinite() );
+ }
+ SMESH_subMeshIteratorPtr smIt = eos->_subMesh->getDependsOnIterator(/*includeSelf=*/false);
+ while ( smIt->more() ) // loop on sub-shapes of the FACE
+ {
+ _EdgesOnShape* eoe = GetShapeEdges( smIt->next()->GetId() );
+ if ( !eoe ) continue;
+
+ vector<_LayerEdge*>& eE = eoe->_edges;
+ for ( size_t iE = 0; iE < eE.size(); ++iE ) // loop on _LayerEdge's on EDGE or VERTEX
+ {
+ if ( eE[iE]->_cosin <= theMinSmoothCosin )
+ continue;
+
+ SMDS_ElemIteratorPtr segIt = eE[iE]->_nodes[0]->GetInverseElementIterator(SMDSAbs_Edge);
+ while ( segIt->more() )
+ {
+ const SMDS_MeshElement* seg = segIt->next();
+ if ( !eos->_subMesh->DependsOn( seg->getshapeId() ))
+ continue;
+ if ( seg->GetNode(0) != eE[iE]->_nodes[0] )
+ continue; // not to check a seg twice
+ for ( size_t iN = 0; iN < eE[iE]->_neibors.size(); ++iN )
+ {
+ _LayerEdge* eN = eE[iE]->_neibors[iN];
+ if ( eN->_nodes[0]->getshapeId() != eos->_shapeID )
+ continue;
+ double dist = SMESH_MeshAlgos::GetDistance( seg, SMESH_TNodeXYZ( eN->_nodes[0] ));
+ double smooLen = getSmoothingThickness( eE[iE]->_cosin, dist );
+ eN->SetSmooLen( Min( smooLen, eN->GetSmooLen() ));
+ eN->Set( _LayerEdge::NEAR_BOUNDARY );
+ }
+ }
+ }
+ }
+ } // if ( eos->ShapeType() == TopAbs_FACE )
+
+ for ( size_t i = 0; i < eos->_edges.size(); ++i )
+ {
+ eos->_edges[i]->_smooFunction = 0;
+ eos->_edges[i]->Set( _LayerEdge::TO_SMOOTH );
+ }
+ bool isCurved = false;
+ for ( size_t i = 0; i < eos->_edges.size(); ++i )
+ {
+ _LayerEdge* edge = eos->_edges[i];
+
+ // get simplices sorted
+ _Simplex::SortSimplices( edge->_simplices );
+
+ // smoothing function
+ edge->ChooseSmooFunction( vertices, _n2eMap );
+
+ // set _curvature
+ double avgNormProj = 0, avgLen = 0;
+ for ( size_t iS = 0; iS < edge->_simplices.size(); ++iS )
+ {
+ _Simplex& s = edge->_simplices[iS];
+
+ gp_XYZ vec = edge->_pos.back() - SMESH_TNodeXYZ( s._nPrev );
+ avgNormProj += edge->_normal * vec;
+ avgLen += vec.Modulus();
+ if ( substituteSrcNodes )
+ {
+ s._nNext = _n2eMap[ s._nNext ]->_nodes.back();
+ s._nPrev = _n2eMap[ s._nPrev ]->_nodes.back();
+ }
+ }
+ avgNormProj /= edge->_simplices.size();
+ avgLen /= edge->_simplices.size();
+ if (( edge->_curvature = _Curvature::New( avgNormProj, avgLen )))
+ {
+ isCurved = true;
+ SMDS_FacePosition* fPos = dynamic_cast<SMDS_FacePosition*>( edge->_nodes[0]->GetPosition() );
+ if ( !fPos )
+ for ( size_t iS = 0; iS < edge->_simplices.size() && !fPos; ++iS )
+ fPos = dynamic_cast<SMDS_FacePosition*>( edge->_simplices[iS]._nPrev->GetPosition() );
+ if ( fPos )
+ edge->_curvature->_uv.SetCoord( fPos->GetUParameter(), fPos->GetVParameter() );
+ }
+ }
+
+ // prepare for putOnOffsetSurface()
+ if (( eos->ShapeType() == TopAbs_FACE ) &&
+ ( isCurved || !eos->_eosConcaVer.empty() ))
+ {
+ eos->_offsetSurf = helper.GetSurface( TopoDS::Face( eos->_shape ));
+ eos->_edgeForOffset = 0;
+
+ double maxCosin = -1;
+ for ( TopExp_Explorer eExp( eos->_shape, TopAbs_EDGE ); eExp.More(); eExp.Next() )
+ {
+ _EdgesOnShape* eoe = GetShapeEdges( eExp.Current() );
+ if ( !eoe || eoe->_edges.empty() ) continue;
+
+ vector<_LayerEdge*>& eE = eoe->_edges;
+ _LayerEdge* e = eE[ eE.size() / 2 ];
+ if ( e->_cosin > maxCosin )
+ {
+ eos->_edgeForOffset = e;
+ maxCosin = e->_cosin;
+ }
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Add faces for smoothing
+ */
+//================================================================================
+
+void _SolidData::AddShapesToSmooth( const set< _EdgesOnShape* >& eosToSmooth,
+ const set< _EdgesOnShape* >* edgesNoAnaSmooth )
+{
+ set< _EdgesOnShape * >::const_iterator eos = eosToSmooth.begin();
+ for ( ; eos != eosToSmooth.end(); ++eos )
+ {
+ if ( !*eos || (*eos)->_toSmooth ) continue;
+
+ (*eos)->_toSmooth = true;
+
+ if ( (*eos)->ShapeType() == TopAbs_FACE )
+ {
+ PrepareEdgesToSmoothOnFace( *eos, /*substituteSrcNodes=*/false );
+ (*eos)->_toSmooth = true;
+ }
+ }
+
+ // avoid _Smoother1D::smoothAnalyticEdge() of edgesNoAnaSmooth
+ if ( edgesNoAnaSmooth )
+ for ( eos = edgesNoAnaSmooth->begin(); eos != edgesNoAnaSmooth->end(); ++eos )
+ {
+ if ( (*eos)->_edgeSmoother )
+ (*eos)->_edgeSmoother->_anaCurve.Nullify();
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Limit _LayerEdge::_maxLen according to local curvature
+ */
+//================================================================================
+
+void _ViscousBuilder::limitMaxLenByCurvature( _SolidData& data, SMESH_MesherHelper& helper )
+{
+ // find intersection of neighbor _LayerEdge's to limit _maxLen
+ // according to local curvature (IPAL52648)
+
+ // This method must be called after findCollisionEdges() where _LayerEdge's
+ // get _lenFactor initialized in the case of eos._hyp.IsOffsetMethod()
+
+ for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
+ {
+ _EdgesOnShape& eosI = data._edgesOnShape[iS];
+ if ( eosI._edges.empty() ) continue;
+ if ( !eosI._hyp.ToSmooth() )
+ {
+ for ( size_t i = 0; i < eosI._edges.size(); ++i )
+ {
+ _LayerEdge* eI = eosI._edges[i];
+ for ( size_t iN = 0; iN < eI->_neibors.size(); ++iN )
+ {
+ _LayerEdge* eN = eI->_neibors[iN];
+ if ( eI->_nodes[0]->GetID() < eN->_nodes[0]->GetID() ) // treat this pair once
+ {
+ _EdgesOnShape* eosN = data.GetShapeEdges( eN );
+ limitMaxLenByCurvature( eI, eN, eosI, *eosN, helper );
+ }
+ }
+ }
+ }
+ else if ( eosI.ShapeType() == TopAbs_EDGE )
+ {
+ const TopoDS_Edge& E = TopoDS::Edge( eosI._shape );
+ if ( SMESH_Algo::IsStraight( E, /*degenResult=*/true )) continue;
+
+ _LayerEdge* e0 = eosI._edges[0];
+ for ( size_t i = 1; i < eosI._edges.size(); ++i )
+ {
+ _LayerEdge* eI = eosI._edges[i];
+ limitMaxLenByCurvature( eI, e0, eosI, eosI, helper );
+ e0 = eI;
+ }
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Limit _LayerEdge::_maxLen according to local curvature
+ */
+//================================================================================
+
+void _ViscousBuilder::limitMaxLenByCurvature( _LayerEdge* e1,
+ _LayerEdge* e2,
+ _EdgesOnShape& eos1,
+ _EdgesOnShape& eos2,
+ SMESH_MesherHelper& helper )
+{
+ gp_XYZ plnNorm = e1->_normal ^ e2->_normal;
+ double norSize = plnNorm.SquareModulus();
+ if ( norSize < std::numeric_limits<double>::min() )
+ return; // parallel normals
+
+ // find closest points of skew _LayerEdge's
+ SMESH_TNodeXYZ src1( e1->_nodes[0] ), src2( e2->_nodes[0] );
+ gp_XYZ dir12 = src2 - src1;
+ gp_XYZ perp1 = e1->_normal ^ plnNorm;
+ gp_XYZ perp2 = e2->_normal ^ plnNorm;
+ double dot1 = perp2 * e1->_normal;
+ double dot2 = perp1 * e2->_normal;
+ double u1 = ( perp2 * dir12 ) / dot1;
+ double u2 = - ( perp1 * dir12 ) / dot2;
+ if ( u1 > 0 && u2 > 0 )
+ {
+ double ovl = ( u1 * e1->_normal * dir12 -
+ u2 * e2->_normal * dir12 ) / dir12.SquareModulus();
+ if ( ovl > theSmoothThickToElemSizeRatio )
+ {
+ e1->_maxLen = Min( e1->_maxLen, 0.75 * u1 / e1->_lenFactor );
+ e2->_maxLen = Min( e2->_maxLen, 0.75 * u2 / e2->_lenFactor );
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Fill data._collisionEdges
+ */
+//================================================================================
+
+void _ViscousBuilder::findCollisionEdges( _SolidData& data, SMESH_MesherHelper& helper )
+{
+ data._collisionEdges.clear();
+
+ // set the full thickness of the layers to LEs
+ for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
+ {
+ _EdgesOnShape& eos = data._edgesOnShape[iS];
+ if ( eos._edges.empty() ) continue;
+ if ( eos.ShapeType() != TopAbs_EDGE && eos.ShapeType() != TopAbs_VERTEX ) continue;
+
+ for ( size_t i = 0; i < eos._edges.size(); ++i )
+ {
+ if ( eos._edges[i]->Is( _LayerEdge::BLOCKED )) continue;
+ double maxLen = eos._edges[i]->_maxLen;
+ eos._edges[i]->_maxLen = Precision::Infinite(); // avoid blocking
+ eos._edges[i]->SetNewLength( 1.5 * maxLen, eos, helper );
+ eos._edges[i]->_maxLen = maxLen;
+ }
+ }
+
+ // make temporary quadrangles got by extrusion of
+ // mesh edges along _LayerEdge._normal's
+
+ vector< const SMDS_MeshElement* > tmpFaces;
+
+ for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
+ {
+ _EdgesOnShape& eos = data._edgesOnShape[ iS ];
+ if ( eos.ShapeType() != TopAbs_EDGE )
+ continue;
+ if ( eos._edges.empty() )
+ {
+ _LayerEdge* edge[2] = { 0, 0 }; // LE of 2 VERTEX'es
+ SMESH_subMeshIteratorPtr smIt = eos._subMesh->getDependsOnIterator(/*includeSelf=*/false);
+ while ( smIt->more() )
+ if ( _EdgesOnShape* eov = data.GetShapeEdges( smIt->next()->GetId() ))
+ if ( eov->_edges.size() == 1 )
+ edge[ bool( edge[0]) ] = eov->_edges[0];
+
+ if ( edge[1] )
+ {
+ _TmpMeshFaceOnEdge* f = new _TmpMeshFaceOnEdge( edge[0], edge[1], --_tmpFaceID );
+ tmpFaces.push_back( f );
+ }
+ }
+ for ( size_t i = 0; i < eos._edges.size(); ++i )
+ {
+ _LayerEdge* edge = eos._edges[i];
+ for ( int j = 0; j < 2; ++j ) // loop on _2NearEdges
+ {
+ const SMDS_MeshNode* src2 = edge->_2neibors->srcNode(j);
+ if ( src2->GetPosition()->GetDim() > 0 &&
+ src2->GetID() < edge->_nodes[0]->GetID() )
+ continue; // avoid using same segment twice
+
+ // a _LayerEdge containg tgt2
+ _LayerEdge* neiborEdge = edge->_2neibors->_edges[j];
+
+ _TmpMeshFaceOnEdge* f = new _TmpMeshFaceOnEdge( edge, neiborEdge, --_tmpFaceID );
+ tmpFaces.push_back( f );
+ }
+ }
+ }
+
+ // Find _LayerEdge's intersecting tmpFaces.
+
+ SMDS_ElemIteratorPtr fIt( new SMDS_ElementVectorIterator( tmpFaces.begin(),
+ tmpFaces.end()));
+ SMESHUtils::Deleter<SMESH_ElementSearcher> searcher
+ ( SMESH_MeshAlgos::GetElementSearcher( *getMeshDS(), fIt ));
+
+ double dist1, dist2, segLen, eps = 0.5;
+ _CollisionEdges collEdges;
+ vector< const SMDS_MeshElement* > suspectFaces;
+ const double angle45 = Cos( 45. * M_PI / 180. );
+
+ for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
+ {
+ _EdgesOnShape& eos = data._edgesOnShape[ iS ];
+ if ( eos.ShapeType() == TopAbs_FACE || !eos._sWOL.IsNull() )
+ continue;
+ // find sub-shapes whose VL can influence VL on eos
+ set< TGeomID > neighborShapes;
+ PShapeIteratorPtr fIt = helper.GetAncestors( eos._shape, *_mesh, TopAbs_FACE );
+ while ( const TopoDS_Shape* face = fIt->next() )
+ {
+ TGeomID faceID = getMeshDS()->ShapeToIndex( *face );
+ if ( _EdgesOnShape* eof = data.GetShapeEdges( faceID ))
+ {
+ SMESH_subMeshIteratorPtr subIt = eof->_subMesh->getDependsOnIterator(/*includeSelf=*/false);
+ while ( subIt->more() )
+ neighborShapes.insert( subIt->next()->GetId() );
+ }
+ }
+ if ( eos.ShapeType() == TopAbs_VERTEX )
+ {
+ PShapeIteratorPtr eIt = helper.GetAncestors( eos._shape, *_mesh, TopAbs_EDGE );
+ while ( const TopoDS_Shape* edge = eIt->next() )
+ neighborShapes.erase( getMeshDS()->ShapeToIndex( *edge ));
+ }
+ // find intersecting _LayerEdge's
+ for ( size_t i = 0; i < eos._edges.size(); ++i )
+ {
+ if ( eos._edges[i]->Is( _LayerEdge::MULTI_NORMAL )) continue;
+ _LayerEdge* edge = eos._edges[i];
+ gp_Ax1 lastSegment = edge->LastSegment( segLen, eos );
+ segLen *= 1.2;
+
+ gp_Vec eSegDir0, eSegDir1;
+ if ( edge->IsOnEdge() )
+ {
+ SMESH_TNodeXYZ eP( edge->_nodes[0] );
+ eSegDir0 = SMESH_TNodeXYZ( edge->_2neibors->srcNode(0) ) - eP;
+ eSegDir1 = SMESH_TNodeXYZ( edge->_2neibors->srcNode(1) ) - eP;
+ }
+ suspectFaces.clear();
+ searcher->GetElementsInSphere( SMESH_TNodeXYZ( edge->_nodes.back()), edge->_len * 2,
+ SMDSAbs_Face, suspectFaces );
+ collEdges._intEdges.clear();
+ for ( size_t j = 0 ; j < suspectFaces.size(); ++j )
+ {
+ const _TmpMeshFaceOnEdge* f = (const _TmpMeshFaceOnEdge*) suspectFaces[j];
+ if ( f->_le1 == edge || f->_le2 == edge ) continue;
+ if ( !neighborShapes.count( f->_le1->_nodes[0]->getshapeId() )) continue;
+ if ( !neighborShapes.count( f->_le2->_nodes[0]->getshapeId() )) continue;
+ if ( edge->IsOnEdge() ) {
+ if ( edge->_2neibors->include( f->_le1 ) ||
+ edge->_2neibors->include( f->_le2 )) continue;
+ }
+ else {
+ if (( f->_le1->IsOnEdge() && f->_le1->_2neibors->include( edge )) ||
+ ( f->_le2->IsOnEdge() && f->_le2->_2neibors->include( edge ))) continue;
+ }
+ dist1 = dist2 = Precision::Infinite();
+ if ( !edge->SegTriaInter( lastSegment, f->_nn[0], f->_nn[1], f->_nn[2], dist1, eps ))
+ dist1 = Precision::Infinite();
+ if ( !edge->SegTriaInter( lastSegment, f->_nn[3], f->_nn[2], f->_nn[0], dist2, eps ))
+ dist2 = Precision::Infinite();
+ if (( dist1 > segLen ) && ( dist2 > segLen ))
+ continue;
+
+ if ( edge->IsOnEdge() )
+ {
+ // skip perpendicular EDGEs
+ gp_Vec fSegDir = SMESH_TNodeXYZ( f->_nn[0] ) - SMESH_TNodeXYZ( f->_nn[3] );
+ bool isParallel = ( isLessAngle( eSegDir0, fSegDir, angle45 ) ||
+ isLessAngle( eSegDir1, fSegDir, angle45 ) ||
+ isLessAngle( eSegDir0, fSegDir.Reversed(), angle45 ) ||
+ isLessAngle( eSegDir1, fSegDir.Reversed(), angle45 ));
+ if ( !isParallel )
+ continue;
+ }
+
+ // either limit inflation of edges or remember them for updating _normal
+ // double dot = edge->_normal * f->GetDir();
+ // if ( dot > 0.1 )
+ {
+ collEdges._intEdges.push_back( f->_le1 );
+ collEdges._intEdges.push_back( f->_le2 );
+ }
+ // else
+ // {
+ // double shortLen = 0.75 * ( Min( dist1, dist2 ) / edge->_lenFactor );
+ // edge->_maxLen = Min( shortLen, edge->_maxLen );
+ // }
+ }
+
+ if ( !collEdges._intEdges.empty() )
+ {
+ collEdges._edge = edge;
+ data._collisionEdges.push_back( collEdges );
+ }
+ }
+ }
+
+ for ( size_t i = 0 ; i < tmpFaces.size(); ++i )
+ delete tmpFaces[i];
+
+ // restore the zero thickness
+ for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
+ {
+ _EdgesOnShape& eos = data._edgesOnShape[iS];
+ if ( eos._edges.empty() ) continue;
+ if ( eos.ShapeType() != TopAbs_EDGE && eos.ShapeType() != TopAbs_VERTEX ) continue;
+
+ for ( size_t i = 0; i < eos._edges.size(); ++i )
+ {
+ eos._edges[i]->InvalidateStep( 1, eos );
+ eos._edges[i]->_len = 0;
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Modify normals of _LayerEdge's on EDGE's to avoid intersection with
+ * _LayerEdge's on neighbor EDGE's
+ */
+//================================================================================
+
+bool _ViscousBuilder::updateNormals( _SolidData& data,
+ SMESH_MesherHelper& helper,
+ int stepNb,
+ double stepSize)
+{
+ updateNormalsOfC1Vertices( data );
+
+ if ( stepNb > 0 && !updateNormalsOfConvexFaces( data, helper, stepNb ))
+ return false;
+
+ // map to store new _normal and _cosin for each intersected edge
+ map< _LayerEdge*, _LayerEdge, _LayerEdgeCmp > edge2newEdge;
+ map< _LayerEdge*, _LayerEdge, _LayerEdgeCmp >::iterator e2neIt;
+ _LayerEdge zeroEdge;
+ zeroEdge._normal.SetCoord( 0,0,0 );
+ zeroEdge._maxLen = Precision::Infinite();
+ zeroEdge._nodes.resize(1); // to init _TmpMeshFaceOnEdge
+
+ set< _EdgesOnShape* > shapesToSmooth, edgesNoAnaSmooth;
+
+ double segLen, dist1, dist2, dist;
+ vector< pair< _LayerEdge*, double > > intEdgesDist;
+ _TmpMeshFaceOnEdge quad( &zeroEdge, &zeroEdge, 0 );
+
+ for ( int iter = 0; iter < 5; ++iter )
+ {
+ edge2newEdge.clear();
+
+ for ( size_t iE = 0; iE < data._collisionEdges.size(); ++iE )
+ {
+ _CollisionEdges& ce = data._collisionEdges[iE];
+ _LayerEdge* edge1 = ce._edge;
+ if ( !edge1 /*|| edge1->Is( _LayerEdge::BLOCKED )*/) continue;
+ _EdgesOnShape* eos1 = data.GetShapeEdges( edge1 );
+ if ( !eos1 ) continue;
+
+ // detect intersections
+ gp_Ax1 lastSeg = edge1->LastSegment( segLen, *eos1 );
+ double testLen = 1.5 * edge1->_maxLen * edge1->_lenFactor;
+ double eps = 0.5;
+ intEdgesDist.clear();
+ double minIntDist = Precision::Infinite();
+ for ( size_t i = 0; i < ce._intEdges.size(); i += 2 )
+ {
+ if ( edge1->Is( _LayerEdge::BLOCKED ) &&
+ ce._intEdges[i ]->Is( _LayerEdge::BLOCKED ) &&
+ ce._intEdges[i+1]->Is( _LayerEdge::BLOCKED ))
+ continue;
+ double dot = edge1->_normal * quad.GetDir( ce._intEdges[i], ce._intEdges[i+1] );
+ double fact = ( 1.1 + dot * dot );
+ SMESH_TNodeXYZ pSrc0( ce.nSrc(i) ), pSrc1( ce.nSrc(i+1) );
+ SMESH_TNodeXYZ pTgt0( ce.nTgt(i) ), pTgt1( ce.nTgt(i+1) );
+ gp_XYZ pLast0 = pSrc0 + ( pTgt0 - pSrc0 ) * fact;
+ gp_XYZ pLast1 = pSrc1 + ( pTgt1 - pSrc1 ) * fact;
+ dist1 = dist2 = Precision::Infinite();
+ if ( !edge1->SegTriaInter( lastSeg, pSrc0, pLast0, pSrc1, dist1, eps ) &&
+ !edge1->SegTriaInter( lastSeg, pSrc1, pLast1, pLast0, dist2, eps ))
+ continue;
+ dist = dist1;
+ if ( dist > testLen || dist <= 0 )
+ {
+ dist = dist2;
+ if ( dist > testLen || dist <= 0 )
+ continue;
+ }
+ // choose a closest edge
+ gp_Pnt intP( lastSeg.Location().XYZ() + lastSeg.Direction().XYZ() * ( dist + segLen ));
+ double d1 = intP.SquareDistance( pSrc0 );
+ double d2 = intP.SquareDistance( pSrc1 );
+ int iClose = i + ( d2 < d1 );
+ _LayerEdge* edge2 = ce._intEdges[iClose];
+ edge2->Unset( _LayerEdge::MARKED );
+
+ // choose a closest edge among neighbors
+ gp_Pnt srcP( SMESH_TNodeXYZ( edge1->_nodes[0] ));
+ d1 = srcP.SquareDistance( SMESH_TNodeXYZ( edge2->_nodes[0] ));
+ for ( size_t j = 0; j < intEdgesDist.size(); ++j )
+ {
+ _LayerEdge * edgeJ = intEdgesDist[j].first;
+ if ( edge2->IsNeiborOnEdge( edgeJ ))
+ {
+ d2 = srcP.SquareDistance( SMESH_TNodeXYZ( edgeJ->_nodes[0] ));
+ ( d1 < d2 ? edgeJ : edge2 )->Set( _LayerEdge::MARKED );
+ }
+ }
+ intEdgesDist.push_back( make_pair( edge2, dist ));
+ // if ( Abs( d2 - d1 ) / Max( d2, d1 ) < 0.5 )
+ // {
+ // iClose = i + !( d2 < d1 );
+ // intEdges.push_back( ce._intEdges[iClose] );
+ // ce._intEdges[iClose]->Unset( _LayerEdge::MARKED );
+ // }
+ minIntDist = Min( edge1->_len * edge1->_lenFactor - segLen + dist, minIntDist );
+ }
+
+ //ce._edge = 0;
+
+ // compute new _normals
+ for ( size_t i = 0; i < intEdgesDist.size(); ++i )
+ {
+ _LayerEdge* edge2 = intEdgesDist[i].first;
+ double distWgt = edge1->_len / intEdgesDist[i].second;
+ // if ( edge1->Is( _LayerEdge::BLOCKED ) &&
+ // edge2->Is( _LayerEdge::BLOCKED )) continue;
+ if ( edge2->Is( _LayerEdge::MARKED )) continue;
+ edge2->Set( _LayerEdge::MARKED );
+
+ // get a new normal
+ gp_XYZ dir1 = edge1->_normal, dir2 = edge2->_normal;
+
+ double cos1 = Abs( edge1->_cosin ), cos2 = Abs( edge2->_cosin );
+ double wgt1 = ( cos1 + 0.001 ) / ( cos1 + cos2 + 0.002 );
+ double wgt2 = ( cos2 + 0.001 ) / ( cos1 + cos2 + 0.002 );
+ // double cos1 = Abs( edge1->_cosin ), cos2 = Abs( edge2->_cosin );
+ // double sgn1 = 0.1 * ( 1 + edge1->_cosin ), sgn2 = 0.1 * ( 1 + edge2->_cosin );
+ // double wgt1 = ( cos1 + sgn1 ) / ( cos1 + cos2 + sgn1 + sgn2 );
+ // double wgt2 = ( cos2 + sgn2 ) / ( cos1 + cos2 + sgn1 + sgn2 );
+ gp_XYZ newNormal = wgt1 * dir1 + wgt2 * dir2;
+ newNormal.Normalize();
+
+ // get new cosin
+ double newCos;
+ double sgn1 = edge1->_cosin / cos1, sgn2 = edge2->_cosin / cos2;
+ if ( cos1 < theMinSmoothCosin )
+ {
+ newCos = cos2 * sgn1;
+ }
+ else if ( cos2 > theMinSmoothCosin ) // both cos1 and cos2 > theMinSmoothCosin
+ {
+ newCos = ( wgt1 * cos1 + wgt2 * cos2 ) * edge1->_cosin / cos1;
+ }
+ else
+ {
+ newCos = edge1->_cosin;
+ }
+
+ e2neIt = edge2newEdge.insert( make_pair( edge1, zeroEdge )).first;
+ e2neIt->second._normal += distWgt * newNormal;
+ e2neIt->second._cosin = newCos;
+ e2neIt->second._maxLen = 0.7 * minIntDist / edge1->_lenFactor;
+ if ( iter > 0 && sgn1 * sgn2 < 0 && edge1->_cosin < 0 )
+ e2neIt->second._normal += dir2;
+ e2neIt = edge2newEdge.insert( make_pair( edge2, zeroEdge )).first;
+ e2neIt->second._normal += distWgt * newNormal;
+ e2neIt->second._cosin = edge2->_cosin;
+ if ( iter > 0 && sgn1 * sgn2 < 0 && edge2->_cosin < 0 )
+ e2neIt->second._normal += dir1;
+ }
+ }
+
+ if ( edge2newEdge.empty() )
+ break; //return true;
+
+ dumpFunction(SMESH_Comment("updateNormals")<< data._index << "_" << stepNb << "_it" << iter);
+
+ // Update data of edges depending on a new _normal
+
+ data.UnmarkEdges();
+ for ( e2neIt = edge2newEdge.begin(); e2neIt != edge2newEdge.end(); ++e2neIt )
+ {
+ _LayerEdge* edge = e2neIt->first;
+ if ( edge->Is( _LayerEdge::BLOCKED )) continue;
+ _LayerEdge& newEdge = e2neIt->second;
+ _EdgesOnShape* eos = data.GetShapeEdges( edge );
+
+ // Check if a new _normal is OK:
+ newEdge._normal.Normalize();
+ if ( !isNewNormalOk( data, *edge, newEdge._normal ))
+ {
+ if ( newEdge._maxLen < edge->_len && iter > 0 ) // limit _maxLen
+ {
+ edge->InvalidateStep( stepNb + 1, *eos, /*restoreLength=*/true );
+ edge->_maxLen = newEdge._maxLen;
+ edge->SetNewLength( newEdge._maxLen, *eos, helper );
+ }
+ continue; // the new _normal is bad
+ }
+ // the new _normal is OK
+
+ // find shapes that need smoothing due to change of _normal
+ if ( edge->_cosin < theMinSmoothCosin &&
+ newEdge._cosin > theMinSmoothCosin )
+ {
+ if ( eos->_sWOL.IsNull() )
+ {
+ SMDS_ElemIteratorPtr fIt = edge->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() )
+ shapesToSmooth.insert( data.GetShapeEdges( fIt->next()->getshapeId() ));
+ }
+ else // edge inflates along a FACE
+ {
+ TopoDS_Shape V = helper.GetSubShapeByNode( edge->_nodes[0], getMeshDS() );
+ PShapeIteratorPtr eIt = helper.GetAncestors( V, *_mesh, TopAbs_EDGE );
+ while ( const TopoDS_Shape* E = eIt->next() )
+ {
+ if ( !helper.IsSubShape( *E, /*FACE=*/eos->_sWOL ))
+ continue;
+ gp_Vec edgeDir = getEdgeDir( TopoDS::Edge( *E ), TopoDS::Vertex( V ));
+ double angle = edgeDir.Angle( newEdge._normal ); // [0,PI]
+ if ( angle < M_PI / 2 )
+ shapesToSmooth.insert( data.GetShapeEdges( *E ));
+ }
+ }
+ }
+
+ double len = edge->_len;
+ edge->InvalidateStep( stepNb + 1, *eos, /*restoreLength=*/true );
+ edge->SetNormal( newEdge._normal );
+ edge->SetCosin( newEdge._cosin );
+ edge->SetNewLength( len, *eos, helper );
+ edge->Set( _LayerEdge::MARKED );
+ edge->Set( _LayerEdge::NORMAL_UPDATED );
+ edgesNoAnaSmooth.insert( eos );
+ }
+
+ // Update normals and other dependent data of not intersecting _LayerEdge's
+ // neighboring the intersecting ones
+
+ for ( e2neIt = edge2newEdge.begin(); e2neIt != edge2newEdge.end(); ++e2neIt )
+ {
+ _LayerEdge* edge1 = e2neIt->first;
+ _EdgesOnShape* eos1 = data.GetShapeEdges( edge1 );
+ if ( !edge1->Is( _LayerEdge::MARKED ))
+ continue;
+
+ if ( edge1->IsOnEdge() )
+ {
+ const SMDS_MeshNode * n1 = edge1->_2neibors->srcNode(0);
+ const SMDS_MeshNode * n2 = edge1->_2neibors->srcNode(1);
+ edge1->SetDataByNeighbors( n1, n2, *eos1, helper );
+ }
+
+ if ( !edge1->_2neibors || !eos1->_sWOL.IsNull() )
+ continue;
+ for ( int j = 0; j < 2; ++j ) // loop on 2 neighbors
+ {
+ _LayerEdge* neighbor = edge1->_2neibors->_edges[j];
+ if ( neighbor->Is( _LayerEdge::MARKED ) /*edge2newEdge.count ( neighbor )*/)
+ continue; // j-th neighbor is also intersected
+ _LayerEdge* prevEdge = edge1;
+ const int nbSteps = 10;
+ for ( int step = nbSteps; step; --step ) // step from edge1 in j-th direction
+ {
+ if ( neighbor->Is( _LayerEdge::BLOCKED ) ||
+ neighbor->Is( _LayerEdge::MARKED ))
+ break;
+ _EdgesOnShape* eos = data.GetShapeEdges( neighbor );
+ if ( !eos ) continue;
+ _LayerEdge* nextEdge = neighbor;
+ if ( neighbor->_2neibors )
+ {
+ int iNext = 0;
+ nextEdge = neighbor->_2neibors->_edges[iNext];
+ if ( nextEdge == prevEdge )
+ nextEdge = neighbor->_2neibors->_edges[ ++iNext ];
+ }
+ double r = double(step-1)/nbSteps/(iter+1);
+ if ( !nextEdge->_2neibors )
+ r = Min( r, 0.5 );
+
+ gp_XYZ newNorm = prevEdge->_normal * r + nextEdge->_normal * (1-r);
+ newNorm.Normalize();
+ if ( !isNewNormalOk( data, *neighbor, newNorm ))
+ break;
+
+ double len = neighbor->_len;
+ neighbor->InvalidateStep( stepNb + 1, *eos, /*restoreLength=*/true );
+ neighbor->SetNormal( newNorm );
+ neighbor->SetCosin( prevEdge->_cosin * r + nextEdge->_cosin * (1-r) );
+ if ( neighbor->_2neibors )
+ neighbor->SetDataByNeighbors( prevEdge->_nodes[0], nextEdge->_nodes[0], *eos, helper );
+ neighbor->SetNewLength( len, *eos, helper );
+ neighbor->Set( _LayerEdge::MARKED );
+ neighbor->Set( _LayerEdge::NORMAL_UPDATED );
+ edgesNoAnaSmooth.insert( eos );
+
+ if ( !neighbor->_2neibors )
+ break; // neighbor is on VERTEX
+
+ // goto the next neighbor
+ prevEdge = neighbor;
+ neighbor = nextEdge;
+ }
+ }
+ }
+ dumpFunctionEnd();
+ } // iterations
+
+ data.AddShapesToSmooth( shapesToSmooth, &edgesNoAnaSmooth );
+
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Check if a new normal is OK
+ */
+//================================================================================
+
+bool _ViscousBuilder::isNewNormalOk( _SolidData& data,
+ _LayerEdge& edge,
+ const gp_XYZ& newNormal)
+{
+ // check a min angle between the newNormal and surrounding faces
+ vector<_Simplex> simplices;
+ SMESH_TNodeXYZ n0( edge._nodes[0] ), n1, n2;
+ _Simplex::GetSimplices( n0._node, simplices, data._ignoreFaceIds, &data );
+ double newMinDot = 1, curMinDot = 1;
+ for ( size_t i = 0; i < simplices.size(); ++i )
+ {
+ n1.Set( simplices[i]._nPrev );
+ n2.Set( simplices[i]._nNext );
+ gp_XYZ normFace = ( n1 - n0 ) ^ ( n2 - n0 );
+ double normLen2 = normFace.SquareModulus();
+ if ( normLen2 < std::numeric_limits<double>::min() )
+ continue;
+ normFace /= Sqrt( normLen2 );
+ newMinDot = Min( newNormal * normFace, newMinDot );
+ curMinDot = Min( edge._normal * normFace, curMinDot );
+ }
+ bool ok = true;
+ if ( newMinDot < 0.5 )
+ {
+ ok = ( newMinDot >= curMinDot * 0.9 );
+ //return ( newMinDot >= ( curMinDot * ( 0.8 + 0.1 * edge.NbSteps() )));
+ // double initMinDot2 = 1. - edge._cosin * edge._cosin;
+ // return ( newMinDot * newMinDot ) >= ( 0.8 * initMinDot2 );
+ }
+
+ return ok;
+}
+
+//================================================================================
+/*!
+ * \brief Modify normals of _LayerEdge's on FACE to reflex smoothing
+ */
+//================================================================================
+
+bool _ViscousBuilder::updateNormalsOfSmoothed( _SolidData& data,
+ SMESH_MesherHelper& helper,
+ const int nbSteps,
+ const double stepSize )
+{
+ if ( data._nbShapesToSmooth == 0 || nbSteps == 0 )
+ return true; // no shapes needing smoothing
+
+ for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
+ {
+ _EdgesOnShape& eos = data._edgesOnShape[ iS ];
+ if ( //!eos._toSmooth || _eosC1 have _toSmooth == false
+ !eos._hyp.ToSmooth() ||
+ eos.ShapeType() != TopAbs_FACE ||
+ eos._edges.empty() )
+ continue;
+
+ bool toSmooth = ( eos._edges[ 0 ]->NbSteps() >= nbSteps+1 );
+ if ( !toSmooth ) continue;
+
+ for ( size_t i = 0; i < eos._edges.size(); ++i )
+ {
+ _LayerEdge* edge = eos._edges[i];
+ if ( !edge->Is( _LayerEdge::SMOOTHED ))
+ continue;
+ if ( edge->Is( _LayerEdge::DIFFICULT ) && nbSteps != 1 )
+ continue;
+
+ const gp_XYZ& pPrev = edge->PrevPos();
+ const gp_XYZ& pLast = edge->_pos.back();
+ gp_XYZ stepVec = pLast - pPrev;
+ double realStepSize = stepVec.Modulus();
+ if ( realStepSize < numeric_limits<double>::min() )
+ continue;
+
+ edge->_lenFactor = realStepSize / stepSize;
+ edge->_normal = stepVec / realStepSize;
+ edge->Set( _LayerEdge::NORMAL_UPDATED );
+ }
+ }
+
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Modify normals of _LayerEdge's on C1 VERTEXes
+ */
+//================================================================================
+
+void _ViscousBuilder::updateNormalsOfC1Vertices( _SolidData& data )
+{
+ for ( size_t iS = 0; iS < data._edgesOnShape.size(); ++iS )
+ {
+ _EdgesOnShape& eov = data._edgesOnShape[ iS ];
+ if ( eov._eosC1.empty() ||
+ eov.ShapeType() != TopAbs_VERTEX ||
+ eov._edges.empty() )
+ continue;
+
+ gp_XYZ newNorm = eov._edges[0]->_normal;
+ double curThick = eov._edges[0]->_len * eov._edges[0]->_lenFactor;
+ bool normChanged = false;
+
+ for ( size_t i = 0; i < eov._eosC1.size(); ++i )
+ {
+ _EdgesOnShape* eoe = eov._eosC1[i];
+ const TopoDS_Edge& e = TopoDS::Edge( eoe->_shape );
+ const double eLen = SMESH_Algo::EdgeLength( e );
+ TopoDS_Shape oppV = SMESH_MesherHelper::IthVertex( 0, e );
+ if ( oppV.IsSame( eov._shape ))
+ oppV = SMESH_MesherHelper::IthVertex( 1, e );
+ _EdgesOnShape* eovOpp = data.GetShapeEdges( oppV );
+ if ( !eovOpp || eovOpp->_edges.empty() ) continue;
+ if ( eov._edges[0]->Is( _LayerEdge::BLOCKED )) continue;
+
+ double curThickOpp = eovOpp->_edges[0]->_len * eovOpp->_edges[0]->_lenFactor;
+ if ( curThickOpp + curThick < eLen )
+ continue;
+
+ double wgt = 2. * curThick / eLen;
+ newNorm += wgt * eovOpp->_edges[0]->_normal;
+ normChanged = true;
+ }
+ if ( normChanged )
+ {
+ eov._edges[0]->SetNormal( newNorm.Normalized() );
+ eov._edges[0]->Set( _LayerEdge::NORMAL_UPDATED );
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Modify normals of _LayerEdge's on _ConvexFace's
+ */
+//================================================================================
+
+bool _ViscousBuilder::updateNormalsOfConvexFaces( _SolidData& data,
+ SMESH_MesherHelper& helper,
+ int stepNb )
+{
+ SMESHDS_Mesh* meshDS = helper.GetMeshDS();
+ bool isOK;
+
+ map< TGeomID, _ConvexFace >::iterator id2face = data._convexFaces.begin();
+ for ( ; id2face != data._convexFaces.end(); ++id2face )
+ {
+ _ConvexFace & convFace = (*id2face).second;
+ if ( convFace._normalsFixed )
+ continue; // already fixed
+ if ( convFace.CheckPrisms() )
+ continue; // nothing to fix
+
+ convFace._normalsFixed = true;
+
+ BRepAdaptor_Surface surface ( convFace._face, false );
+ BRepLProp_SLProps surfProp( surface, 2, 1e-6 );
+
+ // check if the convex FACE is of spherical shape
+
+ Bnd_B3d centersBox; // bbox of centers of curvature of _LayerEdge's on VERTEXes
+ Bnd_B3d nodesBox;
+ gp_Pnt center;
+
+ map< TGeomID, _EdgesOnShape* >::iterator id2eos = convFace._subIdToEOS.begin();
+ for ( ; id2eos != convFace._subIdToEOS.end(); ++id2eos )
+ {
+ _EdgesOnShape& eos = *(id2eos->second);
+ if ( eos.ShapeType() == TopAbs_VERTEX )
+ {
+ _LayerEdge* ledge = eos._edges[ 0 ];
+ if ( convFace.GetCenterOfCurvature( ledge, surfProp, helper, center ))
+ centersBox.Add( center );
+ }
+ for ( size_t i = 0; i < eos._edges.size(); ++i )
+ nodesBox.Add( SMESH_TNodeXYZ( eos._edges[ i ]->_nodes[0] ));
+ }
+ if ( centersBox.IsVoid() )
+ {
+ debugMsg( "Error: centersBox.IsVoid()" );
+ return false;
+ }
+ const bool isSpherical =
+ ( centersBox.SquareExtent() < 1e-6 * nodesBox.SquareExtent() );
+
+ int nbEdges = helper.Count( convFace._face, TopAbs_EDGE, /*ignoreSame=*/false );
+ vector < _CentralCurveOnEdge > centerCurves( nbEdges );
+
+ if ( isSpherical )
+ {
+ // set _LayerEdge::_normal as average of all normals
+
+ // WARNING: different density of nodes on EDGEs is not taken into account that
+ // can lead to an improper new normal
+
+ gp_XYZ avgNormal( 0,0,0 );
+ nbEdges = 0;
+ id2eos = convFace._subIdToEOS.begin();
+ for ( ; id2eos != convFace._subIdToEOS.end(); ++id2eos )
+ {
+ _EdgesOnShape& eos = *(id2eos->second);
+ // set data of _CentralCurveOnEdge
+ if ( eos.ShapeType() == TopAbs_EDGE )
+ {
+ _CentralCurveOnEdge& ceCurve = centerCurves[ nbEdges++ ];
+ ceCurve.SetShapes( TopoDS::Edge( eos._shape ), convFace, data, helper );
+ if ( !eos._sWOL.IsNull() )
+ ceCurve._adjFace.Nullify();
+ else
+ ceCurve._ledges.insert( ceCurve._ledges.end(),
+ eos._edges.begin(), eos._edges.end());
+ }
+ // summarize normals
+ for ( size_t i = 0; i < eos._edges.size(); ++i )
+ avgNormal += eos._edges[ i ]->_normal;
+ }
+ double normSize = avgNormal.SquareModulus();
+ if ( normSize < 1e-200 )
+ {
+ debugMsg( "updateNormalsOfConvexFaces(): zero avgNormal" );
+ return false;
+ }
+ avgNormal /= Sqrt( normSize );
+
+ // compute new _LayerEdge::_cosin on EDGEs
+ double avgCosin = 0;
+ int nbCosin = 0;
+ gp_Vec inFaceDir;
+ for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
+ {
+ _CentralCurveOnEdge& ceCurve = centerCurves[ iE ];
+ if ( ceCurve._adjFace.IsNull() )
+ continue;
+ for ( size_t iLE = 0; iLE < ceCurve._ledges.size(); ++iLE )
+ {
+ const SMDS_MeshNode* node = ceCurve._ledges[ iLE ]->_nodes[0];
+ inFaceDir = getFaceDir( ceCurve._adjFace, ceCurve._edge, node, helper, isOK );
+ if ( isOK )
+ {
+ double angle = inFaceDir.Angle( avgNormal ); // [0,PI]
+ ceCurve._ledges[ iLE ]->_cosin = Cos( angle );
+ avgCosin += ceCurve._ledges[ iLE ]->_cosin;
+ nbCosin++;
+ }
+ }
+ }
+ if ( nbCosin > 0 )
+ avgCosin /= nbCosin;
+
+ // set _LayerEdge::_normal = avgNormal
+ id2eos = convFace._subIdToEOS.begin();
+ for ( ; id2eos != convFace._subIdToEOS.end(); ++id2eos )
+ {
+ _EdgesOnShape& eos = *(id2eos->second);
+ if ( eos.ShapeType() != TopAbs_EDGE )
+ for ( size_t i = 0; i < eos._edges.size(); ++i )
+ eos._edges[ i ]->_cosin = avgCosin;
+
+ for ( size_t i = 0; i < eos._edges.size(); ++i )
+ {
+ eos._edges[ i ]->SetNormal( avgNormal );
+ eos._edges[ i ]->Set( _LayerEdge::NORMAL_UPDATED );
+ }
+ }
+ }
+ else // if ( isSpherical )
+ {
+ // We suppose that centers of curvature at all points of the FACE
+ // lie on some curve, let's call it "central curve". For all _LayerEdge's
+ // having a common center of curvature we define the same new normal
+ // as a sum of normals of _LayerEdge's on EDGEs among them.
+
+ // get all centers of curvature for each EDGE
+
+ helper.SetSubShape( convFace._face );
+ _LayerEdge* vertexLEdges[2], **edgeLEdge, **edgeLEdgeEnd;
+
+ TopExp_Explorer edgeExp( convFace._face, TopAbs_EDGE );
+ for ( int iE = 0; edgeExp.More(); edgeExp.Next(), ++iE )
+ {
+ const TopoDS_Edge& edge = TopoDS::Edge( edgeExp.Current() );
+
+ // set adjacent FACE
+ centerCurves[ iE ].SetShapes( edge, convFace, data, helper );
+
+ // get _LayerEdge's of the EDGE
+ TGeomID edgeID = meshDS->ShapeToIndex( edge );
+ _EdgesOnShape* eos = data.GetShapeEdges( edgeID );
+ if ( !eos || eos->_edges.empty() )
+ {
+ // no _LayerEdge's on EDGE, use _LayerEdge's on VERTEXes
+ for ( int iV = 0; iV < 2; ++iV )
+ {
+ TopoDS_Vertex v = helper.IthVertex( iV, edge );
+ TGeomID vID = meshDS->ShapeToIndex( v );
+ eos = data.GetShapeEdges( vID );
+ vertexLEdges[ iV ] = eos->_edges[ 0 ];
+ }
+ edgeLEdge = &vertexLEdges[0];
+ edgeLEdgeEnd = edgeLEdge + 2;
+
+ centerCurves[ iE ]._adjFace.Nullify();
+ }
+ else
+ {
+ if ( ! eos->_toSmooth )
+ data.SortOnEdge( edge, eos->_edges );
+ edgeLEdge = &eos->_edges[ 0 ];
+ edgeLEdgeEnd = edgeLEdge + eos->_edges.size();
+ vertexLEdges[0] = eos->_edges.front()->_2neibors->_edges[0];
+ vertexLEdges[1] = eos->_edges.back() ->_2neibors->_edges[1];
+
+ if ( ! eos->_sWOL.IsNull() )
+ centerCurves[ iE ]._adjFace.Nullify();
+ }
+
+ // Get curvature centers
+
+ centersBox.Clear();
+
+ if ( edgeLEdge[0]->IsOnEdge() &&
+ convFace.GetCenterOfCurvature( vertexLEdges[0], surfProp, helper, center ))
+ { // 1st VERTEX
+ centerCurves[ iE ].Append( center, vertexLEdges[0] );
+ centersBox.Add( center );
+ }
+ for ( ; edgeLEdge < edgeLEdgeEnd; ++edgeLEdge )
+ if ( convFace.GetCenterOfCurvature( *edgeLEdge, surfProp, helper, center ))
+ { // EDGE or VERTEXes
+ centerCurves[ iE ].Append( center, *edgeLEdge );
+ centersBox.Add( center );
+ }
+ if ( edgeLEdge[-1]->IsOnEdge() &&
+ convFace.GetCenterOfCurvature( vertexLEdges[1], surfProp, helper, center ))
+ { // 2nd VERTEX
+ centerCurves[ iE ].Append( center, vertexLEdges[1] );
+ centersBox.Add( center );
+ }
+ centerCurves[ iE ]._isDegenerated =
+ ( centersBox.IsVoid() || centersBox.SquareExtent() < 1e-6 * nodesBox.SquareExtent() );
+
+ } // loop on EDGES of convFace._face to set up data of centerCurves
+
+ // Compute new normals for _LayerEdge's on EDGEs
+
+ double avgCosin = 0;
+ int nbCosin = 0;
+ gp_Vec inFaceDir;
+ for ( size_t iE1 = 0; iE1 < centerCurves.size(); ++iE1 )
+ {
+ _CentralCurveOnEdge& ceCurve = centerCurves[ iE1 ];
+ if ( ceCurve._isDegenerated )
+ continue;
+ const vector< gp_Pnt >& centers = ceCurve._curvaCenters;
+ vector< gp_XYZ > & newNormals = ceCurve._normals;
+ for ( size_t iC1 = 0; iC1 < centers.size(); ++iC1 )
+ {
+ isOK = false;
+ for ( size_t iE2 = 0; iE2 < centerCurves.size() && !isOK; ++iE2 )
+ {
+ if ( iE1 != iE2 )
+ isOK = centerCurves[ iE2 ].FindNewNormal( centers[ iC1 ], newNormals[ iC1 ]);
+ }
+ if ( isOK && !ceCurve._adjFace.IsNull() )
+ {
+ // compute new _LayerEdge::_cosin
+ const SMDS_MeshNode* node = ceCurve._ledges[ iC1 ]->_nodes[0];
+ inFaceDir = getFaceDir( ceCurve._adjFace, ceCurve._edge, node, helper, isOK );
+ if ( isOK )
+ {
+ double angle = inFaceDir.Angle( newNormals[ iC1 ] ); // [0,PI]
+ ceCurve._ledges[ iC1 ]->_cosin = Cos( angle );
+ avgCosin += ceCurve._ledges[ iC1 ]->_cosin;
+ nbCosin++;
+ }
+ }
+ }
+ }
+ // set new normals to _LayerEdge's of NOT degenerated central curves
+ for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
+ {
+ if ( centerCurves[ iE ]._isDegenerated )
+ continue;
+ for ( size_t iLE = 0; iLE < centerCurves[ iE ]._ledges.size(); ++iLE )
+ {
+ centerCurves[ iE ]._ledges[ iLE ]->SetNormal( centerCurves[ iE ]._normals[ iLE ]);
+ centerCurves[ iE ]._ledges[ iLE ]->Set( _LayerEdge::NORMAL_UPDATED );
+ }
+ }
+ // set new normals to _LayerEdge's of degenerated central curves
+ for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
+ {
+ if ( !centerCurves[ iE ]._isDegenerated ||
+ centerCurves[ iE ]._ledges.size() < 3 )
+ continue;
+ // new normal is an average of new normals at VERTEXes that
+ // was computed on non-degenerated _CentralCurveOnEdge's
+ gp_XYZ newNorm = ( centerCurves[ iE ]._ledges.front()->_normal +
+ centerCurves[ iE ]._ledges.back ()->_normal );
+ double sz = newNorm.Modulus();
+ if ( sz < 1e-200 )
+ continue;
+ newNorm /= sz;
+ double newCosin = ( 0.5 * centerCurves[ iE ]._ledges.front()->_cosin +
+ 0.5 * centerCurves[ iE ]._ledges.back ()->_cosin );
+ for ( size_t iLE = 1, nb = centerCurves[ iE ]._ledges.size() - 1; iLE < nb; ++iLE )
+ {
+ centerCurves[ iE ]._ledges[ iLE ]->SetNormal( newNorm );
+ centerCurves[ iE ]._ledges[ iLE ]->_cosin = newCosin;
+ centerCurves[ iE ]._ledges[ iLE ]->Set( _LayerEdge::NORMAL_UPDATED );
+ }
+ }
+
+ // Find new normals for _LayerEdge's based on FACE
+
+ if ( nbCosin > 0 )
+ avgCosin /= nbCosin;
+ const TGeomID faceID = meshDS->ShapeToIndex( convFace._face );
+ map< TGeomID, _EdgesOnShape* >::iterator id2eos = convFace._subIdToEOS.find( faceID );
+ if ( id2eos != convFace._subIdToEOS.end() )
+ {
+ int iE = 0;
+ gp_XYZ newNorm;
+ _EdgesOnShape& eos = * ( id2eos->second );
+ for ( size_t i = 0; i < eos._edges.size(); ++i )
+ {
+ _LayerEdge* ledge = eos._edges[ i ];
+ if ( !convFace.GetCenterOfCurvature( ledge, surfProp, helper, center ))
+ continue;
+ for ( size_t i = 0; i < centerCurves.size(); ++i, ++iE )
+ {
+ iE = iE % centerCurves.size();
+ if ( centerCurves[ iE ]._isDegenerated )
+ continue;
+ newNorm.SetCoord( 0,0,0 );
+ if ( centerCurves[ iE ].FindNewNormal( center, newNorm ))
+ {
+ ledge->SetNormal( newNorm );
+ ledge->_cosin = avgCosin;
+ ledge->Set( _LayerEdge::NORMAL_UPDATED );
+ break;
+ }
+ }
+ }
+ }
+
+ } // not a quasi-spherical FACE
+
+ // Update _LayerEdge's data according to a new normal
+
+ dumpFunction(SMESH_Comment("updateNormalsOfConvexFaces")<<data._index
+ <<"_F"<<meshDS->ShapeToIndex( convFace._face ));
+
+ id2eos = convFace._subIdToEOS.begin();
+ for ( ; id2eos != convFace._subIdToEOS.end(); ++id2eos )
+ {
+ _EdgesOnShape& eos = * ( id2eos->second );
+ for ( size_t i = 0; i < eos._edges.size(); ++i )
+ {
+ _LayerEdge* & ledge = eos._edges[ i ];
+ double len = ledge->_len;
+ ledge->InvalidateStep( stepNb + 1, eos, /*restoreLength=*/true );
+ ledge->SetCosin( ledge->_cosin );
+ ledge->SetNewLength( len, eos, helper );
+ }
+ if ( eos.ShapeType() != TopAbs_FACE )
+ for ( size_t i = 0; i < eos._edges.size(); ++i )
+ {
+ _LayerEdge* ledge = eos._edges[ i ];
+ for ( size_t iN = 0; iN < ledge->_neibors.size(); ++iN )
+ {
+ _LayerEdge* neibor = ledge->_neibors[iN];
+ if ( neibor->_nodes[0]->GetPosition()->GetDim() == 2 )
+ {
+ neibor->Set( _LayerEdge::NEAR_BOUNDARY );
+ neibor->Set( _LayerEdge::MOVED );
+ neibor->SetSmooLen( neibor->_len );
+ }
+ }
+ }
+ } // loop on sub-shapes of convFace._face
+
+ // Find FACEs adjacent to convFace._face that got necessity to smooth
+ // as a result of normals modification
+
+ set< _EdgesOnShape* > adjFacesToSmooth;
+ for ( size_t iE = 0; iE < centerCurves.size(); ++iE )
+ {
+ if ( centerCurves[ iE ]._adjFace.IsNull() ||
+ centerCurves[ iE ]._adjFaceToSmooth )
+ continue;
+ for ( size_t iLE = 0; iLE < centerCurves[ iE ]._ledges.size(); ++iLE )
+ {
+ if ( centerCurves[ iE ]._ledges[ iLE ]->_cosin > theMinSmoothCosin )
+ {
+ adjFacesToSmooth.insert( data.GetShapeEdges( centerCurves[ iE ]._adjFace ));
+ break;
+ }
+ }
+ }
+ data.AddShapesToSmooth( adjFacesToSmooth );
+
+ dumpFunctionEnd();
+
+
+ } // loop on data._convexFaces
+
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Finds a center of curvature of a surface at a _LayerEdge
+ */
+//================================================================================
+
+bool _ConvexFace::GetCenterOfCurvature( _LayerEdge* ledge,
+ BRepLProp_SLProps& surfProp,
+ SMESH_MesherHelper& helper,
+ gp_Pnt & center ) const
+{
+ gp_XY uv = helper.GetNodeUV( _face, ledge->_nodes[0] );
+ surfProp.SetParameters( uv.X(), uv.Y() );
+ if ( !surfProp.IsCurvatureDefined() )
+ return false;
+
+ const double oriFactor = ( _face.Orientation() == TopAbs_REVERSED ? +1. : -1. );
+ double surfCurvatureMax = surfProp.MaxCurvature() * oriFactor;
+ double surfCurvatureMin = surfProp.MinCurvature() * oriFactor;
+ if ( surfCurvatureMin > surfCurvatureMax )
+ center = surfProp.Value().Translated( surfProp.Normal().XYZ() / surfCurvatureMin * oriFactor );
+ else
+ center = surfProp.Value().Translated( surfProp.Normal().XYZ() / surfCurvatureMax * oriFactor );
+
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Check that prisms are not distorted
+ */
+//================================================================================
+
+bool _ConvexFace::CheckPrisms() const
+{
+ double vol = 0;
+ for ( size_t i = 0; i < _simplexTestEdges.size(); ++i )
+ {
+ const _LayerEdge* edge = _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, vol ))
+ {
+ debugMsg( "Bad simplex of _simplexTestEdges ("
+ << " "<< edge->_nodes[0]->GetID()<< " "<< tgtXYZ._node->GetID()
+ << " "<< edge->_simplices[j]._nPrev->GetID()
+ << " "<< edge->_simplices[j]._nNext->GetID() << " )" );
+ return false;
+ }
+ }
+ return true;
+}
+
+//================================================================================
+/*!
+ * \brief Try to compute a new normal by interpolating normals of _LayerEdge's
+ * stored in this _CentralCurveOnEdge.
+ * \param [in] center - curvature center of a point of another _CentralCurveOnEdge.
+ * \param [in,out] newNormal - current normal at this point, to be redefined
+ * \return bool - true if succeeded.
+ */
+//================================================================================
+
+bool _CentralCurveOnEdge::FindNewNormal( const gp_Pnt& center, gp_XYZ& newNormal )
+{
+ if ( this->_isDegenerated )
+ return false;
+
+ // find two centers the given one lies between
+
+ for ( size_t i = 0, nb = _curvaCenters.size()-1; i < nb; ++i )
+ {
+ double sl2 = 1.001 * _segLength2[ i ];
+
+ double d1 = center.SquareDistance( _curvaCenters[ i ]);
+ if ( d1 > sl2 )
+ continue;
+
+ double d2 = center.SquareDistance( _curvaCenters[ i+1 ]);
+ if ( d2 > sl2 || d2 + d1 < 1e-100 )
+ continue;
+
+ d1 = Sqrt( d1 );
+ d2 = Sqrt( d2 );
+ double r = d1 / ( d1 + d2 );
+ gp_XYZ norm = (( 1. - r ) * _ledges[ i ]->_normal +
+ ( r ) * _ledges[ i+1 ]->_normal );
+ norm.Normalize();
+
+ newNormal += norm;
+ double sz = newNormal.Modulus();
+ if ( sz < 1e-200 )
+ break;
+ newNormal /= sz;
+ return true;
+ }
+ return false;
+}
+
+//================================================================================
+/*!
+ * \brief Set shape members
+ */
+//================================================================================
+
+void _CentralCurveOnEdge::SetShapes( const TopoDS_Edge& edge,
+ const _ConvexFace& convFace,
+ _SolidData& data,
+ SMESH_MesherHelper& helper)
+{
+ _edge = edge;
+
+ PShapeIteratorPtr fIt = helper.GetAncestors( edge, *helper.GetMesh(), TopAbs_FACE );
+ while ( const TopoDS_Shape* F = fIt->next())
+ if ( !convFace._face.IsSame( *F ))
+ {
+ _adjFace = TopoDS::Face( *F );
+ _adjFaceToSmooth = false;
+ // _adjFace already in a smoothing queue ?
+ if ( _EdgesOnShape* eos = data.GetShapeEdges( _adjFace ))
+ _adjFaceToSmooth = eos->_toSmooth;
+ break;
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Looks for intersection of it's last segment with faces
+ * \param distance - returns shortest distance from the last node to intersection
+ */
+//================================================================================
+
+bool _LayerEdge::FindIntersection( SMESH_ElementSearcher& searcher,
+ double & distance,
+ const double& epsilon,
+ _EdgesOnShape& eos,
+ const SMDS_MeshElement** intFace)
+{
+ vector< const SMDS_MeshElement* > suspectFaces;
+ double segLen;
+ gp_Ax1 lastSegment = LastSegment( segLen, eos );
+ searcher.GetElementsNearLine( lastSegment, SMDSAbs_Face, suspectFaces );
+
+ bool segmentIntersected = false;
+ distance = Precision::Infinite();
+ int iFace = -1; // intersected face
+ for ( size_t j = 0 ; j < suspectFaces.size() /*&& !segmentIntersected*/; ++j )
+ {
+ const SMDS_MeshElement* face = suspectFaces[j];
+ if ( face->GetNodeIndex( _nodes.back() ) >= 0 ||
+ face->GetNodeIndex( _nodes[0] ) >= 0 )
+ continue; // face sharing _LayerEdge node
+ const int nbNodes = face->NbCornerNodes();
+ bool intFound = false;
+ double dist;
+ SMDS_MeshElement::iterator nIt = face->begin_nodes();
+ if ( nbNodes == 3 )
+ {
+ intFound = SegTriaInter( lastSegment, *nIt++, *nIt++, *nIt++, dist, epsilon );
+ }
+ else
+ {
+ const SMDS_MeshNode* tria[3];
+ tria[0] = *nIt++;
+ tria[1] = *nIt++;
+ for ( int n2 = 2; n2 < nbNodes && !intFound; ++n2 )
+ {
+ tria[2] = *nIt++;
+ intFound = SegTriaInter(lastSegment, tria[0], tria[1], tria[2], dist, epsilon );
+ tria[1] = tria[2];
+ }
+ }
+ if ( intFound )
+ {
+ if ( dist < segLen*(1.01) && dist > -(_len*_lenFactor-segLen) )
+ segmentIntersected = true;
+ if ( distance > dist )
+ distance = dist, iFace = j;
+ }
+ }
+ if ( intFace ) *intFace = ( iFace != -1 ) ? suspectFaces[iFace] : 0;
+