X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Prism_3D.cxx;h=f5807f09a192f848a28904cad3c1cbd6514297fe;hb=ad3cb4c93852dbc834d7075c087bbc749197454b;hp=0ee4041ce5ef74b936d21ad7b4fc5273cea105ff;hpb=f6b5d2f920970263bd4aa50e6ef7116d5c1b1625;p=modules%2Fsmesh.git diff --git a/src/StdMeshers/StdMeshers_Prism_3D.cxx b/src/StdMeshers/StdMeshers_Prism_3D.cxx index 0ee4041ce..f5807f09a 100644 --- a/src/StdMeshers/StdMeshers_Prism_3D.cxx +++ b/src/StdMeshers/StdMeshers_Prism_3D.cxx @@ -111,6 +111,7 @@ namespace { algo->myProxyMesh.reset( new SMESH_ProxyMesh( *helper->GetMesh() )); algo->myQuadList.clear(); + algo->myHelper = 0; if ( helper ) algo->_quadraticMesh = helper->GetIsQuadratic(); @@ -1182,7 +1183,7 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) // 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 ) { @@ -1193,16 +1194,15 @@ 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 ); @@ -1243,6 +1243,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() ) { @@ -1627,7 +1635,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 @@ -1790,9 +1798,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() )); } } } @@ -2283,7 +2290,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; } @@ -2582,7 +2589,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 //======================================================================= @@ -3055,7 +3062,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; @@ -3416,7 +3423,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(); @@ -3590,10 +3597,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 @@ -4839,7 +4846,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(); @@ -4893,6 +4900,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) { @@ -4917,43 +4925,23 @@ 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; } //================================================================================ @@ -4970,7 +4958,7 @@ bool StdMeshers_Sweeper::ComputeNodesByTrsf( const double tol, 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() ); @@ -4980,7 +4968,11 @@ bool StdMeshers_Sweeper::ComputeNodesByTrsf( const double tol, 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 ); @@ -5006,10 +4998,12 @@ bool StdMeshers_Sweeper::ComputeNodesByTrsf( const double tol, } 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; @@ -5034,6 +5028,22 @@ bool StdMeshers_Sweeper::ComputeNodesByTrsf( const double tol, 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 @@ -5048,10 +5058,12 @@ bool StdMeshers_Sweeper::ComputeNodesByTrsf( const double tol, } 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; @@ -5070,24 +5082,7 @@ bool StdMeshers_Sweeper::ComputeNodesByTrsf( const double tol, (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 ) { @@ -5096,11 +5091,6 @@ bool StdMeshers_Sweeper::ComputeNodesByTrsf( const double tol, 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 { @@ -5111,17 +5101,8 @@ bool StdMeshers_Sweeper::ComputeNodesByTrsf( const double tol, 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 @@ -5153,9 +5134,11 @@ bool StdMeshers_Sweeper::ComputeNodesByTrsf( const double tol, } 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 ) @@ -5186,15 +5169,6 @@ bool StdMeshers_Sweeper::ComputeNodesByTrsf( const double tol, 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 ); @@ -5202,27 +5176,8 @@ bool StdMeshers_Sweeper::ComputeNodesByTrsf( const double tol, } } // 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 ) @@ -5356,7 +5311,7 @@ bool StdMeshers_Sweeper::ComputeNodesOnStraight() return false; // create nodes along a line - SMESH_NodeXYZ botP( botNode ), topP( topNode); + 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 @@ -5454,3 +5409,81 @@ void StdMeshers_Sweeper::prepareTopBotDelaunay() myNodeID2ColID.Bind( botNode->GetID(), i ); } } + +//================================================================================ +/*! + * \brief For each internal node column, find Delaunay triangles including it + * and Barycentric Coordinates within the triangles. Fill in myTopBotTriangles + */ +//================================================================================ + +bool StdMeshers_Sweeper::findDelaunayTriangles() +{ + const SMDS_MeshNode *botNode, *topNode; + const BRepMesh_Triangle *topTria; + TopBotTriangles tbTrias; + bool checkUV = true; + + int nbInternalNodes = myIntColumns.size(); + myTopBotTriangles.resize( nbInternalNodes ); + + myBotDelaunay->InitTraversal( nbInternalNodes ); + + while (( botNode = myBotDelaunay->NextNode( tbTrias.myBotBC, tbTrias.myBotTriaNodes ))) + { + int colID = myNodeID2ColID( botNode->GetID() ); + TNodeColumn* column = myIntColumns[ colID ]; + + // 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(); + + myTopBotTriangles[ colID ] = tbTrias; + } + + if ( myBotDelaunay->NbVisitedNodes() < nbInternalNodes ) + { + myTopBotTriangles.clear(); + return false; + } + + myBotDelaunay.reset(); + myTopDelaunay.reset(); + myNodeID2ColID.Clear(); + + return true; +} + +//================================================================================ +/*! + * \brief Initialize fields + */ +//================================================================================ + +StdMeshers_Sweeper::TopBotTriangles::TopBotTriangles() +{ + 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; +} + +//================================================================================ +/*! + * \brief Set top data equal to bottom data + */ +//================================================================================ + +void StdMeshers_Sweeper::TopBotTriangles::SetTopByBottom() +{ + for ( int i = 0; i < 3; ++i ) + { + myTopBC[i] = myBotBC[i]; + myTopTriaNodes[i] = myBotTriaNodes[0]; + } +}