X-Git-Url: http://git.salome-platform.org/gitweb/?a=blobdiff_plain;f=src%2FStdMeshers%2FStdMeshers_Prism_3D.cxx;h=6ba56c39c926bf4fd20095540b8f854342e29363;hb=4031dc980843986775a48d3f9eae8642be249887;hp=00a53610546f39a9b6cf983e50dba7ed1accbc78;hpb=9a54694a0ab1e5cbc558a35c4606ceea4f7af2ef;p=modules%2Fsmesh.git diff --git a/src/StdMeshers/StdMeshers_Prism_3D.cxx b/src/StdMeshers/StdMeshers_Prism_3D.cxx index 00a536105..6ba56c39c 100644 --- a/src/StdMeshers/StdMeshers_Prism_3D.cxx +++ b/src/StdMeshers/StdMeshers_Prism_3D.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2012 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2013 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 @@ -68,9 +68,15 @@ using namespace std; // gp_Pnt p (xyz); \ // cout << msg << " ("<< p.X() << "; " < & nbUnitePerEdge ) + { + int nbEdges = thePrism.myNbEdgesInWires.front(); // nb outer edges + int nbSides = nbEdges; + + list< TopoDS_Edge >::const_iterator edgeIt = thePrism.myBottomEdges.begin(); + std::advance( edgeIt, nbEdges-1 ); + TopoDS_Edge prevE = *edgeIt; + bool isPrevStraight = SMESH_Algo::isStraight( prevE ); + int iPrev = nbEdges - 1; + + int iUnite = -1; // the first of united EDGEs + + edgeIt = thePrism.myBottomEdges.begin(); + for ( int iE = 0; iE < nbEdges; ++iE, ++edgeIt ) + { + const TopoDS_Edge& curE = *edgeIt; + const bool isCurStraight = SMESH_Algo::isStraight( curE ); + if ( isPrevStraight && isCurStraight && SMESH_Algo::IsContinuous( prevE, curE )) + { + if ( iUnite < 0 ) + iUnite = iPrev; + nbUnitePerEdge[ iUnite ]++; + nbUnitePerEdge[ iE ] = -1; + --nbSides; + } + else + { + iUnite = -1; + } + prevE = curE; + isPrevStraight = isCurStraight; + iPrev = iE; + } + + return nbSides; + } + + void pointsToPython(const std::vector& p) + { +#ifdef _DEBUG_ + for ( int 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; + } +#endif + } } // namespace //======================================================================= @@ -457,21 +518,14 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh if ( nbSolids < 1 ) return true; - Prism_3D::TPrismTopo prism; - - if ( nbSolids == 1 ) - { - return ( initPrism( prism, TopExp_Explorer( theShape, TopAbs_SOLID ).Current() ) && - compute( prism )); - } - TopTools_IndexedDataMapOfShapeListOfShape faceToSolids; TopExp::MapShapesAndAncestors( theShape, TopAbs_FACE, TopAbs_SOLID, faceToSolids ); // look for meshed FACEs ("source" FACEs) that must be prism bottoms - list< TopoDS_Face > meshedFaces;//, notQuadMeshedFaces, notQuadFaces; + list< TopoDS_Face > meshedFaces, notQuadMeshedFaces, notQuadFaces; const bool meshHasQuads = ( theMesh.NbQuadrangles() > 0 ); - for ( int iF = 1; iF < faceToSolids.Extent(); ++iF ) + //StdMeshers_Quadrangle_2D* quadAlgo = TQuadrangleAlgo::instance( this ); + for ( int iF = 1; iF <= faceToSolids.Extent(); ++iF ) { const TopoDS_Face& face = TopoDS::Face( faceToSolids.FindKey( iF )); SMESH_subMesh* faceSM = theMesh.GetSubMesh( face ); @@ -479,17 +533,36 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh { if ( !meshHasQuads || !helper.IsSameElemGeometry( faceSM->GetSubMeshDS(), SMDSGeom_QUADRANGLE ) || - !helper.IsStructured( faceSM )) - // notQuadMeshedFaces are of higher priority + !helper.IsStructured( faceSM ) + ) + notQuadMeshedFaces.push_front( face ); + else if ( myHelper->Count( face, TopAbs_EDGE, /*ignoreSame=*/false ) != 4 ) meshedFaces.push_front( face ); else meshedFaces.push_back( face ); } + // not add not quadrilateral FACE as we can't compute it + // else if ( !quadAlgo->CheckNbEdges( theMesh, face )) + // // not add not quadrilateral FACE as it can be a prism side + // // else if ( myHelper->Count( face, TopAbs_EDGE, /*ignoreSame=*/false ) != 4 ) + // { + // notQuadFaces.push_back( face ); + // } } - //meshedFaces.splice( meshedFaces.begin(), notQuadMeshedFaces ); + // notQuadFaces are of medium priority, put them before ordinary meshed faces + meshedFaces.splice( meshedFaces.begin(), notQuadFaces ); + // notQuadMeshedFaces are of highest priority, put them before notQuadFaces + meshedFaces.splice( meshedFaces.begin(), notQuadMeshedFaces ); - // if ( meshedFaces.empty() ) - // return error( COMPERR_BAD_INPUT_MESH, "No meshed source faces found" ); + Prism_3D::TPrismTopo prism; + + if ( nbSolids == 1 ) + { + if ( !meshedFaces.empty() ) + prism.myBottom = meshedFaces.front(); + return ( initPrism( prism, TopExp_Explorer( theShape, TopAbs_SOLID ).Current() ) && + compute( prism )); + } TopTools_MapOfShape meshedSolids; list< Prism_3D::TPrismTopo > meshedPrism; @@ -619,6 +692,27 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh // TODO. there are other ways to find out the source FACE: // propagation, topological similarity, ect. + // simply try to mesh all not meshed SOLIDs + if ( meshedFaces.empty() ) + { + for ( TopExp_Explorer solid( theShape, TopAbs_SOLID ); solid.More(); solid.Next() ) + { + mySetErrorToSM = false; + prism.Clear(); + if ( !meshedSolids.Contains( solid.Current() ) && + initPrism( prism, solid.Current() )) + { + mySetErrorToSM = true; + if ( !compute( prism )) + return false; + meshedFaces.push_front( prism.myTop ); + meshedFaces.push_front( prism.myBottom ); + meshedPrism.push_back( prism ); + meshedSolids.Add( solid.Current() ); + } + mySetErrorToSM = true; + } + } if ( meshedFaces.empty() ) // set same error to 10 not-computed solids { @@ -633,7 +727,7 @@ bool StdMeshers_Prism_3D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theSh SMESH_subMesh* sm = theMesh.GetSubMesh( solid.Current() ); sm->GetComputeError() = err; } - return false; + return error( err ); } } return true; @@ -666,10 +760,11 @@ 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(); int iE = 0; + double f,l; while ( edge != thePrism.myBottomEdges.end() ) { ++iE; - if ( BRep_Tool::Degenerated( *edge )) + if ( BRep_Tool::Curve( *edge, f,l ).IsNull() ) { edge = thePrism.myBottomEdges.erase( edge ); --iE; @@ -708,7 +803,7 @@ bool StdMeshers_Prism_3D::getWallFaces( Prism_3D::TPrismTopo & thePrism, // ------------------------- // Compose a vector of indixes of right neighbour FACE for each wall FACE - // that is not so evident in case of several WIREs + // 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 ); @@ -950,6 +1045,9 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) if ( !SMESH_Block::ShellPoint( params, myShapeXYZ, coords )) return toSM( error("Can't compute coordinates by normalized parameters")); + // if ( !meshDS->MeshElements( volumeID ) || + // meshDS->MeshElements( volumeID )->NbNodes() == 0 ) + // pointsToPython(myShapeXYZ); SHOWYXZ("TOPFacePoint ",myShapeXYZ[ ID_TOP_FACE]); SHOWYXZ("BOT Node "<< tBotNode.myNode->GetID(),gpXYZ(tBotNode.myNode)); SHOWYXZ("ShellPoint ",coords); @@ -1013,6 +1111,7 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism) { SMESH_Mesh* mesh = myHelper->GetMesh(); SMESHDS_Mesh* meshDS = myHelper->GetMeshDS(); + DBGOUT( endl << "COMPUTE Prism " << meshDS->ShapeToIndex( thePrism.myShape3D )); TProjction1dAlgo* projector1D = TProjction1dAlgo::instance( this ); StdMeshers_Quadrangle_2D* quadAlgo = TQuadrangleAlgo::instance( this, myHelper ); @@ -1078,6 +1177,7 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism) const TopoDS_Edge& srcE = lftSide->Edge(i); SMESH_subMesh* srcSM = mesh->GetSubMesh( srcE ); if ( !srcSM->IsMeshComputed() ) { + DBGOUT( "COMPUTE V edge " << srcSM->GetId() ); srcSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE ); srcSM->ComputeStateEngine ( SMESH_subMesh::COMPUTE ); if ( !srcSM->IsMeshComputed() ) @@ -1112,7 +1212,7 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism) } continue; } - // Compute + // Compute 'vertical projection' if ( nbTgtMeshed == 0 ) { // compute nodes on target VERTEXes @@ -1131,8 +1231,9 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism) } // compute nodes on target EDGEs + DBGOUT( "COMPUTE V edge (proj) " << shapeID( lftSide->Edge(0))); rgtSide->Reverse(); // direct it same as the lftSide - myHelper->SetElementsOnShape( false ); + myHelper->SetElementsOnShape( false ); // myHelper holds the prism shape TopoDS_Edge tgtEdge; for ( size_t iN = 1; iN < srcNodeStr.size()-1; ++iN ) // add nodes { @@ -1143,25 +1244,24 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism) } for ( size_t iN = 1; iN < srcNodeStr.size(); ++iN ) // add segments { - SMDS_MeshElement* newEdge = myHelper->AddEdge( newNodes[ iN-1 ], newNodes[ iN ] ); + // find an EDGE to set a new segment std::pair id2type = myHelper->GetMediumPos( newNodes[ iN-1 ], newNodes[ iN ] ); - if ( id2type.second == TopAbs_EDGE ) - { - meshDS->SetMeshElementOnShape( newEdge, id2type.first ); - } - else // new nodes are on different EDGEs; put one of them on VERTEX + if ( id2type.second != TopAbs_EDGE ) { + // new nodes are on different EDGEs; put one of them on VERTEX const int edgeIndex = rgtSide->EdgeIndex( srcNodeStr[ iN-1 ].normParam ); const double vertexParam = rgtSide->LastParameter( edgeIndex ); const gp_Pnt p = BRep_Tool::Pnt( rgtSide->LastVertex( edgeIndex )); const int isPrev = ( Abs( srcNodeStr[ iN-1 ].normParam - vertexParam ) < Abs( srcNodeStr[ iN ].normParam - vertexParam )); - meshDS->SetMeshElementOnShape( newEdge, newNodes[ iN-(1-isPrev) ]->getshapeId() ); meshDS->UnSetNodeOnShape( newNodes[ iN-isPrev ] ); meshDS->SetNodeOnVertex ( newNodes[ iN-isPrev ], rgtSide->LastVertex( edgeIndex )); - meshDS->MoveNode( newNodes[ iN-isPrev ], p.X(), p.Y(), p.Z() ); + meshDS->MoveNode ( newNodes[ iN-isPrev ], p.X(), p.Y(), p.Z() ); + id2type.first = newNodes[ iN-(1-isPrev) ]->getshapeId(); } + SMDS_MeshElement* newEdge = myHelper->AddEdge( newNodes[ iN-1 ], newNodes[ iN ] ); + meshDS->SetMeshElementOnShape( newEdge, id2type.first ); } myHelper->SetElementsOnShape( true ); for ( int i = 0; i < rgtSide->NbEdges(); ++i ) // update state of sub-meshes @@ -1198,28 +1298,43 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism) // Top EDGEs must be projections from the bottom ones // to compute stuctured quad mesh on wall FACEs // --------------------------------------------------- - const TopoDS_Edge& botE = (*quad)->side[ QUAD_BOTTOM_SIDE ]->Edge(0); - const TopoDS_Edge& topE = (*quad)->side[ QUAD_TOP_SIDE ]->Edge(0); - - projector1D->myHyp.SetSourceEdge( botE ); - - SMESH_subMesh* tgtEdgeSm = mesh->GetSubMesh( topE ); - if ( !tgtEdgeSm->IsMeshComputed() ) { - // compute nodes on VERTEXes - tgtEdgeSm->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE ); - // project segments - projector1D->InitComputeError(); - bool ok = projector1D->Compute( *mesh, topE ); - if ( !ok ) + const TopoDS_Edge& botE = (*quad)->side[ QUAD_BOTTOM_SIDE ]->Edge(0); + const TopoDS_Edge& topE = (*quad)->side[ QUAD_TOP_SIDE ]->Edge(0); + SMESH_subMesh* botSM = mesh->GetSubMesh( botE ); + SMESH_subMesh* topSM = mesh->GetSubMesh( topE ); + SMESH_subMesh* srcSM = botSM; + SMESH_subMesh* tgtSM = topSM; + if ( !srcSM->IsMeshComputed() && topSM->IsMeshComputed() ) + std::swap( srcSM, tgtSM ); + + if ( !srcSM->IsMeshComputed() ) { - SMESH_ComputeErrorPtr err = projector1D->GetComputeError(); - if ( err->IsOK() ) err->myName = COMPERR_ALGO_FAILED; - tgtEdgeSm->GetComputeError() = err; - return false; + DBGOUT( "COMPUTE H edge " << srcSM->GetId()); + srcSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE ); // nodes on VERTEXes + srcSM->ComputeStateEngine( SMESH_subMesh::COMPUTE ); // segments on the EDGE } + srcSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); + + if ( !tgtSM->IsMeshComputed() ) + { + // compute nodes on VERTEXes + tgtSM->ComputeSubMeshStateEngine( SMESH_subMesh::COMPUTE ); + // project segments + DBGOUT( "COMPUTE H edge (proj) " << tgtSM->GetId()); + projector1D->myHyp.SetSourceEdge( TopoDS::Edge( srcSM->GetSubShape() )); + projector1D->InitComputeError(); + bool ok = projector1D->Compute( *mesh, tgtSM->GetSubShape() ); + if ( !ok ) + { + SMESH_ComputeErrorPtr err = projector1D->GetComputeError(); + if ( err->IsOK() ) err->myName = COMPERR_ALGO_FAILED; + tgtSM->GetComputeError() = err; + return false; + } + } + tgtSM->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); } - tgtEdgeSm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); // Compute quad mesh on wall FACEs // ------------------------------- @@ -1234,6 +1349,7 @@ bool StdMeshers_Prism_3D::computeWalls(const Prism_3D::TPrismTopo& thePrism) "Not all edges have valid algorithm and hypothesis")); // mesh the quadAlgo->InitComputeError(); + DBGOUT( "COMPUTE Quad face " << fSM->GetId()); bool ok = quadAlgo->Compute( *mesh, face ); fSM->GetComputeError() = quadAlgo->GetComputeError(); if ( !ok ) @@ -1511,7 +1627,12 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top() SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS(); if ( !botSMDS || botSMDS->NbElements() == 0 ) - return toSM( error(TCom("No elememts on face #") << botSM->GetId() )); + { + _gen->Compute( *myHelper->GetMesh(), botSM->GetSubShape() ); + botSMDS = botSM->GetSubMeshDS(); + if ( !botSMDS || botSMDS->NbElements() == 0 ) + return toSM( error(TCom("No elements on face #") << botSM->GetId() )); + } bool needProject = !topSM->IsMeshComputed(); if ( !needProject && @@ -1646,7 +1767,7 @@ bool StdMeshers_Prism_3D::projectBottomToTop() // if the bottom faces is orienetd OK then top faces must be reversed bool reverseTop = true; if ( myHelper->NbAncestors( botFace, *myBlock.Mesh(), TopAbs_SOLID ) > 1 ) - reverseTop = ! SMESH_Algo::IsReversedSubMesh( TopoDS::Face( botFace ), meshDS ); + reverseTop = ! myHelper->IsReversedSubMesh( TopoDS::Face( botFace )); int iFrw, iRev, *iPtr = &( reverseTop ? iRev : iFrw ); // loop on bottom mesh faces @@ -2056,15 +2177,17 @@ bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism, bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper, const Prism_3D::TPrismTopo& thePrism) { + myHelper = helper; + SMESHDS_Mesh* meshDS = myHelper->GetMeshDS(); + SMESH_Mesh* mesh = myHelper->GetMesh(); + if ( mySide ) { delete mySide; mySide = 0; } vector< TSideFace* > sideFaces( NB_WALL_FACES, 0 ); vector< pair< double, double> > params( NB_WALL_FACES ); - mySide = new TSideFace( sideFaces, params ); + mySide = new TSideFace( *mesh, sideFaces, params ); - myHelper = helper; - SMESHDS_Mesh* meshDS = myHelper->GetMeshDS(); SMESH_Block::init(); myShapeIDMap.Clear(); @@ -2086,9 +2209,16 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper, myParam2ColumnMaps.resize( thePrism.myBottomEdges.size() ); // total nb edges size_t iE, nbEdges = thePrism.myNbEdgesInWires.front(); // nb outer edges - vector< double > edgeLength( nbEdges ); + vector< double > edgeLength( nbEdges ); multimap< double, int > len2edgeMap; + // for each EDGE: either split into several parts, or join with several next EDGEs + vector nbSplitPerEdge( nbEdges, 0 ); + vector nbUnitePerEdge( nbEdges, 0 ); // -1 means "joined to a previous" + + // consider continuous straight EDGEs as one side + const int nbSides = countNbSides( thePrism, nbUnitePerEdge ); + list< TopoDS_Edge >::const_iterator edgeIt = thePrism.myBottomEdges.begin(); for ( iE = 0; iE < nbEdges; ++iE, ++edgeIt ) { @@ -2108,14 +2238,8 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper, edgeLength[ iE ] = SMESH_Algo::EdgeLength( *edgeIt ); - if ( nbEdges < NB_WALL_FACES ) // fill map used to split faces - { - SMESHDS_SubMesh* smDS = meshDS->MeshElements( *edgeIt); - if ( !smDS ) - return error(COMPERR_BAD_INPUT_MESH, TCom("Null submesh on the edge #") - << MeshDS()->ShapeToIndex( *edgeIt )); - len2edgeMap.insert( make_pair( edgeLength[ iE ], iE )); - } + if ( nbSides < NB_WALL_FACES ) // fill map used to split faces + len2edgeMap.insert( make_pair( edgeLength[ iE ], iE )); // sort edges by length } // Load columns of internal edges (forming holes) // and fill map ShapeIndex to TParam2ColumnMap for them @@ -2153,95 +2277,122 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper, // Create 4 wall faces of a block // ------------------------------- - if ( nbEdges <= NB_WALL_FACES ) // ************* Split faces if necessary + if ( nbSides <= NB_WALL_FACES ) // ************* Split faces if necessary { - map< int, int > iE2nbSplit; - if ( nbEdges != NB_WALL_FACES ) // define how to split + if ( nbSides != NB_WALL_FACES ) // define how to split { if ( len2edgeMap.size() != nbEdges ) RETURN_BAD_RESULT("Uniqueness of edge lengths not assured"); - map< double, int >::reverse_iterator maxLen_i = len2edgeMap.rbegin(); - map< double, int >::reverse_iterator midLen_i = ++len2edgeMap.rbegin(); + + multimap< double, int >::reverse_iterator maxLen_i = len2edgeMap.rbegin(); + multimap< double, int >::reverse_iterator midLen_i = ++len2edgeMap.rbegin(); + double maxLen = maxLen_i->first; double midLen = ( len2edgeMap.size() == 1 ) ? 0 : midLen_i->first; switch ( nbEdges ) { case 1: // 0-th edge is split into 4 parts - iE2nbSplit.insert( make_pair( 0, 4 )); break; + nbSplitPerEdge[ 0 ] = 4; + break; case 2: // either the longest edge is split into 3 parts, or both edges into halves if ( maxLen / 3 > midLen / 2 ) { - iE2nbSplit.insert( make_pair( maxLen_i->second, 3 )); + nbSplitPerEdge[ maxLen_i->second ] = 3; } else { - iE2nbSplit.insert( make_pair( maxLen_i->second, 2 )); - iE2nbSplit.insert( make_pair( midLen_i->second, 2 )); + nbSplitPerEdge[ maxLen_i->second ] = 2; + nbSplitPerEdge[ midLen_i->second ] = 2; } break; case 3: - // split longest into halves - iE2nbSplit.insert( make_pair( maxLen_i->second, 2 )); + if ( nbSides == 2 ) + // split longest into 3 parts + nbSplitPerEdge[ maxLen_i->second ] = 3; + else + // split longest into halves + nbSplitPerEdge[ maxLen_i->second ] = 2; } } - // Create TSideFace's - int iSide = 0; - list< TopoDS_Edge >::const_iterator botE = thePrism.myBottomEdges.begin(); - for ( iE = 0; iE < nbEdges; ++iE, ++botE ) + } + else // **************************** Unite faces + { + int nbExraFaces = nbSides - 3; // nb of faces to fuse + for ( iE = 0; iE < nbEdges; ++iE ) { - TFaceQuadStructPtr quad = thePrism.myWallQuads[ iE ].front(); - // split? - map< int, int >::iterator i_nb = iE2nbSplit.find( iE ); - if ( i_nb != iE2nbSplit.end() ) { - // split! - int nbSplit = i_nb->second; - vector< double > params; - splitParams( nbSplit, &myParam2ColumnMaps[ iE ], params ); - const bool isForward = - StdMeshers_PrismAsBlock::IsForwardEdge( myHelper->GetMeshDS(), - myParam2ColumnMaps[iE], - *botE, SMESH_Block::ID_Fx0z ); - for ( int i = 0; i < nbSplit; ++i ) { - double f = ( isForward ? params[ i ] : params[ nbSplit - i-1 ]); - double l = ( isForward ? params[ i+1 ] : params[ nbSplit - i ]); - TSideFace* comp = new TSideFace( myHelper, wallFaceIds[ iSide ], - thePrism.myWallQuads[ iE ], *botE, - &myParam2ColumnMaps[ iE ], f, l ); - mySide->SetComponent( iSide++, comp ); - } + if ( nbUnitePerEdge[ iE ] < 0 ) + continue; + // look for already united faces + for ( int i = iE; i < iE + nbExraFaces; ++i ) + { + if ( nbUnitePerEdge[ i ] > 0 ) // a side including nbUnitePerEdge[i]+1 edge + nbExraFaces += nbUnitePerEdge[ i ]; + nbUnitePerEdge[ i ] = -1; } - else { - TSideFace* comp = new TSideFace( myHelper, wallFaceIds[ iSide ], + nbUnitePerEdge[ iE ] = nbExraFaces; + break; + } + } + + // Create TSideFace's + int iSide = 0; + list< TopoDS_Edge >::const_iterator botE = thePrism.myBottomEdges.begin(); + for ( iE = 0; iE < nbEdges; ++iE, ++botE ) + { + TFaceQuadStructPtr quad = thePrism.myWallQuads[ iE ].front(); + const int nbSplit = nbSplitPerEdge[ iE ]; + const int nbExraFaces = nbUnitePerEdge[ iE ] + 1; + if ( nbSplit > 0 ) // split + { + vector< double > params; + splitParams( nbSplit, &myParam2ColumnMaps[ iE ], params ); + const bool isForward = + StdMeshers_PrismAsBlock::IsForwardEdge( myHelper->GetMeshDS(), + myParam2ColumnMaps[iE], + *botE, SMESH_Block::ID_Fx0z ); + for ( int i = 0; i < nbSplit; ++i ) { + double f = ( isForward ? params[ i ] : params[ nbSplit - i-1 ]); + double l = ( isForward ? params[ i+1 ] : params[ nbSplit - i ]); + TSideFace* comp = new TSideFace( *mesh, wallFaceIds[ iSide ], thePrism.myWallQuads[ iE ], *botE, - &myParam2ColumnMaps[ iE ]); + &myParam2ColumnMaps[ iE ], f, l ); mySide->SetComponent( iSide++, comp ); } } - } - else { // **************************** Unite faces - - // unite first faces - int nbExraFaces = nbEdges - 3; - int iSide = 0, iE; - double u0 = 0, sumLen = 0; - for ( iE = 0; iE < nbExraFaces; ++iE ) - sumLen += edgeLength[ iE ]; - - vector< TSideFace* > components( nbExraFaces ); - vector< pair< double, double> > params( nbExraFaces ); - list< TopoDS_Edge >::const_iterator botE = thePrism.myBottomEdges.begin(); - for ( iE = 0; iE < nbExraFaces; ++iE, ++botE ) + else if ( nbExraFaces > 1 ) // unite { - components[ iE ] = new TSideFace( myHelper, wallFaceIds[ iSide ], - thePrism.myWallQuads[ iE ], *botE, - &myParam2ColumnMaps[ iE ]); - double u1 = u0 + edgeLength[ iE ] / sumLen; - params[ iE ] = make_pair( u0 , u1 ); - u0 = u1; + double u0 = 0, sumLen = 0; + for ( int i = iE; i < iE + nbExraFaces; ++i ) + sumLen += edgeLength[ i ]; + + vector< TSideFace* > components( nbExraFaces ); + vector< pair< double, double> > params( nbExraFaces ); + bool endReached = false; + for ( int i = 0; i < nbExraFaces; ++i, ++botE, ++iE ) + { + if ( iE == nbEdges ) + { + endReached = true; + botE = thePrism.myBottomEdges.begin(); + iE = 0; + } + components[ i ] = new TSideFace( *mesh, wallFaceIds[ iSide ], + thePrism.myWallQuads[ iE ], *botE, + &myParam2ColumnMaps[ iE ]); + double u1 = u0 + edgeLength[ iE ] / sumLen; + params[ i ] = make_pair( u0 , u1 ); + u0 = u1; + } + TSideFace* comp = new TSideFace( *mesh, components, params ); + mySide->SetComponent( iSide++, comp ); + if ( endReached ) + break; + --iE; // for increment in an external loop on iE + --botE; } - mySide->SetComponent( iSide++, new TSideFace( components, params )); - - // fill the rest faces - for ( ; iE < nbEdges; ++iE, ++botE ) + else if ( nbExraFaces < 0 ) // skip already united face + { + } + else // use as is { - TSideFace* comp = new TSideFace( myHelper, wallFaceIds[ iSide ], + TSideFace* comp = new TSideFace( *mesh, wallFaceIds[ iSide ], thePrism.myWallQuads[ iE ], *botE, &myParam2ColumnMaps[ iE ]); mySide->SetComponent( iSide++, comp ); @@ -2357,15 +2508,20 @@ bool StdMeshers_PrismAsBlock::Init(SMESH_MesherHelper* helper, } } -// gp_XYZ testPar(0.25, 0.25, 0), testCoord; -// if ( !FacePoint( ID_BOT_FACE, testPar, testCoord )) -// RETURN_BAD_RESULT("TEST FacePoint() FAILED"); -// SHOWYXZ("IN TEST PARAM" , testPar); -// SHOWYXZ("OUT TEST CORD" , testCoord); -// if ( !ComputeParameters( testCoord, testPar , ID_BOT_FACE)) -// RETURN_BAD_RESULT("TEST ComputeParameters() FAILED"); -// SHOWYXZ("OUT TEST PARAM" , testPar); - + // double _u[]={ 0.1, 0.1, 0.9, 0.9 }; + // double _v[]={ 0.1, 0.9, 0.1, 0.9, }; + // for ( int i = 0; i < 4; ++i ) + // { + // //gp_XYZ testPar(0.25, 0.25, 0), testCoord; + // gp_XYZ testPar(_u[i], _v[i], 0), testCoord; + // if ( !FacePoint( ID_BOT_FACE, testPar, testCoord )) + // RETURN_BAD_RESULT("TEST FacePoint() FAILED"); + // SHOWYXZ("IN TEST PARAM" , testPar); + // SHOWYXZ("OUT TEST CORD" , testCoord); + // if ( !ComputeParameters( testCoord, testPar , ID_BOT_FACE)) + // RETURN_BAD_RESULT("TEST ComputeParameters() FAILED"); + // SHOWYXZ("OUT TEST PARAM" , testPar); + // } return true; } @@ -2519,7 +2675,7 @@ bool StdMeshers_PrismAsBlock::IsForwardEdge(SMESHDS_Mesh* meshDS, */ //================================================================================ -StdMeshers_PrismAsBlock::TSideFace::TSideFace(SMESH_MesherHelper* helper, +StdMeshers_PrismAsBlock::TSideFace::TSideFace(SMESH_Mesh& mesh, const int faceID, const Prism_3D::TQuadList& quadList, const TopoDS_Edge& baseEdge, @@ -2528,20 +2684,22 @@ StdMeshers_PrismAsBlock::TSideFace::TSideFace(SMESH_MesherHelper* helper, const double last): myID( faceID ), myParamToColumnMap( columnsMap ), - myHelper( helper ) + myHelper( mesh ) { myParams.resize( 1 ); myParams[ 0 ] = make_pair( first, last ); mySurface = PSurface( new BRepAdaptor_Surface( quadList.front()->face )); myBaseEdge = baseEdge; - myIsForward = StdMeshers_PrismAsBlock::IsForwardEdge( myHelper->GetMeshDS(), + myIsForward = StdMeshers_PrismAsBlock::IsForwardEdge( myHelper.GetMeshDS(), *myParamToColumnMap, myBaseEdge, myID ); + myHelper.SetSubShape( quadList.front()->face ); + if ( quadList.size() > 1 ) // side is vertically composite { // fill myShapeID2Surf map to enable finding a right surface by any sub-shape ID - SMESHDS_Mesh* meshDS = myHelper->GetMeshDS(); + SMESHDS_Mesh* meshDS = myHelper.GetMeshDS(); TopTools_IndexedDataMapOfShapeListOfShape subToFaces; Prism_3D::TQuadList::const_iterator quad = quadList.begin(); @@ -2566,20 +2724,34 @@ StdMeshers_PrismAsBlock::TSideFace::TSideFace(SMESH_MesherHelper* helper, //================================================================================ /*! - * \brief Constructor of complex side face + * \brief Constructor of a complex side face */ //================================================================================ StdMeshers_PrismAsBlock::TSideFace:: -TSideFace(const vector< TSideFace* >& components, +TSideFace(SMESH_Mesh& mesh, + const vector< TSideFace* >& components, const vector< pair< double, double> > & params) :myID( components[0] ? components[0]->myID : 0 ), myParamToColumnMap( 0 ), myParams( params ), myIsForward( true ), myComponents( components ), - myHelper( components[0] ? components[0]->myHelper : 0 ) -{} + myHelper( mesh ) +{ + if ( myID == ID_Fx1z || myID == ID_F0yz ) + { + // reverse components + std::reverse( myComponents.begin(), myComponents.end() ); + std::reverse( myParams.begin(), myParams.end() ); + for ( size_t i = 0; i < myParams.size(); ++i ) + { + const double f = myParams[i].first; + const double l = myParams[i].second; + myParams[i] = make_pair( 1. - l, 1. - f ); + } + } +} //================================================================================ /*! * \brief Copy constructor @@ -2587,17 +2759,17 @@ TSideFace(const vector< TSideFace* >& components, */ //================================================================================ -StdMeshers_PrismAsBlock::TSideFace::TSideFace( const TSideFace& other ) +StdMeshers_PrismAsBlock::TSideFace::TSideFace( const TSideFace& other ): + myID ( other.myID ), + myParamToColumnMap ( other.myParamToColumnMap ), + mySurface ( other.mySurface ), + myBaseEdge ( other.myBaseEdge ), + myShapeID2Surf ( other.myShapeID2Surf ), + myParams ( other.myParams ), + myIsForward ( other.myIsForward ), + myComponents ( other.myComponents.size() ), + myHelper ( *other.myHelper.GetMesh() ) { - myID = other.myID; - mySurface = other.mySurface; - myBaseEdge = other.myBaseEdge; - myParams = other.myParams; - myIsForward = other.myIsForward; - myHelper = other.myHelper; - myParamToColumnMap = other.myParamToColumnMap; - - myComponents.resize( other.myComponents.size()); for (int i = 0 ; i < myComponents.size(); ++i ) myComponents[ i ] = new TSideFace( *other.myComponents[ i ]); } @@ -2800,16 +2972,16 @@ gp_Pnt StdMeshers_PrismAsBlock::TSideFace::Value(const Standard_Real U, } else { - TopoDS_Shape s = myHelper->GetSubShapeByNode( nn[0], myHelper->GetMeshDS() ); + TopoDS_Shape s = myHelper.GetSubShapeByNode( nn[0], myHelper.GetMeshDS() ); if ( s.ShapeType() != TopAbs_EDGE ) - s = myHelper->GetSubShapeByNode( nn[2], myHelper->GetMeshDS() ); + s = myHelper.GetSubShapeByNode( nn[2], myHelper.GetMeshDS() ); if ( s.ShapeType() == TopAbs_EDGE ) edge = TopoDS::Edge( s ); } if ( !edge.IsNull() ) { - double u1 = myHelper->GetNodeU( edge, nn[0] ); - double u3 = myHelper->GetNodeU( edge, nn[2] ); + double u1 = myHelper.GetNodeU( edge, nn[0] ); + double u3 = myHelper.GetNodeU( edge, nn[2] ); double u = u1 * ( 1 - hR ) + u3 * hR; TopLoc_Location loc; double f,l; Handle(Geom_Curve) curve = BRep_Tool::Curve( edge,loc,f,l ); @@ -2845,10 +3017,10 @@ gp_Pnt StdMeshers_PrismAsBlock::TSideFace::Value(const Standard_Real U, } if ( notFaceID2 ) // no nodes of FACE and nodes are on different FACEs { - SMESHDS_Mesh* meshDS = myHelper->GetMeshDS(); - TopoDS_Shape face = myHelper->GetCommonAncestor( meshDS->IndexToShape( notFaceID1 ), + SMESHDS_Mesh* meshDS = myHelper.GetMeshDS(); + TopoDS_Shape face = myHelper.GetCommonAncestor( meshDS->IndexToShape( notFaceID1 ), meshDS->IndexToShape( notFaceID2 ), - *myHelper->GetMesh(), + *myHelper.GetMesh(), TopAbs_FACE ); if ( face.IsNull() ) throw SALOME_Exception("StdMeshers_PrismAsBlock::TSideFace::Value() face.IsNull()"); @@ -2858,13 +3030,14 @@ gp_Pnt StdMeshers_PrismAsBlock::TSideFace::Value(const Standard_Real U, throw SALOME_Exception("StdMeshers_PrismAsBlock::TSideFace::Value() !mySurface"); } } - - gp_XY uv1 = myHelper->GetNodeUV( mySurface->Face(), nn[0], nn[2]); - gp_XY uv2 = myHelper->GetNodeUV( mySurface->Face(), nn[1], nn[3]); + ((TSideFace*) this)->myHelper.SetSubShape( mySurface->Face() ); + + gp_XY uv1 = myHelper.GetNodeUV( mySurface->Face(), nn[0], nn[2]); + gp_XY uv2 = myHelper.GetNodeUV( mySurface->Face(), nn[1], nn[3]); gp_XY uv12 = uv1 * ( 1 - vR ) + uv2 * vR; - gp_XY uv3 = myHelper->GetNodeUV( mySurface->Face(), nn[2], nn[0]); - gp_XY uv4 = myHelper->GetNodeUV( mySurface->Face(), nn[3], nn[1]); + gp_XY uv3 = myHelper.GetNodeUV( mySurface->Face(), nn[2], nn[0]); + gp_XY uv4 = myHelper.GetNodeUV( mySurface->Face(), nn[3], nn[1]); gp_XY uv34 = uv3 * ( 1 - vR ) + uv4 * vR; gp_XY uv = uv12 * ( 1 - hR ) + uv34 * hR; @@ -2893,7 +3066,7 @@ TopoDS_Edge StdMeshers_PrismAsBlock::TSideFace::GetEdge(const int iEdge) const } TopoDS_Shape edge; const SMDS_MeshNode* node = 0; - SMESHDS_Mesh * meshDS = myHelper->GetMesh()->GetMeshDS(); + SMESHDS_Mesh * meshDS = myHelper.GetMesh()->GetMeshDS(); TNodeColumn* column; switch ( iEdge ) { @@ -2901,7 +3074,7 @@ TopoDS_Edge StdMeshers_PrismAsBlock::TSideFace::GetEdge(const int iEdge) const case BOTTOM_EDGE: column = & (( ++myParamToColumnMap->begin())->second ); node = ( iEdge == TOP_EDGE ) ? column->back() : column->front(); - edge = myHelper->GetSubShapeByNode ( node, meshDS ); + edge = myHelper.GetSubShapeByNode ( node, meshDS ); if ( edge.ShapeType() == TopAbs_VERTEX ) { column = & ( myParamToColumnMap->begin()->second ); node = ( iEdge == TOP_EDGE ) ? column->back() : column->front(); @@ -2916,7 +3089,7 @@ TopoDS_Edge StdMeshers_PrismAsBlock::TSideFace::GetEdge(const int iEdge) const else column = & ( myParamToColumnMap->begin()->second ); if ( column->size() > 0 ) - edge = myHelper->GetSubShapeByNode( (*column)[ 1 ], meshDS ); + edge = myHelper.GetSubShapeByNode( (*column)[ 1 ], meshDS ); if ( edge.IsNull() || edge.ShapeType() == TopAbs_VERTEX ) node = column->front(); break; @@ -2928,10 +3101,10 @@ TopoDS_Edge StdMeshers_PrismAsBlock::TSideFace::GetEdge(const int iEdge) const // find edge by 2 vertices TopoDS_Shape V1 = edge; - TopoDS_Shape V2 = myHelper->GetSubShapeByNode( node, meshDS ); + TopoDS_Shape V2 = myHelper.GetSubShapeByNode( node, meshDS ); if ( !V2.IsNull() && V2.ShapeType() == TopAbs_VERTEX && !V2.IsSame( V1 )) { - TopoDS_Shape ancestor = myHelper->GetCommonAncestor( V1, V2, *myHelper->GetMesh(), TopAbs_EDGE); + TopoDS_Shape ancestor = myHelper.GetCommonAncestor( V1, V2, *myHelper.GetMesh(), TopAbs_EDGE); if ( !ancestor.IsNull() ) return TopoDS::Edge( ancestor ); } @@ -2971,8 +3144,8 @@ int StdMeshers_PrismAsBlock::TSideFace::InsertSubShapes(TBlockShapes& shapeMap) GetColumns(0, col1, col2 ); const SMDS_MeshNode* node0 = col1->second.front(); const SMDS_MeshNode* node1 = col1->second.back(); - TopoDS_Shape v0 = myHelper->GetSubShapeByNode( node0, myHelper->GetMeshDS()); - TopoDS_Shape v1 = myHelper->GetSubShapeByNode( node1, myHelper->GetMeshDS()); + TopoDS_Shape v0 = myHelper.GetSubShapeByNode( node0, myHelper.GetMeshDS()); + TopoDS_Shape v1 = myHelper.GetSubShapeByNode( node1, myHelper.GetMeshDS()); if ( v0.ShapeType() == TopAbs_VERTEX ) { nbInserted += SMESH_Block::Insert( v0, vertIdVec[ 0 ], shapeMap); } @@ -2985,8 +3158,8 @@ int StdMeshers_PrismAsBlock::TSideFace::InsertSubShapes(TBlockShapes& shapeMap) GetColumns(1, col1, col2 ); node0 = col2->second.front(); node1 = col2->second.back(); - v0 = myHelper->GetSubShapeByNode( node0, myHelper->GetMeshDS()); - v1 = myHelper->GetSubShapeByNode( node1, myHelper->GetMeshDS()); + v0 = myHelper.GetSubShapeByNode( node0, myHelper.GetMeshDS()); + v1 = myHelper.GetSubShapeByNode( node1, myHelper.GetMeshDS()); if ( v0.ShapeType() == TopAbs_VERTEX ) { nbInserted += SMESH_Block::Insert( v0, vertIdVec[ 0 ], shapeMap); } @@ -3178,7 +3351,7 @@ gp_Pnt2d StdMeshers_PrismAsBlock::TPCurveOnHorFaceAdaptor::Value(const Standard_ { TParam2ColumnIt u_col1, u_col2; double r = mySide->GetColumns( U, u_col1, u_col2 ); - gp_XY uv1 = mySide->GetNodeUV( myFace, u_col1->second[ myZ ]); - gp_XY uv2 = mySide->GetNodeUV( myFace, u_col2->second[ myZ ]); + gp_XY uv1 = mySide->GetNodeUV( myFace, u_col1->second[ myZ ], u_col2->second[ myZ ]); + gp_XY uv2 = mySide->GetNodeUV( myFace, u_col2->second[ myZ ], u_col1->second[ myZ ]); return uv1 * ( 1 - r ) + uv2 * r; }