X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_ProjectionUtils.cxx;h=cda4346fa753a7b2fff9e103289998ecf62c6a36;hp=6079b0c5442464b005b1fcbf884f3da8e774b355;hb=9ecacf41d1f3bb33c5cee644b5359e164b2302cf;hpb=9c9403c2720f6d575af952664e180ab9e5ed3143 diff --git a/src/StdMeshers/StdMeshers_ProjectionUtils.cxx b/src/StdMeshers/StdMeshers_ProjectionUtils.cxx index 6079b0c54..cda4346fa 100644 --- a/src/StdMeshers/StdMeshers_ProjectionUtils.cxx +++ b/src/StdMeshers/StdMeshers_ProjectionUtils.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2014 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE // // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS @@ -27,30 +27,34 @@ // #include "StdMeshers_ProjectionUtils.hxx" -#include "StdMeshers_ProjectionSource1D.hxx" -#include "StdMeshers_ProjectionSource2D.hxx" -#include "StdMeshers_ProjectionSource3D.hxx" - #include "SMDS_EdgePosition.hxx" +#include "SMDS_FacePosition.hxx" +#include "SMESHDS_Mesh.hxx" #include "SMESH_Algo.hxx" #include "SMESH_Block.hxx" #include "SMESH_Gen.hxx" #include "SMESH_HypoFilter.hxx" #include "SMESH_Hypothesis.hxx" #include "SMESH_Mesh.hxx" +#include "SMESH_MeshAlgos.hxx" #include "SMESH_MesherHelper.hxx" #include "SMESH_subMesh.hxx" #include "SMESH_subMeshEventListener.hxx" -#include "SMESH_MeshAlgos.hxx" +#include "StdMeshers_ProjectionSource1D.hxx" +#include "StdMeshers_ProjectionSource2D.hxx" +#include "StdMeshers_ProjectionSource3D.hxx" #include "utilities.h" #include +#include #include #include #include #include #include +#include +#include #include #include #include @@ -66,6 +70,7 @@ #include #include #include +#include #include #include @@ -76,36 +81,39 @@ using namespace std; #define RETURN_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); return false; } #define CONT_BAD_RESULT(msg) { MESSAGE(")-: Error: " << msg); continue; } #define SHOW_SHAPE(v,msg) \ -// { \ -// if ( (v).IsNull() ) cout << msg << " NULL SHAPE" << endl; \ -// else if ((v).ShapeType() == TopAbs_VERTEX) {\ -// gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( (v) ));\ -// cout<::const_iterator e = l.begin();\ -// for ( int i = 0; e != l.end(); ++e, ++i ) {\ -// cout << i << "V (" << TopExp::FirstVertex( *e, true ).TShape().operator->() << ") "\ -// << i << "E (" << e->TShape().operator->() << "); "; }\ -// cout << endl;\ -// } + // { show_list((msg),(l)); } namespace HERE = StdMeshers_ProjectionUtils; namespace { - static SMESHDS_Mesh* theMeshDS[2] = { 0, 0 }; // used to debug only + static SMESHDS_Mesh* theMeshDS[2] = { 0, 0 }; // used for debug only long shapeIndex(const TopoDS_Shape& S) { if ( theMeshDS[0] && theMeshDS[1] ) return max(theMeshDS[0]->ShapeToIndex(S), theMeshDS[1]->ShapeToIndex(S) ); return long(S.TShape().operator->()); } - + void show_shape( TopoDS_Shape v, const char* msg ) // debug + { + if ( v.IsNull() ) cout << msg << " NULL SHAPE" << endl; + else if (v.ShapeType() == TopAbs_VERTEX) { + gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( v )); + cout<& l ) // debug + { + cout << msg << " "; + list< TopoDS_Edge >::const_iterator e = l.begin(); + for ( int i = 0; e != l.end(); ++e, ++i ) { + cout << i << "V (" << TopExp::FirstVertex( *e, true ).TShape().operator->() << ") " + << i << "E (" << e->TShape().operator->() << "); "; } + cout << endl; + } //================================================================================ /*! * \brief Write shape for debug purposes @@ -118,6 +126,10 @@ namespace { const char* type[] ={"COMPOUND","COMPSOLID","SOLID","SHELL","FACE","WIRE","EDGE","VERTEX"}; BRepTools::Write( shape, SMESH_Comment("/tmp/") << type[shape.ShapeType()] << "_" << shape.TShape().operator->() << ".brep"); + if ( !theMeshDS[0] ) { + show_shape( TopoDS_Shape(), "avoid warning: show_shape() defined but not used"); + show_list( "avoid warning: show_list() defined but not used", list< TopoDS_Edge >() ); + } #endif return false; } @@ -231,6 +243,7 @@ namespace { v2 = SMESH_MesherHelper::IthVertex( 0, *eIt2 ); HERE::InsertAssociation( v1, v2, theMap ); } + theMap.SetAssocType( HERE::TShapeShapeMap::FEW_EF ); return true; } return false; @@ -386,6 +399,7 @@ namespace { } } } + theMap.SetAssocType( HERE::TShapeShapeMap::PROPAGATION ); return true; } @@ -402,9 +416,8 @@ namespace { const gp_Pnt2d& uv, const double& tol2d ) { - TopoDS_Vertex VV[2]; - TopExp::Vertices( edge, VV[0], VV[1], true); - gp_Pnt2d v1UV = BRep_Tool::Parameters( VV[vIndex], face); + TopoDS_Vertex V = SMESH_MesherHelper::IthVertex( vIndex, edge, /*CumOri=*/true ); + gp_Pnt2d v1UV = BRep_Tool::Parameters( V, face); double dist2d = v1UV.Distance( uv ); return dist2d < tol2d; } @@ -430,23 +443,36 @@ namespace { return true; } } + SMESH_MesherHelper helper( mesh ); + helper.SetSubShape( shape ); + TopExp_Explorer expF( shape, TopAbs_FACE ), expE; if ( expF.More() ) { for ( ; expF.More(); expF.Next() ) { TopoDS_Shape wire = StdMeshers_ProjectionUtils::OuterShape( TopoDS::Face( expF.Current() ), TopAbs_WIRE ); for ( expE.Init( wire, TopAbs_EDGE ); expE.More(); expE.Next() ) - if ( !SMESH_MesherHelper::IsClosedEdge( TopoDS::Edge( expE.Current() ))) - allBndEdges.push_back( TopoDS::Edge( expE.Current() )); + if ( ! helper.IsClosedEdge( TopoDS::Edge( expE.Current() ))) + { + if ( helper.IsSeamShape( expE.Current() )) + allBndEdges.push_back( TopoDS::Edge( expE.Current() )); + else + allBndEdges.push_front( TopoDS::Edge( expE.Current() )); + } } } else if ( shape.ShapeType() != TopAbs_EDGE) { // no faces for ( expE.Init( shape, TopAbs_EDGE ); expE.More(); expE.Next() ) - if ( !SMESH_MesherHelper::IsClosedEdge( TopoDS::Edge( expE.Current() ))) - allBndEdges.push_back( TopoDS::Edge( expE.Current() )); + if ( ! helper.IsClosedEdge( TopoDS::Edge( expE.Current() ))) + { + if ( helper.IsSeamShape( expE.Current() )) + allBndEdges.push_back( TopoDS::Edge( expE.Current() )); + else + allBndEdges.push_front( TopoDS::Edge( expE.Current() )); + } } else if ( shape.ShapeType() == TopAbs_EDGE ) { - if ( !SMESH_MesherHelper::IsClosedEdge( TopoDS::Edge( shape ))) + if ( ! helper.IsClosedEdge( TopoDS::Edge( shape ))) allBndEdges.push_back( TopoDS::Edge( shape )); } return !allBndEdges.empty(); @@ -475,9 +501,9 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the { // Structure of this long function is following // 1) Group -> Group projection: theShape1 is a group member, - // theShape2 is another group. We find a group theShape1 is in and recall self. + // theShape2 is another group. We find the group theShape1 is in and recall self. // 2) Accosiate same shapes with different location (partners). - // 3) If vertex association is given, perform accosiation according to shape type: + // 3) If vertex association is given, perform association according to shape type: // switch ( ShapeType ) { // case TopAbs_EDGE: // case ...: @@ -497,7 +523,8 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the // ================================================================================= // 1) Is it the case of associating a group member -> another group? (PAL16202, 16203) // ================================================================================= - if ( theShape1.ShapeType() != theShape2.ShapeType() ) { + if ( theShape1.ShapeType() != theShape2.ShapeType() ) + { TopoDS_Shape group1, group2; if ( theShape1.ShapeType() == TopAbs_COMPOUND ) { group1 = theShape1; @@ -529,11 +556,14 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the TShapePairsList::iterator s1_s2 = shapesQueue.begin(); for ( ; s1_s2 != shapesQueue.end(); ++s1_s2 ) { + if ( theMap.IsBound( s1_s2->first )) // avoid re-binding for a seam edge + continue; // to avoid this: Forward seam -> Reversed seam InsertAssociation( s1_s2->first, s1_s2->second, theMap ); TopoDS_Iterator s1It( s1_s2->first), s2It( s1_s2->second ); for ( ; s1It.More(); s1It.Next(), s2It.Next() ) shapesQueue.push_back( make_pair( s1It.Value(), s2It.Value() )); } + theMap.SetAssocType( TShapeShapeMap::PARTNER ); return true; } @@ -542,11 +572,13 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the //====================================================================== // 3) HAS initial vertex association //====================================================================== + bool isVCloseness = ( theMap._assocType == TShapeShapeMap::CLOSE_VERTEX ); + theMap.SetAssocType( TShapeShapeMap::INIT_VERTEX ); switch ( theShape1.ShapeType() ) { // ---------------------------------------------------------------------- case TopAbs_EDGE: { // TopAbs_EDGE // ---------------------------------------------------------------------- - if ( theMap.Extent() != 2 ) + if ( theMap.Extent() != 1 ) RETURN_BAD_RESULT("Wrong map extent " << theMap.Extent() ); TopoDS_Edge edge1 = TopoDS::Edge( theShape1 ); TopoDS_Edge edge2 = TopoDS::Edge( theShape2 ); @@ -590,7 +622,7 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the } } list< TopoDS_Edge > edges1, edges2; - int nbE = FindFaceAssociation( face1, VV1, face2, VV2, edges1, edges2 ); + int nbE = FindFaceAssociation( face1, VV1, face2, VV2, edges1, edges2, isVCloseness ); if ( !nbE ) RETURN_BAD_RESULT("FindFaceAssociation() failed"); fixAssocByPropagation( nbE, edges1, edges2, theMesh1, theMesh2 ); @@ -654,7 +686,6 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the TopoDS_Shape F1, F2; // get a face sharing edge1 (F1) - TopoDS_Shape FF2[2]; TopTools_ListIteratorOfListOfShape ancestIt1( edgeToFace1.FindFromKey( edge1 )); for ( ; F1.IsNull() && ancestIt1.More(); ancestIt1.Next() ) if ( ancestIt1.Value().ShapeType() == TopAbs_FACE ) @@ -664,6 +695,7 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the RETURN_BAD_RESULT(" Face1 not found"); // get 2 faces sharing edge2 (one of them is F2) + TopoDS_Shape FF2[2]; TopTools_ListIteratorOfListOfShape ancestIt2( edgeToFace2.FindFromKey( edge2 )); for ( int i = 0; FF2[1].IsNull() && ancestIt2.More(); ancestIt2.Next() ) if ( ancestIt2.Value().ShapeType() == TopAbs_FACE ) @@ -694,8 +726,6 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the F2 = FF2[ 1 ]; } - TopTools_MapOfShape boundEdges; - // association of face sub-shapes and neighbour faces list< pair < TopoDS_Face, TopoDS_Edge > > FE1, FE2; list< pair < TopoDS_Face, TopoDS_Edge > >::iterator fe1, fe2; @@ -711,11 +741,11 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the TopExp::Vertices( edge1, VV1[0], VV1[1], true ); TopExp::Vertices( edge2, VV2[0], VV2[1], true ); list< TopoDS_Edge > edges1, edges2; - int nbE = FindFaceAssociation( face1, VV1, face2, VV2, edges1, edges2 ); + int nbE = FindFaceAssociation( face1, VV1, face2, VV2, edges1, edges2, isVCloseness ); if ( !nbE ) RETURN_BAD_RESULT("FindFaceAssociation() failed"); InsertAssociation( face1, face2, theMap ); // assoc faces - MESSAGE("Assoc FACE " << theMesh1->GetMeshDS()->ShapeToIndex( face1 )<< - " to " << theMesh2->GetMeshDS()->ShapeToIndex( face2 )); + // MESSAGE("Assoc FACE " << theMesh1->GetMeshDS()->ShapeToIndex( face1 )<< + // " to " << theMesh2->GetMeshDS()->ShapeToIndex( face2 )); if ( nbE == 2 && (edge1.IsSame( edges1.front())) != (edge2.IsSame( edges2.front()))) { reverseEdges( edges2, nbE ); @@ -724,8 +754,8 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the list< TopoDS_Edge >::iterator eIt2 = edges2.begin(); for ( ; eIt1 != edges1.end(); ++eIt1, ++eIt2 ) { - if ( !boundEdges.Add( *eIt1 )) continue; // already associated - InsertAssociation( *eIt1, *eIt2, theMap ); // assoc edges + if ( !InsertAssociation( *eIt1, *eIt2, theMap )) // assoc edges + continue; // already associated VV1[0] = TopExp::FirstVertex( *eIt1, true ); VV2[0] = TopExp::FirstVertex( *eIt2, true ); InsertAssociation( VV1[0], VV2[0], theMap ); // assoc vertices @@ -734,6 +764,10 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the TopoDS_Face nextFace1 = GetNextFace( edgeToFace1, *eIt1, face1 ); TopoDS_Face nextFace2 = GetNextFace( edgeToFace2, *eIt2, face2 ); if ( !nextFace1.IsNull() && !nextFace2.IsNull() ) { + if ( SMESH_MesherHelper::GetSubShapeOri( nextFace1, *eIt1 ) == eIt1->Orientation() ) + nextFace1.Reverse(); + if ( SMESH_MesherHelper::GetSubShapeOri( nextFace2, *eIt2 ) == eIt2->Orientation() ) + nextFace2.Reverse(); FE1.push_back( make_pair( nextFace1, *eIt1 )); FE2.push_back( make_pair( nextFace2, *eIt2 )); } @@ -815,8 +849,8 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the } // Associate shells // - int nbFaces1 = SMESH_MesherHelper:: Count( shell1, TopAbs_FACE, 0 ); - int nbFaces2 = SMESH_MesherHelper:: Count( shell2, TopAbs_FACE, 0 ); + int nbFaces1 = SMESH_MesherHelper::Count( shell1, TopAbs_FACE, 0 ); + int nbFaces2 = SMESH_MesherHelper::Count( shell2, TopAbs_FACE, 0 ); if ( nbFaces1 != nbFaces2 ) RETURN_BAD_RESULT("Different nb of faces found for shells"); if ( nbFaces1 > 0 ) { @@ -922,14 +956,14 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the v2e[0].UnBind( V[0] ); v2e[1].UnBind( V[1] ); InsertAssociation( e0, e1, theMap ); - MESSAGE("Assoc edge " << theMesh1->GetMeshDS()->ShapeToIndex( e0 )<< - " to " << theMesh2->GetMeshDS()->ShapeToIndex( e1 )); + // MESSAGE("Assoc edge " << theMesh1->GetMeshDS()->ShapeToIndex( e0 )<< + // " to " << theMesh2->GetMeshDS()->ShapeToIndex( e1 )); V[0] = GetNextVertex( e0, V[0] ); V[1] = GetNextVertex( e1, V[1] ); if ( !V[0].IsNull() ) { InsertAssociation( V[0], V[1], theMap ); - MESSAGE("Assoc vertex " << theMesh1->GetMeshDS()->ShapeToIndex( V[0] )<< - " to " << theMesh2->GetMeshDS()->ShapeToIndex( V[1] )); + // MESSAGE("Assoc vertex " << theMesh1->GetMeshDS()->ShapeToIndex( V[0] )<< + // " to " << theMesh2->GetMeshDS()->ShapeToIndex( V[1] )); } } else if ( nbE0 == 2 ) @@ -956,12 +990,12 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the InsertAssociation( e0b, e1b, theMap ); InsertAssociation( e0n, e1n, theMap ); InsertAssociation( v0n, v1n, theMap ); - MESSAGE("Assoc edge " << theMesh1->GetMeshDS()->ShapeToIndex( e0b )<< - " to " << theMesh2->GetMeshDS()->ShapeToIndex( e1b )); - MESSAGE("Assoc edge " << theMesh1->GetMeshDS()->ShapeToIndex( e0n )<< - " to " << theMesh2->GetMeshDS()->ShapeToIndex( e1n )); - MESSAGE("Assoc vertex " << theMesh1->GetMeshDS()->ShapeToIndex( v0n )<< - " to " << theMesh2->GetMeshDS()->ShapeToIndex( v1n )); + // MESSAGE("Assoc edge " << theMesh1->GetMeshDS()->ShapeToIndex( e0b )<< + // " to " << theMesh2->GetMeshDS()->ShapeToIndex( e1b )); + // MESSAGE("Assoc edge " << theMesh1->GetMeshDS()->ShapeToIndex( e0n )<< + // " to " << theMesh2->GetMeshDS()->ShapeToIndex( e1n )); + // MESSAGE("Assoc vertex " << theMesh1->GetMeshDS()->ShapeToIndex( v0n )<< + // " to " << theMesh2->GetMeshDS()->ShapeToIndex( v1n )); v2e[0].UnBind( V[0] ); v2e[1].UnBind( V[1] ); V[0] = v0n; @@ -1007,6 +1041,7 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the InsertAssociation( edge1, prpEdge, theMap ); // insert with a proper orientation } InsertAssociation( theShape1, theShape2, theMap ); + theMap.SetAssocType( TShapeShapeMap::PROPAGATION ); return true; // done } } @@ -1065,6 +1100,11 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the // take care of proper association of propagated edges bool same1 = edge1.IsSame( edges1.front() ); bool same2 = edge2.IsSame( edges2.front() ); + if ( !same1 && !same2 ) + { + same1 = ( edges1.back().Orientation() == edge1.Orientation() ); + same2 = ( edges2.back().Orientation() == edge2.Orientation() ); + } if ( same1 != same2 ) { reverseEdges(edges2, nbE); @@ -1077,11 +1117,12 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the for ( ; eIt1 != edges1.end(); ++eIt1, ++eIt2 ) { InsertAssociation( *eIt1, *eIt2, theMap ); - VV1[0] = TopExp::FirstVertex( *eIt1, true ); - VV2[0] = TopExp::FirstVertex( *eIt2, true ); + VV1[0] = SMESH_MesherHelper::IthVertex( 0, *eIt1, true ); + VV2[0] = SMESH_MesherHelper::IthVertex( 0, *eIt2, true ); InsertAssociation( VV1[0], VV2[0], theMap ); } InsertAssociation( theShape1, theShape2, theMap ); + theMap.SetAssocType( TShapeShapeMap::PROPAGATION ); return true; } } @@ -1146,7 +1187,11 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the if ( !VV1[1].IsNull() ) { InsertAssociation( VV1[0], VV2[0], theMap ); InsertAssociation( VV1[1], VV2[1], theMap ); - return FindSubShapeAssociation( theShape1, theMesh1, theShape2, theMesh2, theMap); + TShapeShapeMap::EAssocType asType = theMap._assocType; + theMap.SetAssocType( TShapeShapeMap::PROPAGATION ); + if ( FindSubShapeAssociation( theShape1, theMesh1, theShape2, theMesh2, theMap )) + return true; + theMap._assocType = asType; } } break; // try by vertex closeness @@ -1169,9 +1214,13 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the RETURN_BAD_RESULT("Different nb of vertices"); } - if ( vMap1.Extent() == 1 ) { + if ( vMap1.Extent() == 1 || vMap2.Extent() == 1 ) { InsertAssociation( vMap1(1), vMap2(1), theMap ); if ( theShape1.ShapeType() == TopAbs_EDGE ) { + if ( vMap1.Extent() == 2 ) + InsertAssociation( vMap1(2), vMap2(1), theMap ); + else if ( vMap2.Extent() == 2 ) + InsertAssociation( vMap2(2), vMap1(1), theMap ); InsertAssociation( theShape1, theShape2, theMap ); return true; } @@ -1198,8 +1247,11 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the { InsertAssociation( VV1[0], VV1[0], theMap ); InsertAssociation( VV1[1], VV1[1], theMap ); - if (FindSubShapeAssociation( theShape1, theMesh1, theShape2, theMesh2, theMap )) + TShapeShapeMap::EAssocType asType = theMap._assocType; + theMap.SetAssocType( TShapeShapeMap::COMMON_VERTEX ); + if ( FindSubShapeAssociation( theShape1, theMesh1, theShape2, theMesh2, theMap )) return true; + theMap._assocType = asType; } } } @@ -1238,8 +1290,9 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the double minDist = std::numeric_limits::max(); for ( int nbChecked=0; edge1 != allBndEdges1.end() && nbChecked++ < 10; ++edge1 ) { - TopExp::Vertices( TopoDS::Edge( edge1->Oriented(TopAbs_FORWARD)), VV1[0], VV1[1]); - if ( VV1[0].IsSame( VV1[1] )) + TopoDS_Vertex edge1VV[2]; + TopExp::Vertices( TopoDS::Edge( edge1->Oriented(TopAbs_FORWARD)), edge1VV[0], edge1VV[1]); + if ( edge1VV[0].IsSame( edge1VV[1] )) continue;//RETURN_BAD_RESULT("Only closed edges"); // find vertices closest to 2 linked vertices of shape 1 @@ -1247,7 +1300,7 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the TopoDS_Vertex edge2VV[2]; for ( int i1 = 0; i1 < 2; ++i1 ) { - gp_Pnt p1 = BRep_Tool::Pnt( VV1[ i1 ]); + gp_Pnt p1 = BRep_Tool::Pnt( edge1VV[ i1 ]); p1.Scale( gc[0], scale ); p1.Translate( vec01 ); if ( !i1 ) { @@ -1283,6 +1336,8 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the } } if ( dist2[0] + dist2[1] < minDist ) { + VV1[0] = edge1VV[0]; + VV1[1] = edge1VV[1]; VV2[0] = edge2VV[0]; VV2[1] = edge2VV[1]; minDist = dist2[0] + dist2[1]; @@ -1290,13 +1345,14 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the break; } } + theMap.SetAssocType( TShapeShapeMap::CLOSE_VERTEX ); InsertAssociation( VV1[ 0 ], VV2[ 0 ], theMap ); InsertAssociation( VV1[ 1 ], VV2[ 1 ], theMap ); - MESSAGE("Initial assoc VERT " << theMesh1->GetMeshDS()->ShapeToIndex( VV1[ 0 ] )<< - " to " << theMesh2->GetMeshDS()->ShapeToIndex( VV2[ 0 ] )<< - "\nand VERT " << theMesh1->GetMeshDS()->ShapeToIndex( VV1[ 1 ] )<< - " to " << theMesh2->GetMeshDS()->ShapeToIndex( VV2[ 1 ] )); + // MESSAGE("Initial assoc VERT " << theMesh1->GetMeshDS()->ShapeToIndex( VV1[ 0 ] )<< + // " to " << theMesh2->GetMeshDS()->ShapeToIndex( VV2[ 0 ] )<< + // "\nand VERT " << theMesh1->GetMeshDS()->ShapeToIndex( VV1[ 1 ] )<< + // " to " << theMesh2->GetMeshDS()->ShapeToIndex( VV2[ 1 ] )); if ( theShape1.ShapeType() == TopAbs_EDGE ) { InsertAssociation( theShape1, theShape2, theMap ); return true; @@ -1314,6 +1370,7 @@ bool StdMeshers_ProjectionUtils::FindSubShapeAssociation(const TopoDS_Shape& the * \param VV2 - vertices of face 2 associated with ones of face 1 * \param edges1 - out list of edges of face 1 * \param edges2 - out list of edges of face 2 + * \param isClosenessAssoc - is association starting by VERTEX closeness * \retval int - nb of edges in an outer wire in a success case, else zero */ //================================================================================ @@ -1323,7 +1380,8 @@ int StdMeshers_ProjectionUtils::FindFaceAssociation(const TopoDS_Face& face1, const TopoDS_Face& face2, TopoDS_Vertex VV2[2], list< TopoDS_Edge > & edges1, - list< TopoDS_Edge > & edges2) + list< TopoDS_Edge > & edges2, + const bool isClosenessAssoc) { bool OK = false; list< int > nbEInW1, nbEInW2; @@ -1352,46 +1410,56 @@ int StdMeshers_ProjectionUtils::FindFaceAssociation(const TopoDS_Face& face1, // Define if we need to reverse one of wires to make edges in lists match each other bool reverse = false; + const bool severalWires = ( nbEInW1.size() > 1 ); - if ( !VV1[1].IsSame( TopExp::LastVertex( edges1.front(), true ))) { + if ( !VV1[1].IsSame( TopExp::LastVertex( edges1.front(), true ))) + { reverse = true; - edgeIt = --edges1.end(); // check if the second vertex belongs to the first or last edge in the wire + edgeIt = --edges1.end(); // pointer to the last edge in the outer wire + if ( severalWires ) { + edgeIt = edges1.begin(); + std::advance( edgeIt, nbEInW1.front()-1 ); + } + if ( TopExp::FirstVertex( *edgeIt ).IsSame( TopExp::LastVertex( *edgeIt )) && + SMESH_Algo::isDegenerated( *edgeIt )) { + --edgeIt; // skip a degenerated edge (test 3D_mesh_Projection_00/A3) + } if ( !VV1[1].IsSame( TopExp::FirstVertex( *edgeIt, true ))) { - bool KO = true; // belongs to none - if ( nbEInW1.size() > 1 ) { // several wires - edgeIt = edges1.begin(); - std::advance( edgeIt, nbEInW1.front()-1 ); - KO = !VV1[1].IsSame( TopExp::FirstVertex( *edgeIt, true )); - } - if ( KO ) - CONT_BAD_RESULT("GetOrderedEdges() failed"); + CONT_BAD_RESULT("GetOrderedEdges() failed"); } } - if ( !VV2[1].IsSame( TopExp::LastVertex( edges2.front(), true ))) { + if ( !VV2[1].IsSame( TopExp::LastVertex( edges2.front(), true ))) + { reverse = !reverse; - edgeIt = --edges2.end(); - // move a degenerated edge from back to front - // http://www.salome-platform.org/forum/forum_11/173031193 - if ( TopExp::FirstVertex( *edgeIt ).IsSame( TopExp::LastVertex( *edgeIt ))) { - edges2.splice( edges2.begin(), edges2, edgeIt ); - edgeIt = --edges2.end(); - } // check if the second vertex belongs to the first or last edge in the wire + edgeIt = --edges2.end(); // pointer to the last edge in the outer wire + if ( severalWires ) { + edgeIt = edges2.begin(); + std::advance( edgeIt, nbEInW2.front()-1 ); + } + if ( TopExp::FirstVertex( *edgeIt ).IsSame( TopExp::LastVertex( *edgeIt )) && + SMESH_Algo::isDegenerated( *edgeIt )) { + --edgeIt; // skip a degenerated edge + } if ( !VV2[1].IsSame( TopExp::FirstVertex( *edgeIt, true ))) { - bool KO = true; // belongs to none - if ( nbEInW2.size() > 1 ) { // several wires - edgeIt = edges2.begin(); - std::advance( edgeIt, nbEInW2.front()-1 ); - KO = !VV2[1].IsSame( TopExp::FirstVertex( *edgeIt, true )); - } - if ( KO ) - CONT_BAD_RESULT("GetOrderedEdges() failed"); + CONT_BAD_RESULT("GetOrderedEdges() failed"); } } if ( reverse ) { reverseEdges( edges2 , nbEInW2.front()); + + if ( SMESH_Algo::isDegenerated( edges2.front() )) + { + // move a degenerated edge to the back of the outer wire + edgeIt = edges2.end(); + if ( severalWires ) { + edgeIt = edges2.begin(); + std::advance( edgeIt, nbEInW2.front() ); + } + edges2.splice( edgeIt, edges2, edges2.begin() ); + } if (( VV1[1].IsSame( TopExp::LastVertex( edges1.front(), true ))) != ( VV2[1].IsSame( TopExp::LastVertex( edges2.front(), true )))) CONT_BAD_RESULT("GetOrderedEdges() failed"); @@ -1399,6 +1467,65 @@ int StdMeshers_ProjectionUtils::FindFaceAssociation(const TopoDS_Face& face1, OK = true; } // loop algos getting an outer wire + + if ( OK && nbEInW1.front() > 4 ) // care of a case where faces are closed (23032) + { + // check if the first edges are seam ones + list< TopoDS_Edge >::iterator revSeam1, revSeam2; + revSeam1 = std::find( ++edges1.begin(), edges1.end(), edges1.front().Reversed()); + revSeam2 = edges2.end(); + if ( revSeam1 != edges1.end() ) + revSeam2 = std::find( ++edges2.begin(), edges2.end(), edges2.front().Reversed()); + if ( revSeam2 != edges2.end() ) // two seams detected + { + bool reverse = + std::distance( edges1.begin(), revSeam1 ) != std::distance( edges2.begin(), revSeam2 ); + if ( !reverse && isClosenessAssoc ) + { + // compare orientations of a non-seam edges using 3D closeness; + // look for a non-seam edges + list< TopoDS_Edge >::iterator edge1 = ++edges1.begin(); + list< TopoDS_Edge >::iterator edge2 = ++edges2.begin(); + for ( ; edge1 != edges1.end(); ++edge1, ++edge2 ) + { + if (( edge1 == revSeam1 ) || + ( SMESH_Algo::isDegenerated( *edge1 )) || + ( std::find( ++edges1.begin(), edges1.end(), edge1->Reversed()) != edges1.end() )) + continue; + gp_Pnt p1 = BRep_Tool::Pnt( VV1[0] ); + gp_Pnt p2 = BRep_Tool::Pnt( VV2[0] ); + gp_Vec vec2to1( p2, p1 ); + + gp_Pnt pp1[2], pp2[2]; + const double r = 0.2345; + double f,l; + Handle(Geom_Curve) C = BRep_Tool::Curve( *edge1, f,l ); + pp1[0] = C->Value( f * r + l * ( 1. - r )); + pp1[1] = C->Value( l * r + f * ( 1. - r )); + if ( edge1->Orientation() == TopAbs_REVERSED ) + std::swap( pp1[0], pp1[1] ); + C = BRep_Tool::Curve( *edge2, f,l ); + if ( C.IsNull() ) return 0; + pp2[0] = C->Value( f * r + l * ( 1. - r )).Translated( vec2to1 ); + pp2[1] = C->Value( l * r + f * ( 1. - r )).Translated( vec2to1 ); + if ( edge2->Orientation() == TopAbs_REVERSED ) + std::swap( pp2[0], pp2[1] ); + + double dist00 = pp1[0].SquareDistance( pp2[0] ); + double dist01 = pp1[0].SquareDistance( pp2[1] ); + reverse = ( dist00 > dist01 ); + break; + } + } + if ( reverse ) // make a seam counterpart be the first + { + list< TopoDS_Edge >::iterator outWireEnd = edges2.begin(); + std::advance( outWireEnd, nbEInW2.front() ); + edges2.splice( outWireEnd, edges2, edges2.begin(), ++revSeam2 ); + reverseEdges( edges2 , nbEInW2.front()); + } + } + } // Try to orient all (if !OK) or only internal wires (issue 0020996) by UV similarity @@ -1407,13 +1534,23 @@ int StdMeshers_ProjectionUtils::FindFaceAssociation(const TopoDS_Face& face1, // Check that Vec(VV1[0],VV1[1]) in 2D on face1 is the same // as Vec(VV2[0],VV2[1]) on face2 double vTol = BRep_Tool::Tolerance( VV1[0] ); - BRepAdaptor_Surface surface1( face1, false ); + BRepAdaptor_Surface surface1( face1, true ); + BRepAdaptor_Surface surface2( face2, true ); + // TODO: use TrsfFinder2D to superpose the faces + gp_Pnt2d v0f1UV( surface1.FirstUParameter(), surface1.FirstVParameter() ); + gp_Pnt2d v0f2UV( surface2.FirstUParameter(), surface2.FirstVParameter() ); + gp_Pnt2d v1f1UV( surface1.LastUParameter(), surface1.LastVParameter() ); + gp_Pnt2d v1f2UV( surface2.LastUParameter(), surface2.LastVParameter() ); double vTolUV = surface1.UResolution( vTol ) + surface1.VResolution( vTol ); // let's be tolerant - gp_Pnt2d v0f1UV = BRep_Tool::Parameters( VV1[0], face1 ); - gp_Pnt2d v0f2UV = BRep_Tool::Parameters( VV2[0], face2 ); - gp_Pnt2d v1f1UV = BRep_Tool::Parameters( VV1[1], face1 ); - gp_Pnt2d v1f2UV = BRep_Tool::Parameters( VV2[1], face2 ); + // VV1[0] = TopExp::FirstVertex( edges1.front(), true ); // ori is important if face is closed + // VV1[1] = TopExp::LastVertex ( edges1.front(), true ); + // VV2[0] = TopExp::FirstVertex( edges2.front(), true ); + // VV2[1] = TopExp::LastVertex ( edges2.front(), true ); + // gp_Pnt2d v0f1UV = BRep_Tool::Parameters( VV1[0], face1 ); + // gp_Pnt2d v0f2UV = BRep_Tool::Parameters( VV2[0], face2 ); + // gp_Pnt2d v1f1UV = BRep_Tool::Parameters( VV1[1], face1 ); + // gp_Pnt2d v1f2UV = BRep_Tool::Parameters( VV2[1], face2 ); gp_Vec2d v01f1Vec( v0f1UV, v1f1UV ); gp_Vec2d v01f2Vec( v0f2UV, v1f2UV ); if ( Abs( v01f1Vec.X()-v01f2Vec.X()) < vTolUV && @@ -1432,7 +1569,6 @@ int StdMeshers_ProjectionUtils::FindFaceAssociation(const TopoDS_Face& face1, list< int >::iterator nbE2, nbE1 = nbEInW1.begin(); list< TopoDS_Edge >::iterator edge2Beg, edge1Beg = edges1.begin(); if ( OK ) std::advance( edge1Beg, *nbE1++ ); - // reach an end of edges of a current wire1 list< TopoDS_Edge >::iterator edge2End, edge1End; // // find corresponding wires of face2 @@ -1442,8 +1578,10 @@ int StdMeshers_ProjectionUtils::FindFaceAssociation(const TopoDS_Face& face1, edge1End = edge1Beg; std::advance( edge1End, *nbE1 ); // UV on face1 to find on face2 - v0f1UV = BRep_Tool::Parameters( TopExp::FirstVertex(*edge1Beg,true), face1 ); - v1f1UV = BRep_Tool::Parameters( TopExp::LastVertex (*edge1Beg,true), face1 ); + TopoDS_Vertex v01 = SMESH_MesherHelper::IthVertex(0,*edge1Beg); + TopoDS_Vertex v11 = SMESH_MesherHelper::IthVertex(1,*edge1Beg); + v0f1UV = BRep_Tool::Parameters( v01, face1 ); + v1f1UV = BRep_Tool::Parameters( v11, face1 ); v0f1UV.ChangeCoord() += dUV; v1f1UV.ChangeCoord() += dUV; // @@ -1458,19 +1596,42 @@ int StdMeshers_ProjectionUtils::FindFaceAssociation(const TopoDS_Face& face1, std::advance( edge2End, *nbE2 ); if ( *nbE1 == *nbE2 && iW2 >= iW1 ) { - // rotate edge2 untill coincidence with edge1 in 2D + // rotate edge2 until coincides with edge1 in 2D int i = *nbE2; - while ( i-- > 0 && !sameVertexUV( *edge2Beg, face2, 0, v0f1UV, vTolUV )) + bool sameUV = false; + while ( !( sameUV = sameVertexUV( *edge2Beg, face2, 0, v0f1UV, vTolUV )) && --i > 0 ) // move edge2Beg to place before edge2End edges2.splice( edge2End, edges2, edge2Beg++ ); - if ( edge2Beg != edges2.end() && - sameVertexUV( *edge2Beg, face2, 0, v0f1UV, vTolUV )) + if ( sameUV ) { if ( iW1 == 0 ) OK = true; // OK is for the first wire + // reverse edges2 if needed - if ( !sameVertexUV( *edge2Beg, face2, 1, v1f1UV, vTolUV )) - reverseEdges( edges2 , *nbE2, std::distance( edges2.begin(),edge2Beg )); + if ( SMESH_MesherHelper::IsClosedEdge( *edge1Beg )) + { + // Commented (so far?) as it's not checked if orientation must be same or reversed + // double f,l; + // Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( *edge1Beg, face1,f,l ); + // if ( edge1Beg->Orientation() == TopAbs_REVERSED ) + // std::swap( f,l ); + // gp_Pnt2d uv1 = dUV + c1->Value( f * 0.8 + l * 0.2 ).XY(); + + // Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( *edge2Beg, face2,f,l ); + // if ( edge2Beg->Orientation() == TopAbs_REVERSED ) + // std::swap( f,l ); + // gp_Pnt2d uv2 = c2->Value( f * 0.8 + l * 0.2 ); + // gp_Pnt2d uv3 = c2->Value( l * 0.8 + f * 0.2 ); + + // if ( uv1.SquareDistance( uv2 ) > uv1.SquareDistance( uv3 )) + // edge2Beg->Reverse(); + } + else + { + if ( !sameVertexUV( *edge2Beg, face2, 1, v1f1UV, vTolUV )) + reverseEdges( edges2 , *nbE2, std::distance( edges2.begin(),edge2Beg )); + } + // put wire2 at a right place within edges2 if ( iW1 != iW2 ) { list< TopoDS_Edge >::iterator place2 = edges2.begin(); @@ -1484,7 +1645,7 @@ int StdMeshers_ProjectionUtils::FindFaceAssociation(const TopoDS_Face& face1, break; } } - // prepare to the next wire loop + // prepare for the next wire loop edge2Beg = edge2End; } edge1Beg = edge1End; @@ -1717,7 +1878,7 @@ StdMeshers_ProjectionUtils::GetPropagationEdge( SMESH_Mesh* aMes int prevChainSize = aChain.Extent(); if ( aChain.Add(anOppE) > prevChainSize ) { // ... anOppE is not in aChain // Add found edge to the chain oriented so that to - // have it co-directed with a forward MainEdge + // have it co-directed with a fromEdge TopAbs_Orientation ori = anE.Orientation(); if ( anOppE.Orientation() == fourEdges[found].Orientation() ) ori = TopAbs::Reverse( ori ); @@ -1780,7 +1941,7 @@ FindMatchingNodesOnFaces( const TopoDS_Face& face1, helper1.SetSubShape( face1 ); helper2.SetSubShape( face2 ); - if ( helper1.HasSeam() != helper2.HasSeam() ) + if ( helper1.HasRealSeam() != helper2.HasRealSeam() ) RETURN_BAD_RESULT("Different faces' geometry"); // Data to call SMESH_MeshEditor::FindMatchingNodes(): @@ -1796,7 +1957,8 @@ FindMatchingNodesOnFaces( const TopoDS_Face& face1, eE.Next(); // edge 1 if ( !assocMap.IsBound( e2, /*is2nd=*/true )) - RETURN_BAD_RESULT("Association not found for edge " << meshDS2->ShapeToIndex( e2 )); + continue; + //RETURN_BAD_RESULT("Association not found for edge " << meshDS2->ShapeToIndex( e2 )); TopoDS_Edge e1 = TopoDS::Edge( assocMap( e2, /*is2nd=*/true )); if ( !helper1.IsSubShape( e1, face1 )) RETURN_BAD_RESULT("Wrong association, edge " << meshDS1->ShapeToIndex( e1 ) << @@ -1840,7 +2002,11 @@ FindMatchingNodesOnFaces( const TopoDS_Face& face1, // get 2 matching vertices TopoDS_Vertex V2 = TopExp::FirstVertex( TopoDS::Edge( edge2 )); if ( !assocMap.IsBound( V2, /*is2nd=*/true )) - RETURN_BAD_RESULT("Association not found for vertex " << meshDS2->ShapeToIndex( V2 )); + { + V2 = TopExp::LastVertex( TopoDS::Edge( edge2 )); + if ( !assocMap.IsBound( V2, /*is2nd=*/true )) + RETURN_BAD_RESULT("Association not found for vertex " << meshDS2->ShapeToIndex( V2 )); + } TopoDS_Vertex V1 = TopoDS::Vertex( assocMap( V2, /*is2nd=*/true )); // nodes on vertices @@ -1903,97 +2069,109 @@ FindMatchingNodesOnFaces( const TopoDS_Face& face1, // 2. face sets - set Elems1, Elems2; - for ( int is2 = 0; is2 < 2; ++is2 ) + int assocRes; + for ( int iAttempt = 0; iAttempt < 2; ++iAttempt ) { - set & elems = is2 ? Elems2 : Elems1; - SMESHDS_SubMesh* sm = is2 ? SM2 : SM1; - SMESH_MesherHelper* helper = is2 ? &helper2 : &helper1; - const TopoDS_Face & face = is2 ? face2 : face1; - SMDS_ElemIteratorPtr eIt = sm->GetElements(); - - if ( !helper->IsRealSeam( is2 ? edge2 : edge1 )) - { - while ( eIt->more() ) elems.insert( eIt->next() ); - } - else + set Elems1, Elems2; + for ( int is2 = 0; is2 < 2; ++is2 ) { - // the only suitable edge is seam, i.e. it is a sphere. - // FindMatchingNodes() will not know which way to go from any edge. - // So we ignore all faces having nodes on edges or vertices except - // one of faces sharing current start nodes - - // find a face to keep - const SMDS_MeshElement* faceToKeep = 0; - const SMDS_MeshNode* vNode = is2 ? vNode2 : vNode1; - const SMDS_MeshNode* eNode = is2 ? eNode2[0] : eNode1[0]; - TIDSortedElemSet inSet, notInSet; - - const SMDS_MeshElement* f1 = - SMESH_MeshAlgos::FindFaceInSet( vNode, eNode, inSet, notInSet ); - if ( !f1 ) RETURN_BAD_RESULT("The first face on seam not found"); - notInSet.insert( f1 ); - - const SMDS_MeshElement* f2 = - SMESH_MeshAlgos::FindFaceInSet( vNode, eNode, inSet, notInSet ); - if ( !f2 ) RETURN_BAD_RESULT("The second face on seam not found"); - - // select a face with less UV of vNode - const SMDS_MeshNode* notSeamNode[2] = {0, 0}; - for ( int iF = 0; iF < 2; ++iF ) { - const SMDS_MeshElement* f = ( iF ? f2 : f1 ); - for ( int i = 0; !notSeamNode[ iF ] && i < f->NbNodes(); ++i ) { - const SMDS_MeshNode* node = f->GetNode( i ); - if ( !helper->IsSeamShape( node->getshapeId() )) - notSeamNode[ iF ] = node; - } + set & elems = is2 ? Elems2 : Elems1; + SMESHDS_SubMesh* sm = is2 ? SM2 : SM1; + SMESH_MesherHelper* helper = is2 ? &helper2 : &helper1; + const TopoDS_Face & face = is2 ? face2 : face1; + SMDS_ElemIteratorPtr eIt = sm->GetElements(); + + if ( !helper->IsRealSeam( is2 ? edge2 : edge1 )) + { + while ( eIt->more() ) elems.insert( elems.end(), eIt->next() ); } - gp_Pnt2d uv1 = helper->GetNodeUV( face, vNode, notSeamNode[0] ); - gp_Pnt2d uv2 = helper->GetNodeUV( face, vNode, notSeamNode[1] ); - if ( uv1.X() + uv1.Y() > uv2.X() + uv2.Y() ) - faceToKeep = f2; else - faceToKeep = f1; - - // fill elem set - elems.insert( faceToKeep ); - while ( eIt->more() ) { - const SMDS_MeshElement* f = eIt->next(); - int nbNodes = f->NbNodes(); - if ( f->IsQuadratic() ) - nbNodes /= 2; - bool onBnd = false; - for ( int i = 0; !onBnd && i < nbNodes; ++i ) { - const SMDS_MeshNode* node = f->GetNode( i ); - onBnd = ( node->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE); + { + // the only suitable edge is seam, i.e. it is a sphere. + // FindMatchingNodes() will not know which way to go from any edge. + // So we ignore all faces having nodes on edges or vertices except + // one of faces sharing current start nodes + + // find a face to keep + const SMDS_MeshElement* faceToKeep = 0; + const SMDS_MeshNode* vNode = is2 ? vNode2 : vNode1; + const SMDS_MeshNode* eNode = is2 ? eNode2[0] : eNode1[0]; + TIDSortedElemSet inSet, notInSet; + + const SMDS_MeshElement* f1 = + SMESH_MeshAlgos::FindFaceInSet( vNode, eNode, inSet, notInSet ); + if ( !f1 ) RETURN_BAD_RESULT("The first face on seam not found"); + notInSet.insert( f1 ); + + const SMDS_MeshElement* f2 = + SMESH_MeshAlgos::FindFaceInSet( vNode, eNode, inSet, notInSet ); + if ( !f2 ) RETURN_BAD_RESULT("The second face on seam not found"); + + // select a face with less UV of vNode + const SMDS_MeshNode* notSeamNode[2] = {0, 0}; + for ( int iF = 0; iF < 2; ++iF ) { + const SMDS_MeshElement* f = ( iF ? f2 : f1 ); + for ( int i = 0; !notSeamNode[ iF ] && i < f->NbNodes(); ++i ) { + const SMDS_MeshNode* node = f->GetNode( i ); + if ( !helper->IsSeamShape( node->getshapeId() )) + notSeamNode[ iF ] = node; + } } - if ( !onBnd ) - elems.insert( f ); - } - // add also faces adjacent to faceToKeep - int nbNodes = faceToKeep->NbNodes(); - if ( faceToKeep->IsQuadratic() ) nbNodes /= 2; - notInSet.insert( f1 ); - notInSet.insert( f2 ); - for ( int i = 0; i < nbNodes; ++i ) { - const SMDS_MeshNode* n1 = faceToKeep->GetNode( i ); - const SMDS_MeshNode* n2 = faceToKeep->GetNode(( i+1 ) % nbNodes ); - f1 = SMESH_MeshAlgos::FindFaceInSet( n1, n2, inSet, notInSet ); - if ( f1 ) - elems.insert( f1 ); - } - } // case on a sphere - } // loop on 2 faces - - // int quadFactor = (*Elems1.begin())->IsQuadratic() ? 2 : 1; + gp_Pnt2d uv1 = helper->GetNodeUV( face, vNode, notSeamNode[0] ); + gp_Pnt2d uv2 = helper->GetNodeUV( face, vNode, notSeamNode[1] ); + if ( uv1.X() + uv1.Y() > uv2.X() + uv2.Y() ) + faceToKeep = f2; + else + faceToKeep = f1; + + // fill elem set + elems.insert( faceToKeep ); + while ( eIt->more() ) { + const SMDS_MeshElement* f = eIt->next(); + int nbNodes = f->NbNodes(); + if ( f->IsQuadratic() ) + nbNodes /= 2; + bool onBnd = false; + for ( int i = 0; !onBnd && i < nbNodes; ++i ) { + const SMDS_MeshNode* node = f->GetNode( i ); + onBnd = ( node->GetPosition()->GetTypeOfPosition() != SMDS_TOP_FACE); + } + if ( !onBnd ) + elems.insert( f ); + } + // add also faces adjacent to faceToKeep + int nbNodes = faceToKeep->NbNodes(); + if ( faceToKeep->IsQuadratic() ) nbNodes /= 2; + notInSet.insert( f1 ); + notInSet.insert( f2 ); + for ( int i = 0; i < nbNodes; ++i ) { + const SMDS_MeshNode* n1 = faceToKeep->GetNode( i ); + const SMDS_MeshNode* n2 = faceToKeep->GetNode(( i+1 ) % nbNodes ); + f1 = SMESH_MeshAlgos::FindFaceInSet( n1, n2, inSet, notInSet ); + if ( f1 ) + elems.insert( f1 ); + } + } // case on a sphere + } // loop on 2 faces + + node1To2Map.clear(); + assocRes = SMESH_MeshEditor::FindMatchingNodes( Elems1, Elems2, + vNode1, vNode2, + eNode1[0], eNode2[0], + node1To2Map); + if (( assocRes != SMESH_MeshEditor::SEW_OK ) && + ( eNode1[1] || eNode2[1] )) // there is another node to try (on a closed EDGE) + { + node1To2Map.clear(); + if ( eNode1[1] ) std::swap( eNode1[0], eNode1[1] ); + else std::swap( eNode2[0], eNode2[1] ); + continue; // one more attempt + } - node1To2Map.clear(); - int res = SMESH_MeshEditor::FindMatchingNodes( Elems1, Elems2, - vNode1, vNode2, - eNode1[0], eNode2[0], - node1To2Map); - if ( res != SMESH_MeshEditor::SEW_OK ) - RETURN_BAD_RESULT("FindMatchingNodes() result " << res ); + break; + } + if ( assocRes != SMESH_MeshEditor::SEW_OK ) + RETURN_BAD_RESULT("FindMatchingNodes() result " << assocRes ); // On a sphere, add matching nodes on the edge @@ -2015,7 +2193,7 @@ FindMatchingNodesOnFaces( const TopoDS_Face& face1, static_cast(node->GetPosition()); pos2nodes.insert( make_pair( pos->GetUParameter(), node )); } - if ( pos2nodes.size() != edgeSM->NbNodes() ) + if ((int) pos2nodes.size() != edgeSM->NbNodes() ) RETURN_BAD_RESULT("Equal params of nodes on edge " << smDS->ShapeToIndex( edge ) << " of face " << is2 ); } @@ -2121,7 +2299,7 @@ bool StdMeshers_ProjectionUtils::MakeComputed(SMESH_subMesh * sm, const int iter mesh->GetHypotheses( shape, hypoFilter, hyps, true, &assignedTo ); if ( nbAlgos > 1 ) // concurrent algos { - list smList; // where an algo is assigned + vector smList; // where an algo is assigned list< TopoDS_Shape >::iterator shapeIt = assignedTo.begin(); for ( ; shapeIt != assignedTo.end(); ++shapeIt ) smList.push_back( mesh->GetSubMesh( *shapeIt )); @@ -2153,7 +2331,7 @@ bool StdMeshers_ProjectionUtils::MakeComputed(SMESH_subMesh * sm, const int iter string algoType = algo->GetName(); if ( algoType.substr(0, 11) != "Projection_") - return gen->Compute( *mesh, shape, /*shapeOnly=*/true ); + return gen->Compute( *mesh, shape, SMESH_Gen::SHAPE_ONLY ); // try to compute source mesh @@ -2194,7 +2372,7 @@ bool StdMeshers_ProjectionUtils::MakeComputed(SMESH_subMesh * sm, const int iter srcMesh = mesh; if ( MakeComputed( srcMesh->GetSubMesh( srcShape ), iterationNb + 1 ) && - gen->Compute( *mesh, shape, /*shapeOnly=*/true )) + gen->Compute( *mesh, shape, SMESH_Gen::SHAPE_ONLY )) return sm->IsMeshComputed(); return false; @@ -2216,7 +2394,7 @@ std::string StdMeshers_ProjectionUtils::SourceNotComputedError( SMESH_subMesh * if ( !sm || sm->GetAlgoState() != SMESH_subMesh::NO_ALGO ) return usualMessage; // algo is OK, anything else is KO. - // Try to find a type of all-dimentional algorithm that would compute the + // Try to find a type of all-dimensional algorithm that would compute the // given sub-mesh if it could be launched before projection const TopoDS_Shape shape = sm->GetSubShape(); const int shapeDim = SMESH_Gen::GetShapeDim( shape ); @@ -2372,10 +2550,570 @@ void StdMeshers_ProjectionUtils::SetEventListener(SMESH_subMesh* subMesh, } else { - subMesh->SetEventListener( getSrcSubMeshListener(), - SMESH_subMeshEventListenerData::MakeData( subMesh ), - srcShapeSM ); + if ( SMESH_subMeshEventListenerData* data = + srcShapeSM->GetEventListenerData( getSrcSubMeshListener() )) + { + bool alreadyIn = + (std::find( data->mySubMeshes.begin(), + data->mySubMeshes.end(), subMesh ) != data->mySubMeshes.end() ); + if ( !alreadyIn ) + data->mySubMeshes.push_back( subMesh ); + } + else + { + subMesh->SetEventListener( getSrcSubMeshListener(), + SMESH_subMeshEventListenerData::MakeData( subMesh ), + srcShapeSM ); + } } } } } + +namespace StdMeshers_ProjectionUtils +{ + + //================================================================================ + /*! + * \brief Computes transformation between two sets of 2D points using + * a least square approximation + * + * See "Surface Mesh Projection For Hexahedral Mesh Generation By Sweeping" + * by X.Roca, J.Sarrate, A.Huerta. (2.2) + */ + //================================================================================ + + bool TrsfFinder2D::Solve( const vector< gp_XY >& srcPnts, + const vector< gp_XY >& tgtPnts ) + { + // find gravity centers + gp_XY srcGC( 0,0 ), tgtGC( 0,0 ); + for ( size_t i = 0; i < srcPnts.size(); ++i ) + { + srcGC += srcPnts[i]; + tgtGC += tgtPnts[i]; + } + srcGC /= srcPnts.size(); + tgtGC /= tgtPnts.size(); + + // find trsf + + math_Matrix mat (1,4,1,4, 0.); + math_Vector vec (1,4, 0.); + + // cout << "m1 = smesh.Mesh('src')" << endl + // << "m2 = smesh.Mesh('tgt')" << endl; + double xx = 0, xy = 0, yy = 0; + for ( size_t i = 0; i < srcPnts.size(); ++i ) + { + gp_XY srcUV = srcPnts[i] - srcGC; + gp_XY tgtUV = tgtPnts[i] - tgtGC; + xx += srcUV.X() * srcUV.X(); + yy += srcUV.Y() * srcUV.Y(); + xy += srcUV.X() * srcUV.Y(); + vec( 1 ) += srcUV.X() * tgtUV.X(); + vec( 2 ) += srcUV.Y() * tgtUV.X(); + vec( 3 ) += srcUV.X() * tgtUV.Y(); + vec( 4 ) += srcUV.Y() * tgtUV.Y(); + // cout << "m1.AddNode( " << srcUV.X() << ", " << srcUV.Y() << ", 0 )" << endl + // << "m2.AddNode( " << tgtUV.X() << ", " << tgtUV.Y() << ", 0 )" << endl; + } + mat( 1,1 ) = mat( 3,3 ) = xx; + mat( 2,2 ) = mat( 4,4 ) = yy; + mat( 1,2 ) = mat( 2,1 ) = mat( 3,4 ) = mat( 4,3 ) = xy; + + math_Gauss solver( mat ); + if ( !solver.IsDone() ) + return false; + solver.Solve( vec ); + if ( vec.Norm2() < gp::Resolution() ) + return false; + // cout << vec( 1 ) << "\t " << vec( 2 ) << endl + // << vec( 3 ) << "\t " << vec( 4 ) << endl; + + _trsf.SetTranslationPart( tgtGC ); + _srcOrig = srcGC; + + gp_Mat2d& M = const_cast< gp_Mat2d& >( _trsf.VectorialPart()); + M( 1,1 ) = vec( 1 ); + M( 2,1 ) = vec( 2 ); // | 1 3 | -- is it correct ???????? + M( 1,2 ) = vec( 3 ); // | 2 4 | + M( 2,2 ) = vec( 4 ); + + return true; + } + + //================================================================================ + /*! + * \brief Transforms a 2D points using a found transformation + */ + //================================================================================ + + gp_XY TrsfFinder2D::Transform( const gp_Pnt2d& srcUV ) const + { + gp_XY uv = srcUV.XY() - _srcOrig ; + _trsf.Transforms( uv ); + return uv; + } + + //================================================================================ + /*! + * \brief Computes transformation between two sets of 3D points using + * a least square approximation + * + * See "Surface Mesh Projection For Hexahedral Mesh Generation By Sweeping" + * by X.Roca, J.Sarrate, A.Huerta. (2.4) + */ + //================================================================================ + + bool TrsfFinder3D::Solve( const vector< gp_XYZ > & srcPnts, + const vector< gp_XYZ > & tgtPnts ) + { + // find gravity center + gp_XYZ srcGC( 0,0,0 ), tgtGC( 0,0,0 ); + for ( size_t i = 0; i < srcPnts.size(); ++i ) + { + srcGC += srcPnts[i]; + tgtGC += tgtPnts[i]; + } + srcGC /= srcPnts.size(); + tgtGC /= tgtPnts.size(); + + gp_XYZ srcOrig = 2 * srcGC - tgtGC; + gp_XYZ tgtOrig = srcGC; + + // find trsf + + math_Matrix mat (1,9,1,9, 0.); + math_Vector vec (1,9, 0.); + + double xx = 0, yy = 0, zz = 0; + double xy = 0, xz = 0, yz = 0; + for ( size_t i = 0; i < srcPnts.size(); ++i ) + { + gp_XYZ src = srcPnts[i] - srcOrig; + gp_XYZ tgt = tgtPnts[i] - tgtOrig; + xx += src.X() * src.X(); + yy += src.Y() * src.Y(); + zz += src.Z() * src.Z(); + xy += src.X() * src.Y(); + xz += src.X() * src.Z(); + yz += src.Y() * src.Z(); + vec( 1 ) += src.X() * tgt.X(); + vec( 2 ) += src.Y() * tgt.X(); + vec( 3 ) += src.Z() * tgt.X(); + vec( 4 ) += src.X() * tgt.Y(); + vec( 5 ) += src.Y() * tgt.Y(); + vec( 6 ) += src.Z() * tgt.Y(); + vec( 7 ) += src.X() * tgt.Z(); + vec( 8 ) += src.Y() * tgt.Z(); + vec( 9 ) += src.Z() * tgt.Z(); + } + mat( 1,1 ) = mat( 4,4 ) = mat( 7,7 ) = xx; + mat( 2,2 ) = mat( 5,5 ) = mat( 8,8 ) = yy; + mat( 3,3 ) = mat( 6,6 ) = mat( 9,9 ) = zz; + mat( 1,2 ) = mat( 2,1 ) = mat( 4,5 ) = mat( 5,4 ) = mat( 7,8 ) = mat( 8,7 ) = xy; + mat( 1,3 ) = mat( 3,1 ) = mat( 4,6 ) = mat( 6,4 ) = mat( 7,9 ) = mat( 9,7 ) = xz; + mat( 2,3 ) = mat( 3,2 ) = mat( 5,6 ) = mat( 6,5 ) = mat( 8,9 ) = mat( 9,8 ) = yz; + + math_Gauss solver( mat ); + if ( !solver.IsDone() ) + return false; + solver.Solve( vec ); + if ( vec.Norm2() < gp::Resolution() ) + return false; + // cout << endl + // << vec( 1 ) << "\t " << vec( 2 ) << "\t " << vec( 3 ) << endl + // << vec( 4 ) << "\t " << vec( 5 ) << "\t " << vec( 6 ) << endl + // << vec( 7 ) << "\t " << vec( 8 ) << "\t " << vec( 9 ) << endl; + + _srcOrig = srcOrig; + _trsf.SetTranslationPart( tgtOrig ); + + gp_Mat& M = const_cast< gp_Mat& >( _trsf.VectorialPart() ); + M.SetRows( gp_XYZ( vec( 1 ), vec( 2 ), vec( 3 )), + gp_XYZ( vec( 4 ), vec( 5 ), vec( 6 )), + gp_XYZ( vec( 7 ), vec( 8 ), vec( 9 ))); + return true; + } + + //================================================================================ + /*! + * \brief Transforms a 3D point using a found transformation + */ + //================================================================================ + + gp_XYZ TrsfFinder3D::Transform( const gp_Pnt& srcP ) const + { + gp_XYZ p = srcP.XYZ() - _srcOrig; + _trsf.Transforms( p ); + return p; + } + + //================================================================================ + /*! + * \brief Transforms a 3D vector using a found transformation + */ + //================================================================================ + + gp_XYZ TrsfFinder3D::TransformVec( const gp_Vec& v ) const + { + return v.XYZ().Multiplied( _trsf.VectorialPart() ); + } + //================================================================================ + /*! + * \brief Inversion + */ + //================================================================================ + + bool TrsfFinder3D::Invert() + { + if (( _trsf.Form() == gp_Translation ) && + ( _srcOrig.X() != 0 || _srcOrig.Y() != 0 || _srcOrig.Z() != 0 )) + { + // seems to be defined via Solve() + gp_XYZ newSrcOrig = _trsf.TranslationPart(); + gp_Mat& M = const_cast< gp_Mat& >( _trsf.VectorialPart() ); + const double D = M.Determinant(); + if ( D < 1e-3 * ( newSrcOrig - _srcOrig ).Modulus() ) + { +#ifdef _DEBUG_ + cerr << "TrsfFinder3D::Invert()" + << "D " << M.Determinant() << " IsSingular " << M.IsSingular() << endl; +#endif + return false; + } + gp_Mat Minv = M.Inverted(); + _trsf.SetTranslationPart( _srcOrig ); + _srcOrig = newSrcOrig; + M = Minv; + } + else + { + _trsf.Invert(); + } + return true; + } + + //================================================================================ + /*! + * \brief Add in-FACE nodes surrounding a given node to a queue + */ + //================================================================================ + + void Morph::AddCloseNodes( const SMDS_MeshNode* srcNode, + const BRepMesh_Triangle* bmTria, + const int srcFaceID, + TNodeTriaList & noTriQueue ) + { + // find in-FACE nodes + SMDS_ElemIteratorPtr elems = srcNode->GetInverseElementIterator(SMDSAbs_Face); + while ( elems->more() ) + { + const SMDS_MeshElement* elem = elems->next(); + if ( elem->getshapeId() == srcFaceID ) + { + for ( int i = 0, nb = elem->NbNodes(); i < nb; ++i ) + { + const SMDS_MeshNode* n = elem->GetNode( i ); + if ( !n->isMarked() /*&& n->getshapeId() == srcFaceID*/ ) + noTriQueue.push_back( make_pair( n, bmTria )); + } + } + } + } + + //================================================================================ + /*! + * \brief Find a delauney triangle containing a given 2D point and return + * barycentric coordinates within the found triangle + */ + //================================================================================ + + const BRepMesh_Triangle* Morph::FindTriangle( const gp_XY& uv, + const BRepMesh_Triangle* bmTria, + double bc[3], + int triaNodes[3] ) + { + int nodeIDs[3]; + gp_XY nodeUVs[3]; + int linkIDs[3]; + Standard_Boolean ori[3]; + + while ( bmTria ) + { + // check bmTria + + _triaDS->ElementNodes( *bmTria, nodeIDs ); + nodeUVs[0] = _triaDS->GetNode( nodeIDs[0] ).Coord(); + nodeUVs[1] = _triaDS->GetNode( nodeIDs[1] ).Coord(); + nodeUVs[2] = _triaDS->GetNode( nodeIDs[2] ).Coord(); + + SMESH_MeshAlgos::GetBarycentricCoords( uv, + nodeUVs[0], nodeUVs[1], nodeUVs[2], + bc[0], bc[1] ); + if ( bc[0] >= 0 && bc[1] >= 0 && bc[0] + bc[1] <= 1 ) + { + bc[2] = 1 - bc[0] - bc[1]; + triaNodes[0] = nodeIDs[0]; + triaNodes[1] = nodeIDs[1]; + triaNodes[2] = nodeIDs[2]; + return bmTria; + } + + // look for a neighbor triangle, which is adjacent to a link intersected + // by a segment( triangle center -> uv ) + + gp_XY gc = ( nodeUVs[0] + nodeUVs[1] + nodeUVs[2] ) / 3.; + gp_XY seg = uv - gc; + + bmTria->Edges( linkIDs, ori ); + int triaID = _triaDS->IndexOf( *bmTria ); + bmTria = 0; + + for ( int i = 0; i < 3; ++i ) + { + const BRepMesh_PairOfIndex & triIDs = _triaDS->ElementsConnectedTo( linkIDs[i] ); + if ( triIDs.Extent() < 2 ) + continue; // no neighbor triangle + + // check if a link intersects gc2uv + const BRepMesh_Edge & link = _triaDS->GetLink( linkIDs[i] ); + const BRepMesh_Vertex & n1 = _triaDS->GetNode( link.FirstNode() ); + const BRepMesh_Vertex & n2 = _triaDS->GetNode( link.LastNode() ); + gp_XY uv1 = n1.Coord(); + gp_XY lin = n2.Coord() - uv1; // link direction + + double crossSegLin = seg ^ lin; + if ( Abs( crossSegLin ) < std::numeric_limits::min() ) + continue; // parallel + + double uSeg = ( uv1 - gc ) ^ lin / crossSegLin; + if ( 0. <= uSeg && uSeg <= 1. ) + { + bmTria = & _triaDS->GetElement( triIDs.Index( 1 + ( triIDs.Index(1) == triaID ))); + break; + } + } + } + return bmTria; + } + + //================================================================================ + /*! + * \brief Return a triangle sharing a given boundary node + * \param [in] iBndNode - index of the boundary node + * \return const BRepMesh_Triangle* - a found triangle + */ + //================================================================================ + + const BRepMesh_Triangle* Morph::GetTriangleNear( int iBndNode ) + { + const BRepMesh::ListOfInteger & linkIds = _triaDS->LinksConnectedTo( iBndNode ); + const BRepMesh_PairOfIndex & triaIds = _triaDS->ElementsConnectedTo( linkIds.First() ); + const BRepMesh_Triangle& tria = _triaDS->GetElement( triaIds.Index(1) ); + return &tria; + } + + //================================================================================ + /*! + * \brief triangulate the srcFace in 2D + * \param [in] srcWires - boundary of the src FACE + */ + //================================================================================ + + Morph::Morph(const TSideVector& srcWires) + { + _srcSubMesh = srcWires[0]->GetMesh()->GetSubMesh( srcWires[0]->Face() ); + + // compute _scale + { + BRepAdaptor_Surface surf( srcWires[0]->Face() ); + const int nbDiv = 100; + const double uRange = surf.LastUParameter() - surf.FirstUParameter(); + const double vRange = surf.LastVParameter() - surf.FirstVParameter(); + const double dU = uRange / nbDiv; + const double dV = vRange / nbDiv; + double u = surf.FirstUParameter(), v = surf.FirstVParameter(); + gp_Pnt p0U = surf.Value( u, v ), p0V = p0U; + double lenU = 0, lenV = 0; + for ( ; u < surf.LastUParameter(); u += dU, v += dV ) + { + gp_Pnt p1U = surf.Value( u, surf.FirstVParameter() ); + lenU += p1U.Distance( p0U ); + p0U = p1U; + gp_Pnt p1V = surf.Value( surf.FirstUParameter(), v ); + lenV += p1V.Distance( p0V ); + p0V = p1V; + } + _scale.SetCoord( lenU / uRange, lenV / vRange ); + } + + // count boundary points + int iP = 1, nbP = 0; + for ( size_t iW = 0; iW < srcWires.size(); ++iW ) + nbP += srcWires[iW]->NbPoints() - 1; // 1st and last points coincide + + _bndSrcNodes.resize( nbP + 1 ); _bndSrcNodes[0] = 0; + + // fill boundary points + BRepMesh::Array1OfVertexOfDelaun srcVert( 1, 1 + nbP ); + BRepMesh_Vertex v( 0, 0, BRepMesh_Frontier ); + for ( size_t iW = 0; iW < srcWires.size(); ++iW ) + { + const UVPtStructVec& srcPnt = srcWires[iW]->GetUVPtStruct(); + for ( int i = 0, nb = srcPnt.size() - 1; i < nb; ++i, ++iP ) + { + _bndSrcNodes[ iP ] = srcPnt[i].node; + srcPnt[i].node->setIsMarked( true ); + + v.ChangeCoord() = srcPnt[i].UV().Multiplied( _scale ); + srcVert( iP ) = v; + } + } + // triangulate the srcFace in 2D + BRepMesh_Delaun delauney( srcVert ); + _triaDS = delauney.Result(); + } + + //================================================================================ + /*! + * \brief Move non-marked target nodes + * \param [in,out] tgtHelper - helper + * \param [in] tgtWires - boundary nodes of the target FACE; must be in the + * same order as the nodes in srcWires given in the constructor + * \param [in] src2tgtNodes - map of src -> tgt nodes + * \param [in] moveAll - to move all nodes; if \c false, move only non-marked nodes + * \return bool - Ok or not + */ + //================================================================================ + + bool Morph::Perform(SMESH_MesherHelper& tgtHelper, + const TSideVector& tgtWires, + Handle(ShapeAnalysis_Surface) tgtSurface, + const TNodeNodeMap& src2tgtNodes, + const bool moveAll) + { + // get tgt boundary points corresponding to _bndSrcNodes + size_t nbP = 0; + for ( size_t iW = 0; iW < tgtWires.size(); ++iW ) + nbP += tgtWires[iW]->NbPoints() - 1; // 1st and last points coincide + if ( nbP != _bndSrcNodes.size() - 1 ) + return false; + + BRepMesh::Array1OfVertexOfDelaun tgtVert( 1, 1 + nbP ); + BRepMesh_Vertex v( 0, 0, BRepMesh_Frontier ); + for ( size_t iW = 0, iP = 1; iW < tgtWires.size(); ++iW ) + { + const UVPtStructVec& tgtPnt = tgtWires[iW]->GetUVPtStruct(); + for ( int i = 0, nb = tgtPnt.size() - 1; i < nb; ++i, ++iP ) + { + v.ChangeCoord() = tgtPnt[i].UV().Multiplied( _scale ); + tgtVert( iP ) = v; + } + } + + const TopoDS_Face& srcFace = TopoDS::Face( _srcSubMesh->GetSubShape() ); + const int srcFaceID = _srcSubMesh->GetId(); + SMESHDS_Mesh* tgtMesh = tgtHelper.GetMeshDS(); + const SMDS_MeshNode *srcNode, *tgtNode; + const BRepMesh_Triangle *bmTria; + + // un-mark internal src nodes; later we will mark moved nodes + int nbSrcNodes = 0; + SMDS_NodeIteratorPtr nIt = _srcSubMesh->GetSubMeshDS()->GetNodes(); + if ( !nIt || !nIt->more() ) return true; + if ( moveAll ) + { + nbSrcNodes = _srcSubMesh->GetSubMeshDS()->NbNodes(); + while ( nIt->more() ) + nIt->next()->setIsMarked( false ); + } + else + { + while ( nIt->more() ) + nbSrcNodes += int( !nIt->next()->isMarked() ); + } + + // Move tgt nodes + + double bc[3]; // barycentric coordinates + int nodeIDs[3]; + bool checkUV = true; + const SMDS_FacePosition* pos; + + // a queue of nodes with starting triangles + TNodeTriaList noTriQueue; + size_t iBndSrcN = 1; + + while ( nbSrcNodes > 0 ) + { + while ( !noTriQueue.empty() ) + { + srcNode = noTriQueue.front().first; + bmTria = noTriQueue.front().second; + noTriQueue.pop_front(); + if ( srcNode->isMarked() ) + continue; + --nbSrcNodes; + srcNode->setIsMarked( true ); + + // find a delauney triangle containing the src node + gp_XY uv = tgtHelper.GetNodeUV( srcFace, srcNode, NULL, &checkUV ); + uv *= _scale; + bmTria = FindTriangle( uv, bmTria, bc, nodeIDs ); + if ( !bmTria ) + continue; + + // compute new coordinates for a corresponding tgt node + gp_XY uvNew( 0., 0. ), nodeUV; + for ( int i = 0; i < 3; ++i ) + uvNew += bc[i] * tgtVert( nodeIDs[i]).Coord(); + uvNew.SetCoord( uvNew.X() / _scale.X(), + uvNew.Y() / _scale.Y() ); + gp_Pnt xyz = tgtSurface->Value( uvNew ); + + // find and move tgt node + TNodeNodeMap::const_iterator n2n = src2tgtNodes.find( srcNode ); + if ( n2n == src2tgtNodes.end() ) continue; + tgtNode = n2n->second; + tgtMesh->MoveNode( tgtNode, xyz.X(), xyz.Y(), xyz.Z() ); + + if (( pos = dynamic_cast< const SMDS_FacePosition* >( tgtNode->GetPosition() ))) + const_cast( pos )->SetParameters( uvNew.X(), uvNew.Y() ); + + AddCloseNodes( srcNode, bmTria, srcFaceID, noTriQueue ); + } + + if ( nbSrcNodes > 0 ) + { + // assure that all src nodes are visited + for ( ; iBndSrcN < _bndSrcNodes.size() && noTriQueue.empty(); ++iBndSrcN ) + { + const BRepMesh_Triangle* tria = GetTriangleNear( iBndSrcN ); + AddCloseNodes( _bndSrcNodes[ iBndSrcN ], tria, srcFaceID, noTriQueue ); + } + if ( noTriQueue.empty() ) + { + SMDS_NodeIteratorPtr nIt = _srcSubMesh->GetSubMeshDS()->GetNodes(); + while ( nIt->more() ) + { + srcNode = nIt->next(); + if ( !srcNode->isMarked() ) + noTriQueue.push_back( make_pair( srcNode, bmTria )); + } + } + } + } + + return true; + + } // Morph::Perform + + gp_XY Morph::GetBndUV(const int iNode) const + { + return _triaDS->GetNode( iNode ).Coord(); + } + + +} // namespace StdMeshers_ProjectionUtils