From 096485adbe5e66ace1182bf6c49dc6853f437d9c Mon Sep 17 00:00:00 2001 From: eap Date: Mon, 17 Apr 2017 15:42:28 +0300 Subject: [PATCH] IPAL52499: Prismatic mesh is not computed on a prismatic shape --- src/SMESH/SMESH_MesherHelper.cxx | 2 +- src/SMESHGUI/SMESHGUI.cxx | 2 +- src/SMESH_I/SMESH_Filter_i.cxx | 2 +- src/SMESH_I/SMESH_Gen_i.cxx | 2 +- src/SMESH_SWIG/smeshBuilder.py | 7 +- src/StdMeshers/StdMeshers_FaceSide.cxx | 16 +- src/StdMeshers/StdMeshers_FaceSide.hxx | 34 +- src/StdMeshers/StdMeshers_Prism_3D.cxx | 348 +++++++++++++++++- src/StdMeshers/StdMeshers_Prism_3D.hxx | 24 +- src/StdMeshers/StdMeshers_ProjectionUtils.cxx | 94 ++--- src/StdMeshers/StdMeshers_ProjectionUtils.hxx | 25 +- src/StdMeshers/StdMeshers_Projection_2D.cxx | 2 +- 12 files changed, 462 insertions(+), 96 deletions(-) diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index e03128dc8..004c856f1 100644 --- a/src/SMESH/SMESH_MesherHelper.cxx +++ b/src/SMESH/SMESH_MesherHelper.cxx @@ -2669,7 +2669,7 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2 } // nb rows of nodes - size_t prevNbRows = theParam2ColumnMap.begin()->second.size(); // current, at least 1 here + size_t prevNbRows = theParam2ColumnMap.begin()->second.size(); // current, at least 1 here size_t expectNbRows = faceSubMesh->NbElements() / ( theParam2ColumnMap.size()-1 ); // to be added // fill theParam2ColumnMap column by column by passing from nodes on diff --git a/src/SMESHGUI/SMESHGUI.cxx b/src/SMESHGUI/SMESHGUI.cxx index 526b6c650..1145efb17 100644 --- a/src/SMESHGUI/SMESHGUI.cxx +++ b/src/SMESHGUI/SMESHGUI.cxx @@ -630,7 +630,7 @@ namespace }; // is typeMsg complete? (compilation failure mains that enum SMDSAbs_EntityType changed) const int nbTypes = sizeof( typeMsg ) / sizeof( const char* ); - int _assert[( nbTypes == SMESH::Entity_Last ) ? 2 : -1 ]; _assert[0]=_assert[1]; + int _assert[( nbTypes == SMESH::Entity_Last ) ? 2 : -1 ]; _assert[0]=_assert[1]=0; QString andStr = " " + QObject::tr("SMESH_AND") + " ", comma(", "); for ( size_t iType = 0; iType < presentNotSupported.size(); ++iType ) { diff --git a/src/SMESH_I/SMESH_Filter_i.cxx b/src/SMESH_I/SMESH_Filter_i.cxx index 739a4d904..9432b9117 100644 --- a/src/SMESH_I/SMESH_Filter_i.cxx +++ b/src/SMESH_I/SMESH_Filter_i.cxx @@ -4109,7 +4109,7 @@ static const char** getFunctNames() #ifdef _DEBUG_ // check if functName is complete, compilation failure means that enum FunctorType changed const int nbFunctors = sizeof(functName) / sizeof(const char*); - int _assert[( nbFunctors == SMESH::FT_Undefined + 1 ) ? 2 : -1 ]; _assert[0]=_assert[1]; + int _assert[( nbFunctors == SMESH::FT_Undefined + 1 ) ? 2 : -1 ]; _assert[0]=_assert[1]=0; #endif return functName; diff --git a/src/SMESH_I/SMESH_Gen_i.cxx b/src/SMESH_I/SMESH_Gen_i.cxx index 44cbff77a..916a1f222 100644 --- a/src/SMESH_I/SMESH_Gen_i.cxx +++ b/src/SMESH_I/SMESH_Gen_i.cxx @@ -2560,7 +2560,7 @@ SMESH_Gen_i::ConcatenateCommon(const SMESH::ListOfIDSources& theMeshesArray, const char* typeNames[] = { "All","Nodes","Edges","Faces","Volumes","0DElems","Balls" }; { // check of typeNames: compilation failure mains that NB_ELEMENT_TYPES changed: const int nbNames = sizeof(typeNames) / sizeof(const char*); - int _assert[( nbNames == SMESH::NB_ELEMENT_TYPES ) ? 2 : -1 ]; _assert[0]=_assert[1]; + int _assert[( nbNames == SMESH::NB_ELEMENT_TYPES ) ? 2 : -1 ]; _assert[0]=_assert[1]=0; } string groupName = "Gr"; SALOMEDS::SObject_wrap aMeshSObj = ObjectToSObject( myCurrentStudy, theMeshesArray[i] ); diff --git a/src/SMESH_SWIG/smeshBuilder.py b/src/SMESH_SWIG/smeshBuilder.py index d87d0d553..900aecf18 100644 --- a/src/SMESH_SWIG/smeshBuilder.py +++ b/src/SMESH_SWIG/smeshBuilder.py @@ -2687,9 +2687,14 @@ class Mesh: ## Return an element based on all given nodes. # @ingroup l1_meshinfo - def FindElementByNodes(self,nodes): + def FindElementByNodes(self, nodes): return self.mesh.FindElementByNodes(nodes) + ## Return elements including all given nodes. + # @ingroup l1_meshinfo + def GetElementsByNodes(self, nodes, elemType=SMESH.ALL): + return self.mesh.GetElementsByNodes( nodes, elemType ) + ## Return true if the given element is a polygon # @ingroup l1_meshinfo def IsPoly(self, id): diff --git a/src/StdMeshers/StdMeshers_FaceSide.cxx b/src/StdMeshers/StdMeshers_FaceSide.cxx index 9f5929e36..1424e4226 100644 --- a/src/StdMeshers/StdMeshers_FaceSide.cxx +++ b/src/StdMeshers/StdMeshers_FaceSide.cxx @@ -88,13 +88,13 @@ StdMeshers_FaceSide::StdMeshers_FaceSide(const TopoDS_Face& theFace, */ //================================================================================ -StdMeshers_FaceSide::StdMeshers_FaceSide(const TopoDS_Face& theFace, - std::list& theEdges, - SMESH_Mesh* theMesh, - const bool theIsForward, - const bool theIgnoreMediumNodes, - SMESH_MesherHelper* theFaceHelper, - SMESH_ProxyMesh::Ptr theProxyMesh) +StdMeshers_FaceSide::StdMeshers_FaceSide(const TopoDS_Face& theFace, + const std::list& theEdges, + SMESH_Mesh* theMesh, + const bool theIsForward, + const bool theIgnoreMediumNodes, + SMESH_MesherHelper* theFaceHelper, + SMESH_ProxyMesh::Ptr theProxyMesh) { int nbEdges = theEdges.size(); myEdge.resize ( nbEdges ); @@ -125,7 +125,7 @@ StdMeshers_FaceSide::StdMeshers_FaceSide(const TopoDS_Face& theFace, SMESHDS_Mesh* meshDS = myProxyMesh->GetMeshDS(); int nbDegen = 0; - std::list::iterator edge = theEdges.begin(); + std::list::const_iterator edge = theEdges.begin(); for ( int index = 0; edge != theEdges.end(); ++index, ++edge ) { int i = theIsForward ? index : nbEdges-index-1; diff --git a/src/StdMeshers/StdMeshers_FaceSide.hxx b/src/StdMeshers/StdMeshers_FaceSide.hxx index 6759f50a9..1964638fe 100644 --- a/src/StdMeshers/StdMeshers_FaceSide.hxx +++ b/src/StdMeshers/StdMeshers_FaceSide.hxx @@ -82,13 +82,13 @@ public: /*! * \brief Wrap several edges. Edges must be properly ordered and oriented. */ - StdMeshers_FaceSide(const TopoDS_Face& theFace, - std::list& theEdges, - SMESH_Mesh* theMesh, - const bool theIsForward, - const bool theIgnoreMediumNodes, - SMESH_MesherHelper* theFaceHelper = NULL, - SMESH_ProxyMesh::Ptr theProxyMesh = SMESH_ProxyMesh::Ptr()); + StdMeshers_FaceSide(const TopoDS_Face& theFace, + const std::list& theEdges, + SMESH_Mesh* theMesh, + const bool theIsForward, + const bool theIgnoreMediumNodes, + SMESH_MesherHelper* theFaceHelper = NULL, + SMESH_ProxyMesh::Ptr theProxyMesh = SMESH_ProxyMesh::Ptr()); /*! * \brief Simulate a side from a vertex using data from other FaceSide */ @@ -120,13 +120,13 @@ public: { return StdMeshers_FaceSidePtr ( new StdMeshers_FaceSide( Face,Edge,Mesh,IsForward,IgnoreMediumNodes,FaceHelper,ProxyMesh )); } - static StdMeshers_FaceSidePtr New (const TopoDS_Face& Face, - std::list& Edges, - SMESH_Mesh* Mesh, - const bool IsForward, - const bool IgnoreMediumNodes, - SMESH_MesherHelper* FaceHelper = NULL, - SMESH_ProxyMesh::Ptr ProxyMesh = SMESH_ProxyMesh::Ptr()) + static StdMeshers_FaceSidePtr New (const TopoDS_Face& Face, + const std::list& Edges, + SMESH_Mesh* Mesh, + const bool IsForward, + const bool IgnoreMediumNodes, + SMESH_MesherHelper* FaceHelper = NULL, + SMESH_ProxyMesh::Ptr ProxyMesh = SMESH_ProxyMesh::Ptr()) { return StdMeshers_FaceSidePtr ( new StdMeshers_FaceSide( Face,Edges,Mesh,IsForward,IgnoreMediumNodes,FaceHelper,ProxyMesh )); } @@ -141,9 +141,11 @@ public: ( new StdMeshers_FaceSide( Side,Node,Pnt2d1,Pnt2d2,C2d,UFirst,ULast )); } static StdMeshers_FaceSidePtr New (UVPtStructVec& theSideNodes, - const TopoDS_Face& theFace = TopoDS_Face()) + const TopoDS_Face& theFace = TopoDS_Face(), + const TopoDS_Edge& theEdge = TopoDS_Edge(), + SMESH_Mesh* theMesh = 0) { - return StdMeshers_FaceSidePtr( new StdMeshers_FaceSide( theSideNodes, theFace )); + return StdMeshers_FaceSidePtr( new StdMeshers_FaceSide( theSideNodes, theFace, theEdge, theMesh )); } /*! diff --git a/src/StdMeshers/StdMeshers_Prism_3D.cxx b/src/StdMeshers/StdMeshers_Prism_3D.cxx index a7dae0283..5e253829a 100644 --- a/src/StdMeshers/StdMeshers_Prism_3D.cxx +++ b/src/StdMeshers/StdMeshers_Prism_3D.cxx @@ -51,6 +51,7 @@ #include #include #include +#include #include #include #include @@ -1173,7 +1174,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; + myUseBlock = false; // is set to true if projection is done using "block approach" myBotToColumnMap.clear(); if ( !assocOrProjBottom2Top( bottomToTopTrsf, thePrism ) ) // it also fills myBotToColumnMap return false; @@ -1181,38 +1182,60 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) // Create nodes inside the block - // use transformation (issue 0020680, IPAL0052499) - StdMeshers_Sweeper sweeper; - double tol; - bool allowHighBndError; - if ( !myUseBlock ) { + // use transformation (issue 0020680, IPAL0052499) or a "straight line" approach + StdMeshers_Sweeper sweeper; + // load boundary nodes into sweeper bool dummy; + const SMDS_MeshNode* prevN0 = 0, *prevN1 = 0; 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 ) + + 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; + + for ( ; u2colIt != u2colEnd; ++u2colIt ) sweeper.myBndColumns.push_back( & u2colIt->second ); } - // load node columns inside the bottom face + // load node columns inside the bottom FACE TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin(); + sweeper.myIntColumns.reserve( myBotToColumnMap.size() ); for ( ; bot_column != myBotToColumnMap.end(); ++bot_column ) sweeper.myIntColumns.push_back( & bot_column->second ); - tol = getSweepTolerance( thePrism ); - allowHighBndError = !isSimpleBottom( thePrism ); + 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 ); + } + else if ( sweeper.CheckSameZ() ) + { + myUseBlock = !sweeper.ComputeNodesOnStraightSameZ( *myHelper ); + } + else + { + myUseBlock = !sweeper.ComputeNodesOnStraight( *myHelper, thePrism.myBottom, thePrism.myTop ); + } + myHelper->SetElementsOnShape( false ); } - if ( !myUseBlock && sweeper.ComputeNodes( *myHelper, tol, allowHighBndError )) - { - } - else // use block approach + if ( myUseBlock ) // use block approach { // loop on nodes inside the bottom face Prism_3D::TNode prevBNode; @@ -2400,7 +2423,7 @@ double StdMeshers_Prism_3D::getSweepTolerance( const Prism_3D::TPrismTopo& thePr bool StdMeshers_Prism_3D::isSimpleBottom( const Prism_3D::TPrismTopo& thePrism ) { - if ( thePrism.myBottomEdges.size() != 4 ) + if ( thePrism.myNbEdgesInWires.front() != 4 ) return false; // analyse angles between edges @@ -2432,6 +2455,43 @@ bool StdMeshers_Prism_3D::isSimpleBottom( const Prism_3D::TPrismTopo& thePrism ) return true; } +//======================================================================= +//function : allVerticalEdgesStraight +//purpose : Defines if all "vertical" EDGEs are straight +//======================================================================= + +bool StdMeshers_Prism_3D::allVerticalEdgesStraight( const Prism_3D::TPrismTopo& thePrism ) +{ + for ( size_t i = 0; i < thePrism.myWallQuads.size(); ++i ) + { + const Prism_3D::TQuadList& quads = thePrism.myWallQuads[i]; + Prism_3D::TQuadList::const_iterator quadIt = quads.begin(); + TopoDS_Edge prevQuadEdge; + for ( ; quadIt != quads.end(); ++quadIt ) + { + StdMeshers_FaceSidePtr rightSide = (*quadIt)->side[ QUAD_RIGHT_SIDE ]; + + if ( !prevQuadEdge.IsNull() && + !SMESH_Algo::IsContinuous( rightSide->Edge( 0 ), prevQuadEdge )) + return false; + + for ( int iE = 0; iE < rightSide->NbEdges(); ++iE ) + { + const TopoDS_Edge & rightE = rightSide->Edge( iE ); + if ( !SMESH_Algo::IsStraight( rightE, /*degenResult=*/true )) + return false; + + if ( iE > 0 && + !SMESH_Algo::IsContinuous( rightSide->Edge( iE-1 ), rightE )) + return false; + + prevQuadEdge = rightE; + } + } + } + return true; +} + //======================================================================= //function : project2dMesh //purpose : Project mesh faces from a source FACE of one prism (theSrcFace) @@ -4732,7 +4792,7 @@ void StdMeshers_Sweeper::applyBoundaryError(const vector< gp_XYZ >& bndPoints, //================================================================================ /*! - * \brief Creates internal nodes of the prism + * \brief Create internal nodes of the prism */ //================================================================================ @@ -5012,3 +5072,257 @@ bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, return true; } + +//================================================================================ +/*! + * \brief Check if all nodes of each layers have same logical Z + */ +//================================================================================ + +bool StdMeshers_Sweeper::CheckSameZ() +{ + myZColumns.resize( myBndColumns.size() ); + fillZColumn( myZColumns[0], *myBndColumns[0] ); + + bool sameZ = true; + const double tol = 0.1 * 1./ myBndColumns[0]->size(); + + // check columns based on VERTEXes + + vector< int > vertexIndex; + vertexIndex.push_back( 0 ); + for ( size_t iC = 1; iC < myBndColumns.size() && sameZ; ++iC ) + { + if ( myBndColumns[iC]->front()->GetPosition()->GetDim() > 0 ) + continue; // not on VERTEX + + vertexIndex.push_back( iC ); + fillZColumn( myZColumns[iC], *myBndColumns[iC] ); + + for ( size_t iZ = 0; iZ < myZColumns[0].size() && sameZ; ++iZ ) + sameZ = ( Abs( myZColumns[0][iZ] - myZColumns[iC][iZ]) < tol ); + } + + // check columns based on EDGEs, one per EDGE + + for ( size_t i = 1; i < vertexIndex.size() && sameZ; ++i ) + { + if ( vertexIndex[i] - vertexIndex[i-1] < 2 ) + continue; + + int iC = ( vertexIndex[i] + vertexIndex[i-1] ) / 2; + fillZColumn( myZColumns[iC], *myBndColumns[iC] ); + + for ( size_t iZ = 0; iZ < myZColumns[0].size() && sameZ; ++iZ ) + sameZ = ( Abs( myZColumns[0][iZ] - myZColumns[iC][iZ]) < tol ); + } + + if ( sameZ ) + { + myZColumns.resize(1); + } + else + { + for ( size_t iC = 1; iC < myBndColumns.size(); ++iC ) + fillZColumn( myZColumns[iC], *myBndColumns[iC] ); + } + + return sameZ; +} + +//================================================================================ +/*! + * \brief Create internal nodes of the prism all located on straight lines with + * the same distribution along the lines. + */ +//================================================================================ + +bool StdMeshers_Sweeper::ComputeNodesOnStraightSameZ( SMESH_MesherHelper& helper ) +{ + TZColumn& z = myZColumns[0]; + + for ( size_t i = 0; i < myIntColumns.size(); ++i ) + { + TNodeColumn& nodes = *myIntColumns[i]; + SMESH_NodeXYZ n0( nodes[0] ), n1( nodes.back() ); + + 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() ); + } + } + + return true; +} + +//================================================================================ +/*! + * \brief Create internal nodes of the prism all located on straight lines with + * different distributions along the lines. + */ +//================================================================================ + +bool StdMeshers_Sweeper::ComputeNodesOnStraight( SMESH_MesherHelper& helper, + const TopoDS_Face& botFace, + const TopoDS_Face& topFace ) +{ + // get data to create a Morph + UVPtStructVec botUV( myBndColumns.size() + 1 ); + UVPtStructVec topUV( myBndColumns.size() + 1 ); + 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] )); + topUV[i].node = nodes.back(); + topUV[i].SetUV( helper.GetNodeUV( topFace, 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() )); + + // 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 ); + + if ( helper.GetIsQuadratic() ) + { + // mark all medium nodes of faces on botFace to avoid their treating + SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( botFace ); + SMDS_ElemIteratorPtr eIt = smDS->GetElements(); + while ( eIt->more() ) + { + const SMDS_MeshElement* e = eIt->next(); + for ( int i = e->NbCornerNodes(), nb = e->NbNodes(); i < nb; ++i ) + e->GetNode( i )->setIsMarked( true ); + } + } + + // map to get a node column by a bottom node + TColStd_DataMapOfIntegerInteger iNode2iCol( myIntColumns.size() ); + + // un-mark nodes to treat (internal bottom nodes); later we will mark treated nodes + for ( size_t i = 0; i < myIntColumns.size(); ++i ) + { + const SMDS_MeshNode* botNode = myIntColumns[i]->front(); + botNode->setIsMarked( false ); + iNode2iCol.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; + + 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 ); + + TNodeColumn* column = myIntColumns[ iNode2iCol( botNode->GetID() )]; + + // 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; + + // 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; + + // 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() ); + } + + // add neighbor nodes to the queue + botDelauney.AddCloseNodes( botNode, botTria, botFaceID, botNoTriQueue ); + } + + 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 )); + } + } + } + } + + return true; +} + +//================================================================================ +/*! + * \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 + } +} diff --git a/src/StdMeshers/StdMeshers_Prism_3D.hxx b/src/StdMeshers/StdMeshers_Prism_3D.hxx index ddfbff96a..48c75ba0e 100644 --- a/src/StdMeshers/StdMeshers_Prism_3D.hxx +++ b/src/StdMeshers/StdMeshers_Prism_3D.hxx @@ -29,15 +29,12 @@ #include "SMESH_StdMeshers.hxx" -#include "SMDS_MeshNode.hxx" -#include "SMDS_TypeOfPosition.hxx" #include "SMESHDS_Mesh.hxx" -#include "SMESH_Algo.hxx" #include "SMESH_Block.hxx" -#include "SMESH_Comment.hxx" #include "SMESH_Mesh.hxx" #include "SMESH_MesherHelper.hxx" #include "SMESH_TypeDefs.hxx" +#include "SMESH_Comment.hxx" #include "SMESH_subMesh.hxx" #include @@ -421,10 +418,21 @@ struct StdMeshers_Sweeper std::vector< TNodeColumn* > myBndColumns; // boundary nodes std::vector< TNodeColumn* > myIntColumns; // internal nodes + typedef std::vector< double > TZColumn; + std::vector< TZColumn > myZColumns; // Z distribution of boundary nodes + bool ComputeNodes( SMESH_MesherHelper& helper, const double tol, const bool allowHighBndError ); + bool CheckSameZ(); + + bool ComputeNodesOnStraightSameZ( SMESH_MesherHelper& helper ); + + bool ComputeNodesOnStraight( SMESH_MesherHelper& helper, + const TopoDS_Face& bottom, + const TopoDS_Face& top); + private: gp_XYZ bndPoint( int iP, int z ) const @@ -446,6 +454,9 @@ private: const double r, std::vector< gp_XYZ >& toIntPoints, std::vector< double >& int2BndDist); + + static void fillZColumn( TZColumn& zColumn, + TNodeColumn& nodes ); }; // =============================================== @@ -545,6 +556,11 @@ public: */ bool isSimpleBottom( const Prism_3D::TPrismTopo& thePrism ); + /*! + * \brief Defines if all "vertical" EDGEs are straight + */ + bool allVerticalEdgesStraight( const Prism_3D::TPrismTopo& thePrism ); + /*! * \brief Project mesh faces from a source FACE of one prism to * a source FACE of another prism diff --git a/src/StdMeshers/StdMeshers_ProjectionUtils.cxx b/src/StdMeshers/StdMeshers_ProjectionUtils.cxx index e7051db9b..bca7292ff 100644 --- a/src/StdMeshers/StdMeshers_ProjectionUtils.cxx +++ b/src/StdMeshers/StdMeshers_ProjectionUtils.cxx @@ -2796,12 +2796,10 @@ namespace StdMeshers_ProjectionUtils */ //================================================================================ - typedef list< pair< const SMDS_MeshNode*, const BRepMesh_Triangle* > > TNodeTriaList; - - void addCloseNodes( const SMDS_MeshNode* srcNode, - const BRepMesh_Triangle* bmTria, - const int srcFaceID, - TNodeTriaList & noTriQueue ) + void Morph::AddCloseNodes( const SMDS_MeshNode* srcNode, + const BRepMesh_Triangle* bmTria, + const int srcFaceID, + TNodeTriaList & noTriQueue ) { // find in-FACE nodes SMDS_ElemIteratorPtr elems = srcNode->GetInverseElementIterator(SMDSAbs_Face); @@ -2813,7 +2811,7 @@ namespace StdMeshers_ProjectionUtils for ( int i = 0, nb = elem->NbNodes(); i < nb; ++i ) { const SMDS_MeshNode* n = elem->GetNode( i ); - if ( !n->isMarked() ) + if ( !n->isMarked() /*&& n->getshapeId() == srcFaceID*/ ) noTriQueue.push_back( make_pair( n, bmTria )); } } @@ -2827,10 +2825,10 @@ namespace StdMeshers_ProjectionUtils */ //================================================================================ - const BRepMesh_Triangle* findTriangle( const gp_XY& uv, - const BRepMesh_Triangle* bmTria, - Handle(BRepMesh_DataStructureOfDelaun)& triaDS, - double bc[3] ) + const BRepMesh_Triangle* Morph::FindTriangle( const gp_XY& uv, + const BRepMesh_Triangle* bmTria, + double bc[3], + int triaNodes[3] ) { int nodeIDs[3]; gp_XY nodeUVs[3]; @@ -2841,10 +2839,10 @@ namespace StdMeshers_ProjectionUtils { // check bmTria - triaDS->ElementNodes( *bmTria, nodeIDs ); - nodeUVs[0] = triaDS->GetNode( nodeIDs[0] ).Coord(); - nodeUVs[1] = triaDS->GetNode( nodeIDs[1] ).Coord(); - nodeUVs[2] = triaDS->GetNode( nodeIDs[2] ).Coord(); + _triaDS->ElementNodes( *bmTria, nodeIDs ); + nodeUVs[0] = _triaDS->GetNode( nodeIDs[0] ).Coord(); + nodeUVs[1] = _triaDS->GetNode( nodeIDs[1] ).Coord(); + nodeUVs[2] = _triaDS->GetNode( nodeIDs[2] ).Coord(); SMESH_MeshAlgos::GetBarycentricCoords( uv, nodeUVs[0], nodeUVs[1], nodeUVs[2], @@ -2852,6 +2850,9 @@ namespace StdMeshers_ProjectionUtils if ( bc[0] >= 0 && bc[1] >= 0 && bc[0] + bc[1] <= 1 ) { bc[2] = 1 - bc[0] - bc[1]; + triaNodes[0] = nodeIDs[0]; + triaNodes[1] = nodeIDs[1]; + triaNodes[2] = nodeIDs[2]; return bmTria; } @@ -2862,19 +2863,19 @@ namespace StdMeshers_ProjectionUtils gp_XY seg = uv - gc; bmTria->Edges( linkIDs, ori ); - int triaID = triaDS->IndexOf( *bmTria ); + int triaID = _triaDS->IndexOf( *bmTria ); bmTria = 0; for ( int i = 0; i < 3; ++i ) { - const BRepMesh_PairOfIndex & triIDs = triaDS->ElementsConnectedTo( linkIDs[i] ); + const BRepMesh_PairOfIndex & triIDs = _triaDS->ElementsConnectedTo( linkIDs[i] ); if ( triIDs.Extent() < 2 ) continue; // no neighbor triangle // check if a link intersects gc2uv - const BRepMesh_Edge & link = triaDS->GetLink( linkIDs[i] ); - const BRepMesh_Vertex & n1 = triaDS->GetNode( link.FirstNode() ); - const BRepMesh_Vertex & n2 = triaDS->GetNode( link.LastNode() ); + const BRepMesh_Edge & link = _triaDS->GetLink( linkIDs[i] ); + const BRepMesh_Vertex & n1 = _triaDS->GetNode( link.FirstNode() ); + const BRepMesh_Vertex & n2 = _triaDS->GetNode( link.LastNode() ); gp_XY uv1 = n1.Coord(); gp_XY lin = n2.Coord() - uv1; // link direction @@ -2885,7 +2886,7 @@ namespace StdMeshers_ProjectionUtils double uSeg = ( uv1 - gc ) ^ lin / crossSegLin; if ( 0. <= uSeg && uSeg <= 1. ) { - bmTria = & triaDS->GetElement( triIDs.Index( 1 + ( triIDs.Index(1) == triaID ))); + bmTria = & _triaDS->GetElement( triIDs.Index( 1 + ( triIDs.Index(1) == triaID ))); break; } } @@ -2893,6 +2894,22 @@ namespace StdMeshers_ProjectionUtils return bmTria; } + //================================================================================ + /*! + * \brief Return a triangle sharing a given boundary node + * \param [in] iBndNode - index of the boundary node + * \return const BRepMesh_Triangle* - a found triangle + */ + //================================================================================ + + const BRepMesh_Triangle* Morph::GetTriangleNear( int iBndNode ) + { + const BRepMesh::ListOfInteger & linkIds = _triaDS->LinksConnectedTo( iBndNode ); + const BRepMesh_PairOfIndex & triaIds = _triaDS->ElementsConnectedTo( linkIds.First() ); + const BRepMesh_Triangle& tria = _triaDS->GetElement( triaIds.Index(1) ); + return &tria; + } + //================================================================================ /*! * \brief triangulate the srcFace in 2D @@ -2997,19 +3014,6 @@ namespace StdMeshers_ProjectionUtils const SMDS_MeshNode *srcNode, *tgtNode; const BRepMesh_Triangle *bmTria; - // initialize a queue of nodes with starting triangles - TNodeTriaList noTriQueue; - size_t iBndSrcN = 1; - for ( ; iBndSrcN < _bndSrcNodes.size() && noTriQueue.empty(); ++iBndSrcN ) - { - // get a triangle - const BRepMesh::ListOfInteger & linkIds = _triaDS->LinksConnectedTo( iBndSrcN ); - const BRepMesh_PairOfIndex & triaIds = _triaDS->ElementsConnectedTo( linkIds.First() ); - const BRepMesh_Triangle& tria = _triaDS->GetElement( triaIds.Index(1) ); - - addCloseNodes( _bndSrcNodes[ iBndSrcN ], &tria, srcFaceID, noTriQueue ); - } - // un-mark internal src nodes; later we will mark moved nodes int nbSrcNodes = 0; SMDS_NodeIteratorPtr nIt = _srcSubMesh->GetSubMeshDS()->GetNodes(); @@ -3033,6 +3037,10 @@ namespace StdMeshers_ProjectionUtils bool checkUV = true; const SMDS_FacePosition* pos; + // a queue of nodes with starting triangles + TNodeTriaList noTriQueue; + size_t iBndSrcN = 1; + while ( nbSrcNodes > 0 ) { while ( !noTriQueue.empty() ) @@ -3048,16 +3056,16 @@ namespace StdMeshers_ProjectionUtils // find a delauney triangle containing the src node gp_XY uv = tgtHelper.GetNodeUV( srcFace, srcNode, NULL, &checkUV ); uv *= _scale; - bmTria = findTriangle( uv, bmTria, _triaDS, bc ); + bmTria = FindTriangle( uv, bmTria, bc, nodeIDs ); if ( !bmTria ) continue; // compute new coordinates for a corresponding tgt node gp_XY uvNew( 0., 0. ), nodeUV; - _triaDS->ElementNodes( *bmTria, nodeIDs ); for ( int i = 0; i < 3; ++i ) uvNew += bc[i] * tgtVert( nodeIDs[i]).Coord(); - uvNew.SetCoord( uvNew.X() / _scale.X(), uvNew.Y() / _scale.Y() ); + uvNew.SetCoord( uvNew.X() / _scale.X(), + uvNew.Y() / _scale.Y() ); gp_Pnt xyz = tgtSurface->Value( uvNew ); // find and move tgt node @@ -3069,7 +3077,7 @@ namespace StdMeshers_ProjectionUtils if (( pos = dynamic_cast< const SMDS_FacePosition* >( tgtNode->GetPosition() ))) const_cast( pos )->SetParameters( uvNew.X(), uvNew.Y() ); - addCloseNodes( srcNode, bmTria, srcFaceID, noTriQueue ); + AddCloseNodes( srcNode, bmTria, srcFaceID, noTriQueue ); } if ( nbSrcNodes > 0 ) @@ -3077,10 +3085,8 @@ namespace StdMeshers_ProjectionUtils // assure that all src nodes are visited for ( ; iBndSrcN < _bndSrcNodes.size() && noTriQueue.empty(); ++iBndSrcN ) { - const BRepMesh::ListOfInteger & linkIds = _triaDS->LinksConnectedTo( iBndSrcN ); - const BRepMesh_PairOfIndex & triaIds = _triaDS->ElementsConnectedTo( linkIds.First() ); - const BRepMesh_Triangle& tria = _triaDS->GetElement( triaIds.Index(1) ); - addCloseNodes( _bndSrcNodes[ iBndSrcN ], &tria, srcFaceID, noTriQueue ); + const BRepMesh_Triangle* tria = GetTriangleNear( iBndSrcN ); + AddCloseNodes( _bndSrcNodes[ iBndSrcN ], tria, srcFaceID, noTriQueue ); } if ( noTriQueue.empty() ) { @@ -3099,8 +3105,8 @@ namespace StdMeshers_ProjectionUtils } // Morph::Perform -gp_XY Morph::GetBndUV(const int iNode) const -{ + gp_XY Morph::GetBndUV(const int iNode) const + { return _triaDS->GetNode( iNode ).Coord(); } diff --git a/src/StdMeshers/StdMeshers_ProjectionUtils.hxx b/src/StdMeshers/StdMeshers_ProjectionUtils.hxx index d44e477f0..6c94e45d0 100644 --- a/src/StdMeshers/StdMeshers_ProjectionUtils.hxx +++ b/src/StdMeshers/StdMeshers_ProjectionUtils.hxx @@ -159,11 +159,34 @@ namespace StdMeshers_ProjectionUtils const TNodeNodeMap& src2tgtNodes, const bool moveAll); - // return source boundary nodes. 0-th node is zero + // find a triangle containing an UV starting from a given triangle; + // return barycentric coordinates of the UV in the found triangle + const BRepMesh_Triangle* FindTriangle( const gp_XY& uv, + const BRepMesh_Triangle* bmTria, + double bc[3], + int triaNodes[3]); + + // return any delauney triangle neighboring a given boundary node + const BRepMesh_Triangle* GetTriangleNear( int iBndNode ); + + // return source boundary nodes. 0-th node is NULL so that indices of + // boundary nodes correspond to indices used by Delauney mesh const std::vector< const SMDS_MeshNode* >& GetBndNodes() const { return _bndSrcNodes; } // return UV of the i-th source boundary node gp_XY GetBndUV(const int iNode) const; + + // return scale factor to convert real UV to/from UV used for Delauney meshing: + // delauney_UV = real_UV * scale + const gp_XY& GetScale() const { return _scale; } + + typedef std::list< std::pair< const SMDS_MeshNode*, const BRepMesh_Triangle* > > TNodeTriaList; + + // add non-marked nodes surrounding a given one to a list + static void AddCloseNodes( const SMDS_MeshNode* node, + const BRepMesh_Triangle* bmTria, + const int faceID, + TNodeTriaList & noTriQueue ); }; /*! diff --git a/src/StdMeshers/StdMeshers_Projection_2D.cxx b/src/StdMeshers/StdMeshers_Projection_2D.cxx index 4063153dc..dc418c4f3 100644 --- a/src/StdMeshers/StdMeshers_Projection_2D.cxx +++ b/src/StdMeshers/StdMeshers_Projection_2D.cxx @@ -1133,7 +1133,7 @@ namespace { { SMESH_subMesh* faceSM = helper.GetMesh()->GetSubMesh( helper.GetSubShape() ); - if ( helper.IsDistorted2D( faceSM, /*checkUV=*/true )) + //if ( helper.IsDistorted2D( faceSM, /*checkUV=*/true )) { SMESH_MeshEditor editor( helper.GetMesh() ); SMESHDS_SubMesh* smDS = faceSM->GetSubMeshDS();