X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Prism_3D.cxx;h=80cdc5edd68542ab14ebe5a018093bfc91dded10;hp=5ae0a4b0bd19826ca57415be8af16fa6c1d3998b;hb=24dd5df5f053455d186962be79786dc9237a1a0e;hpb=87a7f0d0495483874b3fb322ffef2866d6bd30c6 diff --git a/src/StdMeshers/StdMeshers_Prism_3D.cxx b/src/StdMeshers/StdMeshers_Prism_3D.cxx index 5ae0a4b0b..80cdc5edd 100644 --- a/src/StdMeshers/StdMeshers_Prism_3D.cxx +++ b/src/StdMeshers/StdMeshers_Prism_3D.cxx @@ -52,7 +52,6 @@ #include #include #include -#include #include #include #include @@ -112,6 +111,7 @@ namespace { algo->myProxyMesh.reset( new SMESH_ProxyMesh( *helper->GetMesh() )); algo->myQuadList.clear(); + algo->myHelper = 0; if ( helper ) algo->_quadraticMesh = helper->GetIsQuadratic(); @@ -1170,6 +1170,11 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) if ( !assocOrProjBottom2Top( bottomToTopTrsf, thePrism ) ) // it also fills myBotToColumnMap return false; + // If all "vertical" EDGEs are straight, then all nodes of an internal node column + // are located on a line connecting the top node and the bottom node. + bool isStrightColunm = allVerticalEdgesStraight( thePrism ); + if ( isStrightColunm ) + myUseBlock = false; // Create nodes inside the block @@ -1177,10 +1182,13 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) { // use transformation (issue 0020680, IPAL0052499) or a "straight line" approach StdMeshers_Sweeper sweeper; + sweeper.myHelper = myHelper; + sweeper.myBotFace = thePrism.myBottom; + sweeper.myTopFace = thePrism.myTop; // load boundary nodes into sweeper bool dummy; - const SMDS_MeshNode* prevN0 = 0, *prevN1 = 0; + std::set< const SMDS_MeshNode* > usedEndNodes; list< TopoDS_Edge >::const_iterator edge = thePrism.myBottomEdges.begin(); for ( ; edge != thePrism.myBottomEdges.end(); ++edge ) { @@ -1191,37 +1199,33 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) TParam2ColumnMap::iterator u2colIt = u2col->begin(), u2colEnd = u2col->end(); const SMDS_MeshNode* n0 = u2colIt->second[0]; const SMDS_MeshNode* n1 = u2col->rbegin()->second[0]; - if ( n0 == prevN0 || n0 == prevN1 ) ++u2colIt; - if ( n1 == prevN0 || n1 == prevN1 ) --u2colEnd; - prevN0 = n0; prevN1 = n1; + if ( !usedEndNodes.insert ( n0 ).second ) ++u2colIt; + if ( !usedEndNodes.insert ( n1 ).second ) --u2colEnd; for ( ; u2colIt != u2colEnd; ++u2colIt ) sweeper.myBndColumns.push_back( & u2colIt->second ); } // load node columns inside the bottom FACE - TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin(); sweeper.myIntColumns.reserve( myBotToColumnMap.size() ); + TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin(); for ( ; bot_column != myBotToColumnMap.end(); ++bot_column ) sweeper.myIntColumns.push_back( & bot_column->second ); myHelper->SetElementsOnShape( true ); - // If all "vertical" EDGEs are straight, then all nodes of an internal node column - // are located on a line connecting the top node and the bottom node. - bool isStrightColunm = allVerticalEdgesStraight( thePrism ); if ( !isStrightColunm ) { double tol = getSweepTolerance( thePrism ); bool allowHighBndError = !isSimpleBottom( thePrism ); - myUseBlock = !sweeper.ComputeNodes( *myHelper, tol, allowHighBndError ); + myUseBlock = !sweeper.ComputeNodesByTrsf( tol, allowHighBndError ); } else if ( sweeper.CheckSameZ() ) { - myUseBlock = !sweeper.ComputeNodesOnStraightSameZ( *myHelper ); + myUseBlock = !sweeper.ComputeNodesOnStraightSameZ(); } else { - myUseBlock = !sweeper.ComputeNodesOnStraight( *myHelper, thePrism.myBottom, thePrism.myTop ); + myUseBlock = !sweeper.ComputeNodesOnStraight(); } myHelper->SetElementsOnShape( false ); } @@ -1241,6 +1245,14 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) // column nodes; middle part of the column are zero pointers TNodeColumn& column = bot_column->second; + // check if a column is already computed using non-block approach + size_t i; + for ( i = 0; i < column.size(); ++i ) + if ( !column[ i ]) + break; + if ( i == column.size() ) + continue; // all nodes created + gp_XYZ botParams, topParams; if ( !tBotNode.HasParams() ) { @@ -1625,7 +1637,7 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism) TopoDS_Vertex v = myHelper->IthVertex( is2ndV, E ); mesh->GetSubMesh( v )->ComputeStateEngine( SMESH_subMesh::COMPUTE ); const SMDS_MeshNode* n = SMESH_Algo::VertexNode( v, meshDS ); - newNodes[ is2ndV ? 0 : newNodes.size()-1 ] = (SMDS_MeshNode*) n; + newNodes[ is2ndV ? newNodes.size()-1 : 0 ] = (SMDS_MeshNode*) n; } // compute nodes on target EDGEs @@ -1788,9 +1800,8 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism) if ( myHelper->GetIsQuadratic() ) { // fill myHelper with medium nodes built by quadAlgo - SMDS_ElemIteratorPtr fIt = fSM->GetSubMeshDS()->GetElements(); - while ( fIt->more() ) - myHelper->AddTLinks( dynamic_cast( fIt->next() )); + for ( SMDS_ElemIteratorPtr fIt = fSM->GetSubMeshDS()->GetElements(); fIt->more(); ) + myHelper->AddTLinks( SMDS_Mesh::DownCast( fIt->next() )); } } } @@ -2281,7 +2292,7 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf TNode2ColumnMap::iterator bN_col = myBotToColumnMap.insert( make_pair ( bN, TNodeColumn() )).first; TNodeColumn & column = bN_col->second; - column.resize( zSize ); + column.resize( zSize, 0 ); column.front() = botNode; column.back() = topNode; } @@ -2580,7 +2591,7 @@ double StdMeshers_Prism_3D::getSweepTolerance( const Prism_3D::TPrismTopo& thePr //======================================================================= //function : isSimpleQuad -//purpose : check if the bottom FACE is meshable with nice qudrangles, +//purpose : check if the bottom FACE is meshable with nice quadrangles, // if so the block aproach can work rather fast. // This is a temporary mean caused by problems in StdMeshers_Sweeper //======================================================================= @@ -3053,7 +3064,7 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA if ( side._topEdge.IsNull() ) { - // find vertical EDGEs --- EGDEs shared with neighbor side FACEs + // find vertical EDGEs --- EDGEs shared with neighbor side FACEs for ( int is2nd = 0; is2nd < 2 && isOK; ++is2nd ) // 2 adjacent neighbors { int di = is2nd ? 1 : -1; @@ -3414,7 +3425,7 @@ bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism, if ( !botSM ) // find a proper bottom { bool savedSetErrorToSM = mySetErrorToSM; - mySetErrorToSM = false; // ingore errors in initPrism() + mySetErrorToSM = false; // ignore errors in initPrism() // search among meshed FACEs list< SMESH_subMesh* >::iterator sm = meshedSubMesh.begin(); @@ -3588,10 +3599,10 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper, if ( !myHelper->LoadNodeColumns( faceColumns, (*quad)->face, quadBot, meshDS )) return error(COMPERR_BAD_INPUT_MESH, TCom("Can't find regular quadrangle mesh ") << "on a side face #" << MeshDS()->ShapeToIndex( (*quad)->face )); - - if ( !faceColumns.empty() && (int)faceColumns.begin()->second.size() != VerticalSize() ) - return error(COMPERR_BAD_INPUT_MESH, "Different 'vertical' discretization"); } + if ( !faceColumns.empty() && (int)faceColumns.begin()->second.size() != VerticalSize() ) + return error(COMPERR_BAD_INPUT_MESH, "Different 'vertical' discretization"); + // edge columns int id = MeshDS()->ShapeToIndex( *edgeIt ); bool isForward = true; // meaningless for intenal wires @@ -4837,7 +4848,7 @@ TPCurveOnHorFaceAdaptor::TPCurveOnHorFaceAdaptor( const TSideFace* sideFace, } if ( !C2d.IsNull() ) { - double u = static_cast< const SMDS_EdgePosition* >( n->GetPosition() )->GetUParameter(); + double u = SMDS_EdgePositionPtr( n->GetPosition() )->GetUParameter(); if ( f <= u && u <= l ) { uv = C2d->Value( u ).XY(); @@ -4891,6 +4902,7 @@ bool StdMeshers_Sweeper::projectIntPoints(const vector< gp_XYZ >& fromBndPoin const vector< gp_XYZ >& toBndPoints, const vector< gp_XYZ >& fromIntPoints, vector< gp_XYZ >& toIntPoints, + const double r, NSProjUtils::TrsfFinder3D& trsf, vector< gp_XYZ > * bndError) { @@ -4915,60 +4927,40 @@ bool StdMeshers_Sweeper::projectIntPoints(const vector< gp_XYZ >& fromBndPoin (*bndError)[ iP ] = toBndPoints[ iP ] - fromTrsf; } } - return true; -} -//================================================================================ -/*! - * \brief Add boundary error to ineternal points - */ -//================================================================================ - -void StdMeshers_Sweeper::applyBoundaryError(const vector< gp_XYZ >& bndPoints, - const vector< gp_XYZ >& bndError1, - const vector< gp_XYZ >& bndError2, - const double r, - vector< gp_XYZ >& intPoints, - vector< double >& int2BndDist) -{ - // fix each internal point - const double eps = 1e-100; - for ( size_t iP = 0; iP < intPoints.size(); ++iP ) + // apply boundary error + if ( bndError && toIntPoints.size() == myTopBotTriangles.size() ) { - gp_XYZ & intPnt = intPoints[ iP ]; - - // compute distance from intPnt to each boundary node - double int2BndDistSum = 0; - for ( size_t iBnd = 0; iBnd < bndPoints.size(); ++iBnd ) + for ( size_t iP = 0; iP < toIntPoints.size(); ++iP ) { - int2BndDist[ iBnd ] = 1 / (( intPnt - bndPoints[ iBnd ]).SquareModulus() + eps ); - int2BndDistSum += int2BndDist[ iBnd ]; - } - - // apply bndError - for ( size_t iBnd = 0; iBnd < bndPoints.size(); ++iBnd ) - { - intPnt += bndError1[ iBnd ] * ( 1 - r ) * int2BndDist[ iBnd ] / int2BndDistSum; - intPnt += bndError2[ iBnd ] * r * int2BndDist[ iBnd ] / int2BndDistSum; + const TopBotTriangles& tbTrias = myTopBotTriangles[ iP ]; + for ( int i = 0; i < 3; ++i ) // boundary errors at 3 triangle nodes + { + toIntPoints[ iP ] += + ( (*bndError)[ tbTrias.myBotTriaNodes[i] ] * tbTrias.myBotBC[i] * ( 1 - r ) + + (*bndError)[ tbTrias.myTopTriaNodes[i] ] * tbTrias.myTopBC[i] * ( r )); + } } } + + return true; } //================================================================================ /*! - * \brief Create internal nodes of the prism + * \brief Create internal nodes of the prism by computing an affine transformation + * from layer to layer */ //================================================================================ -bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, - const double tol, - const bool allowHighBndError) +bool StdMeshers_Sweeper::ComputeNodesByTrsf( const double tol, + const bool allowHighBndError) { const size_t zSize = myBndColumns[0]->size(); const size_t zSrc = 0, zTgt = zSize-1; if ( zSize < 3 ) return true; - vector< vector< gp_XYZ > > intPntsOfLayer( zSize ); // node coodinates to compute + vector< vector< gp_XYZ > > intPntsOfLayer( zSize ); // node coordinates to compute // set coordinates of src and tgt nodes for ( size_t z = 0; z < intPntsOfLayer.size(); ++z ) intPntsOfLayer[ z ].resize( myIntColumns.size() ); @@ -4978,7 +4970,11 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, intPntsOfLayer[ zTgt ][ iP ] = intPoint( iP, zTgt ); } - // compute coordinates of internal nodes by projecting (transfroming) src and tgt + // for each internal column find boundary nodes whose error to use for correction + prepareTopBotDelaunay(); + bool isErrorCorrectable = findDelaunayTriangles(); + + // compute coordinates of internal nodes by projecting (transforming) src and tgt // nodes towards the central layer vector< NSProjUtils::TrsfFinder3D > trsfOfLayer( zSize ); @@ -5004,10 +5000,12 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, } if (! projectIntPoints( fromSrcBndPnts, toSrcBndPnts, intPntsOfLayer[ zS-1 ], intPntsOfLayer[ zS ], + zS / ( zSize - 1.), trsfOfLayer [ zS-1 ], & bndError[ zS-1 ])) return false; if (! projectIntPoints( fromTgtBndPnts, toTgtBndPnts, intPntsOfLayer[ zT+1 ], intPntsOfLayer[ zT ], + zT / ( zSize - 1.), trsfOfLayer [ zT+1 ], & bndError[ zT+1 ])) return false; @@ -5032,6 +5030,22 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, fromSrcBndPnts.swap( toSrcBndPnts ); } + // Evaluate an error of boundary points + + if ( !isErrorCorrectable && !allowHighBndError ) + { + for ( size_t iP = 0; iP < myBndColumns.size(); ++iP ) + { + double sumError = 0; + for ( size_t z = 1; z < zS; ++z ) // loop on layers + sumError += ( bndError[ z-1 ][ iP ].Modulus() + + bndError[ zSize-z ][ iP ].Modulus() ); + + if ( sumError > tol ) + return false; + } + } + // Compute two projections of internal points to the central layer // in order to evaluate an error of internal points @@ -5046,10 +5060,12 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, } if (! projectIntPoints( fromSrcBndPnts, toSrcBndPnts, intPntsOfLayer[ zS-1 ], centerSrcIntPnts, + zS / ( zSize - 1.), trsfOfLayer [ zS-1 ], & bndError[ zS-1 ])) return false; if (! projectIntPoints( fromTgtBndPnts, toTgtBndPnts, intPntsOfLayer[ zT+1 ], centerTgtIntPnts, + zT / ( zSize - 1.), trsfOfLayer [ zT+1 ], & bndError[ zT+1 ])) return false; @@ -5068,24 +5084,7 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, (intPntsOfLayer[ zS-1 ][ iP ] - centerTgtIntPnts[ iP ]).SquareModulus() < tol*tol; } - // Evaluate an error of boundary points - - bool bndErrorIsSmall = true; - for ( size_t iP = 0; ( iP < myBndColumns.size() && bndErrorIsSmall ); ++iP ) - { - double sumError = 0; - for ( size_t z = 1; z < zS; ++z ) // loop on layers - sumError += ( bndError[ z-1 ][ iP ].Modulus() + - bndError[ zSize-z ][ iP ].Modulus() ); - - bndErrorIsSmall = ( sumError < tol ); - } - - if ( !bndErrorIsSmall && !allowHighBndError ) - return false; - // compute final points on the central layer - std::vector< double > int2BndDist( myBndColumns.size() ); // work array of applyBoundaryError() double r = zS / ( zSize - 1.); if ( zS == zT ) { @@ -5094,11 +5093,6 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, intPntsOfLayer[ zS ][ iP ] = ( 1 - r ) * centerSrcIntPnts[ iP ] + r * centerTgtIntPnts[ iP ]; } - if ( !bndErrorIsSmall ) - { - applyBoundaryError( toSrcBndPnts, bndError[ zS-1 ], bndError[ zS+1 ], r, - intPntsOfLayer[ zS ], int2BndDist ); - } } else { @@ -5109,17 +5103,8 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, intPntsOfLayer[ zT ][ iP ] = r * intPntsOfLayer[ zT ][ iP ] + ( 1 - r ) * centerTgtIntPnts[ iP ]; } - if ( !bndErrorIsSmall ) - { - applyBoundaryError( toSrcBndPnts, bndError[ zS-1 ], bndError[ zS+1 ], r, - intPntsOfLayer[ zS ], int2BndDist ); - applyBoundaryError( toTgtBndPnts, bndError[ zT+1 ], bndError[ zT-1 ], r, - intPntsOfLayer[ zT ], int2BndDist ); - } } - centerIntErrorIsSmall = true; // 3D_mesh_Extrusion_00/A3 - bndErrorIsSmall = true; if ( !centerIntErrorIsSmall ) { // Compensate the central error; continue adding projection @@ -5151,9 +5136,11 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, } projectIntPoints( fromSrcBndPnts, toSrcBndPnts, fromSrcIntPnts, toSrcIntPnts, + zS / ( zSize - 1.), trsfOfLayer[ zS+1 ], & srcBndError ); projectIntPoints( fromTgtBndPnts, toTgtBndPnts, fromTgtIntPnts, toTgtIntPnts, + zT / ( zSize - 1.), trsfOfLayer[ zT-1 ], & tgtBndError ); // if ( zS == zTgt - 1 ) @@ -5184,15 +5171,6 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, zTIntPnts[ iP ] = r * zTIntPnts[ iP ] + ( 1 - r ) * toTgtIntPnts[ iP ]; } - // compensate bnd error - if ( !bndErrorIsSmall ) - { - applyBoundaryError( toSrcBndPnts, srcBndError, bndError[ zS+1 ], r, - intPntsOfLayer[ zS ], int2BndDist ); - applyBoundaryError( toTgtBndPnts, tgtBndError, bndError[ zT-1 ], r, - intPntsOfLayer[ zT ], int2BndDist ); - } - fromSrcBndPnts.swap( toSrcBndPnts ); fromSrcIntPnts.swap( toSrcIntPnts ); fromTgtBndPnts.swap( toTgtBndPnts ); @@ -5200,27 +5178,8 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, } } // if ( !centerIntErrorIsSmall ) - else if ( !bndErrorIsSmall ) - { - zS = zSrc + 1; - zT = zTgt - 1; - for ( ; zS < zT; ++zS, --zT ) // vertical loop on layers - { - for ( size_t iP = 0; iP < myBndColumns.size(); ++iP ) - { - toSrcBndPnts[ iP ] = bndPoint( iP, zS ); - toTgtBndPnts[ iP ] = bndPoint( iP, zT ); - } - // compensate bnd error - applyBoundaryError( toSrcBndPnts, bndError[ zS-1 ], bndError[ zS-1 ], 0.5, - intPntsOfLayer[ zS ], int2BndDist ); - applyBoundaryError( toTgtBndPnts, bndError[ zT+1 ], bndError[ zT+1 ], 0.5, - intPntsOfLayer[ zT ], int2BndDist ); - } - } - // cout << "centerIntErrorIsSmall = " << centerIntErrorIsSmall<< endl; - // cout << "bndErrorIsSmall = " << bndErrorIsSmall<< endl; + //cout << "centerIntErrorIsSmall = " << centerIntErrorIsSmall<< endl; // Create nodes for ( size_t iP = 0; iP < myIntColumns.size(); ++iP ) @@ -5229,7 +5188,7 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, for ( size_t z = zSrc + 1; z < zTgt; ++z ) // vertical loop on layers { const gp_XYZ & xyz = intPntsOfLayer[ z ][ iP ]; - if ( !( nodeCol[ z ] = helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() ))) + if ( !( nodeCol[ z ] = myHelper->AddNode( xyz.X(), xyz.Y(), xyz.Z() ))) return false; } } @@ -5301,7 +5260,7 @@ bool StdMeshers_Sweeper::CheckSameZ() */ //================================================================================ -bool StdMeshers_Sweeper::ComputeNodesOnStraightSameZ( SMESH_MesherHelper& helper ) +bool StdMeshers_Sweeper::ComputeNodesOnStraightSameZ() { TZColumn& z = myZColumns[0]; @@ -5313,7 +5272,7 @@ bool StdMeshers_Sweeper::ComputeNodesOnStraightSameZ( SMESH_MesherHelper& helper for ( size_t iZ = 0; iZ < z.size(); ++iZ ) { gp_XYZ p = n0 * ( 1 - z[iZ] ) + n1 * z[iZ]; - nodes[ iZ+1 ] = helper.AddNode( p.X(), p.Y(), p.Z() ); + nodes[ iZ+1 ] = myHelper->AddNode( p.X(), p.Y(), p.Z() ); } } @@ -5327,38 +5286,110 @@ bool StdMeshers_Sweeper::ComputeNodesOnStraightSameZ( SMESH_MesherHelper& helper */ //================================================================================ -bool StdMeshers_Sweeper::ComputeNodesOnStraight( SMESH_MesherHelper& helper, - const TopoDS_Face& botFace, - const TopoDS_Face& topFace ) +bool StdMeshers_Sweeper::ComputeNodesOnStraight() { - // get data to create a Morph - UVPtStructVec botUV( myBndColumns.size() + 1 ); - UVPtStructVec topUV( myBndColumns.size() + 1 ); + prepareTopBotDelaunay(); + + const SMDS_MeshNode *botNode, *topNode; + const BRepMesh_Triangle *topTria; + double botBC[3], topBC[3]; // barycentric coordinates + int botTriaNodes[3], topTriaNodes[3]; + bool checkUV = true; + + int nbInternalNodes = myIntColumns.size(); + myBotDelaunay->InitTraversal( nbInternalNodes ); + + while (( botNode = myBotDelaunay->NextNode( botBC, botTriaNodes ))) + { + TNodeColumn* column = myIntColumns[ myNodeID2ColID( botNode->GetID() )]; + + // find a Delaunay triangle containing the topNode + topNode = column->back(); + gp_XY topUV = myHelper->GetNodeUV( myTopFace, topNode, NULL, &checkUV ); + // get a starting triangle basing on that top and bot boundary nodes have same index + topTria = myTopDelaunay->GetTriangleNear( botTriaNodes[0] ); + topTria = myTopDelaunay->FindTriangle( topUV, topTria, topBC, topTriaNodes ); + if ( !topTria ) + return false; + + // create nodes along a line + SMESH_NodeXYZ botP( botNode ), topP( topNode ); + for ( size_t iZ = 0; iZ < myZColumns[0].size(); ++iZ ) + { + // use barycentric coordinates as weight of Z of boundary columns + double botZ = 0, topZ = 0; + for ( int i = 0; i < 3; ++i ) + { + botZ += botBC[i] * myZColumns[ botTriaNodes[i] ][ iZ ]; + topZ += topBC[i] * myZColumns[ topTriaNodes[i] ][ iZ ]; + } + double rZ = double( iZ + 1 ) / ( myZColumns[0].size() + 1 ); + double z = botZ * ( 1 - rZ ) + topZ * rZ; + gp_XYZ p = botP * ( 1 - z ) + topP * z; + (*column)[ iZ+1 ] = myHelper->AddNode( p.X(), p.Y(), p.Z() ); + } + } + + return myBotDelaunay->NbVisitedNodes() == nbInternalNodes; +} + +//================================================================================ +/*! + * \brief Compute Z of nodes of a straight column + */ +//================================================================================ + +void StdMeshers_Sweeper::fillZColumn( TZColumn& zColumn, + TNodeColumn& nodes ) +{ + if ( zColumn.size() == nodes.size() - 2 ) + return; + + gp_Pnt p0 = SMESH_NodeXYZ( nodes[0] ); + gp_Vec line( p0, SMESH_NodeXYZ( nodes.back() )); + double len2 = line.SquareMagnitude(); + + zColumn.resize( nodes.size() - 2 ); + for ( size_t i = 0; i < zColumn.size(); ++i ) + { + gp_Vec vec( p0, SMESH_NodeXYZ( nodes[ i+1] )); + zColumn[i] = ( line * vec ) / len2; // param [0,1] on the line + } +} + +//================================================================================ +/*! + * \brief Initialize *Delaunay members + */ +//================================================================================ + +void StdMeshers_Sweeper::prepareTopBotDelaunay() +{ + UVPtStructVec botUV( myBndColumns.size() ); + UVPtStructVec topUV( myBndColumns.size() ); for ( size_t i = 0; i < myBndColumns.size(); ++i ) { TNodeColumn& nodes = *myBndColumns[i]; botUV[i].node = nodes[0]; - botUV[i].SetUV( helper.GetNodeUV( botFace, nodes[0] )); + botUV[i].SetUV( myHelper->GetNodeUV( myBotFace, nodes[0] )); topUV[i].node = nodes.back(); - topUV[i].SetUV( helper.GetNodeUV( topFace, nodes.back() )); + topUV[i].SetUV( myHelper->GetNodeUV( myTopFace, nodes.back() )); botUV[i].node->setIsMarked( true ); } - botUV.back() = botUV[0]; - topUV.back() = topUV[0]; - TopoDS_Edge dummyE; - TSideVector botWires( 1, StdMeshers_FaceSide::New( botUV, botFace, dummyE, helper.GetMesh() )); - TSideVector topWires( 1, StdMeshers_FaceSide::New( topUV, topFace, dummyE, helper.GetMesh() )); + SMESH_Mesh* mesh = myHelper->GetMesh(); + TSideVector botWires( 1, StdMeshers_FaceSide::New( botUV, myBotFace, dummyE, mesh )); + TSideVector topWires( 1, StdMeshers_FaceSide::New( topUV, myTopFace, dummyE, mesh )); - // use Morph to make delauney mesh on the FACEs. Locating of a node within a - // delauney triangle will be used to get a weighted Z. - NSProjUtils::Morph botDelauney( botWires ); - NSProjUtils::Morph topDelauney( topWires ); + // Delaunay mesh on the FACEs. + bool checkUV = false; + myBotDelaunay.reset( new NSProjUtils::Delaunay( botWires, checkUV )); + myTopDelaunay.reset( new NSProjUtils::Delaunay( topWires, checkUV )); - if ( helper.GetIsQuadratic() ) + if ( myHelper->GetIsQuadratic() ) { // mark all medium nodes of faces on botFace to avoid their treating - SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( botFace ); + SMESHDS_SubMesh* smDS = myHelper->GetMeshDS()->MeshElements( myBotFace ); SMDS_ElemIteratorPtr eIt = smDS->GetElements(); while ( eIt->more() ) { @@ -5369,124 +5400,92 @@ bool StdMeshers_Sweeper::ComputeNodesOnStraight( SMESH_MesherHelper& helper, } // map to get a node column by a bottom node - TColStd_DataMapOfIntegerInteger iNode2iCol( myIntColumns.size() ); + myNodeID2ColID.Clear(/*doReleaseMemory=*/false); + myNodeID2ColID.ReSize( myIntColumns.size() ); - // un-mark nodes to treat (internal bottom nodes); later we will mark treated nodes + // un-mark nodes to treat (internal bottom nodes) to be returned by myBotDelaunay for ( size_t i = 0; i < myIntColumns.size(); ++i ) { const SMDS_MeshNode* botNode = myIntColumns[i]->front(); botNode->setIsMarked( false ); - iNode2iCol.Bind( botNode->GetID(), i ); + myNodeID2ColID.Bind( botNode->GetID(), i ); } +} - const int botFaceID = helper.GetMesh()->GetSubMesh( botFace )->GetId(); - const SMDS_MeshNode *botNode, *topNode; - const BRepMesh_Triangle *botTria, *topTria; - double botBC[3], topBC[3]; // barycentric coordinates - int botTriaNodes[3], topTriaNodes[3]; - bool checkUV = true; - - // a queue of bottom nodes with starting delauney triangles - NSProjUtils::Morph::TNodeTriaList botNoTriQueue; +//================================================================================ +/*! + * \brief For each internal node column, find Delaunay triangles including it + * and Barycentric Coordinates within the triangles. Fill in myTopBotTriangles + */ +//================================================================================ - size_t iBndN = 1; // index of a bottom boundary node - int nbNodesToProcess = myIntColumns.size(); - while ( nbNodesToProcess > 0 ) - { - while ( !botNoTriQueue.empty() ) // treat all nodes in the queue - { - botNode = botNoTriQueue.front().first; - botTria = botNoTriQueue.front().second; - botNoTriQueue.pop_front(); - if ( botNode->isMarked() ) - continue; - --nbNodesToProcess; - botNode->setIsMarked( true ); +bool StdMeshers_Sweeper::findDelaunayTriangles() +{ + const SMDS_MeshNode *botNode, *topNode; + const BRepMesh_Triangle *topTria; + TopBotTriangles tbTrias; + bool checkUV = true; - TNodeColumn* column = myIntColumns[ iNode2iCol( botNode->GetID() )]; + int nbInternalNodes = myIntColumns.size(); + myTopBotTriangles.resize( nbInternalNodes ); - // find a delauney triangle containing the botNode - gp_XY botUV = helper.GetNodeUV( botFace, botNode, NULL, &checkUV ); - botUV *= botDelauney.GetScale(); - botTria = botDelauney.FindTriangle( botUV, botTria, botBC, botTriaNodes ); - if ( !botTria ) - return false; + myBotDelaunay->InitTraversal( nbInternalNodes ); - // find a delauney triangle containing the topNode - topNode = column->back(); - gp_XY topUV = helper.GetNodeUV( topFace, topNode, NULL, &checkUV ); - topUV *= topDelauney.GetScale(); - // get a starting triangle basing on that top and bot boundary nodes have same index - topTria = topDelauney.GetTriangleNear( botTriaNodes[0] ); - topTria = topDelauney.FindTriangle( topUV, topTria, topBC, topTriaNodes ); - if ( !topTria ) - return false; + while (( botNode = myBotDelaunay->NextNode( tbTrias.myBotBC, tbTrias.myBotTriaNodes ))) + { + int colID = myNodeID2ColID( botNode->GetID() ); + TNodeColumn* column = myIntColumns[ colID ]; - // create nodes along a line - SMESH_NodeXYZ botP( botNode ), topP( topNode); - for ( size_t iZ = 0; iZ < myZColumns[0].size(); ++iZ ) - { - // use barycentric coordinates as weight of Z of boundary columns - double botZ = 0, topZ = 0; - for ( int i = 0; i < 3; ++i ) - { - botZ += botBC[i] * myZColumns[ botTriaNodes[i]-1 ][ iZ ]; - topZ += topBC[i] * myZColumns[ topTriaNodes[i]-1 ][ iZ ]; - } - double rZ = double( iZ + 1 ) / ( myZColumns[0].size() + 1 ); - double z = botZ * ( 1 - rZ ) + topZ * rZ; - gp_XYZ p = botP * ( 1 - z ) + topP * z; - (*column)[ iZ+1 ] = helper.AddNode( p.X(), p.Y(), p.Z() ); - } + // find a Delaunay triangle containing the topNode + topNode = column->back(); + gp_XY topUV = myHelper->GetNodeUV( myTopFace, topNode, NULL, &checkUV ); + // get a starting triangle basing on that top and bot boundary nodes have same index + topTria = myTopDelaunay->GetTriangleNear( tbTrias.myBotTriaNodes[0] ); + topTria = myTopDelaunay->FindTriangle( topUV, topTria, + tbTrias.myTopBC, tbTrias.myTopTriaNodes ); + if ( !topTria ) + tbTrias.SetTopByBottom(); - // add neighbor nodes to the queue - botDelauney.AddCloseNodes( botNode, botTria, botFaceID, botNoTriQueue ); - } + myTopBotTriangles[ colID ] = tbTrias; + } - if ( nbNodesToProcess > 0 ) // fill the queue - { - // assure that all bot nodes are visited - for ( ; iBndN-1 < myBndColumns.size() && botNoTriQueue.empty(); ++iBndN ) - { - botTria = botDelauney.GetTriangleNear( iBndN ); - const SMDS_MeshNode* bndNode = botDelauney.GetBndNodes()[ iBndN ]; - botDelauney.AddCloseNodes( bndNode, botTria, botFaceID, botNoTriQueue ); - } - if ( botNoTriQueue.empty() ) - { - for ( size_t i = 0; i < myIntColumns.size(); ++i ) - { - botNode = myIntColumns[i]->front(); - if ( !botNode->isMarked() ) - botNoTriQueue.push_back( make_pair( botNode, botTria )); - } - } - } + if ( myBotDelaunay->NbVisitedNodes() < nbInternalNodes ) + { + myTopBotTriangles.clear(); + return false; } + myBotDelaunay.reset(); + myTopDelaunay.reset(); + myNodeID2ColID.Clear(); + return true; } //================================================================================ /*! - * \brief Compute Z of nodes of a straight column + * \brief Initialize fields */ //================================================================================ -void StdMeshers_Sweeper::fillZColumn( TZColumn& zColumn, - TNodeColumn& nodes ) +StdMeshers_Sweeper::TopBotTriangles::TopBotTriangles() { - if ( zColumn.size() == nodes.size() - 2 ) - return; + myBotBC[0] = myBotBC[1] = myBotBC[2] = myTopBC[0] = myTopBC[1] = myTopBC[2] = 0.; + myBotTriaNodes[0] = myBotTriaNodes[1] = myBotTriaNodes[2] = 0; + myTopTriaNodes[0] = myTopTriaNodes[1] = myTopTriaNodes[2] = 0; +} - gp_Pnt p0 = SMESH_NodeXYZ( nodes[0] ); - gp_Vec line( p0, SMESH_NodeXYZ( nodes.back() )); - double len2 = line.SquareMagnitude(); +//================================================================================ +/*! + * \brief Set top data equal to bottom data + */ +//================================================================================ - zColumn.resize( nodes.size() - 2 ); - for ( size_t i = 0; i < zColumn.size(); ++i ) +void StdMeshers_Sweeper::TopBotTriangles::SetTopByBottom() +{ + for ( int i = 0; i < 3; ++i ) { - gp_Vec vec( p0, SMESH_NodeXYZ( nodes[ i+1] )); - zColumn[i] = ( line * vec ) / len2; // param [0,1] on the line + myTopBC[i] = myBotBC[i]; + myTopTriaNodes[i] = myBotTriaNodes[0]; } }