X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Prism_3D.cxx;h=8b7347d7b9702cc27cd8b02b523e7da47039e776;hp=62f538d84d72e927508e3add571619d4ba4d557f;hb=3da8fefe9c957f4538e9eacf013ce678df4d6c91;hpb=17d0c28bc8421608e7ededd4db97a0e0486dcf73 diff --git a/src/StdMeshers/StdMeshers_Prism_3D.cxx b/src/StdMeshers/StdMeshers_Prism_3D.cxx index 62f538d84..8b7347d7b 100644 --- a/src/StdMeshers/StdMeshers_Prism_3D.cxx +++ b/src/StdMeshers/StdMeshers_Prism_3D.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2014 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2015 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 @@ -49,6 +49,7 @@ #include #include #include +#include #include #include #include @@ -61,6 +62,7 @@ #include #include +#include using namespace std; @@ -70,8 +72,7 @@ using namespace std; #ifdef _DEBUG_ #define DBGOUT(msg) //cout << msg << endl; #define SHOWYXZ(msg, xyz) \ - // { gp_Pnt p (xyz); \ - // cout << msg << " ("<< p.X() << "; " <= columns.size() ) + if ( xColumn < 0 || xColumn >= (int) columns.size() ) { // select a column for X dir double maxDist = 0; - for ( int i = 0; i < columns.size(); ++i ) + for ( size_t i = 0; i < columns.size(); ++i ) { double dist = ( O - gpXYZ((*columns[ i ])[ z ])).SquareModulus(); if ( dist > maxDist ) @@ -407,8 +408,8 @@ namespace { if ( nbQuads > 0 ) toRemove = helper->IsStructured( faceSm ); else - toRemove = quadAlgo->CheckNbEdges( *helper->GetMesh(), - faceSm->GetSubShape() ); + toRemove = ( quadAlgo->CheckNbEdges( *helper->GetMesh(), + faceSm->GetSubShape() ) != NULL ); nbRemoved += toRemove; if ( toRemove ) smIt = notQuadSubMesh.erase( smIt ); @@ -454,9 +455,9 @@ namespace { std::advance( edgeIt, nbEdges-1 ); TopoDS_Edge prevE = *edgeIt; // bool isPrevStraight = SMESH_Algo::IsStraight( prevE ); - int iPrev = nbEdges - 1; + // int iPrev = nbEdges - 1; - int iUnite = -1; // the first of united EDGEs + // int iUnite = -1; // the first of united EDGEs // analyse angles between EDGEs int nbCorners = 0; @@ -524,7 +525,7 @@ namespace { void pointsToPython(const std::vector& p) { #ifdef _DEBUG_ - for ( int i = SMESH_Block::ID_V000; i < p.size(); ++i ) + for ( size_t i = SMESH_Block::ID_V000; i < p.size(); ++i ) { cout << "mesh.AddNode( " << p[i].X() << ", "<< p[i].Y() << ", "<< p[i].Z() << ") # " << i <<" " ; SMESH_Block::DumpShapeID( i, cout ) << endl; @@ -661,12 +662,14 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh Prism_3D::TPrismTopo prism; myPropagChains = 0; + bool selectBottom = meshedFaces.empty(); if ( nbSolids == 1 ) { + TopoDS_Shape solid = TopExp_Explorer( theShape, TopAbs_SOLID ).Current(); if ( !meshedFaces.empty() ) prism.myBottom = meshedFaces.front(); - return ( initPrism( prism, TopExp_Explorer( theShape, TopAbs_SOLID ).Current() ) && + return ( initPrism( prism, solid, selectBottom ) && compute( prism )); } @@ -709,7 +712,7 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh { prism.Clear(); prism.myBottom = face; - if ( !initPrism( prism, solid ) || + if ( !initPrism( prism, solid, selectBottom ) || !compute( prism )) return false; @@ -762,15 +765,24 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh while ( const TopoDS_Shape* f = faceIt->next() ) { const TopoDS_Face& candidateF = TopoDS::Face( *f ); + if ( candidateF.IsSame( wFace )) continue; + // select a source FACE: prismIt->myBottom or prismIt->myTop + TopoDS_Face sourceF = prismIt->myBottom; + for ( TopExp_Explorer v( prismIt->myTop, TopAbs_VERTEX ); v.More(); v.Next() ) + if ( myHelper->IsSubShape( v.Current(), candidateF )) { + sourceF = prismIt->myTop; + break; + } prism.Clear(); - prism.myBottom = candidateF; + prism.myBottom = candidateF; mySetErrorToSM = false; if ( !myHelper->IsSubShape( candidateF, prismIt->myShape3D ) && - myHelper->IsSubShape( candidateF, solid ) && + myHelper ->IsSubShape( candidateF, solid ) && !myHelper->GetMesh()->GetSubMesh( candidateF )->IsMeshComputed() && - initPrism( prism, solid ) && + initPrism( prism, solid, /*selectBottom=*/false ) && !myHelper->GetMesh()->GetSubMesh( prism.myTop )->IsMeshComputed() && - project2dMesh( prismIt->myBottom, candidateF)) + !myHelper->GetMesh()->GetSubMesh( prism.myBottom )->IsMeshComputed() && + project2dMesh( sourceF, prism.myBottom )) { mySetErrorToSM = true; if ( !compute( prism )) @@ -780,6 +792,7 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh { meshedFaces.push_front( prism.myTop ); meshedFaces.push_front( prism.myBottom ); + selectBottom = false; } meshedPrism.push_back( prism ); meshedSolids.Add( solid ); @@ -802,6 +815,7 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh if ( meshedFaces.empty() ) { meshedFaces.splice( meshedFaces.end(), suspectSourceFaces ); + selectBottom = true; } // find FACEs with local 1D hyps, which has to be computed by now, @@ -822,6 +836,7 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh { if ( !meshedFaces.empty() ) meshedFaces.pop_back(); meshedFaces.push_back( face ); // lower priority + selectBottom = true; prevNbFaces = nbFaces; } } @@ -836,6 +851,7 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh faceSM->ComputeStateEngine( SMESH_subMesh::COMPUTE ); if ( !faceSM->IsEmpty() ) { meshedFaces.push_front( face ); // higher priority + selectBottom = true; break; } else { @@ -867,6 +883,7 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh meshedFaces.push_front( prism.myBottom ); meshedPrism.push_back( prism ); meshedSolids.Add( solid.Current() ); + selectBottom = true; } mySetErrorToSM = true; } @@ -888,7 +905,7 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh return error( err ); } } - return true; + return error( COMPERR_OK ); } //================================================================================ @@ -917,6 +934,7 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism, list< TopoDS_Edge >::iterator edge = thePrism.myBottomEdges.begin(); std::list< int >::iterator nbE = thePrism.myNbEdgesInWires.begin(); + std::list< int > nbQuadsPerWire; int iE = 0; double f,l; while ( edge != thePrism.myBottomEdges.end() ) @@ -960,6 +978,8 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism, { iE = 0; ++nbE; + int nbQuadPrev = std::accumulate( nbQuadsPerWire.begin(), nbQuadsPerWire.end(), 0 ); + nbQuadsPerWire.push_back( thePrism.myWallQuads.size() - nbQuadPrev ); } } @@ -971,12 +991,14 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism, // that is not so evident in case of several WIREs in the bottom FACE thePrism.myRightQuadIndex.clear(); for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i ) - thePrism.myRightQuadIndex.push_back( i+1 ); - list< int >::iterator nbEinW = thePrism.myNbEdgesInWires.begin(); - for ( int iLeft = 0; nbEinW != thePrism.myNbEdgesInWires.end(); ++nbEinW ) { - thePrism.myRightQuadIndex[ iLeft + *nbEinW - 1 ] = iLeft; // 1st EDGE index of a current WIRE - iLeft += *nbEinW; + thePrism.myRightQuadIndex.push_back( i+1 ); // OK for all but the last EDGE of a WIRE + } + list< int >::iterator nbQinW = nbQuadsPerWire.begin(); + for ( int iLeft = 0; nbQinW != nbQuadsPerWire.end(); ++nbQinW ) + { + thePrism.myRightQuadIndex[ iLeft + *nbQinW - 1 ] = iLeft; // for the last EDGE of a WIRE + iLeft += *nbQinW; } while ( totalNbFaces - faceMap.Extent() > 2 ) @@ -1056,7 +1078,7 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism, { // now only top and bottom FACEs are not in the faceMap faceMap.Add( thePrism.myBottom ); - for ( TopExp_Explorer f( thePrism.myShape3D, TopAbs_FACE );f.More(); f.Next() ) + for ( TopExp_Explorer f( thePrism.myShape3D, TopAbs_FACE ); f.More(); f.Next() ) if ( !faceMap.Contains( f.Current() )) { thePrism.myTop = TopoDS::Face( f.Current() ); break; @@ -1071,7 +1093,7 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism, StdMeshers_FaceSidePtr topSide = thePrism.myWallQuads[i].back()->side[ QUAD_TOP_SIDE ]; const TopoDS_Edge & topE = topSide->Edge( 0 ); if ( !myHelper->IsSubShape( topE, thePrism.myTop )) - return toSM( error( TCom("Wrong source face (#") << shapeID( thePrism.myBottom ))); + return toSM( error( TCom("Wrong source face: #") << shapeID( thePrism.myBottom ))); } return true; @@ -1088,13 +1110,23 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) if ( _computeCanceled ) return toSM( error( SMESH_ComputeError::New(COMPERR_CANCELED))); + // Assure the bottom is meshed + SMESH_subMesh * botSM = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom ); + if (( botSM->IsEmpty() ) && + ( ! botSM->GetAlgo() || + ! _gen->Compute( *botSM->GetFather(), botSM->GetSubShape(), /*shapeOnly=*/true ))) + return error( COMPERR_BAD_INPUT_MESH, + TCom( "No mesher defined to compute the face #") + << shapeID( thePrism.myBottom )); + // Make all side FACEs of thePrism meshed with quads if ( !computeWalls( thePrism )) return false; // Analyse mesh and geometry to find all block sub-shapes and submeshes - // (after fixing IPAL52499 myBlock is used only as a holder of boundary nodes - // and location of internal nodes is computed by StdMeshers_Sweeper) + // (after fixing IPAL52499 myBlock is used as a holder of boundary nodes + // and for 2D projection in hard cases where StdMeshers_Projection_2D fails; + // location of internal nodes is usually computed by StdMeshers_Sweeper) if ( !myBlock.Init( myHelper, thePrism )) return toSM( error( myBlock.GetError())); @@ -1123,6 +1155,7 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) // Projections on the top and bottom faces are taken from nodes existing // on these faces; find correspondence between bottom and top nodes + myUseBlock = false; myBotToColumnMap.clear(); if ( !assocOrProjBottom2Top( bottomToTopTrsf, thePrism ) ) // it also fills myBotToColumnMap return false; @@ -1132,27 +1165,33 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) // use transformation (issue 0020680, IPAL0052499) StdMeshers_Sweeper sweeper; + double tol; + bool allowHighBndError; - // load boundary nodes - bool dummy; - list< TopoDS_Edge >::const_iterator edge = thePrism.myBottomEdges.begin(); - for ( ; edge != thePrism.myBottomEdges.end(); ++edge ) + if ( !myUseBlock ) { - int edgeID = meshDS->ShapeToIndex( *edge ); - TParam2ColumnMap* u2col = const_cast - ( myBlock.GetParam2ColumnMap( edgeID, dummy )); - TParam2ColumnMap::iterator u2colIt = u2col->begin(); - for ( ; u2colIt != u2col->end(); ++u2colIt ) - sweeper.myBndColumns.push_back( & u2colIt->second ); - } - // load node columns inside the bottom face - TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin(); - for ( ; bot_column != myBotToColumnMap.end(); ++bot_column ) - sweeper.myIntColumns.push_back( & bot_column->second ); + // load boundary nodes into sweeper + bool dummy; + list< TopoDS_Edge >::const_iterator edge = thePrism.myBottomEdges.begin(); + for ( ; edge != thePrism.myBottomEdges.end(); ++edge ) + { + int edgeID = meshDS->ShapeToIndex( *edge ); + TParam2ColumnMap* u2col = const_cast + ( myBlock.GetParam2ColumnMap( edgeID, dummy )); + TParam2ColumnMap::iterator u2colIt = u2col->begin(); + for ( ; u2colIt != u2col->end(); ++u2colIt ) + sweeper.myBndColumns.push_back( & u2colIt->second ); + } + // load node columns inside the bottom face + TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin(); + for ( ; bot_column != myBotToColumnMap.end(); ++bot_column ) + sweeper.myIntColumns.push_back( & bot_column->second ); - const double tol = getSweepTolerance( thePrism ); + tol = getSweepTolerance( thePrism ); + allowHighBndError = !isSimpleBottom( thePrism ); + } - if ( sweeper.ComputeNodes( *myHelper, tol )) + if ( !myUseBlock && sweeper.ComputeNodes( *myHelper, tol, allowHighBndError )) { } else // use block approach @@ -1253,6 +1292,7 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) if ( !smDS ) return toSM( error(COMPERR_BAD_INPUT_MESH, "Null submesh")); // loop on bottom mesh faces + vector< const TNodeColumn* > columns; SMDS_ElemIteratorPtr faceIt = smDS->GetElements(); while ( faceIt->more() ) { @@ -1262,7 +1302,7 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) // find node columns for each node int nbNodes = face->NbCornerNodes(); - vector< const TNodeColumn* > columns( nbNodes ); + columns.resize( nbNodes ); for ( int i = 0; i < nbNodes; ++i ) { const SMDS_MeshNode* n = face->GetNode( i ); @@ -1279,14 +1319,25 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) } } // create prisms - AddPrisms( columns, myHelper ); + if ( !AddPrisms( columns, myHelper )) + return toSM( error("Different 'vertical' discretization")); } // loop on bottom mesh faces // clear data myBotToColumnMap.clear(); myBlock.Clear(); - + + // update state of sub-meshes (mostly in order to erase improper errors) + SMESH_subMesh* sm = myHelper->GetMesh()->GetSubMesh( thePrism.myShape3D ); + SMESH_subMeshIteratorPtr smIt = sm->getDependsOnIterator(/*includeSelf=*/false); + while ( smIt->more() ) + { + sm = smIt->next(); + sm->GetComputeError().reset(); + sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); + } + return true; } @@ -1528,7 +1579,7 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism) for ( ; quad != thePrism.myWallQuads[iW].end(); ++quad ) { const TopoDS_Face& face = (*quad)->face; - SMESH_subMesh* fSM = mesh->GetSubMesh( face ); + SMESH_subMesh* fSM = mesh->GetSubMesh( face ); if ( ! fSM->IsMeshComputed() ) { // Top EDGEs must be projections from the bottom ones @@ -1767,17 +1818,20 @@ bool StdMeshers_Prism_3D::Evaluate(SMESH_Mesh& theMesh, */ //================================================================================ -void StdMeshers_Prism_3D::AddPrisms( vector & columns, +bool StdMeshers_Prism_3D::AddPrisms( vector & columns, SMESH_MesherHelper* helper) { - int nbNodes = columns.size(); - int nbZ = columns[0]->size(); - if ( nbZ < 2 ) return; + size_t nbNodes = columns.size(); + size_t nbZ = columns[0]->size(); + if ( nbZ < 2 ) return false; + for ( size_t i = 1; i < nbNodes; ++i ) + if ( columns[i]->size() != nbZ ) + return false; // find out orientation bool isForward = true; SMDS_VolumeTool vTool; - int z = 1; + size_t z = 1; switch ( nbNodes ) { case 3: { SMDS_VolumeOfNodes tmpPenta ( (*columns[0])[z-1], // bottom @@ -1863,7 +1917,7 @@ void StdMeshers_Prism_3D::AddPrisms( vector & columns, vector nodes( 2*nbNodes + 4*nbNodes); for ( z = 1; z < nbZ; ++z ) { - for ( int i = 0; i < nbNodes; ++i ) { + for ( size_t i = 0; i < nbNodes; ++i ) { nodes[ i ] = (*columns[ i ])[z+iBase1]; // bottom or top nodes[ 2*nbNodes-i-1 ] = (*columns[ i ])[z+iBase2]; // top or bottom // side @@ -1877,6 +1931,8 @@ void StdMeshers_Prism_3D::AddPrisms( vector & columns, } } // switch ( nbNodes ) + + return true; } //================================================================================ @@ -1932,14 +1988,14 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf n2nMapPtr = & TProjction2dAlgo::instance( this )->GetNodesMap(); } - if ( !n2nMapPtr || n2nMapPtr->size() < botSMDS->NbNodes() ) + if ( !n2nMapPtr || (int) n2nMapPtr->size() < botSMDS->NbNodes() ) { // associate top and bottom faces NSProjUtils::TShapeShapeMap shape2ShapeMap; const bool sameTopo = NSProjUtils::FindSubShapeAssociation( thePrism.myBottom, myHelper->GetMesh(), - thePrism.myTop, myHelper->GetMesh(), - shape2ShapeMap); + thePrism.myTop, myHelper->GetMesh(), + shape2ShapeMap); if ( !sameTopo ) for ( size_t iQ = 0; iQ < thePrism.myWallQuads.size(); ++iQ ) { @@ -1951,10 +2007,10 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf for ( int iE = 0; iE < botSide->NbEdges(); ++iE ) { NSProjUtils::InsertAssociation( botSide->Edge( iE ), - topSide->Edge( iE ), shape2ShapeMap ); + topSide->Edge( iE ), shape2ShapeMap ); NSProjUtils::InsertAssociation( myHelper->IthVertex( 0, botSide->Edge( iE )), - myHelper->IthVertex( 0, topSide->Edge( iE )), - shape2ShapeMap ); + myHelper->IthVertex( 0, topSide->Edge( iE )), + shape2ShapeMap ); } } else @@ -1969,7 +2025,7 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf vt.IsSame( sideT->LastVertex() )) { NSProjUtils::InsertAssociation( botSide->Edge( 0 ), - topSide->Edge( 0 ), shape2ShapeMap ); + topSide->Edge( 0 ), shape2ShapeMap ); NSProjUtils::InsertAssociation( vb, vt, shape2ShapeMap ); } vb = myHelper->IthVertex( 1, botSide->Edge( botSide->NbEdges()-1 )); @@ -1980,8 +2036,8 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf vt.IsSame( sideT->LastVertex() )) { NSProjUtils::InsertAssociation( botSide->Edge( botSide->NbEdges()-1 ), - topSide->Edge( topSide->NbEdges()-1 ), - shape2ShapeMap ); + topSide->Edge( topSide->NbEdges()-1 ), + shape2ShapeMap ); NSProjUtils::InsertAssociation( vb, vt, shape2ShapeMap ); } } @@ -1990,8 +2046,8 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf // Find matching nodes of top and bottom faces n2nMapPtr = & n2nMap; if ( ! NSProjUtils::FindMatchingNodesOnFaces( thePrism.myBottom, myHelper->GetMesh(), - thePrism.myTop, myHelper->GetMesh(), - shape2ShapeMap, n2nMap )) + thePrism.myTop, myHelper->GetMesh(), + shape2ShapeMap, n2nMap )) { if ( sameTopo ) return toSM( error(TCom("Mesh on faces #") << botSM->GetId() @@ -2014,7 +2070,7 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf continue; // wall columns are contained in myBlock // create node column Prism_3D::TNode bN( botNode ); - TNode2ColumnMap::iterator bN_col = + TNode2ColumnMap::iterator bN_col = myBotToColumnMap.insert( make_pair ( bN, TNodeColumn() )).first; TNodeColumn & column = bN_col->second; column.resize( zSize ); @@ -2038,6 +2094,11 @@ bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottom { return true; } + NSProjUtils::TNodeNodeMap& n2nMap = + (NSProjUtils::TNodeNodeMap&) TProjction2dAlgo::instance( this )->GetNodesMap(); + n2nMap.clear(); + + myUseBlock = true; SMESHDS_Mesh* meshDS = myHelper->GetMeshDS(); SMESH_subMesh * botSM = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom ); @@ -2047,7 +2108,13 @@ bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottom SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS(); if ( topSMDS && topSMDS->NbElements() > 0 ) - topSM->ComputeStateEngine( SMESH_subMesh::CLEAN ); + { + //topSM->ComputeStateEngine( SMESH_subMesh::CLEAN ); -- avoid propagation of events + for ( SMDS_ElemIteratorPtr eIt = topSMDS->GetElements(); eIt->more(); ) + meshDS->RemoveFreeElement( eIt->next(), topSMDS, /*fromGroups=*/false ); + for ( SMDS_NodeIteratorPtr nIt = topSMDS->GetNodes(); nIt->more(); ) + meshDS->RemoveFreeNode( nIt->next(), topSMDS, /*fromGroups=*/false ); + } const TopoDS_Face& botFace = thePrism.myBottom; // oriented within const TopoDS_Face& topFace = thePrism.myTop; // the 3D SHAPE @@ -2119,6 +2186,8 @@ bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottom column.front() = botNode; column.back() = topNode; + n2nMap.insert( n2nMap.end(), make_pair( botNode, topNode )); + if ( _computeCanceled ) return toSM( error( SMESH_ComputeError::New(COMPERR_CANCELED))); } @@ -2168,11 +2237,11 @@ bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottom case 3: { newFace = myHelper->AddFace(nodes[0], nodes[1], nodes[2]); break; - } + } case 4: { newFace = myHelper->AddFace( nodes[0], nodes[1], nodes[2], nodes[3] ); break; - } + } default: newFace = meshDS->AddPolygonalFace( nodes ); } @@ -2180,7 +2249,42 @@ bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottom meshDS->SetMeshElementOnShape( newFace, topFaceID ); } - myHelper->SetElementsOnShape( oldSetElemsOnShape ); + myHelper->SetElementsOnShape( oldSetElemsOnShape ); + + // Check the projected mesh + + if ( thePrism.myNbEdgesInWires.size() > 1 && // there are holes + topHelper.IsDistorted2D( topSM, /*checkUV=*/false )) + { + SMESH_MeshEditor editor( topHelper.GetMesh() ); + + // smooth in 2D or 3D? + TopLoc_Location loc; + Handle(Geom_Surface) surface = BRep_Tool::Surface( topFace, loc ); + bool isPlanar = GeomLib_IsPlanarSurface( surface ).IsPlanar(); + + bool isFixed = false; + set fixedNodes; + for ( int iAttemp = 0; !isFixed && iAttemp < 10; ++iAttemp ) + { + TIDSortedElemSet faces; + for ( faceIt = topSMDS->GetElements(); faceIt->more(); ) + faces.insert( faces.end(), faceIt->next() ); + + SMESH_MeshEditor::SmoothMethod algo = + iAttemp ? SMESH_MeshEditor::CENTROIDAL : SMESH_MeshEditor::LAPLACIAN; + + // smoothing + editor.Smooth( faces, fixedNodes, algo, /*nbIterations=*/ 10, + /*theTgtAspectRatio=*/1.0, /*the2D=*/!isPlanar); + + isFixed = !topHelper.IsDistorted2D( topSM, /*checkUV=*/true ); + } + if ( !isFixed ) + return toSM( error( TCom("Projection from face #") << botSM->GetId() + << " to face #" << topSM->GetId() + << " failed: inverted elements created")); + } return true; } @@ -2256,6 +2360,44 @@ double StdMeshers_Prism_3D::getSweepTolerance( const Prism_3D::TPrismTopo& thePr return 0.1 * Sqrt ( minDist ); } +//======================================================================= +//function : isSimpleQuad +//purpose : check if the bottom FACE is meshable with nice qudrangles, +// if so the block aproach can work rather fast. +// This is a temporary mean caused by problems in StdMeshers_Sweeper +//======================================================================= + +bool StdMeshers_Prism_3D::isSimpleBottom( const Prism_3D::TPrismTopo& thePrism ) +{ + // analyse angles between edges + double nbConcaveAng = 0, nbConvexAng = 0; + TopoDS_Face reverseBottom = TopoDS::Face( thePrism.myBottom.Reversed() ); // see initPrism() + TopoDS_Vertex commonV; + const list< TopoDS_Edge >& botEdges = thePrism.myBottomEdges; + list< TopoDS_Edge >::const_iterator edge = botEdges.begin(); + while ( edge != botEdges.end() ) + { + if ( SMESH_Algo::isDegenerated( *edge )) + return false; + TopoDS_Edge e1 = *edge++; + TopoDS_Edge e2 = ( edge == botEdges.end() ? botEdges.front() : *edge ); + if ( ! TopExp::CommonVertex( e1, e2, commonV )) + { + e2 = botEdges.front(); + if ( ! TopExp::CommonVertex( e1, e2, commonV )) + break; + } + double angle = myHelper->GetAngle( e1, e2, reverseBottom, commonV ); + if ( angle < -5 * M_PI/180 ) + if ( ++nbConcaveAng > 1 ) + return false; + if ( angle > 85 * M_PI/180 ) + if ( ++nbConvexAng > 4 ) + return false; + } + return true; +} + //======================================================================= //function : project2dMesh //purpose : Project mesh faces from a source FACE of one prism (theSrcFace) @@ -2270,6 +2412,15 @@ bool StdMeshers_Prism_3D::project2dMesh(const TopoDS_Face& theSrcFace, bool ok = projector2D->Compute( *myHelper->GetMesh(), theTgtFace ); SMESH_subMesh* tgtSM = myHelper->GetMesh()->GetSubMesh( theTgtFace ); + if ( !ok && tgtSM->GetSubMeshDS() ) { + //tgtSM->ComputeStateEngine( SMESH_subMesh::CLEAN ); -- avoid propagation of events + SMESHDS_Mesh* meshDS = myHelper->GetMeshDS(); + SMESHDS_SubMesh* tgtSMDS = tgtSM->GetSubMeshDS(); + for ( SMDS_ElemIteratorPtr eIt = tgtSMDS->GetElements(); eIt->more(); ) + meshDS->RemoveFreeElement( eIt->next(), tgtSMDS, /*fromGroups=*/false ); + for ( SMDS_NodeIteratorPtr nIt = tgtSMDS->GetNodes(); nIt->more(); ) + meshDS->RemoveFreeNode( nIt->next(), tgtSMDS, /*fromGroups=*/false ); + } tgtSM->ComputeStateEngine ( SMESH_subMesh::CHECK_COMPUTE_STATE ); tgtSM->ComputeSubMeshStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); @@ -2361,17 +2512,20 @@ namespace // utils used by StdMeshers_Prism_3D::IsApplicable() } EdgeWithNeighbors() {} }; - struct PrismSide + // PrismSide contains all FACEs linking a bottom EDGE with a top one. + struct PrismSide { - TopoDS_Face _face; - TopTools_IndexedMapOfShape *_faces; // pointer because its copy constructor is private - TopoDS_Edge _topEdge; - vector< EdgeWithNeighbors >*_edges; - int _iBotEdge; - vector< bool > _isCheckedEdge; + TopoDS_Face _face; // a currently treated upper FACE + TopTools_IndexedMapOfShape *_faces; // all FACEs (pointer because of a private copy constructor) + TopoDS_Edge _topEdge; // a current top EDGE + vector< EdgeWithNeighbors >*_edges; // all EDGEs of _face + int _iBotEdge; // index of _topEdge within _edges + vector< bool > _isCheckedEdge; // mark EDGEs whose two owner FACEs found int _nbCheckedEdges; // nb of EDGEs whose location is defined - PrismSide *_leftSide; + PrismSide *_leftSide; // neighbor sides PrismSide *_rightSide; + void SetExcluded() { _leftSide = _rightSide = NULL; } + bool IsExcluded() const { return !_leftSide; } const TopoDS_Edge& Edge( int i ) const { return (*_edges)[ i ]._edge; @@ -2410,7 +2564,7 @@ namespace // utils used by StdMeshers_Prism_3D::IsApplicable() for ( iE = 0; iE < *nbE; ++e, ++iE ) if ( SMESH_Algo::isDegenerated( *e )) { - ee.erase( e ); + e = --ee.erase( e ); --(*nbE); --iE; } @@ -2433,16 +2587,34 @@ namespace // utils used by StdMeshers_Prism_3D::IsApplicable() /*! * \brief Return another faces sharing an edge */ - const TopoDS_Shape & getAnotherFace( const TopoDS_Face& face, - const TopoDS_Edge& edge, - TopTools_IndexedDataMapOfShapeListOfShape& facesOfEdge) + const TopoDS_Face & getAnotherFace( const TopoDS_Face& face, + const TopoDS_Edge& edge, + TopTools_IndexedDataMapOfShapeListOfShape& facesOfEdge) { TopTools_ListIteratorOfListOfShape faceIt( facesOfEdge.FindFromKey( edge )); for ( ; faceIt.More(); faceIt.Next() ) if ( !face.IsSame( faceIt.Value() )) - return faceIt.Value(); + return TopoDS::Face( faceIt.Value() ); return face; } + + //-------------------------------------------------------------------------------- + /*! + * \brief Return number of faces sharing given edges + */ + int nbAdjacentFaces( const std::vector< EdgeWithNeighbors >& edges, + const TopTools_IndexedDataMapOfShapeListOfShape& facesOfEdge ) + { + TopTools_MapOfShape adjFaces; + + for ( size_t i = 0; i < edges.size(); ++i ) + { + TopTools_ListIteratorOfListOfShape faceIt( facesOfEdge.FindFromKey( edges[i]._edge )); + for ( ; faceIt.More(); faceIt.Next() ) + adjFaces.Add( faceIt.Value() ); + } + return adjFaces.Extent(); + } } //================================================================================ @@ -2508,11 +2680,13 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA typedef vector< EdgeWithNeighbors > TEdgeWithNeighborsVec; vector< TEdgeWithNeighborsVec > faceEdgesVec( allFaces.Extent() + 1 ); - TopTools_IndexedMapOfShape* facesOfSide = new TopTools_IndexedMapOfShape[ faceEdgesVec.size() ]; + const size_t nbEdgesMax = facesOfEdge.Extent() * 2; // there can be seam EDGEs + TopTools_IndexedMapOfShape* facesOfSide = new TopTools_IndexedMapOfShape[ nbEdgesMax ]; SMESHUtils::ArrayDeleter delFacesOfSide( facesOfSide ); // try to use each face as a bottom one bool prismDetected = false; + vector< PrismSide > sides; for ( int iF = 1; iF < allFaces.Extent() && !prismDetected; ++iF ) { const TopoDS_Face& botF = TopoDS::Face( allFaces( iF )); @@ -2525,11 +2699,12 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA continue; // all faces are adjacent to botF - no top FACE // init data of side FACEs - vector< PrismSide > sides( botEdges.size() ); - for ( int iS = 0; iS < botEdges.size(); ++iS ) + sides.clear(); + sides.resize( botEdges.size() ); + for ( size_t iS = 0; iS < botEdges.size(); ++iS ) { - sides[ iS ]._topEdge = botEdges[ iS ]._edge; - sides[ iS ]._face = botF; + sides[ iS ]._topEdge = botEdges[ iS ]._edge; + sides[ iS ]._face = botF; sides[ iS ]._leftSide = & sides[ botEdges[ iS ]._iR ]; sides[ iS ]._rightSide = & sides[ botEdges[ iS ]._iL ]; sides[ iS ]._faces = & facesOfSide[ iS ]; @@ -2561,8 +2736,8 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA if ( side._isCheckedEdge[ iE ] ) continue; const TopoDS_Edge& vertE = side.Edge( iE ); const TopoDS_Shape& neighborF = getAnotherFace( side._face, vertE, facesOfEdge ); - bool isEdgeShared = adjSide->IsSideFace( neighborF ); - if ( isEdgeShared ) + bool isEdgeShared = adjSide->IsSideFace( neighborF ); + if ( isEdgeShared ) // vertE is shared with adjSide { isAdvanced = true; side._isCheckedEdge[ iE ] = true; @@ -2641,6 +2816,10 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA side._isCheckedEdge[ side._iBotEdge ] = true; side._nbCheckedEdges = 1; // bottom EDGE is known } + else // probably a triangular top face found + { + side._face.Nullify(); + } side._topEdge.Nullify(); isOK = ( !side._edges->empty() || side._faces->Extent() > 1 ); @@ -2673,7 +2852,7 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA const TopoDS_Shape& topFace = sides[0]._faces->FindKey( nbFaces ); size_t iS; for ( iS = 1; iS < sides.size(); ++iS ) - if ( !sides[ iS ]._faces->Contains( topFace )) + if ( ! sides[ iS ]._faces->Contains( topFace )) break; prismDetected = ( iS == sides.size() ); } @@ -2785,16 +2964,18 @@ void StdMeshers_PrismAsBlock::Clear() //======================================================================= bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism, - const TopoDS_Shape& shape3D) + const TopoDS_Shape& theShape3D, + const bool selectBottom) { - myHelper->SetSubShape( shape3D ); + myHelper->SetSubShape( theShape3D ); - SMESH_subMesh* mainSubMesh = myHelper->GetMesh()->GetSubMeshContaining( shape3D ); + SMESH_subMesh* mainSubMesh = myHelper->GetMesh()->GetSubMeshContaining( theShape3D ); if ( !mainSubMesh ) return toSM( error(COMPERR_BAD_INPUT_MESH,"Null submesh of shape3D")); // detect not-quad FACE sub-meshes of the 3D SHAPE list< SMESH_subMesh* > notQuadGeomSubMesh; list< SMESH_subMesh* > notQuadElemSubMesh; + list< SMESH_subMesh* > meshedSubMesh; int nbFaces = 0; // SMESH_subMesh* anyFaceSM = 0; @@ -2816,10 +2997,14 @@ bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism, if ( nbWires != 1 || nbEdgesInWires.front() != 4 ) notQuadGeomSubMesh.push_back( sm ); - // look for not quadrangle mesh elements - if ( SMESHDS_SubMesh* smDS = sm->GetSubMeshDS() ) - if ( !myHelper->IsSameElemGeometry( smDS, SMDSGeom_QUADRANGLE )) + // look for a not structured sub-mesh + if ( !sm->IsEmpty() ) + { + meshedSubMesh.push_back( sm ); + if ( !myHelper->IsSameElemGeometry( sm->GetSubMeshDS(), SMDSGeom_QUADRANGLE ) || + !myHelper->IsStructured ( sm )) notQuadElemSubMesh.push_back( sm ); + } } int nbNotQuadMeshed = notQuadElemSubMesh.size(); @@ -2884,33 +3069,50 @@ bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism, // use thePrism.myBottom if ( !thePrism.myBottom.IsNull() ) { - if ( botSM ) { + if ( botSM ) { // <-- not quad geom or mesh on botSM if ( ! botSM->GetSubShape().IsSame( thePrism.myBottom )) { std::swap( botSM, topSM ); - if ( !botSM || ! botSM->GetSubShape().IsSame( thePrism.myBottom )) - return toSM( error( COMPERR_BAD_INPUT_MESH, - "Incompatible non-structured sub-meshes")); + if ( !botSM || ! botSM->GetSubShape().IsSame( thePrism.myBottom )) { + if ( !selectBottom ) + return toSM( error( COMPERR_BAD_INPUT_MESH, + "Incompatible non-structured sub-meshes")); + std::swap( botSM, topSM ); + thePrism.myBottom = TopoDS::Face( botSM->GetSubShape() ); + } } } - else { + else if ( !selectBottom ) { botSM = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom ); } } - else if ( !botSM ) // find a proper bottom + if ( !botSM ) // find a proper bottom { - // composite walls or not prism shape - for ( TopExp_Explorer f( shape3D, TopAbs_FACE ); f.More(); f.Next() ) + bool savedSetErrorToSM = mySetErrorToSM; + mySetErrorToSM = false; // ingore errors in initPrism() + + // search among meshed FACEs + list< SMESH_subMesh* >::iterator sm = meshedSubMesh.begin(); + for ( ; !botSM && sm != meshedSubMesh.end(); ++sm ) + { + thePrism.Clear(); + botSM = *sm; + thePrism.myBottom = TopoDS::Face( botSM->GetSubShape() ); + if ( !initPrism( thePrism, theShape3D, /*selectBottom=*/false )) + botSM = NULL; + } + // search among all FACEs + for ( TopExp_Explorer f( theShape3D, TopAbs_FACE ); !botSM && f.More(); f.Next() ) { int minNbFaces = 2 + myHelper->Count( f.Current(), TopAbs_EDGE, false); - if ( nbFaces >= minNbFaces) - { - thePrism.Clear(); - thePrism.myBottom = TopoDS::Face( f.Current() ); - if ( initPrism( thePrism, shape3D )) - return true; - } - return toSM( error( COMPERR_BAD_SHAPE )); + if ( nbFaces < minNbFaces) continue; + thePrism.Clear(); + thePrism.myBottom = TopoDS::Face( f.Current() ); + botSM = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom ); + if ( !initPrism( thePrism, theShape3D, /*selectBottom=*/false )) + botSM = NULL; } + mySetErrorToSM = savedSetErrorToSM; + return botSM ? true : toSM( error( COMPERR_BAD_SHAPE )); } // find vertex 000 - the one with smallest coordinates (for easy DEBUG :-) @@ -2929,11 +3131,11 @@ bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism, } } - thePrism.myShape3D = shape3D; + thePrism.myShape3D = theShape3D; if ( thePrism.myBottom.IsNull() ) thePrism.myBottom = TopoDS::Face( botSM->GetSubShape() ); - thePrism.myBottom.Orientation( myHelper->GetSubShapeOri( shape3D, thePrism.myBottom )); - thePrism.myTop. Orientation( myHelper->GetSubShapeOri( shape3D, thePrism.myTop )); + thePrism.myBottom.Orientation( myHelper->GetSubShapeOri( theShape3D, thePrism.myBottom )); + thePrism.myTop. Orientation( myHelper->GetSubShapeOri( theShape3D, thePrism.myTop )); // Get ordered bottom edges TopoDS_Face reverseBottom = // to have order of top EDGEs as in the top FACE @@ -2943,7 +3145,7 @@ bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism, thePrism.myNbEdgesInWires, V000 ); // Get Wall faces corresponding to the ordered bottom edges and the top FACE - if ( !getWallFaces( thePrism, nbFaces )) + if ( !getWallFaces( thePrism, nbFaces )) // it also sets thePrism.myTop return false; //toSM( error(COMPERR_BAD_SHAPE, "Can't find side faces")); if ( topSM ) @@ -3124,7 +3326,7 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper, if ( nbUnitePerEdge[ iE ] < 0 ) continue; // look for already united faces - for ( int i = iE; i < iE + nbExraFaces; ++i ) + for ( size_t i = iE; i < iE + nbExraFaces; ++i ) { if ( nbUnitePerEdge[ i ] > 0 ) // a side including nbUnitePerEdge[i]+1 edge nbExraFaces += nbUnitePerEdge[ i ]; @@ -3163,7 +3365,7 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper, else if ( nbExraFaces > 1 ) // unite { double u0 = 0, sumLen = 0; - for ( int i = iE; i < iE + nbExraFaces; ++i ) + for ( size_t i = iE; i < iE + nbExraFaces; ++i ) sumLen += edgeLength[ i ]; vector< TSideFace* > components( nbExraFaces ); @@ -3314,10 +3516,8 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper, } } -// #define SHOWYXZ(msg, xyz) { \ -// gp_Pnt p (xyz); \ -// cout << msg << " ("<< p.X() << "; " < & double tol2; { Bnd_B3d bndBox; - for ( int i = 0; i < columns.size(); ++i ) + for ( size_t i = 0; i < columns.size(); ++i ) bndBox.Add( gpXYZ( columns[i]->front() )); tol2 = bndBox.SquareExtent() * 1e-5; } @@ -3432,7 +3632,7 @@ bool StdMeshers_PrismAsBlock::GetLayersTransformation(vector & //t.SetScaleFactor( distZ/dist0 ); - it does not work properly, wrong base point // check a transformation - for ( int i = 0; i < columns.size(); ++i ) + for ( size_t i = 0; i < columns.size(); ++i ) { gp_Pnt p0 = gpXYZ( (*columns[i])[0] ); gp_Pnt pz = gpXYZ( (*columns[i])[z] ); @@ -3636,7 +3836,7 @@ StdMeshers_PrismAsBlock::TSideFace::TSideFace( const TSideFace& other ): myComponents ( other.myComponents.size() ), myHelper ( *other.myHelper.GetMesh() ) { - for (int i = 0 ; i < myComponents.size(); ++i ) + for ( size_t i = 0 ; i < myComponents.size(); ++i ) myComponents[ i ] = new TSideFace( *other.myComponents[ i ]); } @@ -3648,7 +3848,7 @@ StdMeshers_PrismAsBlock::TSideFace::TSideFace( const TSideFace& other ): StdMeshers_PrismAsBlock::TSideFace::~TSideFace() { - for (int i = 0 ; i < myComponents.size(); ++i ) + for ( size_t i = 0 ; i < myComponents.size(); ++i ) if ( myComponents[ i ] ) delete myComponents[ i ]; } @@ -3748,7 +3948,7 @@ StdMeshers_PrismAsBlock::TSideFace::GetComponent(const double U,double & localU) if ( myComponents.empty() ) return const_cast( this ); - int i; + size_t i; for ( i = 0; i < myComponents.size(); ++i ) if ( U < myParams[ i ].second ) break; @@ -3887,7 +4087,7 @@ gp_Pnt StdMeshers_PrismAsBlock::TSideFace::Value(const Standard_Real U, TopoDS_Shape s = myHelper.GetSubShapeByNode( nn[0], myHelper.GetMeshDS() ); if ( s.ShapeType() != TopAbs_EDGE ) s = myHelper.GetSubShapeByNode( nn[2], myHelper.GetMeshDS() ); - if ( s.ShapeType() == TopAbs_EDGE ) + if ( !s.IsNull() && s.ShapeType() == TopAbs_EDGE ) edge = TopoDS::Edge( s ); } if ( !edge.IsNull() ) @@ -4187,9 +4387,9 @@ gp_Pnt StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::Value(const Standard_Real void StdMeshers_PrismAsBlock::TVerticalEdgeAdaptor::dumpNodes(int nbNodes) const { #ifdef _DEBUG_ - for ( int i = 0; i < nbNodes && i < myNodeColumn->size(); ++i ) + for ( int i = 0; i < nbNodes && i < (int)myNodeColumn->size(); ++i ) cout << (*myNodeColumn)[i]->GetID() << " "; - if ( nbNodes < myNodeColumn->size() ) + if ( nbNodes < (int) myNodeColumn->size() ) cout << myNodeColumn->back()->GetID(); #endif } @@ -4430,7 +4630,8 @@ void StdMeshers_Sweeper::applyBoundaryError(const vector< gp_XYZ >& bndPoints, //================================================================================ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, - const double tol) + const double tol, + const bool allowHighBndError) { const size_t zSize = myBndColumns[0]->size(); const size_t zSrc = 0, zTgt = zSize-1; @@ -4549,6 +4750,9 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, 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.);