X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MesherHelper.cxx;h=1bc159afa2d0c95e7811ad26858037e75ba86147;hp=0556659dfe412cea035503dd439117b6c4d4529b;hb=ebb518296392b51b6ca9d80037a12703a160518f;hpb=8f3082ad8e9e635572a99753b7d0a5e86c95adc4 diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index 0556659df..1bc159afa 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() ); @@ -2571,7 +2581,7 @@ bool SMESH_MesherHelper::IsStructured( SMESH_subMesh* faceSM ) faceAnalyser.SetSubShape( faceSM->GetSubShape() ); // rotate edges to get the first node being at corner - // (in principle it's not necessary but so far none SALOME algo can make + // (in principle it's not necessary because so far none SALOME algo can make // such a structured mesh that all corner nodes are not on VERTEXes) bool isCorner = false; int nbRemainEdges = nbEdgesInWires.front(); @@ -2733,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; @@ -2775,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 ); @@ -3207,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 )) @@ -3247,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; @@ -3267,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()); @@ -3749,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 " << @@ -4740,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; } @@ -4768,12 +4824,14 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, // mesure chain length and compute link position along the chain double chainLen = 0; vector< double > linkPos; + TChain savedChain; // backup MSGBEG( "Link medium nodes: "); TChain::iterator link0 = chain.begin(), link1 = chain.begin(), link2; for ( ++link1; link1 != chain.end(); ++link1, ++link0 ) { MSGBEG( (*link0)->_mediumNode->GetID() << "-" <<(*link1)->_mediumNode->GetID()<<" "); double len = ((*link0)->MiddlePnt() - (*link1)->MiddlePnt()).Modulus(); while ( len < numeric_limits::min() ) { // remove degenerated link + if ( savedChain.empty() ) savedChain = chain; link1 = chain.erase( link1 ); if ( link1 == chain.end() ) break; @@ -4783,9 +4841,16 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, linkPos.push_back( chainLen ); } MSG(""); - if ( linkPos.size() < 2 ) - continue; - + if ( linkPos.size() <= 2 && savedChain.size() > 2 ) { + //continue; + linkPos.clear(); + chainLen = 0; + chain = savedChain; + for ( link1 = chain.begin(); link1 != chain.end(); ++link1 ) { + chainLen += 1; + linkPos.push_back( chainLen ); + } + } gp_Vec move0 = chain.front()->_nodeMove; gp_Vec move1 = chain.back ()->_nodeMove; @@ -4864,9 +4929,14 @@ void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& compError, // transform to global gp_Vec x01( (*link0)->MiddlePnt(), (*link1)->MiddlePnt() ); gp_Vec x12( (*link1)->MiddlePnt(), (*link2)->MiddlePnt() ); - gp_Vec x = x01.Normalized() + x12.Normalized(); - trsf.SetTransformation( gp_Ax3( gp::Origin(), link1->Normal(), x), gp_Ax3() ); + try { + gp_Vec x = x01.Normalized() + x12.Normalized(); + trsf.SetTransformation( gp_Ax3( gp::Origin(), link1->Normal(), x), gp_Ax3() ); + } catch ( Standard_Failure ) { + trsf.Invert(); + } move.Transform(trsf); + (*link1)->Move( move, /*sum=*/false, /*is2dFixed=*/false ); } else { // compute 3D displacement by 2D one @@ -4891,8 +4961,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() << @@ -4911,11 +4981,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