+ }
+ avgNormProj /= edge->_simplices.size();
+ avgLen /= edge->_simplices.size();
+ if (( edge->_curvature = _Curvature::New( avgNormProj, avgLen )))
+ {
+ edge->Set( _LayerEdge::SMOOTHED_C1 );
+ isCurved = true;
+ SMDS_FacePositionPtr fPos = edge->_nodes[0]->GetPosition();
+ if ( !fPos )
+ for ( size_t iS = 0; iS < edge->_simplices.size() && !fPos; ++iS )
+ fPos = 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, eosI._hyp.ToSmooth() );
+ }
+ }
+ }
+ }
+ 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, eosI._hyp.ToSmooth() );
+ e0 = eI;
+ }
+ }
+ }
+}
+
+//================================================================================
+/*!
+ * \brief Limit _LayerEdge::_maxLen according to local curvature
+ */
+//================================================================================
+
+void _ViscousBuilder::limitMaxLenByCurvature( _LayerEdge* e1,
+ _LayerEdge* e2,
+ _EdgesOnShape& eos1,
+ _EdgesOnShape& eos2,
+ const bool isSmoothable )
+{
+ if (( e1->_nodes[0]->GetPosition()->GetDim() !=
+ e2->_nodes[0]->GetPosition()->GetDim() ) &&
+ ( e1->_cosin < 0.75 ))
+ return; // angle > 90 deg at e1
+
+ 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 )
+ {
+ const double coef = 0.75;
+ e1->SetMaxLen( Min( e1->_maxLen, coef * u1 / e1->_lenFactor ));
+ e2->SetMaxLen( Min( e2->_maxLen, coef * 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;
+ if ( !eos._sWOL.IsNull() ) continue; // PAL23566
+
+ 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 containing 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->n(0), f->n(1), f->n(2), dist1, eps ))
+ dist1 = Precision::Infinite();
+ if ( !edge->SegTriaInter( lastSegment, f->n(3), f->n(2), f->n(0), dist2, eps ))
+ dist2 = Precision::Infinite();
+ if (( dist1 > segLen ) && ( dist2 > segLen ))
+ continue;
+
+ if ( edge->IsOnEdge() )
+ {
+ // skip perpendicular EDGEs
+ gp_Vec fSegDir = SMESH_TNodeXYZ( f->n(0) ) - SMESH_TNodeXYZ( f->n(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->SetMaxLen( 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 Find _LayerEdge's located on boundary of a convex FACE whose normal
+ * will be updated at each inflation step
+ */
+//================================================================================
+
+void _ViscousBuilder::findEdgesToUpdateNormalNearConvexFace( _ConvexFace & convFace,
+ _SolidData& data,
+ SMESH_MesherHelper& helper )
+{
+ const TGeomID convFaceID = getMeshDS()->ShapeToIndex( convFace._face );
+ const double preci = BRep_Tool::Tolerance( convFace._face );
+ Handle(ShapeAnalysis_Surface) surface = helper.GetSurface( convFace._face );
+
+ bool edgesToUpdateFound = false;
+
+ map< TGeomID, _EdgesOnShape* >::iterator id2eos = convFace._subIdToEOS.begin();
+ for ( ; id2eos != convFace._subIdToEOS.end(); ++id2eos )
+ {
+ _EdgesOnShape& eos = * id2eos->second;
+ if ( !eos._sWOL.IsNull() ) continue;
+ if ( !eos._hyp.ToSmooth() ) continue;
+ for ( size_t i = 0; i < eos._edges.size(); ++i )
+ {
+ _LayerEdge* ledge = eos._edges[ i ];
+ if ( ledge->Is( _LayerEdge::UPD_NORMAL_CONV )) continue; // already checked
+ if ( ledge->Is( _LayerEdge::MULTI_NORMAL )) continue; // not inflatable
+
+ gp_XYZ tgtPos = ( SMESH_NodeXYZ( ledge->_nodes[0] ) +
+ ledge->_normal * ledge->_lenFactor * ledge->_maxLen );
+
+ // the normal must be updated if distance from tgtPos to surface is less than
+ // target thickness
+
+ // find an initial UV for search of a projection of tgtPos to surface
+ const SMDS_MeshNode* nodeInFace = 0;
+ SMDS_ElemIteratorPtr fIt = ledge->_nodes[0]->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() && !nodeInFace )