X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MesherHelper.cxx;h=22b28da7d571a4156a4b60a7bb07b2d21c9dfb43;hp=f1b39ed75c15c29d1e7e3a6112c40f514bdd2820;hb=0b1674e3557c3f2d5028b6f24e44a16b6e19bd2d;hpb=929ba51821ae4494080b840a7c30e2db0f078f5d diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index f1b39ed75..22b28da7d 100644 --- a/src/SMESH/SMESH_MesherHelper.cxx +++ b/src/SMESH/SMESH_MesherHelper.cxx @@ -1,4 +1,4 @@ -// Copyright (C) 2007-2014 CEA/DEN, EDF R&D, OPEN CASCADE +// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE // // Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS @@ -70,7 +70,7 @@ using namespace std; namespace { - gp_XYZ XYZ(const SMDS_MeshNode* n) { return gp_XYZ(n->X(), n->Y(), n->Z()); } + inline SMESH_TNodeXYZ XYZ(const SMDS_MeshNode* n) { return SMESH_TNodeXYZ(n); } enum { U_periodic = 1, V_periodic = 2 }; } @@ -302,10 +302,20 @@ void SMESH_MesherHelper::SetSubShape(const TopoDS_Shape& aSh) isSeam = ( Abs( uv1.Coord(2) - myPar1[1] ) < Precision::PConfusion() || Abs( uv1.Coord(2) - myPar2[1] ) < Precision::PConfusion() ); } + if ( isSeam ) // vertices are on period boundary, check a middle point (23032) + { + double f,l, r = 0.2345; + Handle(Geom2d_Curve) C2d = BRep_Tool::CurveOnSurface( edge, face, f, l ); + uv2 = C2d->Value( f * r + l * ( 1.-r )); + if ( du < Precision::PConfusion() ) + isSeam = ( Abs( uv1.Coord(1) - uv2.Coord(1) ) < Precision::PConfusion() ); + else + isSeam = ( Abs( uv1.Coord(2) - uv2.Coord(2) ) < Precision::PConfusion() ); + } } if ( isSeam ) { - // store seam shape indices, negative if shape encounters twice + // store seam shape indices, negative if shape encounters twice ('real seam') mySeamShapeIds.insert( IsSeamShape( edgeID ) ? -edgeID : edgeID ); for ( TopExp_Explorer v( edge, TopAbs_VERTEX ); v.More(); v.Next() ) { int vertexID = meshDS->ShapeToIndex( v.Current() ); @@ -579,7 +589,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, const SMDS_FacePosition* fpos = static_cast( Pos ); uv.SetCoord( fpos->GetUParameter(), fpos->GetVParameter() ); if ( check ) - uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*getFaceMaxTol( F )); + uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 2.*getFaceMaxTol( F )); // 2. from 22830 } else if ( Pos->GetTypeOfPosition() == SMDS_TOP_EDGE ) { @@ -595,7 +605,7 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, if ( validU ) uv = C2d->Value( u ); else uv.SetCoord( Precision::Infinite(),0.); if ( check || !validU ) - uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*getFaceMaxTol( F ),/*force=*/ !validU ); + uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 2.*getFaceMaxTol( F ),/*force=*/ !validU ); // for a node on a seam EDGE select one of UVs on 2 pcurves if ( n2 && IsSeamShape( edgeID )) @@ -691,10 +701,10 @@ gp_XY SMESH_MesherHelper::GetNodeUV(const TopoDS_Face& F, } else { - uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 10*getFaceMaxTol( F )); + uvOK = CheckNodeUV( F, n, uv.ChangeCoord(), 2.*getFaceMaxTol( F )); } - if ( check ) + if ( check && !uvOK ) *check = uvOK; return uv.XY(); @@ -1415,21 +1425,20 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetCentralNode(const SMDS_MeshNode* n1, TopoDS_Face F; gp_XY uvAvg; - bool badTria=false; if ( shapeType == TopAbs_FACE ) { F = TopoDS::Face( meshDS->IndexToShape( faceID )); - bool check; - gp_XY uv1 = GetNodeUV( F, n1, n23, &check ); - gp_XY uv2 = GetNodeUV( F, n2, n31, &check ); - gp_XY uv3 = GetNodeUV( F, n3, n12, &check ); - gp_XY uv12 = GetNodeUV( F, n12, n3, &check ); - gp_XY uv23 = GetNodeUV( F, n23, n1, &check ); - gp_XY uv31 = GetNodeUV( F, n31, n2, &check ); + bool checkOK = true, badTria = false; + gp_XY uv1 = GetNodeUV( F, n1, n23, &checkOK ); + gp_XY uv2 = GetNodeUV( F, n2, n31, &checkOK ); + gp_XY uv3 = GetNodeUV( F, n3, n12, &checkOK ); + gp_XY uv12 = GetNodeUV( F, n12, n3, &checkOK ); + gp_XY uv23 = GetNodeUV( F, n23, n1, &checkOK ); + gp_XY uv31 = GetNodeUV( F, n31, n2, &checkOK ); uvAvg = GetCenterUV( uv1,uv2,uv3, uv12,uv23,uv31, &badTria ); - if ( badTria ) - force3d = false; + if ( badTria || !checkOK ) + force3d = true; } // Create a central node @@ -1499,7 +1508,7 @@ const SMDS_MeshNode* SMESH_MesherHelper::GetMediumNode(const SMDS_MeshNode* n1, int faceID = -1, edgeID = -1; TopoDS_Edge E; double u [2]; TopoDS_Face F; gp_XY uv[2]; - bool uvOK[2] = { false, false }; + bool uvOK[2] = { true, true }; const bool useCurSubShape = ( !myShape.IsNull() && myShape.ShapeType() == TopAbs_EDGE ); pair pos = GetMediumPos( n1, n2, useCurSubShape, expectedSupport ); @@ -2734,38 +2743,80 @@ bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace) if ( !aSubMeshDSFace ) return isReversed; - // find an element with a good normal - gp_Vec Ne; - bool normalOK = false; - gp_XY uv; + // find an element on a bounday of theFace SMDS_ElemIteratorPtr iteratorElem = aSubMeshDSFace->GetElements(); - while ( !normalOK && iteratorElem->more() ) // loop on elements on theFace + const SMDS_MeshNode* nn[2]; + while ( iteratorElem->more() ) // loop on elements on theFace { const SMDS_MeshElement* elem = iteratorElem->next(); - if ( elem && elem->NbCornerNodes() > 2 ) + if ( ! elem ) continue; + + // look for 2 nodes on EDGE + int nbNodes = elem->NbCornerNodes(); + nn[0] = elem->GetNode( nbNodes-1 ); + for ( int iN = 0; iN < nbNodes; ++iN ) { - SMESH_TNodeXYZ nPnt[3]; - SMDS_ElemIteratorPtr nodesIt = elem->nodesIterator(); - int iNodeOnFace = 0, iPosDim = SMDS_TOP_VERTEX; - for ( int iN = 0; nodesIt->more() && iN < 3; ++iN) // loop on nodes + nn[1] = elem->GetNode( iN ); + if ( nn[0]->GetPosition()->GetDim() < 2 && + nn[1]->GetPosition()->GetDim() < 2 ) { - nPnt[ iN ] = nodesIt->next(); - if ( nPnt[ iN ]._node->GetPosition()->GetTypeOfPosition() > iPosDim ) + TopoDS_Shape s0 = GetSubShapeByNode( nn[0], GetMeshDS() ); + TopoDS_Shape s1 = GetSubShapeByNode( nn[1], GetMeshDS() ); + TopoDS_Shape E = GetCommonAncestor( s0, s1, *myMesh, TopAbs_EDGE ); + if ( !E.IsNull() && !s0.IsSame( s1 )) { - iNodeOnFace = iN; - iPosDim = nPnt[ iN ]._node->GetPosition()->GetTypeOfPosition(); + // is E seam edge? + int nb = 0; + for ( TopExp_Explorer exp( theFace, TopAbs_EDGE ); exp.More(); exp.Next() ) + if ( E.IsSame( exp.Current() )) { + ++nb; + E = exp.Current(); // to know orientation + } + if ( nb == 1 ) + { + bool ok = true; + double u0 = GetNodeU( TopoDS::Edge( E ), nn[0], nn[1], &ok ); + double u1 = GetNodeU( TopoDS::Edge( E ), nn[1], nn[0], &ok ); + if ( ok ) + { + isReversed = ( u0 > u1 ); + if ( E.Orientation() == TopAbs_REVERSED ) + isReversed = !isReversed; + return isReversed; + } + } } } - // compute normal - gp_Vec v01( nPnt[0], nPnt[1] ), v02( nPnt[0], nPnt[2] ); - if ( v01.SquareMagnitude() > RealSmall() && - v02.SquareMagnitude() > RealSmall() ) + nn[0] = nn[1]; + } + } + + // find an element with a good normal + gp_Vec Ne; + bool normalOK = false; + gp_XY uv; + iteratorElem = aSubMeshDSFace->GetElements(); + while ( !normalOK && iteratorElem->more() ) // loop on elements on theFace + { + const SMDS_MeshElement* elem = iteratorElem->next(); + if ( ! SMESH_MeshAlgos::FaceNormal( elem, const_cast( Ne.XYZ() ), /*normalized=*/0 )) + continue; + normalOK = true; + + // get UV of a node inside theFACE + SMDS_ElemIteratorPtr nodesIt = elem->nodesIterator(); + const SMDS_MeshNode* nInFace = 0; + int iPosDim = SMDS_TOP_VERTEX; + while ( nodesIt->more() ) // loop on nodes + { + const SMDS_MeshNode* n = static_cast( nodesIt->next() ); + if ( n->GetPosition()->GetTypeOfPosition() >= iPosDim ) { - Ne = v01 ^ v02; - if (( normalOK = ( Ne.SquareMagnitude() > RealSmall() ))) - uv = GetNodeUV( theFace, nPnt[iNodeOnFace]._node, 0, &normalOK ); + nInFace = n; + iPosDim = n->GetPosition()->GetTypeOfPosition(); } } + uv = GetNodeUV( theFace, nInFace, 0, &normalOK ); } if ( !normalOK ) return isReversed; @@ -2776,11 +2827,8 @@ bool SMESH_MesherHelper::IsReversedSubMesh (const TopoDS_Face& theFace) // if ( surf.IsNull() || surf->Continuity() < GeomAbs_C1 ) // some surfaces not detected as GeomAbs_C1 are nevertheless correct for meshing if ( surf.IsNull() || surf->Continuity() < GeomAbs_C0 ) - { - if (!surf.IsNull()) - MESSAGE("surf->Continuity() < GeomAbs_C1 " << (surf->Continuity() < GeomAbs_C1)); - return isReversed; - } + return isReversed; + gp_Vec d1u, d1v; gp_Pnt p; surf->D1( uv.X(), uv.Y(), p, d1u, d1v ); gp_Vec Nf = (d1u ^ d1v).Transformed( loc ); @@ -3208,6 +3256,11 @@ TopoDS_Shape SMESH_MesherHelper::GetCommonAncestor(const TopoDS_Shape& shape1, TopoDS_Shape commonAnc; if ( !shape1.IsNull() && !shape2.IsNull() ) { + if ( shape1.ShapeType() == ancestorType && IsSubShape( shape2, shape1 )) + return shape1; + if ( shape2.ShapeType() == ancestorType && IsSubShape( shape1, shape2 )) + return shape2; + PShapeIteratorPtr ancIt = GetAncestors( shape1, mesh, ancestorType ); while ( const TopoDS_Shape* anc = ancIt->next() ) if ( IsSubShape( shape2, *anc )) @@ -3248,12 +3301,13 @@ namespace { // Structures used by FixQuadraticElements() mutable vector _faces; mutable gp_Vec _nodeMove; mutable int _nbMoves; + mutable bool _is2dFixed; // is moved along surface or in 3D QLink(const SMDS_MeshNode* n1, const SMDS_MeshNode* n2, const SMDS_MeshNode* nm): SMESH_TLink( n1,n2 ), _mediumNode(nm), _nodeMove(0,0,0), _nbMoves(0) { _faces.reserve(4); - //if ( MediumPos() != SMDS_TOP_3DSPACE ) - _nodeMove = MediumPnt() - MiddlePnt(); + _nodeMove = MediumPnt() - MiddlePnt(); + _is2dFixed = ( MediumPos() != SMDS_TOP_FACE ); } void SetContinuesFaces() const; const QFace* GetContinuesFace( const QFace* face ) const; @@ -3268,10 +3322,11 @@ namespace { // Structures used by FixQuadraticElements() const SMDS_MeshNode* EndPosNode(SMDS_TypeOfPosition pos) const { return EndPos(0) == pos ? node1() : EndPos(1) == pos ? node2() : 0; } - void Move(const gp_Vec& move, bool sum=false) const - { _nodeMove += move; _nbMoves += sum ? (_nbMoves==0) : 1; } + void Move(const gp_Vec& move, bool sum=false, bool is2dFixed=false) const + { _nodeMove += move; _nbMoves += sum ? (_nbMoves==0) : 1; _is2dFixed |= is2dFixed; } gp_XYZ Move() const { return _nodeMove.XYZ() / _nbMoves; } bool IsMoved() const { return (_nbMoves > 0 /*&& !IsStraight()*/); } + bool IsFixedOnSurface() const { return _is2dFixed; } bool IsStraight() const { return isStraightLink( (XYZ(node1())-XYZ(node2())).SquareModulus(), _nodeMove.SquareMagnitude()); @@ -3750,7 +3805,7 @@ namespace { // Structures used by FixQuadraticElements() double r = thePrevLen / fullLen; gp_Vec move = linkNorm * refProj * ( 1 - r ); - theLink->Move( move, true ); + theLink->Move( move, /*sum=*/true ); MSG(string(theStep,'.')<<" Move "<< theLink->_mediumNode->GetID()<< " by " << refProj * ( 1 - r ) << " following " << @@ -4741,7 +4796,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, } else if ( error == ERR_TRI ) { // chain contains continues triangles TSplitTriaResult res = splitTrianglesIntoChains( rawChain, chains, pos ); - if ( res != _OK ) { // not quadrangles split into triangles + if ( res != _OK ) { // not 'quadrangles split into triangles' in chain fixTriaNearBoundary( rawChain, *this ); break; } @@ -4868,6 +4923,7 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, gp_Vec x = x01.Normalized() + x12.Normalized(); trsf.SetTransformation( gp_Ax3( gp::Origin(), link1->Normal(), x), gp_Ax3() ); move.Transform(trsf); + (*link1)->Move( move, /*sum=*/false, /*is2dFixed=*/false ); } else { // compute 3D displacement by 2D one @@ -4892,8 +4948,8 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, "newUV: "<Move( move, /*sum=*/false, /*is2dFixed=*/true ); } - (*link1)->Move( move ); MSG( "Move " << (*link1)->_mediumNode->GetID() << " following " << chain.front()->_mediumNode->GetID() <<"-" << chain.back ()->_mediumNode->GetID() << @@ -4912,11 +4968,28 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, const bool toFixCentralNodes = ( myMesh->NbBiQuadQuadrangles() + myMesh->NbBiQuadTriangles() + myMesh->NbTriQuadraticHexas() ); + double distXYZ[4]; + faceHlp.ToFixNodeParameters( true ); for ( pLink = links.begin(); pLink != links.end(); ++pLink ) { if ( pLink->IsMoved() ) { gp_Pnt p = pLink->MiddlePnt() + pLink->Move(); + + // put on surface nodes on FACE but moved in 3D (23050) + if ( !pLink->IsFixedOnSurface() ) + { + faceHlp.SetSubShape( pLink->_mediumNode->getshapeId() ); + if ( faceHlp.GetSubShape().ShapeType() == TopAbs_FACE ) + { + const_cast( pLink->_mediumNode )->setXYZ( p.X(), p.Y(), p.Z()); + p.Coord( distXYZ[1], distXYZ[2], distXYZ[3] ); + gp_XY uv( Precision::Infinite(), 0 ); + if ( faceHlp.CheckNodeUV( TopoDS::Face( faceHlp.GetSubShape() ), pLink->_mediumNode, + uv, /*tol=*/pLink->Move().Modulus(), /*force=*/true, distXYZ )) + p.SetCoord( distXYZ[1], distXYZ[2], distXYZ[3] ); + } + } GetMeshDS()->MoveNode( pLink->_mediumNode, p.X(), p.Y(), p.Z()); // collect bi-quadratic elements @@ -4990,17 +5063,28 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, // nodes nodes.assign( tria->begin_nodes(), tria->end_nodes() ); // UV + bool uvOK = true, badTria; for ( int i = 0; i < 6; ++i ) { - uv[ i ] = GetNodeUV( F, nodes[i], nodes[(i+1)%3], &checkUV ); + uv[ i ] = GetNodeUV( F, nodes[i], nodes[(i+1)%3], &uvOK ); // as this method is used after mesh generation, UV of nodes is not // updated according to bending links, so we update if ( nodes[i]->GetPosition()->GetTypeOfPosition() == SMDS_TOP_FACE ) CheckNodeUV( F, nodes[i], uv[ i ], 2*tol, /*force=*/true ); } // move the central node - gp_XY uvCent = GetCenterUV( uv[0], uv[1], uv[2], uv[3], uv[4], uv[5] ); - gp_Pnt p = surf->Value( uvCent.X(), uvCent.Y() ).Transformed( loc ); + gp_Pnt p; + if ( !uvOK || badTria ) + { + p = ( SMESH_TNodeXYZ( nodes[3] ) + + SMESH_TNodeXYZ( nodes[4] ) + + SMESH_TNodeXYZ( nodes[5] )) / 3; + } + else + { + gp_XY uvCent = GetCenterUV( uv[0], uv[1], uv[2], uv[3], uv[4], uv[5], &badTria ); + p = surf->Value( uvCent.X(), uvCent.Y() ).Transformed( loc ); + } GetMeshDS()->MoveNode( tria->GetNode(6), p.X(), p.Y(), p.Z() ); } }