X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Projection_2D.cxx;h=494989655917c4b505ad4df7df72da684b394554;hp=4ef3b3f671d6ed10964de5512a9f23231ff6bfe3;hb=HEAD;hpb=665d037f93371114bf4b00bf11b0f95be418fb77 diff --git a/src/StdMeshers/StdMeshers_Projection_2D.cxx b/src/StdMeshers/StdMeshers_Projection_2D.cxx index 4ef3b3f67..b73b3a0fe 100644 --- a/src/StdMeshers/StdMeshers_Projection_2D.cxx +++ b/src/StdMeshers/StdMeshers_Projection_2D.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2020 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2024 CEA, EDF, OPEN CASCADE // // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS @@ -28,25 +28,30 @@ // #include "StdMeshers_Projection_2D.hxx" +#include "StdMeshers_FaceSide.hxx" +#include "StdMeshers_NumberOfSegments.hxx" #include "StdMeshers_ProjectionSource2D.hxx" #include "StdMeshers_ProjectionUtils.hxx" -#include "StdMeshers_FaceSide.hxx" - -#include "SMDS_EdgePosition.hxx" -#include "SMDS_FacePosition.hxx" -#include "SMESHDS_Hypothesis.hxx" -#include "SMESHDS_Mesh.hxx" -#include "SMESHDS_SubMesh.hxx" -#include "SMESH_Block.hxx" -#include "SMESH_Comment.hxx" -#include "SMESH_Gen.hxx" -#include "SMESH_Mesh.hxx" -#include "SMESH_MeshAlgos.hxx" -#include "SMESH_MeshEditor.hxx" -#include "SMESH_MesherHelper.hxx" -#include "SMESH_Pattern.hxx" -#include "SMESH_subMesh.hxx" -#include "SMESH_subMeshEventListener.hxx" +#include "StdMeshers_Quadrangle_2D.hxx" +#include "StdMeshers_Regular_1D.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include @@ -55,8 +60,21 @@ #include #include #include +#include #include +#include +#include + +#include + +#if OCC_VERSION_LARGE < 0x07070000 +#include +#include +#endif + #include +#include +#include #include #include #include @@ -81,9 +99,15 @@ using namespace std; namespace TAssocTool = StdMeshers_ProjectionUtils; //typedef StdMeshers_ProjectionUtils TAssocTool; +// allow range iteration on NCollection_IndexedMap +template < class IMAP > +typename IMAP::const_iterator begin( const IMAP & m ) { return m.cbegin(); } +template < class IMAP > +typename IMAP::const_iterator end( const IMAP & m ) { return m.cend(); } + //======================================================================= //function : StdMeshers_Projection_2D -//purpose : +//purpose : //======================================================================= StdMeshers_Projection_2D::StdMeshers_Projection_2D(int hypId, SMESH_Gen* gen) @@ -105,7 +129,7 @@ StdMeshers_Projection_2D::~StdMeshers_Projection_2D() //======================================================================= //function : CheckHypothesis -//purpose : +//purpose : //======================================================================= bool StdMeshers_Projection_2D::CheckHypothesis(SMESH_Mesh& theMesh, @@ -301,7 +325,7 @@ namespace { break; } case TopAbs_EDGE: { - + // Get submeshes of sub-vertices const map< int, SMESH_subMesh * >& subSM = sm->DependsOn(); if ( subSM.size() != 2 ) @@ -688,7 +712,7 @@ namespace { while ( elemIt->more() ) // loop on all mesh faces on srcFace { const SMDS_MeshElement* elem = elemIt->next(); - const int nbN = elem->NbCornerNodes(); + const int nbN = elem->NbCornerNodes(); tgtNodes.resize( nbN ); helper->SetElementsOnShape( false ); for ( int i = 0; i < nbN; ++i ) // loop on nodes of the source element @@ -762,7 +786,7 @@ namespace { // check node positions - if ( !tgtFace.IsPartner( srcFace ) ) + // if ( !tgtFace.IsPartner( srcFace ) ) for NETGEN 6 which sets wrong UV { helper->ToFixNodeParameters( true ); @@ -967,29 +991,382 @@ namespace { } // loop on all mesh faces on srcFace return true; + + } // projectBy2DSimilarity() + + //================================================================================ + /*! + * \brief Quadrangle algorithm computing structured triangle mesh + */ + //================================================================================ + + struct QuadAlgo : public StdMeshers_Quadrangle_2D + { + QuadAlgo( int hypId, SMESH_Gen* gen ): StdMeshers_Quadrangle_2D( hypId, gen ) {} + + bool Compute( SMESH_MesherHelper & theHelper, const StdMeshers_FaceSidePtr& theWire ) + { + SMESH_Mesh& theMesh = *theHelper.GetMesh(); + + // set sides of a quad FACE + FaceQuadStruct::Ptr quad( new FaceQuadStruct ); + quad->side.reserve( 4 ); + quad->face = theWire->Face(); + for ( int i = 0; i < 4; ++i ) + quad->side.push_back + ( StdMeshers_FaceSide::New( quad->face, theWire->Edge(i), &theMesh, i < QUAD_TOP_SIDE, + /*skipMedium=*/true, theWire->FaceHelper() )); + if ( !setNormalizedGrid( quad )) + return false; + + // make internal nodes + SMESHDS_Mesh * meshDS = theMesh.GetMeshDS(); + int geomFaceID = meshDS->ShapeToIndex( quad->face ); + Handle(Geom_Surface) S = BRep_Tool::Surface( quad->face ); + for ( int i = 1; i < quad->iSize - 1; i++) + for ( int j = 1; j < quad->jSize - 1; j++) + { + UVPtStruct& uvPnt = quad->UVPt( i, j ); + gp_Pnt P = S->Value( uvPnt.u, uvPnt.v ); + uvPnt.node = meshDS->AddNode(P.X(), P.Y(), P.Z()); + meshDS->SetNodeOnFace( uvPnt.node, geomFaceID, uvPnt.u, uvPnt.v ); + } + + // make triangles + for ( int i = 0; i < quad->iSize-1; i++) { + for ( int j = 0; j < quad->jSize-1; j++) + { + const SMDS_MeshNode* a = quad->uv_grid[ j * quad->iSize + i ].node; + const SMDS_MeshNode* b = quad->uv_grid[ j * quad->iSize + i + 1].node; + const SMDS_MeshNode* c = quad->uv_grid[(j + 1) * quad->iSize + i + 1].node; + const SMDS_MeshNode* d = quad->uv_grid[(j + 1) * quad->iSize + i ].node; + theHelper.AddFace(a, b, c); + theHelper.AddFace(a, c, d); + } + } + return true; + } + }; + + //================================================================================ + /*! + * \brief Local coordinate system of a triangle. Return barycentric coordinates of a point + */ + //================================================================================ + + struct TriaCoordSys + { + gp_Pnt myO; //!< origin + gp_Vec myX; //!< X axis + gp_Vec myY; //!< Y axis + gp_XY myUV1; //!< UV of 2nd node in this CS + gp_XY myUV2; //!< UV of 3d node in this CS + + void Init( const gp_Pnt p1, const gp_Pnt p2, const gp_Pnt p3 ) + { + myO = p1; + + myX = gp_Vec( p1, p2 ); + myUV1.SetCoord( myX.Magnitude(), 0 ); + myX /= myUV1.X(); + + gp_Vec v13( p1, p3 ); + myY = myX.CrossCrossed( myX, v13 ); + myY.Normalize(); + myUV2.SetCoord( v13 * myX, v13 * myY ); + + return; + } + + void GetBaryCoords( const gp_Pnt p, double& bc1, double& bc2, double& bc3 ) const + { + gp_Vec op( myO, p ); + gp_XY uv( op * myX, op * myY ); + + SMESH_MeshAlgos::GetBarycentricCoords( uv, + gp::Origin2d().XY(), myUV1, myUV2, + bc1, bc2 ); + bc3 = 1 - bc1 - bc2; + } + }; + + + //================================================================================ + /*! + * \brief Structured 2D mesh of a quadrilateral FACE; is used in projectQuads() + */ + //================================================================================ + + struct QuadMesh : public SMESH_SequentialMesh + { + ObjectPool< TriaCoordSys > _traiLCSPool; + SMESH_ElementSearcher* _elemSearcher; + SMESH_Gen _sgen; + SMESH_MesherHelper _helper; + + QuadMesh(const TopoDS_Face& face): + _elemSearcher( nullptr ), _helper( *this ) + { + _meshDS = new SMESHDS_Mesh( 0, true ); + _gen = &_sgen; + ShapeToMesh( face ); + } + ~QuadMesh() { delete _elemSearcher; } + + // -------------------------------------------------------------------------------- + /*! + * \brief Compute quadrangle mesh and prepare for face search + */ + bool Compute( const TSideVector& wires, int nbSeg1, int nbSeg2, bool isSourceMesh ) + { + if ( wires.size() > 1 || wires[0]->NbEdges() != 4 ) + return false; + + // compute quadrangle mesh + + SMESH_Hypothesis* algo1D = new StdMeshers_Regular_1D( _sgen.GetANewId(), &_sgen ); + AddHypothesis( GetShapeToMesh(), algo1D->GetID() ); + + StdMeshers_NumberOfSegments * nbHyp1, *nbHyp2; + nbHyp1 = new StdMeshers_NumberOfSegments( _sgen.GetANewId(), &_sgen ); + nbHyp1->SetNumberOfSegments( nbSeg1 ); + AddHypothesis( wires[0]->Edge(0), nbHyp1->GetID() ); + AddHypothesis( wires[0]->Edge(2), nbHyp1->GetID() ); + + nbHyp2 = new StdMeshers_NumberOfSegments( _sgen.GetANewId(), &_sgen ); + nbHyp2->SetNumberOfSegments( nbSeg2 ); + AddHypothesis( wires[0]->Edge(1), nbHyp2->GetID() ); + AddHypothesis( wires[0]->Edge(3), nbHyp2->GetID() ); + + if ( !_sgen.Compute( *this, GetShapeToMesh(), SMESH_Gen::SHAPE_ONLY_UPWARD )) + return false; + + QuadAlgo algo2D( _sgen.GetANewId(), &_sgen ); + if ( !algo2D.Compute( _helper, wires[0] )) + return false; + + // remove edges + // for ( SMDS_ElemIteratorPtr eIt = _meshDS->elementsIterator( SMDSAbs_Edge ); eIt->more(); ) + // _meshDS->RemoveFreeElement( eIt->next(), /*sm=*/0, /*groups=*/false ); + + // _meshDS->Modified(); // setMyModified(); + // _meshDS->CompactMesh(); + + // create TriaCoordSys for every triangle + if ( isSourceMesh ) + { + for ( SMDS_ElemIteratorPtr fIt = _meshDS->elementsIterator( SMDSAbs_Face ); fIt->more(); ) + { + const SMDS_MeshElement* tria = fIt->next(); + TriaCoordSys* triaLCS = _traiLCSPool.getNew(); + triaLCS->Init( SMESH_NodeXYZ( tria->GetNode( 0 )), + SMESH_NodeXYZ( tria->GetNode( 1 )), + SMESH_NodeXYZ( tria->GetNode( 2 ))); + // int i= tria->GetID() - NbEdges() - 1; + // cout << "ID from TRIA " << i << " - poolSize " << _traiLCSPool.nbElements() << + // ( _traiLCSPool[i]!= triaLCS ? " KO" : "" ) << endl; + } + _elemSearcher = SMESH_MeshAlgos::GetElementSearcher( *_meshDS ); + } + return true; + } + // -------------------------------------------------------------------------------- + /*! + * \brief Find a source triangle including a point and return its barycentric coordinates + */ + const SMDS_MeshElement* FindFaceByPoint( const gp_Pnt p, + double & bc1, double & bc2, double & bc3 ) + { + const SMDS_MeshElement* tria = nullptr; + gp_XYZ projPnt = _elemSearcher->Project( p, SMDSAbs_Face, &tria ); + + int lcsID = tria->GetID() - NbEdges() - 1; + const TriaCoordSys* triaLCS = _traiLCSPool[ lcsID ]; + triaLCS->GetBaryCoords( projPnt, bc1, bc2, bc3 ); + + return tria; + } + // -------------------------------------------------------------------------------- + /*! + * \brief Return a point lying on a corresponding target triangle + */ + gp_Pnt GetPoint( const SMDS_MeshElement* srcTria, double & bc1, double & bc2, double & bc3 ) + { + const SMDS_MeshElement* tgtTria = _meshDS->FindElement( srcTria->GetID() ); + gp_Pnt p = ( bc1 * SMESH_NodeXYZ( tgtTria->GetNode(0) ) + + bc2 * SMESH_NodeXYZ( tgtTria->GetNode(1) ) + + bc3 * SMESH_NodeXYZ( tgtTria->GetNode(2) ) ); + return p; + } + // -------------------------------------------------------------------------------- + /*! + * \brief Return an UV of point lying on a corresponding target triangle + */ + gp_XY GetUV( const SMDS_MeshElement* srcTria, + double & bc1, double & bc2, double & bc3 ) + { + const SMDS_MeshElement* tgtTria = _meshDS->FindElement( srcTria->GetID() ); + TopoDS_Shape tgtShape = GetShapeToMesh(); + const TopoDS_Face& face = TopoDS::Face( tgtShape ); + + gp_XY p = ( bc1 * _helper.GetNodeUV( face, tgtTria->GetNode(0) ) + + bc2 * _helper.GetNodeUV( face, tgtTria->GetNode(1) ) + + bc3 * _helper.GetNodeUV( face, tgtTria->GetNode(2) ) ); + return p; + } + }; + + //================================================================================ + /*! + * \brief Calculate average size of faces + * Actually calculate average of min and max face size + */ + //================================================================================ + + double calcAverageFaceSize( SMESHDS_SubMesh* sm ) + { + double minLen = Precision::Infinite(), maxLen = 0; + for ( SMDS_ElemIteratorPtr fIt = sm->GetElements(); fIt->more(); ) + { + const SMDS_MeshElement* face = fIt->next(); + int nbNodes = face->NbCornerNodes(); + gp_XYZ pPrev = SMESH_NodeXYZ( face->GetNode( nbNodes - 1 )); + for ( int i = 0; i < nbNodes; ++i ) + { + SMESH_NodeXYZ p( face->GetNode( i )); + double len = ( p - pPrev ).SquareModulus(); + minLen = Min( len, minLen ); + maxLen = Max( len, maxLen ); + pPrev = p; + } + } + return 0.5 * ( Sqrt( minLen ) + Sqrt( maxLen )); } //================================================================================ /*! - * \brief Perform projection in case of quadrilateral faces + * \brief Perform projection from a quadrilateral FACE to another quadrilateral one */ //================================================================================ - bool projectQuads(const TopoDS_Face& /*tgtFace*/, - const TopoDS_Face& /*srcFace*/, - const TSideVector& /*tgtWires*/, - const TSideVector& /*srcWires*/, - const TAssocTool::TShapeShapeMap& /*shape2ShapeMap*/, - TAssocTool::TNodeNodeMap& /*src2tgtNodes*/, - const bool /*is1DComputed*/) + bool projectQuads(const TopoDS_Face& tgtFace, + const TopoDS_Face& srcFace, + const TSideVector& tgtWires, + const TSideVector& srcWires, + const TAssocTool::TShapeShapeMap& shape2ShapeMap, + TAssocTool::TNodeNodeMap& src2tgtNodes, + const bool is1DComputed) { - // SMESH_Mesh * tgtMesh = tgtWires[0]->GetMesh(); - // SMESH_Mesh * srcMesh = srcWires[0]->GetMesh(); - // //SMESHDS_Mesh * tgtMeshDS = tgtMesh->GetMeshDS(); - // SMESHDS_Mesh * srcMeshDS = srcMesh->GetMeshDS(); + SMESH_Mesh * tgtMesh = tgtWires[0]->GetMesh(); + SMESH_Mesh * srcMesh = srcWires[0]->GetMesh(); + SMESHDS_Mesh * tgtMeshDS = tgtMesh->GetMeshDS(); + SMESHDS_Mesh * srcMeshDS = srcMesh->GetMeshDS(); + + if ( srcWires.size() != 1 || // requirements below can be weaken + SMESH_MesherHelper::Count( tgtFace, TopAbs_EDGE, /*ignoreSame=*/true) != 4 || + SMESH_MesherHelper::Count( srcFace, TopAbs_EDGE, /*ignoreSame=*/true) != 4 ) + return false; + + // make auxiliary structured meshes that will be used to get corresponding + // points on the target FACE + QuadMesh srcQuadMesh( srcFace ), tgtQuadMesh( tgtFace ); + double avgSize = calcAverageFaceSize( srcMeshDS->MeshElements( srcFace )); + int nbSeg1 = (int) Max( 2., Max( srcWires[0]->EdgeLength(0), + srcWires[0]->EdgeLength(2)) / avgSize ); + int nbSeg2 = (int) Max( 2., Max( srcWires[0]->EdgeLength(1), + srcWires[0]->EdgeLength(3)) / avgSize ); + if ( !srcQuadMesh.Compute( srcWires, nbSeg1, nbSeg2, /*isSrc=*/true ) || + !tgtQuadMesh.Compute( tgtWires, nbSeg1, nbSeg2, /*isSrc=*/false )) + return false; + + // Make new faces + + // prepare the helper to adding quadratic elements if necessary + SMESH_MesherHelper* helper = tgtWires[0]->FaceHelper(); + helper->IsQuadraticSubMesh( tgtFace ); + + SMESHDS_SubMesh* srcSubDS = srcMeshDS->MeshElements( srcFace ); + if ( !is1DComputed && srcSubDS->NbElements() ) + helper->SetIsQuadratic( srcSubDS->GetElements()->next()->IsQuadratic() ); + + SMESH_MesherHelper* srcHelper = srcWires[0]->FaceHelper(); + SMESH_MesherHelper edgeHelper( *tgtMesh ); + edgeHelper.ToFixNodeParameters( true ); + + const SMDS_MeshNode* nullNode = 0; + TAssocTool::TNodeNodeMap::iterator srcN_tgtN; + + SMDS_ElemIteratorPtr elemIt = srcSubDS->GetElements(); + vector< const SMDS_MeshNode* > tgtNodes; + while ( elemIt->more() ) // loop on all mesh faces on srcFace + { + const SMDS_MeshElement* elem = elemIt->next(); + const int nbN = elem->NbCornerNodes(); + tgtNodes.resize( nbN ); + helper->SetElementsOnShape( false ); + for ( int i = 0; i < nbN; ++i ) // loop on nodes of the source element + { + const SMDS_MeshNode* srcNode = elem->GetNode(i); + srcN_tgtN = src2tgtNodes.insert( make_pair( srcNode, nullNode )).first; + if ( srcN_tgtN->second == nullNode ) + { + // create a new node + gp_Pnt srcP = SMESH_TNodeXYZ( srcNode ); + double bc[3]; + const SMDS_MeshElement* auxF = srcQuadMesh.FindFaceByPoint( srcP, bc[0], bc[1], bc[2] ); + gp_Pnt tgtP = tgtQuadMesh.GetPoint( auxF, bc[0], bc[1], bc[2] ); + SMDS_MeshNode* n = helper->AddNode( tgtP.X(), tgtP.Y(), tgtP.Z() ); + srcN_tgtN->second = n; + switch ( srcNode->GetPosition()->GetTypeOfPosition() ) + { + case SMDS_TOP_FACE: + { + gp_XY tgtUV = tgtQuadMesh.GetUV( auxF, bc[0], bc[1], bc[2] ); + tgtMeshDS->SetNodeOnFace( n, helper->GetSubShapeID(), tgtUV.X(), tgtUV.Y() ); + break; + } + case SMDS_TOP_EDGE: + { + const TopoDS_Edge& srcE = TopoDS::Edge( srcMeshDS->IndexToShape( srcNode->GetShapeID())); + const TopoDS_Edge& tgtE = TopoDS::Edge( shape2ShapeMap( srcE, /*isSrc=*/true )); + double srcU = srcHelper->GetNodeU( srcE, srcNode ); + tgtMeshDS->SetNodeOnEdge( n, tgtE, srcU ); + edgeHelper.SetSubShape( tgtE ); + double tol = BRep_Tool::MaxTolerance( tgtE, TopAbs_VERTEX ), distXYZ[4]; + /*isOk = */edgeHelper.CheckNodeU( tgtE, n, srcU, 2 * tol, /*force=*/true, distXYZ ); + //if ( isOk ) + tgtMeshDS->MoveNode( n, distXYZ[1], distXYZ[2], distXYZ[3] ); + SMDS_EdgePositionPtr( n->GetPosition() )->SetUParameter( srcU ); + break; + } + case SMDS_TOP_VERTEX: + { + const TopoDS_Shape & srcV = srcMeshDS->IndexToShape( srcNode->getshapeId() ); + const TopoDS_Shape & tgtV = shape2ShapeMap( srcV, /*isSrc=*/true ); + tgtP = BRep_Tool::Pnt( TopoDS::Vertex( tgtV )); + tgtMeshDS->MoveNode( n, tgtP.X(), tgtP.Y(), tgtP.Z() ); + tgtMeshDS->SetNodeOnVertex( n, TopoDS::Vertex( tgtV )); + break; + } + default:; + } + } + tgtNodes[i] = srcN_tgtN->second; + } + // create a new face + helper->SetElementsOnShape( true ); + switch ( nbN ) + { + case 3: helper->AddFace(tgtNodes[0], tgtNodes[1], tgtNodes[2]); break; + case 4: helper->AddFace(tgtNodes[0], tgtNodes[1], tgtNodes[2], tgtNodes[3]); break; + default: helper->AddPolygonalFace( tgtNodes ); + } + } // // loop on all mesh faces on srcFace + + return true; + + // below is projection of a structured source mesh - // if ( srcWires[0]->NbEdges() != 4 ) - // return false; // if ( !is1DComputed ) // return false; // for ( int iE = 0; iE < 4; ++iE ) @@ -1052,7 +1429,7 @@ namespace { // const SMDS_MeshElement* elem = elemIt->next(); // TFaceConn& tgtNodes = newFacesVec[ iFaceSrc++ ]; - // const int nbN = elem->NbCornerNodes(); + // const int nbN = elem->NbCornerNodes(); // tgtNodes.resize( nbN ); // for ( int i = 0; i < nbN; ++i ) // loop on nodes of the source element // { @@ -1066,7 +1443,7 @@ namespace { // { // tgtNodeOrXY.first = srcN_tgtN->second; // tgt node exists // } - // else + // else // { // // find XY of src node within the quadrilateral srcFace // if ( !block.ComputeParameters( SMESH_TNodeXYZ( srcNode ), @@ -1203,7 +1580,7 @@ namespace { tgtNbEW.front() != 4 || srcNbEW.front() != 4 ) return; // not quads - int srcNbSeg[4]; + smIdType srcNbSeg[4]; list< TopoDS_Edge >::iterator edgeS = srcEdges.begin(), edgeT = tgtEdges.begin(); for ( int i = 0; edgeS != srcEdges.end(); ++i, ++edgeS ) if ( SMESHDS_SubMesh* sm = srcMesh->GetMeshDS()->MeshElements( *edgeS )) @@ -1224,6 +1601,423 @@ namespace { tgtHelper.IthVertex( 1,*edgeS ), assocMap ); } + //================================================================================ + /*! + * \brief Find sub-shape association such that corresponding VERTEXes of + * two corresponding FACEs lie on lines parallel to thePiercingLine + */ + //================================================================================ + + bool findSubShapeAssociationByPiercing( const TopoDS_Face& theTgtFace, + SMESH_Mesh * /*theTgtMesh*/, + const TopoDS_Shape& theSrcShape, + SMESH_Mesh* theSrcMesh, + TAssocTool::TShapeShapeMap& theShape2ShapeMap, + Handle(Geom_Line) & thePiercingLine ) + { + list< TopoDS_Edge > tgtEdges, srcEdges; + list< int > tgtNbEW, srcNbEW; + int tgtNbW = SMESH_Block::GetOrderedEdges( TopoDS::Face( theTgtFace ), tgtEdges, tgtNbEW ); + + TopTools_IndexedMapOfShape tgtVV, srcVV; + for ( const TopoDS_Edge& tgtEdge : tgtEdges ) + tgtVV.Add( SMESH_MesherHelper::IthVertex( 0, tgtEdge )); + // if ( tgtVV.Size() < 2 ) + // return false; + + const int nbVV = tgtVV.Size(); + const gp_Pnt tgtP0 = BRep_Tool::Pnt( TopoDS::Vertex( tgtVV( 1 ))); + double minVertexDist = Precision::Infinite(), assocTol; + gp_Lin piercingLine; + TopoDS_Face assocSrcFace; + double tol; + + for ( TopExp_Explorer faceExp( theSrcShape, TopAbs_FACE ); faceExp.More(); faceExp.Next()) + { + const TopoDS_Face& srcFace = TopoDS::Face( faceExp.Current() ); + + int srcNbW = SMESH_Block::GetOrderedEdges( srcFace, srcEdges, srcNbEW ); + if ( tgtNbW != srcNbW ) + continue; + + srcVV.Clear( false ); + for ( const TopoDS_Edge& srcEdge : srcEdges ) + srcVV.Add( SMESH_MesherHelper::IthVertex( 0, srcEdge )); + if ( srcVV.Extent() != tgtVV.Extent() ) + continue; + + // make srcFace computed + SMESH_subMesh* srcFaceSM = theSrcMesh->GetSubMesh( srcFace ); + if ( !TAssocTool::MakeComputed( srcFaceSM )) + continue; + + // compute tolerance + double sumLen = 0, nbEdges = 0; + for ( const TopoDS_Edge& srcEdge : srcEdges ) + { + SMESH_subMesh* srcSM = theSrcMesh->GetSubMesh( srcEdge ); + if ( !srcSM->GetSubMeshDS() ) + continue; + SMDS_ElemIteratorPtr edgeIt = srcSM->GetSubMeshDS()->GetElements(); + while ( edgeIt->more() ) + { + const SMDS_MeshElement* edge = edgeIt->next(); + sumLen += SMESH_NodeXYZ( edge->GetNode( 0 )).Distance( edge->GetNode( 1 )); + nbEdges += 1; + } + } + if ( nbEdges == 0 ) + continue; + + tol = 0.1 * sumLen / nbEdges; + + // try to find corresponding VERTEXes + + gp_Lin line; + double vertexDist; + for ( int iSrcV0 = 1; iSrcV0 <= srcVV.Size(); ++iSrcV0 ) + { + const gp_Pnt srcP0 = BRep_Tool::Pnt( TopoDS::Vertex( srcVV( iSrcV0 ))); + try { + line.SetDirection( gp_Vec( srcP0, tgtP0 )); + } + catch (...) { + continue; + } + bool correspond; + for ( int iDir : { -1, 1 }) // move connected VERTEX forward and backward + { + correspond = true; + vertexDist = 0; + int iTgtV = 0, iSrcV = iSrcV0 - 1; + for ( int i = 1; i < tgtVV.Size() && correspond; ++i ) + { + iTgtV = ( iTgtV + 1 ) % nbVV; + iSrcV = ( iSrcV + iDir + nbVV ) % nbVV; + gp_Pnt tgtP = BRep_Tool::Pnt( TopoDS::Vertex( tgtVV( iTgtV + 1 ))); + gp_Pnt srcP = BRep_Tool::Pnt( TopoDS::Vertex( srcVV( iSrcV + 1 ))); + line.SetLocation( tgtP ); + correspond = ( line.SquareDistance( srcP ) < tol * tol ); + vertexDist += tgtP.SquareDistance( srcP ); + } + if ( correspond ) + break; + } + if ( correspond ) + { + if ( vertexDist < minVertexDist ) + { + minVertexDist = vertexDist; + piercingLine = line; + assocSrcFace = srcFace; + assocTol = tol; + } + break; + } + } + continue; + + } // loop on src FACEs + + if ( Precision::IsInfinite( minVertexDist )) + return false; // no correspondence found + + thePiercingLine = new Geom_Line( piercingLine ); + + // fill theShape2ShapeMap + + TAssocTool::InsertAssociation( theTgtFace, assocSrcFace, theShape2ShapeMap ); + + for ( const TopoDS_Shape& tgtV : tgtVV ) // fill theShape2ShapeMap with VERTEXes + { + gp_Pnt tgtP = BRep_Tool::Pnt( TopoDS::Vertex( tgtV )); + piercingLine.SetLocation( tgtP ); + bool found = false; + for ( const TopoDS_Shape& srcV : srcVV ) + { + gp_Pnt srcP = BRep_Tool::Pnt( TopoDS::Vertex( srcV )); + if ( piercingLine.SquareDistance( srcP ) < assocTol * assocTol ) + { + TAssocTool::InsertAssociation( tgtV, srcV, theShape2ShapeMap ); + found = true; + break; + } + } + if ( !found ) + return false; + } + + TopoDS_Vertex vvT[2], vvS[2], vvMapped[2]; + for ( const TopoDS_Edge& tgtEdge : tgtEdges ) // fill theShape2ShapeMap with EDGEs + { + if ( SMESH_Algo::isDegenerated( tgtEdge )) + continue; + + TopExp::Vertices( tgtEdge, vvT[0], vvT[1], true ); + if ( !theShape2ShapeMap.IsBound( vvT[0] ) || + !theShape2ShapeMap.IsBound( vvT[1] )) + return false; + + vvMapped[0] = TopoDS::Vertex( theShape2ShapeMap( vvT[0] )); + vvMapped[1] = TopoDS::Vertex( theShape2ShapeMap( vvT[1] )); + + bool found = false; + for ( TopExp_Explorer eExp( assocSrcFace, TopAbs_EDGE ); eExp.More(); eExp.Next()) + { + TopoDS_Edge srcEdge = TopoDS::Edge( eExp.Current() ); + TopExp::Vertices( srcEdge, vvS[0], vvS[1], true ); + found = (( vvMapped[0].IsSame( vvS[0] ) && vvMapped[1].IsSame( vvS[1] )) || + ( vvMapped[0].IsSame( vvS[1] ) && vvMapped[1].IsSame( vvS[0] ))); + + if ( found && nbVV < 3 ) + { + BRepAdaptor_Curve tgtCurve( tgtEdge ); + gp_Pnt tgtP = tgtCurve.Value( 0.5 * ( tgtCurve.FirstParameter() + + tgtCurve.LastParameter() )); + thePiercingLine->SetLocation( tgtP ); + + double f,l; + Handle(Geom_Curve) srcCurve = BRep_Tool::Curve( srcEdge, f,l ); + if ( srcCurve.IsNull() ) + { + found = false; + continue; + } + GeomAPI_ExtremaCurveCurve extrema( thePiercingLine, srcCurve ); + if ( !extrema.Extrema().IsDone() || + extrema.Extrema().IsParallel() || + extrema.NbExtrema() == 0 || + extrema.LowerDistance() > tol ) + found = false; + } + if ( found ) + { + if ( !vvMapped[0].IsSame( vvS[0] )) + srcEdge.Reverse(); + TAssocTool::InsertAssociation( tgtEdge, srcEdge, theShape2ShapeMap ); + break; + } + } + if ( !found ) + return false; + } + + return true; + + } // findSubShapeAssociationByPiercing() + + //================================================================================ + /*! + * \brief Project by piercing theTgtFace by lines parallel to thePiercingLine + */ + //================================================================================ + + bool projectByPiercing(Handle(Geom_Line) thePiercingLine, + const TopoDS_Face& theTgtFace, + const TopoDS_Face& theSrcFace, + const TSideVector& theTgtWires, + const TSideVector& theSrcWires, + const TAssocTool::TShapeShapeMap& theShape2ShapeMap, + TAssocTool::TNodeNodeMap& theSrc2tgtNodes, + const bool theIs1DComputed) + { + SMESH_Mesh * tgtMesh = theTgtWires[0]->GetMesh(); + SMESH_Mesh * srcMesh = theSrcWires[0]->GetMesh(); + + if ( thePiercingLine.IsNull() ) + { + // try to set thePiercingLine by VERTEX association of theShape2ShapeMap + + const double tol = 0.1 * theSrcWires[0]->Length() / theSrcWires[0]->NbSegments(); + + for ( TopExp_Explorer vExp( theTgtFace, TopAbs_VERTEX ); vExp.More(); vExp.Next() ) + { + const TopoDS_Vertex & tgtV = TopoDS::Vertex( vExp.Current() ); + const TopoDS_Vertex & srcV = TopoDS::Vertex( theShape2ShapeMap( tgtV )); + gp_Pnt tgtP = BRep_Tool::Pnt( tgtV ); + gp_Pnt srcP = BRep_Tool::Pnt( srcV ); + if ( thePiercingLine.IsNull() ) // set thePiercingLine + { + gp_Lin line; + try { + line.SetDirection( gp_Vec( srcP, tgtP )); + line.SetLocation( tgtP ); + thePiercingLine = new Geom_Line( line ); + } + catch ( ... ) + { + continue; + } + } + else // check thePiercingLine + { + thePiercingLine->SetLocation( tgtP ); + if ( thePiercingLine->Lin().SquareDistance( srcP ) > tol * tol ) + return false; + } + } + + for ( TopExp_Explorer eExp( theTgtFace, TopAbs_EDGE ); eExp.More(); eExp.Next() ) + { + BRepAdaptor_Curve tgtCurve( TopoDS::Edge( eExp.Current() )); + gp_Pnt tgtP = tgtCurve.Value( 0.5 * ( tgtCurve.FirstParameter() + + tgtCurve.LastParameter() )); + thePiercingLine->SetLocation( tgtP ); + + double f,l; + TopoDS_Edge srcEdge = TopoDS::Edge( theShape2ShapeMap( eExp.Current() )); + Handle(Geom_Curve) srcCurve = BRep_Tool::Curve( srcEdge, f,l ); + if ( srcCurve.IsNull() ) + continue; + GeomAPI_ExtremaCurveCurve extrema( thePiercingLine, srcCurve, + -Precision::Infinite(), Precision::Infinite(), f, l ); + if ( !extrema.Extrema().IsDone() || + extrema.Extrema().IsParallel() || + extrema.NbExtrema() == 0 || + extrema.LowerDistance() > tol ) + return false; + } + } // if ( thePiercingLine.IsNull() ) + + SMESHDS_SubMesh* srcSubDS = srcMesh->GetMeshDS()->MeshElements( theSrcFace ); + + SMESH_MesherHelper* helper = theTgtWires[0]->FaceHelper(); + if ( theIs1DComputed ) + helper->IsQuadraticSubMesh( theTgtFace ); + else + helper->SetIsQuadratic( srcSubDS->GetElements()->next()->IsQuadratic() ); + helper->SetElementsOnShape( true ); + SMESHDS_Mesh* tgtMeshDS = tgtMesh->GetMeshDS(); + + Handle(Geom_Surface) tgtSurface = BRep_Tool::Surface( theTgtFace ); +#if OCC_VERSION_LARGE < 0x07070000 + Handle(GeomAdaptor_HSurface) tgtSurfAdaptor = new GeomAdaptor_HSurface( tgtSurface ); + Handle(GeomAdaptor_HCurve) piercingCurve = new GeomAdaptor_HCurve( thePiercingLine ); +#else + Handle(GeomAdaptor_Surface) tgtSurfAdaptor = new GeomAdaptor_Surface( tgtSurface ); + Handle(GeomAdaptor_Curve) piercingCurve = new GeomAdaptor_Curve( thePiercingLine ); +#endif + IntCurveSurface_HInter intersect; + + SMESH_MesherHelper* srcHelper = theSrcWires[0]->FaceHelper(); + + const SMDS_MeshNode* nullNode = 0; + TAssocTool::TNodeNodeMap::iterator srcN_tgtN; + vector< const SMDS_MeshNode* > tgtNodes; + + SMDS_ElemIteratorPtr elemIt = srcSubDS->GetElements(); + while ( elemIt->more() ) // loop on all mesh faces on srcFace + { + const SMDS_MeshElement* elem = elemIt->next(); + const int nbN = elem->NbCornerNodes(); + tgtNodes.resize( nbN ); + for ( int i = 0; i < nbN; ++i ) // loop on nodes of the source element + { + const SMDS_MeshNode* srcNode = elem->GetNode(i); + srcN_tgtN = theSrc2tgtNodes.insert( make_pair( srcNode, nullNode )).first; + if ( srcN_tgtN->second == nullNode ) + { + // create a new node + thePiercingLine->SetLocation( SMESH_NodeXYZ( srcNode )); + intersect.Perform( piercingCurve, tgtSurfAdaptor ); + bool pierced = ( intersect.IsDone() && intersect.NbPoints() > 0 ); + double U, V; + const SMDS_MeshNode* n = nullNode; + if ( pierced ) + { + double W, minW = Precision::Infinite(); + gp_Pnt tgtP; + for ( int iInt = 1; iInt <= intersect.NbPoints(); ++iInt ) + { + W = intersect.Point( iInt ).W(); + if ( 0 < W && W < minW ) + { + U = intersect.Point( iInt ).U(); + V = intersect.Point( iInt ).V(); + tgtP = intersect.Point( iInt ).Pnt(); + minW = W; + } + } + n = tgtMeshDS->AddNode( tgtP.X(), tgtP.Y(), tgtP.Z() ); + } + + SMDS_TypeOfPosition shapeType = srcNode->GetPosition()->GetTypeOfPosition(); + TopoDS_Shape srcShape; + if ( shapeType != SMDS_TOP_FACE ) + { + srcShape = srcHelper->GetSubShapeByNode( srcNode, srcHelper->GetMeshDS() ); + if ( !theShape2ShapeMap.IsBound( srcShape, /*isSrc=*/true )) + { + if ( n ) // INTERNAL shape w/o corresponding target shape (3D_mesh_Extrusion_02/E0) + shapeType = SMDS_TOP_FACE; + else + return false; + } + } + + switch ( shapeType ) + { + case SMDS_TOP_FACE: { + if ( !n ) + return false; + tgtMeshDS->SetNodeOnFace( n, helper->GetSubShapeID(), U, V ); + break; + } + case SMDS_TOP_EDGE: { + TopoDS_Edge tgtEdge = TopoDS::Edge( theShape2ShapeMap( srcShape, /*isSrc=*/true )); + if ( n ) + { + U = Precision::Infinite(); + helper->CheckNodeU( tgtEdge, n, U, Precision::PConfusion()); + } + else + { + Handle(Geom_Curve) tgtCurve = BRep_Tool::Curve( tgtEdge, U,V ); + if ( tgtCurve.IsNull() ) + return false; + GeomAPI_ExtremaCurveCurve extrema( thePiercingLine, tgtCurve ); + if ( !extrema.Extrema().IsDone() || + extrema.Extrema().IsParallel() || + extrema.NbExtrema() == 0 ) + return false; + gp_Pnt pOnLine, pOnEdge; + extrema.NearestPoints( pOnLine, pOnEdge ); + extrema.LowerDistanceParameters( V, U ); + n = tgtMeshDS->AddNode( pOnEdge.X(), pOnEdge.Y(), pOnEdge.Z() ); + } + tgtMeshDS->SetNodeOnEdge( n, tgtEdge, U ); + break; + } + case SMDS_TOP_VERTEX: { + TopoDS_Shape tgtV = theShape2ShapeMap( srcShape, /*isSrc=*/true ); + if ( !n ) + { + gp_Pnt tgtP = BRep_Tool::Pnt( TopoDS::Vertex( tgtV )); + n = tgtMeshDS->AddNode( tgtP.X(), tgtP.Y(), tgtP.Z() ); + } + tgtMeshDS->SetNodeOnVertex( n, TopoDS::Vertex( tgtV )); + break; + } + default:; + } + srcN_tgtN->second = n; + } + tgtNodes[i] = srcN_tgtN->second; + } + // create a new face (with reversed orientation) + switch ( nbN ) + { + case 3: helper->AddFace(tgtNodes[0], tgtNodes[2], tgtNodes[1]); break; + case 4: helper->AddFace(tgtNodes[0], tgtNodes[3], tgtNodes[2], tgtNodes[1]); break; + } + } // loop on all mesh faces on srcFace + + return true; + + } // projectByPiercing() + + + } // namespace @@ -1260,20 +2054,29 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& TAssocTool::InitVertexAssociation( _sourceHypo, shape2ShapeMap ); if ( shape2ShapeMap.IsEmpty() ) initAssoc4Quad2Closed( tgtFace, helper, srcShape, srcMesh, shape2ShapeMap ); + + Handle(Geom_Line) piercingLine; + bool piercingTried = false; + if ( !TAssocTool::FindSubShapeAssociation( tgtFace, tgtMesh, srcShape, srcMesh, shape2ShapeMap) || !shape2ShapeMap.IsBound( tgtFace )) { - if ( srcShape.ShapeType() == TopAbs_FACE ) + piercingTried = true; + if ( !findSubShapeAssociationByPiercing( tgtFace, tgtMesh, srcShape, srcMesh, + shape2ShapeMap, piercingLine )) { - int nbE1 = helper.Count( tgtFace, TopAbs_EDGE, /*ignoreSame=*/true ); - int nbE2 = helper.Count( srcShape, TopAbs_EDGE, /*ignoreSame=*/true ); - if ( nbE1 != nbE2 ) - return error(COMPERR_BAD_SHAPE, - SMESH_Comment("Different number of edges in source and target faces: ") - << nbE2 << " and " << nbE1 ); + if ( srcShape.ShapeType() == TopAbs_FACE ) + { + int nbE1 = helper.Count( tgtFace, TopAbs_EDGE, /*ignoreSame=*/true ); + int nbE2 = helper.Count( srcShape, TopAbs_EDGE, /*ignoreSame=*/true ); + if ( nbE1 != nbE2 ) + return error(COMPERR_BAD_SHAPE, + SMESH_Comment("Different number of edges in source and target faces: ") + << nbE2 << " and " << nbE1 ); + } + return error(COMPERR_BAD_SHAPE,"Topology of source and target faces seems different" ); } - return error(COMPERR_BAD_SHAPE,"Topology of source and target faces seems different" ); } TopoDS_Face srcFace = TopoDS::Face( shape2ShapeMap( tgtFace ).Oriented(TopAbs_FORWARD)); @@ -1310,6 +2113,13 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& bool projDone = false; + if ( !projDone && !piercingLine.IsNull() ) + { + // project by piercing tgtFace by lines parallel to piercingLine + projDone = projectByPiercing( piercingLine, tgtFace, srcFace, tgtWires, srcWires, + shape2ShapeMap, _src2tgtNodes, is1DComputed ); + piercingTried = true; + } if ( !projDone ) { // try to project from the same face with different location @@ -1325,10 +2135,15 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& if ( !projDone ) { // projection in case of quadrilateral faces - // NOT IMPLEMENTED, returns false projDone = projectQuads( tgtFace, srcFace, tgtWires, srcWires, shape2ShapeMap, _src2tgtNodes, is1DComputed); } + if ( !projDone && !piercingTried ) + { + // project by piercing tgtFace by lines parallel to piercingLine + projDone = projectByPiercing( piercingLine, tgtFace, srcFace, tgtWires, srcWires, + shape2ShapeMap, _src2tgtNodes, is1DComputed ); + } // it will remove mesh built on edges and vertices in failure case MeshCleaner cleaner( tgtSubMesh ); @@ -1337,7 +2152,7 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& { _src2tgtNodes.clear(); // -------------------- - // Prepare to mapping + // Prepare to mapping // -------------------- // Check if node projection to a face is needed @@ -1444,7 +2259,7 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& // compare nb nodes on srcEdge1 and srcEdge2 if ( srcEdge2 != srcEdges.end() ) { - int nbN1 = 0, nbN2 = 0; + smIdType nbN1 = 0, nbN2 = 0; if ( SMESHDS_SubMesh* sm = srcMesh->GetMeshDS()->MeshElements( srcEdge1 )) nbN1 = sm->NbNodes(); if ( SMESHDS_SubMesh* sm = srcMesh->GetMeshDS()->MeshElements( *srcEdge2 )) @@ -1607,7 +2422,7 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& if ( u2nodesMaps[ NEW_NODES ].size() > 0 && u2nodesMaps[ OLD_NODES ].size() > 0 ) { - u_oldNode = u2nodesMaps[ OLD_NODES ].begin(); + u_oldNode = u2nodesMaps[ OLD_NODES ].begin(); newEnd = u2nodesMaps[ OLD_NODES ].end(); for ( ; u_oldNode != newEnd; ++u_oldNode ) SMESH_Algo::addBadInputElement( u_oldNode->second ); @@ -1640,7 +2455,7 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& // Make groups of nodes to merge - u_oldNode = u2nodesMaps[ OLD_NODES ].begin(); + u_oldNode = u2nodesMaps[ OLD_NODES ].begin(); u_newNode = u2nodesMaps[ NEW_NODES ].begin(); newEnd = u2nodesMaps[ NEW_NODES ].end(); u_newOnSeam = u2nodesOnSeam.begin(); @@ -1666,9 +2481,9 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& // Merge SMESH_MeshEditor editor( tgtMesh ); - int nbFaceBeforeMerge = tgtSubMesh->GetSubMeshDS()->NbElements(); + smIdType nbFaceBeforeMerge = tgtSubMesh->GetSubMeshDS()->NbElements(); editor.MergeNodes( groupsOfNodes ); - int nbFaceAtferMerge = tgtSubMesh->GetSubMeshDS()->NbElements(); + smIdType nbFaceAtferMerge = tgtSubMesh->GetSubMeshDS()->NbElements(); if ( nbFaceBeforeMerge != nbFaceAtferMerge && !helper.HasDegeneratedEdges() ) return error(COMPERR_BAD_INPUT_MESH, "Probably invalid node parameters on geom faces"); @@ -1704,8 +2519,10 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& TAssocTool::Morph morph( srcWires ); morph.Perform( helper, tgtWires, helper.GetSurface( tgtFace ), _src2tgtNodes, /*moveAll=*/true ); + if(SALOME::VerbosityActivated()) + cout << "StdMeshers_Projection_2D: Projection mesh IsDistorted2D() ==> do morph" << endl; - if ( !fixDistortedFaces( helper, tgtWires )) + if ( !fixDistortedFaces( helper, tgtWires )) // smooth and check return error("Invalid mesh generated"); } // --------------------------- @@ -1761,7 +2578,7 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& //======================================================================= //function : Evaluate -//purpose : +//purpose : //======================================================================= bool StdMeshers_Projection_2D::Evaluate(SMESH_Mesh& theMesh, @@ -1796,7 +2613,7 @@ bool StdMeshers_Projection_2D::Evaluate(SMESH_Mesh& theMesh, // Assure that mesh on a source Face is computed/evaluated // ------------------------------------------------------- - std::vector aVec; + std::vector aVec; SMESH_subMesh* srcSubMesh = srcMesh->GetSubMesh( srcFace ); if ( srcSubMesh->IsMeshComputed() )