From b385a0679d4079d84272d9fe9006c2736beeb907 Mon Sep 17 00:00:00 2001 From: eap Date: Wed, 24 Sep 2014 17:38:27 +0400 Subject: [PATCH] 52499: Prismatic mesh is not computed on a prismatic shape 22701: EDF SMESH: Crash when creating mesh 1) Search of transformation using a least-square approximation 2) Compute internal prism nodes using the found transformation --- src/SMESHUtils/SMESH_TypeDefs.hxx | 3 +- src/StdMeshers/StdMeshers_Prism_3D.cxx | 693 ++++++++++++++---- src/StdMeshers/StdMeshers_Prism_3D.hxx | 49 +- src/StdMeshers/StdMeshers_ProjectionUtils.cxx | 229 +++++- src/StdMeshers/StdMeshers_ProjectionUtils.hxx | 55 +- src/StdMeshers/StdMeshers_Projection_2D.cxx | 547 ++++++++------ src/StdMeshers/StdMeshers_Projection_2D.hxx | 8 +- src/StdMeshers/StdMeshers_Projection_3D.cxx | 2 +- src/StdMeshers/StdMeshers_RadialPrism_3D.cxx | 2 +- src/StdMeshers/StdMeshers_ViscousLayers.cxx | 2 +- 10 files changed, 1232 insertions(+), 358 deletions(-) diff --git a/src/SMESHUtils/SMESH_TypeDefs.hxx b/src/SMESHUtils/SMESH_TypeDefs.hxx index 9fec2e5bd..f5dedaec0 100644 --- a/src/SMESHUtils/SMESH_TypeDefs.hxx +++ b/src/SMESHUtils/SMESH_TypeDefs.hxx @@ -43,7 +43,8 @@ typedef std::map, TIDCompare > TElemOfElemListMap; typedef std::map, TIDCompare > TElemOfNodeListMap; -typedef std::map TNodeNodeMap; +typedef std::map TNodeNodeMap; //!< Set of elements sorted by ID, to be used to assure predictability of edition typedef std::set< const SMDS_MeshElement*, TIDCompare > TIDSortedElemSet; diff --git a/src/StdMeshers/StdMeshers_Prism_3D.cxx b/src/StdMeshers/StdMeshers_Prism_3D.cxx index afac9bebe..eb43e8895 100644 --- a/src/StdMeshers/StdMeshers_Prism_3D.cxx +++ b/src/StdMeshers/StdMeshers_Prism_3D.cxx @@ -77,7 +77,7 @@ using namespace std; #define SHOWYXZ(msg, xyz) #endif -namespace TAssocTool = StdMeshers_ProjectionUtils; +namespace NSProjUtils = StdMeshers_ProjectionUtils; typedef SMESH_Comment TCom; @@ -157,6 +157,10 @@ namespace { fatherAlgo->GetGen() ); return algo; } + const NSProjUtils::TNodeNodeMap& GetNodesMap() + { + return _src2tgtNodes; + } }; //======================================================================= /*! @@ -539,7 +543,7 @@ StdMeshers_Prism_3D::StdMeshers_Prism_3D(int hypId, int studyId, SMESH_Gen* gen) { _name = "Prism_3D"; _shapeType = (1 << TopAbs_SOLID); // 1 bit per shape type - _onlyUnaryInput = false; // accept all SOLIDs at once + _onlyUnaryInput = false; // mesh all SOLIDs at once _requireDiscreteBoundary = false; // mesh FACEs and EDGEs by myself _supportSubmeshes = true; // "source" FACE must be meshed by other algo _neededLowerHyps[ 1 ] = true; // suppress warning on hiding a global 1D algo @@ -578,12 +582,12 @@ bool StdMeshers_Prism_3D::CheckHypothesis(SMESH_Mesh& a for ( ; exp.More(); exp.Next() ) { ++nbFace; const TopoDS_Shape& face = exp.Current(); - nbEdge = TAssocTool::Count( face, TopAbs_EDGE, 0 ); - nbWire = TAssocTool::Count( face, TopAbs_WIRE, 0 ); + nbEdge = NSProjUtils::Count( face, TopAbs_EDGE, 0 ); + nbWire = NSProjUtils::Count( face, TopAbs_WIRE, 0 ); if ( nbEdge!= 4 || nbWire!= 1 ) { if ( !notQuadFaces.empty() ) { - if ( TAssocTool::Count( notQuadFaces.back(), TopAbs_EDGE, 0 ) != nbEdge || - TAssocTool::Count( notQuadFaces.back(), TopAbs_WIRE, 0 ) != nbWire ) + if ( NSProjUtils::Count( notQuadFaces.back(), TopAbs_EDGE, 0 ) != nbEdge || + NSProjUtils::Count( notQuadFaces.back(), TopAbs_WIRE, 0 ) != nbWire ) RETURN_BAD_RESULT("Different not quad faces"); } notQuadFaces.push_back( face ); @@ -595,7 +599,7 @@ bool StdMeshers_Prism_3D::CheckHypothesis(SMESH_Mesh& a RETURN_BAD_RESULT("Bad nb not quad faces: " << notQuadFaces.size()); // check total nb faces - nbEdge = TAssocTool::Count( notQuadFaces.back(), TopAbs_EDGE, 0 ); + nbEdge = NSProjUtils::Count( notQuadFaces.back(), TopAbs_EDGE, 0 ); if ( nbFace != nbEdge + 2 ) RETURN_BAD_RESULT("Bad nb of faces: " << nbFace << " but must be " << nbEdge + 2); } @@ -1057,6 +1061,8 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& 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) if ( !myBlock.Init( myHelper, thePrism )) return toSM( error( myBlock.GetError())); @@ -1067,10 +1073,10 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) // Try to get gp_Trsf to get all nodes from bottom ones vector trsf; gp_Trsf bottomToTopTrsf; - if ( !myBlock.GetLayersTransformation( trsf, thePrism )) - trsf.clear(); - else if ( !trsf.empty() ) - bottomToTopTrsf = trsf.back(); + // if ( !myBlock.GetLayersTransformation( trsf, thePrism )) + // trsf.clear(); + // else if ( !trsf.empty() ) + // bottomToTopTrsf = trsf.back(); // To compute coordinates of a node inside a block, it is necessary to know // 1. normalized parameters of the node by which @@ -1092,31 +1098,30 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) // Create nodes inside the block - // try to use transformation (issue 0020680) - if ( !trsf.empty() ) + // use transformation (issue 0020680, IPAL0052499) + StdMeshers_Sweeper sweeper; + + // load boundary nodes + bool dummy; + list< TopoDS_Edge >::const_iterator edge = thePrism.myBottomEdges.begin(); + for ( ; edge != thePrism.myBottomEdges.end(); ++edge ) { - // loop on nodes inside the bottom face - TNode2ColumnMap::iterator bot_column = myBotToColumnMap.begin(); - for ( ; bot_column != myBotToColumnMap.end(); ++bot_column ) - { - const Prism_3D::TNode& tBotNode = bot_column->first; // bottom TNode - if ( tBotNode.GetPositionType() != SMDS_TOP_FACE ) - continue; // node is not inside face + 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 ); - // column nodes; middle part of the column are zero pointers - TNodeColumn& column = bot_column->second; - TNodeColumn::iterator columnNodes = column.begin(); - for ( int z = 0; columnNodes != column.end(); ++columnNodes, ++z) - { - const SMDS_MeshNode* & node = *columnNodes; - if ( node ) continue; // skip bottom or top node + const double tol = getSweepTolerance( thePrism ); - gp_XYZ coords = tBotNode.GetCoords(); - trsf[z-1].Transforms( coords ); - node = meshDS->AddNode( coords.X(), coords.Y(), coords.Z() ); - meshDS->SetNodeInVolume( node, volumeID ); - } - } // loop on bottom nodes + if ( sweeper.ComputeNodes( *myHelper, tol )) + { } else // use block approach { @@ -1127,7 +1132,7 @@ bool StdMeshers_Prism_3D::compute(const Prism_3D::TPrismTopo& thePrism) { const Prism_3D::TNode& tBotNode = bot_column->first; // bottom TNode if ( tBotNode.GetPositionType() != SMDS_TOP_FACE ) - continue; // node is not inside the FACE + continue; // node is not inside the FACE // column nodes; middle part of the column are zero pointers TNodeColumn& column = bot_column->second; @@ -1853,8 +1858,8 @@ void StdMeshers_Prism_3D::AddPrisms( vector & columns, bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf, const Prism_3D::TPrismTopo& thePrism) { - SMESH_subMesh * botSM = myBlock.SubMesh( ID_BOT_FACE ); - SMESH_subMesh * topSM = myBlock.SubMesh( ID_TOP_FACE ); + SMESH_subMesh * botSM = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom ); + SMESH_subMesh * topSM = myHelper->GetMesh()->GetSubMesh( thePrism.myTop ); SMESHDS_SubMesh * botSMDS = botSM->GetSubMeshDS(); SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS(); @@ -1885,86 +1890,90 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf <<" and #"<< topSM->GetId() << " seems different" )); ///RETURN_BAD_RESULT("Need to project but not allowed"); + NSProjUtils::TNodeNodeMap n2nMap; + const NSProjUtils::TNodeNodeMap* n2nMapPtr = & n2nMap; if ( needProject ) { - return projectBottomToTop( bottomToTopTrsf ); + if ( !projectBottomToTop( bottomToTopTrsf, thePrism )) + return false; + n2nMapPtr = & TProjction2dAlgo::instance( this )->GetNodesMap(); } - TopoDS_Face botFace = TopoDS::Face( myBlock.Shape( ID_BOT_FACE )); - TopoDS_Face topFace = TopoDS::Face( myBlock.Shape( ID_TOP_FACE )); - // associate top and bottom faces - TAssocTool::TShapeShapeMap shape2ShapeMap; - const bool sameTopo = - TAssocTool::FindSubShapeAssociation( botFace, myBlock.Mesh(), - topFace, myBlock.Mesh(), - shape2ShapeMap); - if ( !sameTopo ) - for ( size_t iQ = 0; iQ < thePrism.myWallQuads.size(); ++iQ ) - { - const Prism_3D::TQuadList& quadList = thePrism.myWallQuads[iQ]; - StdMeshers_FaceSidePtr botSide = quadList.front()->side[ QUAD_BOTTOM_SIDE ]; - StdMeshers_FaceSidePtr topSide = quadList.back ()->side[ QUAD_TOP_SIDE ]; - if ( botSide->NbEdges() == topSide->NbEdges() ) - { - for ( int iE = 0; iE < botSide->NbEdges(); ++iE ) - { - TAssocTool::InsertAssociation( botSide->Edge( iE ), - topSide->Edge( iE ), shape2ShapeMap ); - TAssocTool::InsertAssociation( myHelper->IthVertex( 0, botSide->Edge( iE )), - myHelper->IthVertex( 0, topSide->Edge( iE )), - shape2ShapeMap ); - } - } - else + if ( !n2nMapPtr || 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); + if ( !sameTopo ) + for ( size_t iQ = 0; iQ < thePrism.myWallQuads.size(); ++iQ ) { - TopoDS_Vertex vb, vt; - StdMeshers_FaceSidePtr sideB, sideT; - vb = myHelper->IthVertex( 0, botSide->Edge( 0 )); - vt = myHelper->IthVertex( 0, topSide->Edge( 0 )); - sideB = quadList.front()->side[ QUAD_LEFT_SIDE ]; - sideT = quadList.back ()->side[ QUAD_LEFT_SIDE ]; - if ( vb.IsSame( sideB->FirstVertex() ) && - vt.IsSame( sideT->LastVertex() )) + const Prism_3D::TQuadList& quadList = thePrism.myWallQuads[iQ]; + StdMeshers_FaceSidePtr botSide = quadList.front()->side[ QUAD_BOTTOM_SIDE ]; + StdMeshers_FaceSidePtr topSide = quadList.back ()->side[ QUAD_TOP_SIDE ]; + if ( botSide->NbEdges() == topSide->NbEdges() ) { - TAssocTool::InsertAssociation( botSide->Edge( 0 ), - topSide->Edge( 0 ), shape2ShapeMap ); - TAssocTool::InsertAssociation( vb, vt, shape2ShapeMap ); + for ( int iE = 0; iE < botSide->NbEdges(); ++iE ) + { + NSProjUtils::InsertAssociation( botSide->Edge( iE ), + topSide->Edge( iE ), shape2ShapeMap ); + NSProjUtils::InsertAssociation( myHelper->IthVertex( 0, botSide->Edge( iE )), + myHelper->IthVertex( 0, topSide->Edge( iE )), + shape2ShapeMap ); + } } - vb = myHelper->IthVertex( 1, botSide->Edge( botSide->NbEdges()-1 )); - vt = myHelper->IthVertex( 1, topSide->Edge( topSide->NbEdges()-1 )); - sideB = quadList.front()->side[ QUAD_RIGHT_SIDE ]; - sideT = quadList.back ()->side[ QUAD_RIGHT_SIDE ]; - if ( vb.IsSame( sideB->FirstVertex() ) && - vt.IsSame( sideT->LastVertex() )) + else { - TAssocTool::InsertAssociation( botSide->Edge( botSide->NbEdges()-1 ), - topSide->Edge( topSide->NbEdges()-1 ), - shape2ShapeMap ); - TAssocTool::InsertAssociation( vb, vt, shape2ShapeMap ); + TopoDS_Vertex vb, vt; + StdMeshers_FaceSidePtr sideB, sideT; + vb = myHelper->IthVertex( 0, botSide->Edge( 0 )); + vt = myHelper->IthVertex( 0, topSide->Edge( 0 )); + sideB = quadList.front()->side[ QUAD_LEFT_SIDE ]; + sideT = quadList.back ()->side[ QUAD_LEFT_SIDE ]; + if ( vb.IsSame( sideB->FirstVertex() ) && + vt.IsSame( sideT->LastVertex() )) + { + NSProjUtils::InsertAssociation( botSide->Edge( 0 ), + topSide->Edge( 0 ), shape2ShapeMap ); + NSProjUtils::InsertAssociation( vb, vt, shape2ShapeMap ); + } + vb = myHelper->IthVertex( 1, botSide->Edge( botSide->NbEdges()-1 )); + vt = myHelper->IthVertex( 1, topSide->Edge( topSide->NbEdges()-1 )); + sideB = quadList.front()->side[ QUAD_RIGHT_SIDE ]; + sideT = quadList.back ()->side[ QUAD_RIGHT_SIDE ]; + if ( vb.IsSame( sideB->FirstVertex() ) && + vt.IsSame( sideT->LastVertex() )) + { + NSProjUtils::InsertAssociation( botSide->Edge( botSide->NbEdges()-1 ), + topSide->Edge( topSide->NbEdges()-1 ), + shape2ShapeMap ); + NSProjUtils::InsertAssociation( vb, vt, shape2ShapeMap ); + } } } - } - // Find matching nodes of top and bottom faces - TNodeNodeMap n2nMap; - if ( ! TAssocTool::FindMatchingNodesOnFaces( botFace, myBlock.Mesh(), - topFace, myBlock.Mesh(), - shape2ShapeMap, n2nMap )) - { - if ( sameTopo ) - return toSM( error(TCom("Mesh on faces #") << botSM->GetId() - <<" and #"<< topSM->GetId() << " seems different" )); - else - return toSM( error(TCom("Topology of faces #") << botSM->GetId() - <<" and #"<< topSM->GetId() << " seems different" )); + // Find matching nodes of top and bottom faces + n2nMapPtr = & n2nMap; + if ( ! NSProjUtils::FindMatchingNodesOnFaces( thePrism.myBottom, myHelper->GetMesh(), + thePrism.myTop, myHelper->GetMesh(), + shape2ShapeMap, n2nMap )) + { + if ( sameTopo ) + return toSM( error(TCom("Mesh on faces #") << botSM->GetId() + <<" and #"<< topSM->GetId() << " seems different" )); + else + return toSM( error(TCom("Topology of faces #") << botSM->GetId() + <<" and #"<< topSM->GetId() << " seems different" )); + } } // Fill myBotToColumnMap int zSize = myBlock.VerticalSize(); - //TNode prevTNode; - TNodeNodeMap::iterator bN_tN = n2nMap.begin(); - for ( ; bN_tN != n2nMap.end(); ++bN_tN ) + TNodeNodeMap::const_iterator bN_tN = n2nMapPtr->begin(); + for ( ; bN_tN != n2nMapPtr->end(); ++bN_tN ) { const SMDS_MeshNode* botNode = bN_tN->first; const SMDS_MeshNode* topNode = bN_tN->second; @@ -1984,17 +1993,22 @@ bool StdMeshers_Prism_3D::assocOrProjBottom2Top( const gp_Trsf & bottomToTopTrsf //================================================================================ /*! - * \brief Remove quadrangles from the top face and - * create triangles there by projection from the bottom + * \brief Remove faces from the top face and re-create them by projection from the bottom * \retval bool - a success or not */ //================================================================================ -bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottomToTopTrsf ) +bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottomToTopTrsf, + const Prism_3D::TPrismTopo& thePrism ) { - SMESHDS_Mesh* meshDS = myBlock.MeshDS(); - SMESH_subMesh * botSM = myBlock.SubMesh( ID_BOT_FACE ); - SMESH_subMesh * topSM = myBlock.SubMesh( ID_TOP_FACE ); + if ( project2dMesh( thePrism.myBottom, thePrism.myTop )) + { + return true; + } + + SMESHDS_Mesh* meshDS = myHelper->GetMeshDS(); + SMESH_subMesh * botSM = myHelper->GetMesh()->GetSubMesh( thePrism.myBottom ); + SMESH_subMesh * topSM = myHelper->GetMesh()->GetSubMesh( thePrism.myTop ); SMESHDS_SubMesh * botSMDS = botSM->GetSubMeshDS(); SMESHDS_SubMesh * topSMDS = topSM->GetSubMeshDS(); @@ -2002,9 +2016,9 @@ bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottomToTopTrsf ) if ( topSMDS && topSMDS->NbElements() > 0 ) topSM->ComputeStateEngine( SMESH_subMesh::CLEAN ); - const TopoDS_Face& botFace = TopoDS::Face( myBlock.Shape( ID_BOT_FACE )); // oriented within - const TopoDS_Face& topFace = TopoDS::Face( myBlock.Shape( ID_TOP_FACE )); // the 3D SHAPE - int topFaceID = meshDS->ShapeToIndex( topFace ); + const TopoDS_Face& botFace = thePrism.myBottom; // oriented within + const TopoDS_Face& topFace = thePrism.myTop; // the 3D SHAPE + int topFaceID = meshDS->ShapeToIndex( thePrism.myTop ); SMESH_MesherHelper botHelper( *myHelper->GetMesh() ); botHelper.SetSubShape( botFace ); @@ -2071,6 +2085,9 @@ bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottomToTopTrsf ) column.resize( zSize ); column.front() = botNode; column.back() = topNode; + + if ( _computeCanceled ) + return toSM( error( SMESH_ComputeError::New(COMPERR_CANCELED))); } // Create top faces @@ -2135,6 +2152,77 @@ bool StdMeshers_Prism_3D::projectBottomToTop( const gp_Trsf & bottomToTopTrsf ) return true; } +//======================================================================= +//function : getSweepTolerance +//purpose : Compute tolerance to pass to StdMeshers_Sweeper +//======================================================================= + +double StdMeshers_Prism_3D::getSweepTolerance( const Prism_3D::TPrismTopo& thePrism ) +{ + SMESHDS_Mesh* meshDS = myHelper->GetMeshDS(); + SMESHDS_SubMesh * sm[2] = { meshDS->MeshElements( thePrism.myBottom ), + meshDS->MeshElements( thePrism.myTop ) }; + double minDist = 1e100; + + vector< SMESH_TNodeXYZ > nodes; + for ( int iSM = 0; iSM < 2; ++iSM ) + { + if ( !sm[ iSM ]) continue; + + SMDS_ElemIteratorPtr fIt = sm[ iSM ]->GetElements(); + while ( fIt->more() ) + { + const SMDS_MeshElement* face = fIt->next(); + const int nbNodes = face->NbCornerNodes(); + SMDS_ElemIteratorPtr nIt = face->nodesIterator(); + + nodes.resize( nbNodes + 1 ); + for ( int iN = 0; iN < nbNodes; ++iN ) + nodes[ iN ] = nIt->next(); + nodes.back() = nodes[0]; + + // loop on links + double dist2; + for ( int iN = 0; iN < nbNodes; ++iN ) + { + if ( nodes[ iN ]._node->GetPosition()->GetDim() < 2 && + nodes[ iN+1 ]._node->GetPosition()->GetDim() < 2 ) + { + // it's a boundary link; measure distance of other + // nodes to this link + gp_XYZ linkDir = nodes[ iN ] - nodes[ iN+1 ]; + double linkLen = linkDir.Modulus(); + bool isDegen = ( linkLen < numeric_limits::min() ); + if ( !isDegen ) linkDir /= linkLen; + for ( int iN2 = 0; iN2 < nbNodes; ++iN2 ) // loop on other nodes + { + if ( nodes[ iN2 ] == nodes[ iN ] || + nodes[ iN2 ] == nodes[ iN+1 ]) continue; + if ( isDegen ) + { + dist2 = ( nodes[ iN ] - nodes[ iN2 ]).SquareModulus(); + } + else + { + dist2 = linkDir.CrossSquareMagnitude( nodes[ iN ] - nodes[ iN2 ]); + } + if ( dist2 > numeric_limits::min() ) + minDist = Min ( minDist, dist2 ); + } + } + // measure length link + else if ( nodes[ iN ]._node < nodes[ iN+1 ]._node ) // not to measure same link twice + { + dist2 = ( nodes[ iN ] - nodes[ iN+1 ]).SquareModulus(); + if ( dist2 > numeric_limits::min() ) + minDist = Min ( minDist, dist2 ); + } + } + } + } + return 0.1 * Sqrt ( minDist ); +} + //======================================================================= //function : project2dMesh //purpose : Project mesh faces from a source FACE of one prism (theSrcFace) @@ -2261,6 +2349,12 @@ namespace // utils used by StdMeshers_Prism_3D::IsApplicable() if ( E.IsSame( Edge( i ))) return i; return -1; } + bool IsSideFace( const TopoDS_Shape& face ) const + { + if ( _faces->Contains( face )) // avoid returning true for a prism top FACE + return ( !_face.IsNull() || !( face.IsSame( _faces->FindKey( _faces->Extent() )))); + return false; + } }; //-------------------------------------------------------------------------------- /*! @@ -2392,12 +2486,11 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA TEdgeWithNeighborsVec& botEdges = faceEdgesVec[ iF ]; if ( botEdges.empty() ) - { if ( !getEdges( botF, botEdges, /*noHoles=*/false )) break; - if ( allFaces.Extent()-1 <= (int) botEdges.size() ) - continue; // all faces are adjacent to botF - no top FACE - } + if ( allFaces.Extent()-1 <= (int) botEdges.size() ) + 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 ) @@ -2411,8 +2504,8 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA } bool isOK = true; // ok for a current botF - bool isAdvanced = true; - int nbFoundSideFaces = 0; + bool isAdvanced = true; // is new data found in a current loop + int nbFoundSideFaces = 0; for ( int iLoop = 0; isOK && isAdvanced; ++iLoop ) { isAdvanced = false; @@ -2420,7 +2513,8 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA { PrismSide& side = sides[ iS ]; if ( side._face.IsNull() ) - continue; + continue; // probably the prism top face is the last of side._faces + if ( side._topEdge.IsNull() ) { // find vertical EDGEs --- EGDEs shared with neighbor side FACEs @@ -2434,7 +2528,7 @@ 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->_faces->Contains( neighborF ); + bool isEdgeShared = adjSide->IsSideFace( neighborF ); if ( isEdgeShared ) { isAdvanced = true; @@ -2480,13 +2574,13 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA { if ( side._leftSide->_faces->Contains( f )) { - stop = true; + stop = true; // probably f is the prism top face side._leftSide->_face.Nullify(); side._leftSide->_topEdge.Nullify(); } if ( side._rightSide->_faces->Contains( f )) { - stop = true; + stop = true; // probably f is the prism top face side._rightSide->_face.Nullify(); side._rightSide->_topEdge.Nullify(); } @@ -2497,8 +2591,8 @@ bool StdMeshers_Prism_3D::IsApplicable(const TopoDS_Shape & shape, bool toCheckA side._topEdge.Nullify(); continue; } - side._face = TopoDS::Face( f ); - int faceID = allFaces.FindIndex( side._face ); + side._face = TopoDS::Face( f ); + int faceID = allFaces.FindIndex( side._face ); side._edges = & faceEdgesVec[ faceID ]; if ( side._edges->empty() ) if ( !getEdges( side._face, * side._edges, /*noHoles=*/true )) @@ -2805,8 +2899,9 @@ bool StdMeshers_Prism_3D::initPrism(Prism_3D::TPrismTopo& thePrism, thePrism.myShape3D = shape3D; if ( thePrism.myBottom.IsNull() ) thePrism.myBottom = TopoDS::Face( botSM->GetSubShape() ); - thePrism.myBottom.Orientation( myHelper->GetSubShapeOri( shape3D, - thePrism.myBottom )); + thePrism.myBottom.Orientation( myHelper->GetSubShapeOri( shape3D, thePrism.myBottom )); + thePrism.myTop. Orientation( myHelper->GetSubShapeOri( shape3D, thePrism.myTop )); + // Get ordered bottom edges TopoDS_Face reverseBottom = // to have order of top EDGEs as in the top FACE TopoDS::Face( thePrism.myBottom.Reversed() ); @@ -4221,3 +4316,355 @@ gp_Pnt2d StdMeshers_PrismAsBlock::TPCurveOnHorFaceAdaptor::Value(const Standard_ double r = ( U - i1->first ) / ( i2->first - i1->first ); return i1->second * ( 1 - r ) + i2->second * r; } + +//================================================================================ +/*! + * \brief Projects internal nodes using transformation found by boundary nodes + */ +//================================================================================ + +bool StdMeshers_Sweeper::projectIntPoints(const vector< gp_XYZ >& fromBndPoints, + const vector< gp_XYZ >& toBndPoints, + const vector< gp_XYZ >& fromIntPoints, + vector< gp_XYZ >& toIntPoints, + NSProjUtils::TrsfFinder3D& trsf, + vector< gp_XYZ > * bndError) +{ + // find transformation + if ( trsf.IsIdentity() && !trsf.Solve( fromBndPoints, toBndPoints )) + return false; + + // compute internal points using the found trsf + for ( size_t iP = 0; iP < fromIntPoints.size(); ++iP ) + { + toIntPoints[ iP ] = trsf.Transform( fromIntPoints[ iP ]); + } + + // compute boundary error + if ( bndError ) + { + bndError->resize( fromBndPoints.size() ); + gp_XYZ fromTrsf; + for ( size_t iP = 0; iP < fromBndPoints.size(); ++iP ) + { + fromTrsf = trsf.Transform( fromBndPoints[ iP ] ); + (*bndError)[ iP ] = toBndPoints[ iP ] - fromTrsf; + } + } + return true; +} + +//================================================================================ +/*! + * \brief Add boundary error to ineternal points + */ +//================================================================================ + +void StdMeshers_Sweeper::applyBoundaryError(const vector< gp_XYZ >& bndPoints, + const vector< gp_XYZ >& bndError1, + const vector< gp_XYZ >& bndError2, + const double r, + vector< gp_XYZ >& intPoints, + vector< double >& int2BndDist) +{ + // fix each internal point + const double eps = 1e-100; + for ( size_t iP = 0; iP < intPoints.size(); ++iP ) + { + gp_XYZ & intPnt = intPoints[ iP ]; + + // compute distance from intPnt to each boundary node + double int2BndDistSum = 0; + for ( size_t iBnd = 0; iBnd < bndPoints.size(); ++iBnd ) + { + int2BndDist[ iBnd ] = 1 / (( intPnt - bndPoints[ iBnd ]).SquareModulus() + eps ); + int2BndDistSum += int2BndDist[ iBnd ]; + } + + // apply bndError + for ( size_t iBnd = 0; iBnd < bndPoints.size(); ++iBnd ) + { + intPnt += bndError1[ iBnd ] * ( 1 - r ) * int2BndDist[ iBnd ] / int2BndDistSum; + intPnt += bndError2[ iBnd ] * r * int2BndDist[ iBnd ] / int2BndDistSum; + } + } +} + +//================================================================================ +/*! + * \brief Creates internal nodes of the prism + */ +//================================================================================ + +bool StdMeshers_Sweeper::ComputeNodes( SMESH_MesherHelper& helper, + const double tol) +{ + const size_t zSize = myBndColumns[0]->size(); + const size_t zSrc = 0, zTgt = zSize-1; + if ( zSize < 3 ) return true; + + vector< vector< gp_XYZ > > intPntsOfLayer( zSize ); // node coodinates to compute + // set coordinates of src and tgt nodes + for ( size_t z = 0; z < intPntsOfLayer.size(); ++z ) + intPntsOfLayer[ z ].resize( myIntColumns.size() ); + for ( size_t iP = 0; iP < myIntColumns.size(); ++iP ) + { + intPntsOfLayer[ zSrc ][ iP ] = intPoint( iP, zSrc ); + intPntsOfLayer[ zTgt ][ iP ] = intPoint( iP, zTgt ); + } + + // compute coordinates of internal nodes by projecting (transfroming) src and tgt + // nodes towards the central layer + + vector< NSProjUtils::TrsfFinder3D > trsfOfLayer( zSize ); + vector< vector< gp_XYZ > > bndError( zSize ); + + // boundary points used to compute an affine transformation from a layer to a next one + vector< gp_XYZ > fromSrcBndPnts( myBndColumns.size() ), fromTgtBndPnts( myBndColumns.size() ); + vector< gp_XYZ > toSrcBndPnts ( myBndColumns.size() ), toTgtBndPnts ( myBndColumns.size() ); + for ( size_t iP = 0; iP < myBndColumns.size(); ++iP ) + { + fromSrcBndPnts[ iP ] = bndPoint( iP, zSrc ); + fromTgtBndPnts[ iP ] = bndPoint( iP, zTgt ); + } + + size_t zS = zSrc + 1; + size_t zT = zTgt - 1; + for ( ; zS < zT; ++zS, --zT ) // vertical loop on layers + { + for ( size_t iP = 0; iP < myBndColumns.size(); ++iP ) + { + toSrcBndPnts[ iP ] = bndPoint( iP, zS ); + toTgtBndPnts[ iP ] = bndPoint( iP, zT ); + } + if (! projectIntPoints( fromSrcBndPnts, toSrcBndPnts, + intPntsOfLayer[ zS-1 ], intPntsOfLayer[ zS ], + trsfOfLayer [ zS-1 ], & bndError[ zS-1 ])) + return false; + if (! projectIntPoints( fromTgtBndPnts, toTgtBndPnts, + intPntsOfLayer[ zT+1 ], intPntsOfLayer[ zT ], + trsfOfLayer [ zT+1 ], & bndError[ zT+1 ])) + return false; + + // if ( zT == zTgt - 1 ) + // { + // for ( size_t iP = 0; iP < myBndColumns.size(); ++iP ) + // { + // gp_XYZ fromTrsf = trsfOfLayer [ zT+1].Transform( fromTgtBndPnts[ iP ] ); + // cout << "mesh.AddNode( " + // << fromTrsf.X() << ", " + // << fromTrsf.Y() << ", " + // << fromTrsf.Z() << ") " << endl; + // } + // for ( size_t iP = 0; iP < myIntColumns.size(); ++iP ) + // cout << "mesh.AddNode( " + // << intPntsOfLayer[ zT ][ iP ].X() << ", " + // << intPntsOfLayer[ zT ][ iP ].Y() << ", " + // << intPntsOfLayer[ zT ][ iP ].Z() << ") " << endl; + // } + + fromTgtBndPnts.swap( toTgtBndPnts ); + fromSrcBndPnts.swap( toSrcBndPnts ); + } + + // Compute two projections of internal points to the central layer + // in order to evaluate an error of internal points + + bool centerIntErrorIsSmall; + vector< gp_XYZ > centerSrcIntPnts( myIntColumns.size() ); + vector< gp_XYZ > centerTgtIntPnts( myIntColumns.size() ); + + for ( size_t iP = 0; iP < myBndColumns.size(); ++iP ) + { + toSrcBndPnts[ iP ] = bndPoint( iP, zS ); + toTgtBndPnts[ iP ] = bndPoint( iP, zT ); + } + if (! projectIntPoints( fromSrcBndPnts, toSrcBndPnts, + intPntsOfLayer[ zS-1 ], centerSrcIntPnts, + trsfOfLayer [ zS-1 ], & bndError[ zS-1 ])) + return false; + if (! projectIntPoints( fromTgtBndPnts, toTgtBndPnts, + intPntsOfLayer[ zT+1 ], centerTgtIntPnts, + trsfOfLayer [ zT+1 ], & bndError[ zT+1 ])) + return false; + + // evaluate an error of internal points on the central layer + centerIntErrorIsSmall = true; + if ( zS == zT ) // odd zSize + { + for ( size_t iP = 0; ( iP < myIntColumns.size() && centerIntErrorIsSmall ); ++iP ) + centerIntErrorIsSmall = + (centerSrcIntPnts[ iP ] - centerTgtIntPnts[ iP ]).SquareModulus() < tol*tol; + } + else // even zSize + { + for ( size_t iP = 0; ( iP < myIntColumns.size() && centerIntErrorIsSmall ); ++iP ) + centerIntErrorIsSmall = + (intPntsOfLayer[ zS-1 ][ iP ] - centerTgtIntPnts[ iP ]).SquareModulus() < tol*tol; + } + + // Evaluate an error of boundary points + + bool bndErrorIsSmall = true; + for ( size_t iP = 0; ( iP < myBndColumns.size() && bndErrorIsSmall ); ++iP ) + { + double sumError = 0; + for ( size_t z = 1; z < zS; ++z ) // loop on layers + sumError += ( bndError[ z-1 ][ iP ].Modulus() + + bndError[ zSize-z ][ iP ].Modulus() ); + + bndErrorIsSmall = ( sumError < tol ); + } + + // compute final points on the central layer + std::vector< double > int2BndDist( myBndColumns.size() ); // work array of applyBoundaryError() + double r = zS / ( zSize - 1.); + if ( zS == zT ) + { + for ( size_t iP = 0; iP < myIntColumns.size(); ++iP ) + { + intPntsOfLayer[ zS ][ iP ] = + ( 1 - r ) * centerSrcIntPnts[ iP ] + r * centerTgtIntPnts[ iP ]; + } + if ( !bndErrorIsSmall ) + { + applyBoundaryError( toSrcBndPnts, bndError[ zS-1 ], bndError[ zS+1 ], r, + intPntsOfLayer[ zS ], int2BndDist ); + } + } + else + { + for ( size_t iP = 0; iP < myIntColumns.size(); ++iP ) + { + intPntsOfLayer[ zS ][ iP ] = + r * intPntsOfLayer[ zS ][ iP ] + ( 1 - r ) * centerSrcIntPnts[ iP ]; + intPntsOfLayer[ zT ][ iP ] = + r * intPntsOfLayer[ zT ][ iP ] + ( 1 - r ) * centerTgtIntPnts[ iP ]; + } + if ( !bndErrorIsSmall ) + { + applyBoundaryError( toSrcBndPnts, bndError[ zS-1 ], bndError[ zS+1 ], r, + intPntsOfLayer[ zS ], int2BndDist ); + applyBoundaryError( toTgtBndPnts, bndError[ zT+1 ], bndError[ zT-1 ], r, + intPntsOfLayer[ zT ], int2BndDist ); + } + } + + //centerIntErrorIsSmall = true; + //bndErrorIsSmall = true; + if ( !centerIntErrorIsSmall ) + { + // Compensate the central error; continue adding projection + // by going from central layer to the source and target ones + + vector< gp_XYZ >& fromSrcIntPnts = centerSrcIntPnts; + vector< gp_XYZ >& fromTgtIntPnts = centerTgtIntPnts; + vector< gp_XYZ > toSrcIntPnts( myIntColumns.size() ); + vector< gp_XYZ > toTgtIntPnts( myIntColumns.size() ); + vector< gp_XYZ > srcBndError( myBndColumns.size() ); + vector< gp_XYZ > tgtBndError( myBndColumns.size() ); + + fromTgtBndPnts.swap( toTgtBndPnts ); + fromSrcBndPnts.swap( toSrcBndPnts ); + + for ( ++zS, --zT; zS < zTgt; ++zS, --zT ) // vertical loop on layers + { + // invert transformation + if ( !trsfOfLayer[ zS+1 ].Invert() ) + trsfOfLayer[ zS+1 ] = NSProjUtils::TrsfFinder3D(); // to recompute + if ( !trsfOfLayer[ zT-1 ].Invert() ) + trsfOfLayer[ zT-1 ] = NSProjUtils::TrsfFinder3D(); + + // project internal nodes and compute bnd error + for ( size_t iP = 0; iP < myBndColumns.size(); ++iP ) + { + toSrcBndPnts[ iP ] = bndPoint( iP, zS ); + toTgtBndPnts[ iP ] = bndPoint( iP, zT ); + } + projectIntPoints( fromSrcBndPnts, toSrcBndPnts, + fromSrcIntPnts, toSrcIntPnts, + trsfOfLayer[ zS+1 ], & srcBndError ); + projectIntPoints( fromTgtBndPnts, toTgtBndPnts, + fromTgtIntPnts, toTgtIntPnts, + trsfOfLayer[ zT-1 ], & tgtBndError ); + + // if ( zS == zTgt - 1 ) + // { + // cout << "mesh2 = smesh.Mesh()" << endl; + // for ( size_t iP = 0; iP < myBndColumns.size(); ++iP ) + // { + // gp_XYZ fromTrsf = trsfOfLayer [ zS+1].Transform( fromSrcBndPnts[ iP ] ); + // cout << "mesh2.AddNode( " + // << fromTrsf.X() << ", " + // << fromTrsf.Y() << ", " + // << fromTrsf.Z() << ") " << endl; + // } + // for ( size_t iP = 0; iP < myIntColumns.size(); ++iP ) + // cout << "mesh2.AddNode( " + // << toSrcIntPnts[ iP ].X() << ", " + // << toSrcIntPnts[ iP ].Y() << ", " + // << toSrcIntPnts[ iP ].Z() << ") " << endl; + // } + + // sum up 2 projections + r = zS / ( zSize - 1.); + vector< gp_XYZ >& zSIntPnts = intPntsOfLayer[ zS ]; + vector< gp_XYZ >& zTIntPnts = intPntsOfLayer[ zT ]; + for ( size_t iP = 0; iP < myIntColumns.size(); ++iP ) + { + zSIntPnts[ iP ] = r * zSIntPnts[ iP ] + ( 1 - r ) * toSrcIntPnts[ iP ]; + zTIntPnts[ iP ] = r * zTIntPnts[ iP ] + ( 1 - r ) * toTgtIntPnts[ iP ]; + } + + // compensate bnd error + if ( !bndErrorIsSmall ) + { + applyBoundaryError( toSrcBndPnts, srcBndError, bndError[ zS+1 ], r, + intPntsOfLayer[ zS ], int2BndDist ); + applyBoundaryError( toTgtBndPnts, tgtBndError, bndError[ zT-1 ], r, + intPntsOfLayer[ zT ], int2BndDist ); + } + + fromSrcBndPnts.swap( toSrcBndPnts ); + fromSrcIntPnts.swap( toSrcIntPnts ); + fromTgtBndPnts.swap( toTgtBndPnts ); + fromTgtIntPnts.swap( toTgtIntPnts ); + } + } // if ( !centerIntErrorIsSmall ) + + else if ( !bndErrorIsSmall ) + { + zS = zSrc + 1; + zT = zTgt - 1; + for ( ; zS < zT; ++zS, --zT ) // vertical loop on layers + { + for ( size_t iP = 0; iP < myBndColumns.size(); ++iP ) + { + toSrcBndPnts[ iP ] = bndPoint( iP, zS ); + toTgtBndPnts[ iP ] = bndPoint( iP, zT ); + } + // compensate bnd error + applyBoundaryError( toSrcBndPnts, bndError[ zS-1 ], bndError[ zS-1 ], 0.5, + intPntsOfLayer[ zS ], int2BndDist ); + applyBoundaryError( toTgtBndPnts, bndError[ zT+1 ], bndError[ zT+1 ], 0.5, + intPntsOfLayer[ zT ], int2BndDist ); + } + } + + // cout << "centerIntErrorIsSmall = " << centerIntErrorIsSmall<< endl; + // cout << "bndErrorIsSmall = " << bndErrorIsSmall<< endl; + + // Create nodes + for ( size_t iP = 0; iP < myIntColumns.size(); ++iP ) + { + vector< const SMDS_MeshNode* > & nodeCol = *myIntColumns[ iP ]; + for ( size_t z = zSrc + 1; z < zTgt; ++z ) // vertical loop on layers + { + const gp_XYZ & xyz = intPntsOfLayer[ z ][ iP ]; + if ( !( nodeCol[ z ] = helper.AddNode( xyz.X(), xyz.Y(), xyz.Z() ))) + return false; + } + } + + return true; +} diff --git a/src/StdMeshers/StdMeshers_Prism_3D.hxx b/src/StdMeshers/StdMeshers_Prism_3D.hxx index f2929c31a..6c7c52fe1 100644 --- a/src/StdMeshers/StdMeshers_Prism_3D.hxx +++ b/src/StdMeshers/StdMeshers_Prism_3D.hxx @@ -56,6 +56,10 @@ namespace Prism_3D struct TNode; struct TPrismTopo; } +namespace StdMeshers_ProjectionUtils +{ + class TrsfFinder3D; +} class SMESHDS_SubMesh; class TopoDS_Edge; @@ -398,7 +402,42 @@ private: }; // class StdMeshers_PrismAsBlock -// ============================================= +// =============================================== +/*! + * \brief Tool building internal nodes in a prism + */ +struct StdMeshers_Sweeper +{ + std::vector< TNodeColumn* > myBndColumns; // boundary nodes + std::vector< TNodeColumn* > myIntColumns; // internal nodes + + bool ComputeNodes( SMESH_MesherHelper& helper, + const double tol ); + +private: + + gp_XYZ bndPoint( int iP, int z ) const + { return SMESH_TNodeXYZ( (*myBndColumns[ iP ])[ z ]); } + + gp_XYZ intPoint( int iP, int z ) const + { return SMESH_TNodeXYZ( (*myIntColumns[ iP ])[ z ]); } + + static bool projectIntPoints(const std::vector< gp_XYZ >& fromBndPoints, + const std::vector< gp_XYZ >& toBndPoints, + const std::vector< gp_XYZ >& fromIntPoints, + std::vector< gp_XYZ >& toIntPoints, + StdMeshers_ProjectionUtils::TrsfFinder3D& trsf, + std::vector< gp_XYZ > * bndError); + + static void applyBoundaryError(const std::vector< gp_XYZ >& bndPoints, + const std::vector< gp_XYZ >& bndError1, + const std::vector< gp_XYZ >& bndError2, + const double r, + std::vector< gp_XYZ >& toIntPoints, + std::vector< double >& int2BndDist); +}; + +// =============================================== /*! * \brief Algo building prisms on a prism shape */ @@ -480,7 +519,13 @@ public: * create triangles there by projection from the bottom * \retval bool - a success or not */ - bool projectBottomToTop( const gp_Trsf & bottomToTopTrsf ); + bool projectBottomToTop( const gp_Trsf & bottomToTopTrsf, + const Prism_3D::TPrismTopo& thePrism ); + + /*! + * \brief Compute tolerance to pass to StdMeshers_Sweeper + */ + double getSweepTolerance( const Prism_3D::TPrismTopo& thePrism ); /*! * \brief Project mesh faces from a source FACE of one prism to diff --git a/src/StdMeshers/StdMeshers_ProjectionUtils.cxx b/src/StdMeshers/StdMeshers_ProjectionUtils.cxx index 92474f205..df96a9c53 100644 --- a/src/StdMeshers/StdMeshers_ProjectionUtils.cxx +++ b/src/StdMeshers/StdMeshers_ProjectionUtils.cxx @@ -66,6 +66,7 @@ #include #include #include +#include #include #include @@ -98,7 +99,7 @@ namespace HERE = StdMeshers_ProjectionUtils; namespace { - static SMESHDS_Mesh* theMeshDS[2] = { 0, 0 }; // used to debug only + static SMESHDS_Mesh* theMeshDS[2] = { 0, 0 }; // used for debug only long shapeIndex(const TopoDS_Shape& S) { if ( theMeshDS[0] && theMeshDS[1] ) @@ -2402,3 +2403,229 @@ void StdMeshers_ProjectionUtils::SetEventListener(SMESH_subMesh* subMesh, } } } + +namespace StdMeshers_ProjectionUtils +{ + + //================================================================================ + /*! + * \brief Computes transformation beween two sets of 2D points using + * a least square approximation + * + * See "Surface Mesh Projection For Hexahedral Mesh Generation By Sweeping" + * by X.Roca, J.Sarrate, A.Huerta. (2.2) + */ + //================================================================================ + + bool TrsfFinder2D::Solve( const vector< gp_XY >& srcPnts, + const vector< gp_XY >& tgtPnts ) + { + // find gravity centers + gp_XY srcGC( 0,0 ), tgtGC( 0,0 ); + for ( size_t i = 0; i < srcPnts.size(); ++i ) + { + srcGC += srcPnts[i]; + tgtGC += tgtPnts[i]; + } + srcGC /= srcPnts.size(); + tgtGC /= tgtPnts.size(); + + // find trsf + + math_Matrix mat (1,4,1,4, 0.); + math_Vector vec (1,4, 0.); + + // cout << "m1 = smesh.Mesh('src')" << endl + // << "m2 = smesh.Mesh('tgt')" << endl; + double xx = 0, xy = 0, yy = 0; + for ( size_t i = 0; i < srcPnts.size(); ++i ) + { + gp_XY srcUV = srcPnts[i] - srcGC; + gp_XY tgtUV = tgtPnts[i] - tgtGC; + xx += srcUV.X() * srcUV.X(); + yy += srcUV.Y() * srcUV.Y(); + xy += srcUV.X() * srcUV.Y(); + vec( 1 ) += srcUV.X() * tgtUV.X(); + vec( 2 ) += srcUV.Y() * tgtUV.X(); + vec( 3 ) += srcUV.X() * tgtUV.Y(); + vec( 4 ) += srcUV.Y() * tgtUV.Y(); + // cout << "m1.AddNode( " << srcUV.X() << ", " << srcUV.Y() << ", 0 )" << endl + // << "m2.AddNode( " << tgtUV.X() << ", " << tgtUV.Y() << ", 0 )" << endl; + } + mat( 1,1 ) = mat( 3,3 ) = xx; + mat( 2,2 ) = mat( 4,4 ) = yy; + mat( 1,2 ) = mat( 2,1 ) = mat( 3,4 ) = mat( 4,3 ) = xy; + + math_Gauss solver( mat ); + if ( !solver.IsDone() ) + return false; + solver.Solve( vec ); + if ( vec.Norm2() < gp::Resolution() ) + return false; + // cout << vec( 1 ) << "\t " << vec( 2 ) << endl + // << vec( 3 ) << "\t " << vec( 4 ) << endl; + + _trsf.SetTranslation( tgtGC ); + _srcOrig = srcGC; + + gp_Mat2d& M = const_cast< gp_Mat2d& >( _trsf.HVectorialPart()); + M( 1,1 ) = vec( 1 ); + M( 2,1 ) = vec( 2 ); + M( 1,2 ) = vec( 3 ); + M( 2,2 ) = vec( 4 ); + + return true; + } + + //================================================================================ + /*! + * \brief Transforms a 2D points using a found transformation + */ + //================================================================================ + + gp_XY TrsfFinder2D::Transform( const gp_Pnt2d& srcUV ) const + { + gp_XY uv = srcUV.XY() - _srcOrig ; + _trsf.Transforms( uv ); + return uv; + } + + //================================================================================ + /*! + * \brief Computes transformation beween two sets of 3D points using + * a least square approximation + * + * See "Surface Mesh Projection For Hexahedral Mesh Generation By Sweeping" + * by X.Roca, J.Sarrate, A.Huerta. (2.4) + */ + //================================================================================ + + bool TrsfFinder3D::Solve( const vector< gp_XYZ > & srcPnts, + const vector< gp_XYZ > & tgtPnts ) + { + // find gravity center + gp_XYZ srcGC( 0,0,0 ), tgtGC( 0,0,0 ); + for ( size_t i = 0; i < srcPnts.size(); ++i ) + { + srcGC += srcPnts[i]; + tgtGC += tgtPnts[i]; + } + srcGC /= srcPnts.size(); + tgtGC /= tgtPnts.size(); + + gp_XYZ srcOrig = 2 * srcGC - tgtGC; + gp_XYZ tgtOrig = srcGC; + + // find trsf + + math_Matrix mat (1,9,1,9, 0.); + math_Vector vec (1,9, 0.); + + double xx = 0, yy = 0, zz = 0; + double xy = 0, xz = 0, yz = 0; + for ( size_t i = 0; i < srcPnts.size(); ++i ) + { + gp_XYZ src = srcPnts[i] - srcOrig; + gp_XYZ tgt = tgtPnts[i] - tgtOrig; + xx += src.X() * src.X(); + yy += src.Y() * src.Y(); + zz += src.Z() * src.Z(); + xy += src.X() * src.Y(); + xz += src.X() * src.Z(); + yz += src.Y() * src.Z(); + vec( 1 ) += src.X() * tgt.X(); + vec( 2 ) += src.Y() * tgt.X(); + vec( 3 ) += src.Z() * tgt.X(); + vec( 4 ) += src.X() * tgt.Y(); + vec( 5 ) += src.Y() * tgt.Y(); + vec( 6 ) += src.Z() * tgt.Y(); + vec( 7 ) += src.X() * tgt.Z(); + vec( 8 ) += src.Y() * tgt.Z(); + vec( 9 ) += src.Z() * tgt.Z(); + } + mat( 1,1 ) = mat( 4,4 ) = mat( 7,7 ) = xx; + mat( 2,2 ) = mat( 5,5 ) = mat( 8,8 ) = yy; + mat( 3,3 ) = mat( 6,6 ) = mat( 9,9 ) = zz; + mat( 1,2 ) = mat( 2,1 ) = mat( 4,5 ) = mat( 5,4 ) = mat( 7,8 ) = mat( 8,7 ) = xy; + mat( 1,3 ) = mat( 3,1 ) = mat( 4,6 ) = mat( 6,4 ) = mat( 7,9 ) = mat( 9,7 ) = xz; + mat( 2,3 ) = mat( 3,2 ) = mat( 5,6 ) = mat( 6,5 ) = mat( 8,9 ) = mat( 9,8 ) = yz; + + math_Gauss solver( mat ); + if ( !solver.IsDone() ) + return false; + solver.Solve( vec ); + if ( vec.Norm2() < gp::Resolution() ) + return false; + // cout << endl + // << vec( 1 ) << "\t " << vec( 2 ) << "\t " << vec( 3 ) << endl + // << vec( 4 ) << "\t " << vec( 5 ) << "\t " << vec( 6 ) << endl + // << vec( 7 ) << "\t " << vec( 8 ) << "\t " << vec( 9 ) << endl; + + _srcOrig = srcOrig; + _trsf.SetTranslation( tgtOrig ); + + gp_Mat& M = const_cast< gp_Mat& >( _trsf.HVectorialPart() ); + M.SetRows( gp_XYZ( vec( 1 ), vec( 2 ), vec( 3 )), + gp_XYZ( vec( 4 ), vec( 5 ), vec( 6 )), + gp_XYZ( vec( 7 ), vec( 8 ), vec( 9 ))); + return true; + } + + //================================================================================ + /*! + * \brief Transforms a 3D point using a found transformation + */ + //================================================================================ + + gp_XYZ TrsfFinder3D::Transform( const gp_Pnt& srcP ) const + { + gp_XYZ p = srcP.XYZ() - _srcOrig; + _trsf.Transforms( p ); + return p; + } + + //================================================================================ + /*! + * \brief Transforms a 3D vector using a found transformation + */ + //================================================================================ + + gp_XYZ TrsfFinder3D::TransformVec( const gp_Vec& v ) const + { + return v.XYZ().Multiplied( _trsf.HVectorialPart() ); + } + //================================================================================ + /*! + * \brief Inversion + */ + //================================================================================ + + bool TrsfFinder3D::Invert() + { + if (( _trsf.Form() == gp_Translation ) && + ( _srcOrig.X() != 0 || _srcOrig.Y() != 0 || _srcOrig.Z() != 0 )) + { + // seems to be defined via Solve() + gp_XYZ newSrcOrig = _trsf.TranslationPart(); + gp_Mat& M = const_cast< gp_Mat& >( _trsf.HVectorialPart() ); + const double D = M.Determinant(); + if ( D < 1e-3 * ( newSrcOrig - _srcOrig ).Modulus() ) + { +#ifdef _DEBUG_ + cerr << "TrsfFinder3D::Invert()" + << "D " << M.Determinant() << " IsSingular " << M.IsSingular() << endl; +#endif + return false; + } + gp_Mat Minv = M.Inverted(); + _trsf.SetTranslation( _srcOrig ); + _srcOrig = newSrcOrig; + M = Minv; + } + else + { + _trsf.Invert(); + } + return true; + } +} diff --git a/src/StdMeshers/StdMeshers_ProjectionUtils.hxx b/src/StdMeshers/StdMeshers_ProjectionUtils.hxx index 31afb072f..e7ea4423d 100644 --- a/src/StdMeshers/StdMeshers_ProjectionUtils.hxx +++ b/src/StdMeshers/StdMeshers_ProjectionUtils.hxx @@ -30,10 +30,14 @@ #include "SMESH_StdMeshers.hxx" +#include "SMDS_MeshElement.hxx" + #include #include -#include #include +#include +#include +#include #include #include @@ -77,7 +81,54 @@ namespace StdMeshers_ProjectionUtils { typedef StdMeshers_ShapeShapeBiDirectionMap TShapeShapeMap; typedef TopTools_IndexedDataMapOfShapeListOfShape TAncestorMap; - typedef std::map TNodeNodeMap; + typedef std::map TNodeNodeMap; + + + /*! + * \brief Finds transformation beween two sets of 2D points using + * a least square approximation + */ + class TrsfFinder2D + { + gp_Trsf2d _trsf; + gp_XY _srcOrig; + public: + TrsfFinder2D(): _srcOrig(0,0) {} + + void Set( const gp_Trsf2d& t ) { _trsf = t; } // it's an alternative to Solve() + + bool Solve( const std::vector< gp_XY >& srcPnts, + const std::vector< gp_XY >& tgtPnts ); + + gp_XY Transform( const gp_Pnt2d& srcUV ) const; + + bool IsIdentity() const { return ( _trsf.Form() == gp_Identity ); } + }; + /*! + * \brief Finds transformation beween two sets of 3D points using + * a least square approximation + */ + class TrsfFinder3D + { + gp_Trsf _trsf; + gp_XYZ _srcOrig; + public: + TrsfFinder3D(): _srcOrig(0,0,0) {} + + void Set( const gp_Trsf& t ) { _trsf = t; } // it's an alternative to Solve() + + bool Solve( const std::vector< gp_XYZ > & srcPnts, + const std::vector< gp_XYZ > & tgtPnts ); + + gp_XYZ Transform( const gp_Pnt& srcP ) const; + + gp_XYZ TransformVec( const gp_Vec& v ) const; + + bool IsIdentity() const { return ( _trsf.Form() == gp_Identity ); } + + bool Invert(); + }; /*! * \brief Looks for association of all sub-shapes of two shapes diff --git a/src/StdMeshers/StdMeshers_Projection_2D.cxx b/src/StdMeshers/StdMeshers_Projection_2D.cxx index b0118c719..38832904d 100644 --- a/src/StdMeshers/StdMeshers_Projection_2D.cxx +++ b/src/StdMeshers/StdMeshers_Projection_2D.cxx @@ -47,8 +47,10 @@ #include "utilities.h" +#include #include #include +#include #include #include #include @@ -372,125 +374,197 @@ namespace { } // bool getBoundaryNodes() + //================================================================================ + /*! + * \brief Check if two consecutive EDGEs are connected in 2D + * \param [in] E1 - a well oriented non-seam EDGE + * \param [in] E2 - a possibly well oriented seam EDGE + * \param [in] F - a FACE + * \return bool - result + */ + //================================================================================ + + bool are2dConnected( const TopoDS_Edge & E1, + const TopoDS_Edge & E2, + const TopoDS_Face & F ) + { + double f,l; + Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( E1, F, f, l ); + gp_Pnt2d uvLast1 = c1->Value( E1.Orientation() == TopAbs_REVERSED ? f : l ); + + Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( E2, F, f, l ); + gp_Pnt2d uvFirst2 = c2->Value( f ); + gp_Pnt2d uvLast2 = c2->Value( l ); + double tol2 = 1e-5 * uvLast2.SquareDistance( uvFirst2 ); + + return (( uvLast1.SquareDistance( uvFirst2 ) < tol2 ) || + ( uvLast1.SquareDistance( uvLast2 ) < tol2 )); + } + + //================================================================================ + /*! + * \brief Compose TSideVector for both FACEs keeping matching order of EDGEs + */ + //================================================================================ + + bool getWires(const TopoDS_Face& tgtFace, + const TopoDS_Face& srcFace, + SMESH_Mesh * tgtMesh, + SMESH_Mesh * srcMesh, + const TAssocTool::TShapeShapeMap& shape2ShapeMap, + TSideVector& srcWires, + TSideVector& tgtWires, + const bool is1DComputed) + { + // get ordered src EDGEs + TError err; + srcWires = StdMeshers_FaceSide::GetFaceWires( srcFace, *srcMesh,/*skipMediumNodes=*/0, err); + if ( err && !err->IsOK() ) + return false; + + // make corresponding sequence of tgt EDGEs + tgtWires.resize( srcWires.size() ); + for ( size_t iW = 0; iW < srcWires.size(); ++iW ) + { + list< TopoDS_Edge > tgtEdges; + StdMeshers_FaceSidePtr srcWire = srcWires[iW]; + TopTools_IndexedMapOfShape edgeMap; // to detect seam edges + for ( int iE = 0; iE < srcWire->NbEdges(); ++iE ) + { + TopoDS_Edge E = TopoDS::Edge( shape2ShapeMap( srcWire->Edge( iE ), /*isSrc=*/true)); + // reverse a seam edge encountered for the second time + const int index = edgeMap.Add( E ); + if ( index < edgeMap.Extent() ) // E is a seam + { + // check which of edges to reverse, E or one already being in tgtEdges + if ( are2dConnected( tgtEdges.back(), E, tgtFace )) + { + list< TopoDS_Edge >::iterator eIt = tgtEdges.begin(); + std::advance( eIt, index-1 ); + eIt->Reverse(); + } + else + { + E.Reverse(); + } + } + if ( srcWire->NbEdges() == 1 && tgtMesh == srcMesh ) // circle + { + // try to verify ori by propagation + pair nE = StdMeshers_ProjectionUtils::GetPropagationEdge + ( srcMesh, E, srcWire->Edge( iE )); + if ( !nE.second.IsNull() ) + E = nE.second; + } + tgtEdges.push_back( E ); + } + tgtWires[ iW ].reset( new StdMeshers_FaceSide( tgtFace, tgtEdges, tgtMesh, + /*theIsForward = */ true, + /*theIgnoreMediumNodes = */false)); + if ( is1DComputed && + srcWires[iW]->GetUVPtStruct().size() != + tgtWires[iW]->GetUVPtStruct().size()) + return false; + } + return true; + } + //================================================================================ /*! * \brief Preform projection in case if tgtFace.IsPartner( srcFace ) and in case - * if projection by transformation is possible + * if projection by 3D transformation is possible */ //================================================================================ - bool projectPartner(const TopoDS_Face& tgtFace, - const TopoDS_Face& srcFace, - SMESH_Mesh * tgtMesh, - SMESH_Mesh * srcMesh, - const TAssocTool::TShapeShapeMap& shape2ShapeMap) + bool projectPartner(const TopoDS_Face& tgtFace, + const TopoDS_Face& srcFace, + const TSideVector& tgtWires, + const TSideVector& srcWires, + const TAssocTool::TShapeShapeMap& shape2ShapeMap, + TAssocTool::TNodeNodeMap& src2tgtNodes) { + SMESH_Mesh * tgtMesh = tgtWires[0]->GetMesh(); + SMESH_Mesh * srcMesh = srcWires[0]->GetMesh(); SMESHDS_Mesh* tgtMeshDS = tgtMesh->GetMeshDS(); SMESHDS_Mesh* srcMeshDS = srcMesh->GetMeshDS(); + SMESH_MesherHelper helper( *tgtMesh ); const double tol = 1.e-7 * srcMeshDS->getMaxDim(); - gp_Trsf trsf; // transformation to get location of target nodes from source ones + // transformation to get location of target nodes from source ones + StdMeshers_ProjectionUtils::TrsfFinder3D trsf; if ( tgtFace.IsPartner( srcFace )) { gp_Trsf srcTrsf = srcFace.Location(); gp_Trsf tgtTrsf = tgtFace.Location(); - trsf = srcTrsf.Inverted() * tgtTrsf; + trsf.Set( srcTrsf.Inverted() * tgtTrsf ); } else { - // Try to find the transformation - - // make any local coord systems of src and tgt faces - vector srcPP, tgtPP; // 3 points on face boundaries to make axes of CS - int tgtNbVert = SMESH_MesherHelper::Count( tgtFace, TopAbs_VERTEX, /*ignoreSame=*/true ); - int srcNbVert = SMESH_MesherHelper::Count( srcFace, TopAbs_VERTEX, /*ignoreSame=*/true ); - SMESH_subMesh * srcSM = srcMesh->GetSubMesh( srcFace ); - SMESH_subMeshIteratorPtr smIt = srcSM->getDependsOnIterator(/*includeSelf=*/false,false); - srcSM = smIt->next(); // sm of a vertex - while ( smIt->more() && srcPP.size() < 3 ) + // Try to find the 3D transformation + + const int totNbSeg = 50; + vector< gp_XYZ > srcPnts, tgtPnts; + srcPnts.reserve( totNbSeg ); + tgtPnts.reserve( totNbSeg ); + for ( size_t iW = 0; iW < srcWires.size(); ++iW ) { - srcSM = smIt->next(); - SMESHDS_SubMesh* srcSmds = srcSM->GetSubMeshDS(); - if ( !srcSmds ) continue; - SMDS_NodeIteratorPtr nIt = srcSmds->GetNodes(); - while ( nIt->more() ) + const double minSegLen = srcWires[iW]->Length() / totNbSeg; + for ( int iE = 0; iE < srcWires[iW]->NbEdges(); ++iE ) { - SMESH_TNodeXYZ p ( nIt->next()); - bool pOK = false; - switch ( srcPP.size() ) - { - case 0: pOK = true; break; - - case 1: pOK = ( srcPP[0].SquareDistance( p ) > 10*tol ); break; - - case 2: - { - gp_Vec p0p1( srcPP[0], srcPP[1] ), p0p( srcPP[0], p ); - // pOK = !p0p1.IsParallel( p0p, tol ); - pOK = !p0p1.IsParallel( p0p, 3.14/20 ); // angle min 18 degrees - break; - } - } - if ( !pOK ) - continue; - - // find corresponding point on target shape - pOK = false; - gp_Pnt tgtP; - const TopoDS_Shape& tgtShape = shape2ShapeMap( srcSM->GetSubShape(), /*isSrc=*/true ); - if ( tgtShape.ShapeType() == TopAbs_VERTEX ) + int nbSeg = Max( 1, int( srcWires[iW]->EdgeLength( iE ) / minSegLen )); + double srcU = srcWires[iW]->FirstParameter( iE ); + double tgtU = tgtWires[iW]->FirstParameter( iE ); + double srcDu = ( srcWires[iW]->LastParameter( iE )- srcU ) / nbSeg; + double tgtDu = ( tgtWires[iW]->LastParameter( iE )- tgtU ) / nbSeg; + for ( size_t i = 0; i < nbSeg; ++i ) { - tgtP = BRep_Tool::Pnt( TopoDS::Vertex( tgtShape )); - if ( srcNbVert == tgtNbVert || tgtPP.empty() ) - pOK = true; - else - pOK = (( tgtP.Distance( tgtPP[0] ) > tol*tol ) && - ( tgtPP.size() == 1 || tgtP.Distance( tgtPP[1] ) > tol*tol )); - //cout << "V - nS " << p._node->GetID() << " - nT " << SMESH_Algo::VertexNode(TopoDS::Vertex( tgtShape),tgtMeshDS)->GetID() << endl; + srcPnts.push_back( srcWires[iW]->Value3d( srcU ).XYZ() ); + tgtPnts.push_back( tgtWires[iW]->Value3d( tgtU ).XYZ() ); + srcU += srcDu; + tgtU += tgtDu; } - else if ( tgtPP.size() > 0 ) - { - if ( SMESHDS_SubMesh* tgtSmds = tgtMeshDS->MeshElements( tgtShape )) - { - double srcDist = srcPP[0].Distance( p ); - double eTol = BRep_Tool::Tolerance( TopoDS::Edge( tgtShape )); - if (eTol < tol) eTol = tol; - SMDS_NodeIteratorPtr nItT = tgtSmds->GetNodes(); - while ( nItT->more() && !pOK ) - { - const SMDS_MeshNode* n = nItT->next(); - tgtP = SMESH_TNodeXYZ( n ); - pOK = ( fabs( srcDist - tgtPP[0].Distance( tgtP )) < 2*eTol ); - //cout << "E - nS " << p._node->GetID() << " - nT " << n->GetID()<< " OK - " << pOK<< " " << fabs( srcDist - tgtPP[0].Distance( tgtP ))<< " tol " << eTol<< endl; - } - } - } - if ( !pOK ) - continue; - - srcPP.push_back( p ); - tgtPP.push_back( tgtP ); } } - if ( srcPP.size() != 3 ) + if ( !trsf.Solve( srcPnts, tgtPnts )) return false; - // make transformation - gp_Trsf fromTgtCS, toSrcCS; // from/to global CS - gp_Ax2 srcCS( srcPP[0], gp_Vec( srcPP[0], srcPP[1] ), gp_Vec( srcPP[0], srcPP[2])); - gp_Ax2 tgtCS( tgtPP[0], gp_Vec( tgtPP[0], tgtPP[1] ), gp_Vec( tgtPP[0], tgtPP[2])); - toSrcCS .SetTransformation( gp_Ax3( srcCS )); - fromTgtCS.SetTransformation( gp_Ax3( tgtCS )); - fromTgtCS.Invert(); + // check trsf - trsf = fromTgtCS * toSrcCS; + bool trsfIsOK = true; + const int nbTestPnt = 20; + const size_t iStep = Max( 1, int( srcPnts.size() / nbTestPnt )); + // check boundary + for ( size_t i = 0; ( i < srcPnts.size() && trsfIsOK ); i += iStep ) + { + gp_Pnt trsfTgt = trsf.Transform( srcPnts[i] ); + trsfIsOK = ( trsfTgt.SquareDistance( tgtPnts[i] ) < tol*tol ); + } + // check an in-FACE point + if ( trsfIsOK ) + { + BRepAdaptor_Surface srcSurf( srcFace ); + gp_Pnt srcP = + srcSurf.Value( 0.5 * ( srcSurf.FirstUParameter() + srcSurf.LastUParameter() ), + 0.5 * ( srcSurf.FirstVParameter() + srcSurf.LastVParameter() )); + gp_Pnt tgtTrsfP = trsf.Transform( srcP ); + TopLoc_Location loc; + GeomAPI_ProjectPointOnSurf& proj = helper.GetProjector( tgtFace, loc, 0.1*tol ); + if ( !loc.IsIdentity() ) + tgtTrsfP.Transform( loc.Transformation().Inverted() ); + proj.Perform( tgtTrsfP ); + trsfIsOK = ( proj.IsDone() && + proj.NbPoints() > 0 && + proj.LowerDistance() < tol ); + } + if ( !trsfIsOK ) + return false; } // Fill map of src to tgt nodes with nodes on edges - map src2tgtNodes; - map::iterator srcN_tgtN; + src2tgtNodes.clear(); + TAssocTool::TNodeNodeMap::iterator srcN_tgtN; bool tgtEdgesMeshed = false; for ( TopExp_Explorer srcExp( srcFace, TopAbs_EDGE); srcExp.More(); srcExp.Next() ) @@ -518,31 +592,6 @@ namespace { ) return false; - if ( !tgtEdge.IsPartner( srcEdge )) - { - if ( tgtNodes.empty() ) - return false; - // check that transformation is OK by three nodes - gp_Pnt p0S = SMESH_TNodeXYZ( (srcNodes.begin()) ->second); - gp_Pnt p1S = SMESH_TNodeXYZ( (srcNodes.rbegin()) ->second); - gp_Pnt p2S = SMESH_TNodeXYZ( (++srcNodes.begin())->second); - - gp_Pnt p0T = SMESH_TNodeXYZ( (tgtNodes.begin()) ->second); - gp_Pnt p1T = SMESH_TNodeXYZ( (tgtNodes.rbegin()) ->second); - gp_Pnt p2T = SMESH_TNodeXYZ( (++tgtNodes.begin())->second); - - // transform source points, they must coincide with target ones - if ( p0T.SquareDistance( p0S.Transformed( trsf )) > tol || - p1T.SquareDistance( p1S.Transformed( trsf )) > tol || - p2T.SquareDistance( p2S.Transformed( trsf )) > tol ) - { - //cout << "KO trsf, 3 dist: " - //<< p0T.SquareDistance( p0S.Transformed( trsf ))<< ", " - //<< p1T.SquareDistance( p1S.Transformed( trsf ))<< ", " - //<< p2T.SquareDistance( p2S.Transformed( trsf ))<< ", "<::iterator u_tn = tgtNodes.begin(); @@ -563,16 +612,6 @@ namespace { if ( !tgtN || tgtV.ShapeType() != TopAbs_VERTEX ) return false; - if ( !tgtV.IsPartner( srcV )) - { - // check that transformation is OK by three nodes - gp_Pnt p0S = SMESH_TNodeXYZ( srcN ); - gp_Pnt p0T = SMESH_TNodeXYZ( tgtN ); - if ( p0T.SquareDistance( p0S.Transformed( trsf )) > tol ) - { - return false; - } - } src2tgtNodes.insert( make_pair( srcN, tgtN )); } @@ -580,7 +619,6 @@ namespace { // Make new faces // prepare the helper to adding quadratic elements if necessary - SMESH_MesherHelper helper( *tgtMesh ); helper.SetSubShape( tgtFace ); helper.IsQuadraticSubMesh( tgtFace ); @@ -594,7 +632,7 @@ namespace { const SMDS_MeshNode* nullNode = 0; // indices of nodes to create properly oriented faces - bool isReverse = ( trsf.Form() != gp_Identity ); + bool isReverse = ( !trsf.IsIdentity() ); int tri1 = 1, tri2 = 2, quad1 = 1, quad3 = 3; if ( isReverse ) std::swap( tri1, tri2 ), std::swap( quad1, quad3 ); @@ -614,7 +652,7 @@ namespace { if ( srcN_tgtN->second == nullNode ) { // create a new node - gp_Pnt tgtP = gp_Pnt( SMESH_TNodeXYZ( srcNode )).Transformed( trsf ); + gp_Pnt tgtP = trsf.Transform( SMESH_TNodeXYZ( srcNode )); SMDS_MeshNode* n = helper.AddNode( tgtP.X(), tgtP.Y(), tgtP.Z() ); srcN_tgtN->second = n; switch ( srcNode->GetPosition()->GetTypeOfPosition() ) @@ -698,95 +736,29 @@ namespace { } // bool projectPartner() - //================================================================================ - /*! - * \brief Check if two consecutive EDGEs are connected in 2D - * \param [in] E1 - a well oriented non-seam EDGE - * \param [in] E2 - a possibly well oriented seam EDGE - * \param [in] F - a FACE - * \return bool - result - */ - //================================================================================ - - bool are2dConnected( const TopoDS_Edge & E1, - const TopoDS_Edge & E2, - const TopoDS_Face & F ) - { - double f,l; - Handle(Geom2d_Curve) c1 = BRep_Tool::CurveOnSurface( E1, F, f, l ); - gp_Pnt2d uvLast1 = c1->Value( E1.Orientation() == TopAbs_REVERSED ? f : l ); - - Handle(Geom2d_Curve) c2 = BRep_Tool::CurveOnSurface( E2, F, f, l ); - gp_Pnt2d uvFirst2 = c2->Value( f ); - gp_Pnt2d uvLast2 = c2->Value( l ); - double tol2 = 1e-5 * uvLast2.SquareDistance( uvFirst2 ); - - return (( uvLast1.SquareDistance( uvFirst2 ) < tol2 ) || - ( uvLast1.SquareDistance( uvLast2 ) < tol2 )); - } - //================================================================================ /*! * \brief Preform projection in case if the faces are similar in 2D space */ //================================================================================ - bool projectBy2DSimilarity(const TopoDS_Face& tgtFace, - const TopoDS_Face& srcFace, - SMESH_Mesh * tgtMesh, - SMESH_Mesh * srcMesh, - const TAssocTool::TShapeShapeMap& shape2ShapeMap, - const bool is1DComputed) + bool projectBy2DSimilarity(const TopoDS_Face& tgtFace, + const TopoDS_Face& srcFace, + const TSideVector& tgtWires, + const TSideVector& srcWires, + const TAssocTool::TShapeShapeMap& shape2ShapeMap, + const bool is1DComputed, + TAssocTool::TNodeNodeMap& src2tgtNodes) { - // 1) Preparation + SMESH_Mesh * tgtMesh = tgtWires[0]->GetMesh(); + SMESH_Mesh * srcMesh = srcWires[0]->GetMesh(); - // get ordered src EDGEs - TError err; - TSideVector srcWires = - StdMeshers_FaceSide::GetFaceWires( srcFace, *srcMesh,/*ignoreMediumNodes = */false, err); - if ( err && !err->IsOK() ) - return false; + // WARNING: we can have problems if the FACE is symmetrical in 2D, + // then the projection can be mirrored relating to what is expected - // make corresponding sequence of tgt EDGEs - TSideVector tgtWires( srcWires.size() ); - for ( size_t iW = 0; iW < srcWires.size(); ++iW ) - { - list< TopoDS_Edge > tgtEdges; - StdMeshers_FaceSidePtr srcWire = srcWires[iW]; - TopTools_IndexedMapOfShape edgeMap; // to detect seam edges - for ( int iE = 0; iE < srcWire->NbEdges(); ++iE ) - { - TopoDS_Edge E = TopoDS::Edge( shape2ShapeMap( srcWire->Edge( iE ), /*isSrc=*/true)); - // reverse a seam edge encountered for the second time - const int index = edgeMap.Add( E ); - if ( index < edgeMap.Extent() ) // E is a seam - { - // check which of edges to reverse, E or one already being in tgtEdges - if ( are2dConnected( tgtEdges.back(), E, tgtFace )) - { - list< TopoDS_Edge >::iterator eIt = tgtEdges.begin(); - std::advance( eIt, index-1 ); - eIt->Reverse(); - } - else - { - E.Reverse(); - } - } - tgtEdges.push_back( E ); - } - tgtWires[ iW ].reset( new StdMeshers_FaceSide( tgtFace, tgtEdges, tgtMesh, - /*theIsForward = */ true, - /*theIgnoreMediumNodes = */false)); - if ( is1DComputed && - srcWires[iW]->GetUVPtStruct().size() != - tgtWires[iW]->GetUVPtStruct().size()) - return false; - } - - // 2) Find transformation + // 1) Find 2D transformation - gp_Trsf2d trsf; + StdMeshers_ProjectionUtils::TrsfFinder2D trsf; { // get 2 pairs of corresponding UVs gp_Pnt2d srcP0 = srcWires[0]->Value2d(0.0); @@ -801,34 +773,74 @@ namespace { toSrcCS .SetTransformation( srcCS ); fromTgtCS.SetTransformation( tgtCS ); fromTgtCS.Invert(); - - trsf = fromTgtCS * toSrcCS; + trsf.Set( fromTgtCS * toSrcCS ); // check transformation + bool trsfIsOK = true; const double tol = 1e-5 * gp_Vec2d( srcP0, srcP1 ).Magnitude(); - for ( double u = 0.12; u < 1.; u += 0.1 ) + for ( double u = 0.12; ( u < 1. && trsfIsOK ); u += 0.1 ) + { + gp_Pnt2d srcUV = srcWires[0]->Value2d( u ); + gp_Pnt2d tgtUV = tgtWires[0]->Value2d( u ); + gp_Pnt2d tgtUV2 = trsf.Transform( srcUV ); + trsfIsOK = ( tgtUV.Distance( tgtUV2 ) < tol ); + } + + // Find trsf using a least-square approximation + if ( !trsfIsOK ) { - gp_Pnt2d srcUV = srcWires[0]->Value2d( u ); - gp_Pnt2d tgtUV = tgtWires[0]->Value2d( u ); - gp_Pnt2d tgtUV2 = srcUV.Transformed( trsf ); - if ( tgtUV.Distance( tgtUV2 ) > tol ) + // find trsf + const int totNbSeg = 50; + vector< gp_XY > srcPnts, tgtPnts; + srcPnts.resize( totNbSeg ); + tgtPnts.resize( totNbSeg ); + for ( size_t iW = 0; iW < srcWires.size(); ++iW ) + { + const double minSegLen = srcWires[iW]->Length() / totNbSeg; + for ( int iE = 0; iE < srcWires[iW]->NbEdges(); ++iE ) + { + int nbSeg = Max( 1, int( srcWires[iW]->EdgeLength( iE ) / minSegLen )); + double srcU = srcWires[iW]->FirstParameter( iE ); + double tgtU = tgtWires[iW]->FirstParameter( iE ); + double srcDu = ( srcWires[iW]->LastParameter( iE )- srcU ) / nbSeg; + double tgtDu = ( tgtWires[iW]->LastParameter( iE )- tgtU ) / nbSeg; + for ( size_t i = 0; i < nbSeg; ++i, srcU += srcDu, tgtU += tgtDu ) + { + srcPnts.push_back( srcWires[iW]->Value2d( srcU ).XY() ); + tgtPnts.push_back( tgtWires[iW]->Value2d( tgtU ).XY() ); + } + } + } + if ( !trsf.Solve( srcPnts, tgtPnts )) + return false; + + // check trsf + + trsfIsOK = true; + const int nbTestPnt = 10; + const size_t iStep = Max( 1, int( srcPnts.size() / nbTestPnt )); + for ( size_t i = 0; ( i < srcPnts.size() && trsfIsOK ); i += iStep ) + { + gp_Pnt2d trsfTgt = trsf.Transform( srcPnts[i] ); + trsfIsOK = ( trsfTgt.Distance( tgtPnts[i] ) < tol ); + } + if ( !trsfIsOK ) return false; } - } + } // "Find transformation" block - // 3) Projection + // 2) Projection - typedef map TN2NMap; - TN2NMap src2tgtNodes; - TN2NMap::iterator srcN_tgtN; + src2tgtNodes.clear(); + TAssocTool::TNodeNodeMap::iterator srcN_tgtN; // fill src2tgtNodes in with nodes on EDGEs - for ( unsigned iW = 0; iW < srcWires.size(); ++iW ) + for ( size_t iW = 0; iW < srcWires.size(); ++iW ) if ( is1DComputed ) { const vector& srcUVs = srcWires[iW]->GetUVPtStruct(); const vector& tgtUVs = tgtWires[iW]->GetUVPtStruct(); - for ( unsigned i = 0; i < srcUVs.size(); ++i ) + for ( size_t i = 0; i < srcUVs.size(); ++i ) src2tgtNodes.insert( make_pair( srcUVs[i].node, tgtUVs[i].node )); } else @@ -880,8 +892,8 @@ namespace { // create a new node gp_Pnt2d srcUV = srcHelper.GetNodeUV( srcFace, srcNode, elem->GetNode( helper.WrapIndex(i+1,nbN)), &uvOK); - gp_Pnt2d tgtUV = srcUV.Transformed( trsf ); - gp_Pnt tgtP = tgtSurface->Value( tgtUV.X(), tgtUV.Y() ); + gp_Pnt2d tgtUV = trsf.Transform( srcUV ); + gp_Pnt tgtP = tgtSurface->Value( tgtUV.X(), tgtUV.Y() ); SMDS_MeshNode* n = tgtMeshDS->AddNode( tgtP.X(), tgtP.Y(), tgtP.Z() ); switch ( srcNode->GetPosition()->GetTypeOfPosition() ) { @@ -920,6 +932,74 @@ namespace { } // bool projectBy2DSimilarity(...) + //================================================================================ + /*! + * \brief Fix bad faces by smoothing + */ + //================================================================================ + + void fixDistortedFaces( SMESH_MesherHelper& helper ) + { + // Detect bad faces + + bool haveBadFaces = false; + + const TopoDS_Face& F = TopoDS::Face( helper.GetSubShape() ); + SMESHDS_SubMesh* smDS = helper.GetMeshDS()->MeshElements( F ); + if ( !smDS || smDS->NbElements() == 0 ) return; + + SMDS_ElemIteratorPtr faceIt = smDS->GetElements(); + double prevArea2D = 0; + vector< const SMDS_MeshNode* > nodes; + vector< gp_XY > uv; + while ( faceIt->more() && !haveBadFaces ) + { + const SMDS_MeshElement* face = faceIt->next(); + + // get nodes + nodes.resize( face->NbCornerNodes() ); + SMDS_MeshElement::iterator n = face->begin_nodes(); + for ( size_t i = 0; i < nodes.size(); ++n, ++i ) + nodes[ i ] = *n; + + // get UVs + const SMDS_MeshNode* inFaceNode = 0; + if ( helper.HasSeam() ) + for ( size_t i = 0; ( i < nodes.size() && !inFaceNode ); ++i ) + if ( !helper.IsSeamShape( nodes[ i ]->getshapeId() )) + inFaceNode = nodes[ i ]; + + uv.resize( nodes.size() ); + for ( size_t i = 0; i < nodes.size(); ++i ) + uv[ i ] = helper.GetNodeUV( F, nodes[ i ], inFaceNode ); + + // compare orientation of triangles + for ( int iT = 0, nbT = nodes.size()-2; iT < nbT; ++iT ) + { + gp_XY v1 = uv[ iT+1 ] - uv[ 0 ]; + gp_XY v2 = uv[ iT+2 ] - uv[ 0 ]; + double area2D = v2 ^ v1; + if (( haveBadFaces = ( area2D * prevArea2D < 0 ))) + break; + prevArea2D = area2D; + } + } + + // Fix faces + + if ( haveBadFaces ) + { + SMESH_MeshEditor editor( helper.GetMesh() ); + + TIDSortedElemSet faces; + for ( faceIt = smDS->GetElements(); faceIt->more(); ) + faces.insert( faces.end(), faceIt->next() ); + + set fixedNodes; + editor.Smooth( faces, fixedNodes, SMESH_MeshEditor::CENTROIDAL, 5 ); + } + } + } // namespace @@ -930,6 +1010,8 @@ namespace { bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& theShape) { + _src2tgtNodes.clear(); + MESSAGE("Projection_2D Compute"); if ( !_sourceHypo ) return false; @@ -1001,17 +1083,25 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& is1DComputed = sm->IsMeshComputed(); } + // get ordered src and tgt EDGEs + TSideVector srcWires, tgtWires; + if ( !getWires( tgtFace, srcFace, tgtMesh, srcMesh, + shape2ShapeMap, srcWires, tgtWires, is1DComputed )) + return false; + bool done = false; - if ( !done ) + if ( !done ) { // try to project from the same face with different location - done = projectPartner( tgtFace, srcFace, tgtMesh, srcMesh, shape2ShapeMap ); + done = projectPartner( tgtFace, srcFace, tgtWires, srcWires, + shape2ShapeMap, _src2tgtNodes ); } if ( !done ) { // projection in case if the faces are similar in 2D space - done = projectBy2DSimilarity( tgtFace, srcFace, tgtMesh, srcMesh, shape2ShapeMap, is1DComputed); + done = projectBy2DSimilarity( tgtFace, srcFace, tgtWires, srcWires, + shape2ShapeMap, is1DComputed, _src2tgtNodes); } SMESH_MesherHelper helper( theMesh ); @@ -1019,6 +1109,7 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& if ( !done ) { + _src2tgtNodes.clear(); // -------------------- // Prepare to mapping // -------------------- @@ -1100,7 +1191,7 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& return error(COMPERR_BAD_INPUT_MESH,"Can't load mesh pattern from the source face"); // -------------------- - // Perform 2D mapping + // Perform 2D mapping // -------------------- // Compute mesh on a target face @@ -1278,6 +1369,14 @@ bool StdMeshers_Projection_2D::Compute(SMESH_Mesh& theMesh, const TopoDS_Shape& if ( nbFaceBeforeMerge != nbFaceAtferMerge && !helper.HasDegeneratedEdges() ) return error(COMPERR_BAD_INPUT_MESH, "Probably invalid node parameters on geom faces"); + + // ---------------------------------------------------------------- + // The mapper can create distorted faces by placing nodes out of the FACE + // boundary -- fix bad faces by smoothing + // ---------------------------------------------------------------- + + fixDistortedFaces( helper ); + // ---------------------------------------------------------------- // The mapper can't create quadratic elements, so convert if needed // ---------------------------------------------------------------- diff --git a/src/StdMeshers/StdMeshers_Projection_2D.hxx b/src/StdMeshers/StdMeshers_Projection_2D.hxx index 57e34dbc4..878d21849 100644 --- a/src/StdMeshers/StdMeshers_Projection_2D.hxx +++ b/src/StdMeshers/StdMeshers_Projection_2D.hxx @@ -30,6 +30,7 @@ #include "SMESH_StdMeshers.hxx" #include "SMESH_Algo.hxx" +#include "StdMeshers_ProjectionUtils.hxx" class StdMeshers_ProjectionSource2D; @@ -59,9 +60,12 @@ public: */ virtual void SetEventListener(SMESH_subMesh* whenSetToSubMesh); -protected: - const StdMeshers_ProjectionSource2D* _sourceHypo; + protected: + + const StdMeshers_ProjectionSource2D* _sourceHypo; + + StdMeshers_ProjectionUtils::TNodeNodeMap _src2tgtNodes; }; diff --git a/src/StdMeshers/StdMeshers_Projection_3D.cxx b/src/StdMeshers/StdMeshers_Projection_3D.cxx index 3079c7bd8..93969dfeb 100644 --- a/src/StdMeshers/StdMeshers_Projection_3D.cxx +++ b/src/StdMeshers/StdMeshers_Projection_3D.cxx @@ -307,7 +307,7 @@ bool StdMeshers_Projection_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& aS } } // Find matching nodes of tgt and src faces - TNodeNodeMap faceMatchingNodes; + TAssocTool::TNodeNodeMap faceMatchingNodes; if ( ! TAssocTool::FindMatchingNodesOnFaces( srcFace, srcMesh, tgtFace, tgtMesh, shape2ShapeMap, faceMatchingNodes )) return error(COMPERR_BAD_INPUT_MESH,SMESH_Comment("Mesh on faces #") diff --git a/src/StdMeshers/StdMeshers_RadialPrism_3D.cxx b/src/StdMeshers/StdMeshers_RadialPrism_3D.cxx index 17411417e..a06cf8608 100644 --- a/src/StdMeshers/StdMeshers_RadialPrism_3D.cxx +++ b/src/StdMeshers/StdMeshers_RadialPrism_3D.cxx @@ -239,7 +239,7 @@ bool StdMeshers_RadialPrism_3D::Compute(SMESH_Mesh& aMesh, const TopoDS_Shape& a } // Find matching nodes of in and out faces - TNodeNodeMap nodeIn2OutMap; + ProjectionUtils::TNodeNodeMap nodeIn2OutMap; if ( ! ProjectionUtils::FindMatchingNodesOnFaces( inFace, &aMesh, outFace, &aMesh, shape2ShapeMap, nodeIn2OutMap )) return error(COMPERR_BAD_INPUT_MESH,SMESH_Comment("Mesh on faces #") diff --git a/src/StdMeshers/StdMeshers_ViscousLayers.cxx b/src/StdMeshers/StdMeshers_ViscousLayers.cxx index c04a6326f..eeec08ff2 100644 --- a/src/StdMeshers/StdMeshers_ViscousLayers.cxx +++ b/src/StdMeshers/StdMeshers_ViscousLayers.cxx @@ -101,7 +101,7 @@ namespace VISCOUS_3D const double theSmoothThickToElemSizeRatio = 0.3; // what part of thickness is allowed till intersection - // defined by SALOME_TESTS/Grids/smesh/viscous_layers_00/A5 + // (defined by SALOME_TESTS/Grids/smesh/viscous_layers_00/A5) const double theThickToIntersection = 1.5; bool needSmoothing( double cosin, double tgtThick, double elemSize ) -- 2.30.2