X-Git-Url: http://git.salome-platform.org/gitweb/?p=modules%2Fsmesh.git;a=blobdiff_plain;f=src%2FSMESH%2FSMESH_MesherHelper.cxx;h=ffca47e6aa48d5dc3a864f35b634f2bbde2707d7;hp=efaf7a853362f95f3f2dc082aec184ab173b4e68;hb=f5016d85b7b4b88623723027a1585c6414c4dc66;hpb=e4f02cdb389c8e4170ac26760a3f0257a009fd3b diff --git a/src/SMESH/SMESH_MesherHelper.cxx b/src/SMESH/SMESH_MesherHelper.cxx index efaf7a853..ffca47e6a 100644 --- a/src/SMESH/SMESH_MesherHelper.cxx +++ b/src/SMESH/SMESH_MesherHelper.cxx @@ -26,13 +26,17 @@ // #include "SMESH_MesherHelper.hxx" -#include "SMDS_FacePosition.hxx" #include "SMDS_EdgePosition.hxx" +#include "SMDS_FaceOfNodes.hxx" +#include "SMDS_FacePosition.hxx" +#include "SMDS_IteratorOnIterators.hxx" #include "SMDS_VolumeTool.hxx" -#include "SMESH_subMesh.hxx" #include "SMESH_ProxyMesh.hxx" +#include "SMESH_subMesh.hxx" +#include #include +#include #include #include #include @@ -403,7 +407,7 @@ void SMESH_MesherHelper::AddTLinks(const SMDS_MeshVolume* volume) { int iN1 = iNodes[i++]; int iN12 = iNodes[i++]; - int iN2 = iNodes[i++]; + int iN2 = iNodes[i]; if ( iN1 > iN2 ) std::swap( iN1, iN2 ); int linkID = iN1 * vTool.NbNodes() + iN2; pair< set::iterator, bool > it_isNew = addedLinks.insert( linkID ); @@ -1675,6 +1679,24 @@ SMESH_MesherHelper::AddPolyhedralVolume (const std::vector return elem; } +namespace +{ + //================================================================================ + /*! + * \brief Check if a node belongs to any face of sub-mesh + */ + //================================================================================ + + bool isNodeInSubMesh( const SMDS_MeshNode* n, const SMESHDS_SubMesh* sm ) + { + SMDS_ElemIteratorPtr fIt = n->GetInverseElementIterator( SMDSAbs_Face ); + while ( fIt->more() ) + if ( sm->Contains( fIt->next() )) + return true; + return false; + } +} + //======================================================================= //function : LoadNodeColumns //purpose : Load nodes bound to face into a map of node columns @@ -1746,13 +1768,36 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2 SMESH_Algo::GetSortedNodesOnEdge( theMesh, *edge,/*noMedium=*/true, sortedBaseNodes); if ( sortedBaseNodes.empty() ) continue; + map< double, const SMDS_MeshNode*>::iterator u_n = sortedBaseNodes.begin(); + if ( theProxyMesh ) // from sortedBaseNodes remove nodes not shared by faces of faceSubMesh + { + const SMDS_MeshNode* n1 = sortedBaseNodes.begin()->second; + const SMDS_MeshNode* n2 = sortedBaseNodes.rbegin()->second; + bool allNodesAreProxy = ( n1 != theProxyMesh->GetProxyNode( n1 ) && + n2 != theProxyMesh->GetProxyNode( n2 )); + if ( allNodesAreProxy ) + for ( u_n = sortedBaseNodes.begin(); u_n != sortedBaseNodes.end(); u_n++ ) + u_n->second = theProxyMesh->GetProxyNode( u_n->second ); + + if ( u_n = sortedBaseNodes.begin(), !isNodeInSubMesh( u_n->second, faceSubMesh )) + { + while ( ++u_n != sortedBaseNodes.end() && !isNodeInSubMesh( u_n->second, faceSubMesh )); + sortedBaseNodes.erase( sortedBaseNodes.begin(), u_n ); + } + else if ( u_n = --sortedBaseNodes.end(), !isNodeInSubMesh( u_n->second, faceSubMesh )) + { + while ( u_n != sortedBaseNodes.begin() && !isNodeInSubMesh( (--u_n)->second, faceSubMesh )); + sortedBaseNodes.erase( ++u_n, sortedBaseNodes.end() ); + } + if ( sortedBaseNodes.empty() ) continue; + } + double f, l; BRep_Tool::Range( *edge, f, l ); if ( edge->Orientation() == TopAbs_REVERSED ) std::swap( f, l ); const double coeff = 1. / ( l - f ) * length[iE] / fullLen; const double prevPar = theParam2ColumnMap.empty() ? 0 : theParam2ColumnMap.rbegin()->first; - map< double, const SMDS_MeshNode*>::iterator u_n = sortedBaseNodes.begin(); - for ( ; u_n != sortedBaseNodes.end(); u_n++ ) + for ( u_n = sortedBaseNodes.begin(); u_n != sortedBaseNodes.end(); u_n++ ) { double par = prevPar + coeff * ( u_n->first - f ); TParam2ColumnMap::iterator u2nn = @@ -1760,21 +1805,16 @@ bool SMESH_MesherHelper::LoadNodeColumns(TParam2ColumnMap & theParam2 u2nn->second.push_back( u_n->second ); } } - TParam2ColumnMap::iterator par_nVec_2, par_nVec_1 = theParam2ColumnMap.begin(); - if ( theProxyMesh ) - { - for ( ; par_nVec_1 != theParam2ColumnMap.end(); ++par_nVec_1 ) - { - const SMDS_MeshNode* & n = par_nVec_1->second[0]; - n = theProxyMesh->GetProxyNode( n ); - } - } + if ( theParam2ColumnMap.empty() ) + return false; + int nbRows = 1 + faceSubMesh->NbElements() / ( theParam2ColumnMap.size()-1 ); // fill theParam2ColumnMap column by column by passing from nodes on // theBaseEdge up via mesh faces on theFace + TParam2ColumnMap::iterator par_nVec_1, par_nVec_2; par_nVec_2 = theParam2ColumnMap.begin(); par_nVec_1 = par_nVec_2++; TIDSortedElemSet emptySet, avoidSet; @@ -1947,6 +1987,30 @@ TopoDS_Vertex SMESH_MesherHelper::IthVertex( const bool is2nd, return ( vIt.More() ? TopoDS::Vertex(vIt.Value()) : TopoDS_Vertex() ); } +//================================================================================ +/*! + * \brief Return type of shape contained in a group + * \param group - a shape of type TopAbs_COMPOUND + * \param avoidCompound - not to return TopAbs_COMPOUND + */ +//================================================================================ + +TopAbs_ShapeEnum SMESH_MesherHelper::GetGroupType(const TopoDS_Shape& group, + const bool avoidCompound) +{ + if ( !group.IsNull() ) + { + if ( group.ShapeType() != TopAbs_COMPOUND ) + return group.ShapeType(); + + // iterate on a compound + TopoDS_Iterator it( group ); + if ( it.More() ) + return avoidCompound ? GetGroupType( it.Value() ) : it.Value().ShapeType(); + } + return TopAbs_SHAPE; +} + //======================================================================= //function : IsQuadraticMesh //purpose : Check mesh without geometry for: if all elements on this shape are quadratic, @@ -1961,6 +2025,8 @@ SMESH_MesherHelper:: MType SMESH_MesherHelper::IsQuadraticMesh() int NbFacesAndEdges=0; //All faces and edges NbAllEdgsAndFaces = myMesh->NbEdges() + myMesh->NbFaces(); + if ( NbAllEdgsAndFaces == 0 ) + return SMESH_MesherHelper::LINEAR; //Quadratic faces and edges NbQuadFacesAndEdgs = myMesh->NbEdges(ORDER_QUADRATIC) + myMesh->NbFaces(ORDER_QUADRATIC); @@ -3021,18 +3087,349 @@ namespace { // Structures used by FixQuadraticElements() return _OK; } + + //================================================================================ + /*! + * \brief Place medium nodes at the link middle for elements whose corner nodes + * are out of geometrical boundary to prevent distorting elements. + * Issue 0020982, note 0013990 + */ + //================================================================================ + + void force3DOutOfBoundary( SMESH_MesherHelper& theHelper, + SMESH_ComputeErrorPtr& theError) + { + SMESHDS_Mesh* meshDS = theHelper.GetMeshDS(); + TopoDS_Shape shape = theHelper.GetSubShape().Oriented( TopAbs_FORWARD ); + if ( shape.IsNull() ) return; + + if ( !theError ) theError = SMESH_ComputeError::New(); + + gp_XYZ faceNorm; + + if ( shape.ShapeType() == TopAbs_FACE ) // 2D + { + if ( theHelper.GetMesh()->NbTriangles( ORDER_QUADRATIC ) < 1 ) return; + + SMESHDS_SubMesh* faceSM = meshDS->MeshElements( shape ); + if ( !faceSM ) return; + + const TopoDS_Face& face = TopoDS::Face( shape ); + Handle(Geom_Surface) surface = BRep_Tool::Surface( face ); + + TopExp_Explorer edgeIt( face, TopAbs_EDGE ); + for ( ; edgeIt.More(); edgeIt.Next() ) // loop on EDGEs of a FACE + { + // check if the EDGE needs checking + const TopoDS_Edge& edge = TopoDS::Edge( edgeIt.Current() ); + if ( BRep_Tool::Degenerated( edge ) ) + continue; + if ( theHelper.IsRealSeam( edge ) && + edge.Orientation() == TopAbs_REVERSED ) + continue; + + SMESHDS_SubMesh* edgeSM = meshDS->MeshElements( edge ); + if ( !edgeSM ) continue; + + double f,l; + Handle(Geom2d_Curve) pcurve = BRep_Tool::CurveOnSurface( edge, face, f, l ); + BRepAdaptor_Curve curve3D( edge ); + switch ( curve3D.GetType() ) { + case GeomAbs_Line: continue; + case GeomAbs_Circle: + case GeomAbs_Ellipse: + case GeomAbs_Hyperbola: + case GeomAbs_Parabola: + try + { + gp_Vec D1, D2, Du1, Dv1; gp_Pnt p; + curve3D.D2( 0.5 * ( f + l ), p, D1, D2 ); + gp_Pnt2d uv = pcurve->Value( 0.5 * ( f + l ) ); + surface->D1( uv.X(), uv.Y(), p, Du1, Dv1 ); + gp_Vec fNorm = Du1 ^ Dv1; + if ( fNorm.IsParallel( D2, M_PI * 25./180. )) + continue; // face is normal to the curve3D + + gp_Vec curvNorm = fNorm ^ D1; + if ( edge.Orientation() == TopAbs_REVERSED ) curvNorm.Reverse(); + if ( curvNorm * D2 > 0 ) + continue; // convex edge + } + catch ( Standard_Failure ) + { + continue; + } + } + // get nodes shared by faces that may be distorted + SMDS_NodeIteratorPtr nodeIt; + if ( edgeSM->NbNodes() > 0 ) { + nodeIt = edgeSM->GetNodes(); + } + else { + SMESHDS_SubMesh* vertexSM = meshDS->MeshElements( theHelper.IthVertex( 0, edge )); + if ( !vertexSM ) + vertexSM = meshDS->MeshElements( theHelper.IthVertex( 1, edge )); + if ( !vertexSM ) continue; + nodeIt = vertexSM->GetNodes(); + } + + // find suspicious faces + TIDSortedElemSet checkedFaces; + vector< const SMDS_MeshNode* > nOnEdge( 2 ); + const SMDS_MeshNode* nOnFace; + while ( nodeIt->more() ) + { + const SMDS_MeshNode* n = nodeIt->next(); + SMDS_ElemIteratorPtr faceIt = n->GetInverseElementIterator( SMDSAbs_Face ); + while ( faceIt->more() ) + { + const SMDS_MeshElement* f = faceIt->next(); + if ( !faceSM->Contains( f ) || + f->NbNodes() != 6 || // check quadratic triangles only + !checkedFaces.insert( f ).second ) + continue; + + // get nodes on EDGE and on FACE of a suspicious face + nOnEdge.clear(); nOnFace = 0; + SMDS_MeshElement::iterator triNode = f->begin_nodes(); + for ( int nbN = 0; nbN < 3; ++triNode, ++nbN ) + { + n = *triNode; + if ( n->GetPosition()->GetDim() == 2 ) + nOnFace = n; + else + nOnEdge.push_back( n ); + } + + // check if nOnFace is inside the FACE + if ( nOnFace && nOnEdge.size() == 2 ) + { + theHelper.AddTLinks( static_cast< const SMDS_MeshFace* > ( f )); + if ( !SMESH_Algo::FaceNormal( f, faceNorm, /*normalized=*/false )) + continue; + gp_XYZ edgeDir = SMESH_TNodeXYZ( nOnEdge[0] ) - SMESH_TNodeXYZ( nOnEdge[1] ); + gp_XYZ edgeNorm = faceNorm ^ edgeDir; + n = theHelper.GetMediumNode( nOnEdge[0], nOnEdge[1], true ); + gp_XYZ pN0 = SMESH_TNodeXYZ( nOnEdge[0] ); + gp_XYZ pMedium = SMESH_TNodeXYZ( n ); // on-edge node location + gp_XYZ pFaceN = SMESH_TNodeXYZ( nOnFace ); // on-face node location + double hMedium = edgeNorm * gp_Vec( pN0, pMedium ).XYZ(); + double hFace = edgeNorm * gp_Vec( pN0, pFaceN ).XYZ(); + if ( Abs( hMedium ) > Abs( hFace * 0.6 )) + { + // nOnFace is out of FACE, move a medium on-edge node to the middle + gp_XYZ pMid3D = 0.5 * ( pN0 + SMESH_TNodeXYZ( nOnEdge[1] )); + meshDS->MoveNode( n, pMid3D.X(), pMid3D.Y(), pMid3D.Z() ); + MSG( "move OUT of face " << n ); + theError->myBadElements.push_back( f ); + } + } + } + } + } + if ( !theError->myBadElements.empty() ) + theError->myName = EDITERR_NO_MEDIUM_ON_GEOM; + return; + + } // 2D ============================================================================== + + if ( shape.ShapeType() == TopAbs_SOLID ) // 3D + { + if ( theHelper.GetMesh()->NbTetras ( ORDER_QUADRATIC ) < 1 && + theHelper.GetMesh()->NbPyramids( ORDER_QUADRATIC ) < 1 ) return; + + SMESHDS_SubMesh* solidSM = meshDS->MeshElements( shape ); + if ( !solidSM ) return; + + // check if the SOLID is bound by concave FACEs + vector< TopoDS_Face > concaveFaces; + TopExp_Explorer faceIt( shape, TopAbs_FACE ); + for ( ; faceIt.More(); faceIt.Next() ) // loop on FACEs of a SOLID + { + const TopoDS_Face& face = TopoDS::Face( faceIt.Current() ); + if ( !meshDS->MeshElements( face )) continue; + + BRepAdaptor_Surface surface( face ); + switch ( surface.GetType() ) { + case GeomAbs_Plane: continue; + case GeomAbs_Cylinder: + case GeomAbs_Cone: + case GeomAbs_Sphere: + try + { + double u = 0.5 * ( surface.FirstUParameter() + surface.LastUParameter() ); + double v = 0.5 * ( surface.FirstVParameter() + surface.LastVParameter() ); + gp_Vec Du1, Dv1, Du2, Dv2, Duv2; gp_Pnt p; + surface.D2( u,v, p, Du1, Dv1, Du2, Dv2, Duv2 ); + gp_Vec fNorm = Du1 ^ Dv1; + if ( face.Orientation() == TopAbs_REVERSED ) fNorm.Reverse(); + bool concaveU = ( fNorm * Du2 > 1e-100 ); + bool concaveV = ( fNorm * Dv2 > 1e-100 ); + if ( concaveU || concaveV ) + concaveFaces.push_back( face ); + } + catch ( Standard_Failure ) + { + concaveFaces.push_back( face ); + } + } + } + if ( concaveFaces.empty() ) + return; + + // fix 2D mesh on the SOLID + for ( faceIt.ReInit(); faceIt.More(); faceIt.Next() ) // loop on FACEs of a SOLID + { + SMESH_MesherHelper faceHelper( *theHelper.GetMesh() ); + faceHelper.SetSubShape( faceIt.Current() ); + force3DOutOfBoundary( faceHelper, theError ); + } + + // get an iterator over faces on concaveFaces + vector< SMDS_ElemIteratorPtr > faceIterVec( concaveFaces.size() ); + for ( size_t i = 0; i < concaveFaces.size(); ++i ) + faceIterVec[i] = meshDS->MeshElements( concaveFaces[i] )->GetElements(); + typedef SMDS_IteratorOnIterators + < const SMDS_MeshElement*, vector< SMDS_ElemIteratorPtr > > TIterOnIter; + SMDS_ElemIteratorPtr faceIter( new TIterOnIter( faceIterVec )); + + // a seacher to check if a volume is close to a concave face + std::auto_ptr< SMESH_ElementSearcher > faceSearcher + ( SMESH_MeshEditor( theHelper.GetMesh() ).GetElementSearcher( faceIter )); + + // classifier + //BRepClass3d_SolidClassifier solidClassifier( shape ); + + TIDSortedElemSet checkedVols, movedNodes; + for ( faceIt.ReInit(); faceIt.More(); faceIt.Next() ) // loop on FACEs of a SOLID + { + const TopoDS_Shape& face = faceIt.Current(); + SMESHDS_SubMesh* faceSM = meshDS->MeshElements( face ); + if ( !faceSM ) continue; + + // get nodes shared by volumes (tet and pyra) on the FACE that may be distorted + SMDS_NodeIteratorPtr nodeIt; + if ( faceSM->NbNodes() > 0 ) { + nodeIt = faceSM->GetNodes(); + } + else { + TopExp_Explorer vertex( face, TopAbs_VERTEX ); + SMESHDS_SubMesh* vertexSM = meshDS->MeshElements( vertex.Current() ); + if ( !vertexSM ) continue; + nodeIt = vertexSM->GetNodes(); + } + + // find suspicious volumes adjacent to the FACE + vector< const SMDS_MeshNode* > nOnFace( 4 ); + const SMDS_MeshNode* nInSolid; + //vector< const SMDS_MeshElement* > intersectedFaces; + while ( nodeIt->more() ) + { + const SMDS_MeshNode* n = nodeIt->next(); + SMDS_ElemIteratorPtr volIt = n->GetInverseElementIterator( SMDSAbs_Volume ); + while ( volIt->more() ) + { + const SMDS_MeshElement* vol = volIt->next(); + int nbN = vol->NbCornerNodes(); + if ( ( nbN != 4 && nbN != 5 ) || + !solidSM->Contains( vol ) || + !checkedVols.insert( vol ).second ) + continue; + + // get nodes on FACE and in SOLID of a suspicious volume + nOnFace.clear(); nInSolid = 0; + SMDS_MeshElement::iterator volNode = vol->begin_nodes(); + for ( int nb = nbN; nb > 0; ++volNode, --nb ) + { + n = *volNode; + if ( n->GetPosition()->GetDim() == 3 ) + nInSolid = n; + else + nOnFace.push_back( n ); + } + if ( !nInSolid || nOnFace.size() != nbN - 1 ) + continue; + + // get size of the vol + SMESH_TNodeXYZ pInSolid( nInSolid ), pOnFace0( nOnFace[0] ); + double volLength = pInSolid.SquareDistance( nOnFace[0] ); + for ( size_t i = 1; i < nOnFace.size(); ++i ) + { + volLength = Max( volLength, pOnFace0.SquareDistance( nOnFace[i] )); + } + + // check if vol is close to concaveFaces + const SMDS_MeshElement* closeFace = + faceSearcher->FindClosestTo( pInSolid, SMDSAbs_Face ); + if ( !closeFace || + pInSolid.SquareDistance( closeFace->GetNode(0) ) > 4 * volLength ) + continue; + + // check if vol is distorted, i.e. a medium node is much closer + // to nInSolid than the link middle + bool isDistorted = false; + SMDS_FaceOfNodes onFaceTria( nOnFace[0], nOnFace[1], nOnFace[2] ); + if ( !SMESH_Algo::FaceNormal( &onFaceTria, faceNorm, /*normalized=*/false )) + continue; + theHelper.AddTLinks( static_cast< const SMDS_MeshVolume* > ( vol )); + vector< pair< SMESH_TLink, const SMDS_MeshNode* > > links; + for ( size_t i = 0; i < nOnFace.size(); ++i ) // loop on links between nOnFace + for ( size_t j = i+1; j < nOnFace.size(); ++j ) + { + SMESH_TLink link( nOnFace[i], nOnFace[j] ); + TLinkNodeMap::const_iterator linkIt = + theHelper.GetTLinkNodeMap().find( link ); + if ( linkIt != theHelper.GetTLinkNodeMap().end() ) + { + links.push_back( make_pair( linkIt->first, linkIt->second )); + if ( !isDistorted ) { + // compare projections of nInSolid and nMedium to face normal + gp_Pnt pMedium = SMESH_TNodeXYZ( linkIt->second ); + double hMedium = faceNorm * gp_Vec( pOnFace0, pMedium ).XYZ(); + double hVol = faceNorm * gp_Vec( pOnFace0, pInSolid ).XYZ(); + isDistorted = ( Abs( hMedium ) > Abs( hVol * 0.5 )); + } + } + } + // move medium nodes to link middle + if ( isDistorted ) + { + for ( size_t i = 0; i < links.size(); ++i ) + { + const SMDS_MeshNode* nMedium = links[i].second; + if ( movedNodes.insert( nMedium ).second ) + { + gp_Pnt pMid3D = 0.5 * ( SMESH_TNodeXYZ( links[i].first.node1() ) + + SMESH_TNodeXYZ( links[i].first.node2() )); + meshDS->MoveNode( nMedium, pMid3D.X(), pMid3D.Y(), pMid3D.Z() ); + MSG( "move OUT of solid " << nMedium ); + } + } + theError->myBadElements.push_back( vol ); + } + } // loop on volumes sharing a node on FACE + } // loop on nodes on FACE + } // loop on FACEs of a SOLID + + if ( !theError->myBadElements.empty() ) + theError->myName = EDITERR_NO_MEDIUM_ON_GEOM; + } // 3D case + } + } //namespace //======================================================================= /*! * \brief Move medium nodes of faces and volumes to fix distorted elements + * \param error - container of fixed distorted elements * \param volumeOnly - to fix nodes on faces or not, if the shape is solid * * Issue 0020307: EDF 992 SMESH : Linea/Quadratic with Medium Node on Geometry */ //======================================================================= -void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) +void SMESH_MesherHelper::FixQuadraticElements(SMESH_ComputeErrorPtr& error, + bool volumeOnly) { // setenv NO_FixQuadraticElements to know if FixQuadraticElements() is guilty of bad conversion if ( getenv("NO_FixQuadraticElements") ) @@ -3065,7 +3462,7 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) #endif SMESH_MesherHelper h(*myMesh); h.SetSubShape( s.Current() ); - h.FixQuadraticElements(false); + h.FixQuadraticElements( error, false ); } } // fix nodes on geom faces @@ -3076,10 +3473,13 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) MSG("FIX FACE " << nbfaces-- << " #" << GetMeshDS()->ShapeToIndex(fIt.Key())); SMESH_MesherHelper h(*myMesh); h.SetSubShape( fIt.Key() ); - h.FixQuadraticElements(true); + h.FixQuadraticElements( error, true); h.ToFixNodeParameters(true); } //perf_print_all_meters(1); + if ( error && error->myName == EDITERR_NO_MEDIUM_ON_GEOM ) + error->myComment = "during conversion to quadratic, " + "some medium nodes were not placed on geometry to avoid distorting elements"; return; } @@ -3118,6 +3518,11 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) TIDSortedNodeSet apexOfPyramid; const int apexIndex = 4; + // Issue 0020982 + // Move medium nodes to the link middle for elements whose corner nodes + // are out of geometrical boundary to fix distorted elements. + force3DOutOfBoundary( *this, error ); + if ( elemType == SMDSAbs_Volume ) { while ( elemIt->more() ) // loop on volumes @@ -3392,109 +3797,73 @@ void SMESH_MesherHelper::FixQuadraticElements(bool volumeOnly) } // loop on chains of links } // loop on 2 directions of propagation from quadrangle } // loop on faces - } + } // fix faces and/or volumes // 4. Move nodes // ------------- -// vector vols( 100 ); -// vector volSize( 100 ); -// int nbVols; -// bool ok; for ( pLink = links.begin(); pLink != links.end(); ++pLink ) { if ( pLink->IsMoved() ) { gp_Pnt p = pLink->MiddlePnt() + pLink->Move(); GetMeshDS()->MoveNode( pLink->_mediumNode, p.X(), p.Y(), p.Z()); - // -// gp_Pnt pNew = pLink->MiddlePnt() + pLink->Move(); -// if ( pLink->MediumPos() != SMDS_TOP_3DSPACE ) -// { -// // avoid making distorted volumes near boundary -// SMDS_ElemIteratorPtr volIt = -// (*pLink)._mediumNode->GetInverseElementIterator( SMDSAbs_Volume ); -// for ( nbVols = 0; volIt->more() && volTool.Set( volIt->next() ); ++nbVols ) -// { -// vols [ nbVols ] = volTool.Element(); -// volSize[ nbVols ] = volTool.GetSize(); -// } -// gp_Pnt pOld = pLink->MediumPnt(); -// const_cast( pLink->_mediumNode )->setXYZ( pNew.X(), pNew.Y(), pNew.Z() ); -// ok = true; -// while ( nbVols-- && ok ) -// { -// volTool.Set( vols[ nbVols ]); -// ok = ( volSize[ nbVols ] * volTool.GetSize() > 1e-20 ); -// } -// if ( !ok ) -// { -// const_cast( pLink->_mediumNode )->setXYZ( pOld.X(), pOld.Y(), pOld.Z() ); -// MSG( "Do NOT move \t" << pLink->_mediumNode->GetID() -// << " because of distortion of volume " << vols[ nbVols+1 ]->GetID()); -// continue; -// } -// } -// GetMeshDS()->MoveNode( pLink->_mediumNode, pNew.X(), pNew.Y(), pNew.Z() ); - } - } - - //return; - - // issue 0020982 - // Move the apex of pyramid together with the most curved link - - TIDSortedNodeSet::iterator apexIt = apexOfPyramid.begin(); - for ( ; apexIt != apexOfPyramid.end(); ++apexIt ) - { - SMESH_TNodeXYZ apex = *apexIt; - - gp_Vec maxMove( 0,0,0 ); - double maxMoveSize2 = 0; - - // shift of node index to get medium nodes between the base nodes - const int base2MediumShift = 5; - - // find maximal movement of medium node - SMDS_ElemIteratorPtr volIt = apex._node->GetInverseElementIterator( SMDSAbs_Volume ); - vector< const SMDS_MeshElement* > pyramids; - while ( volIt->more() ) - { - const SMDS_MeshElement* pyram = volIt->next(); - if ( pyram->GetEntityType() != SMDSEntity_Quad_Pyramid ) continue; - pyramids.push_back( pyram ); - - for ( int iBase = 0; iBase < apexIndex; ++iBase ) - { - SMESH_TNodeXYZ medium = pyram->GetNode( iBase + base2MediumShift ); - if ( medium._node->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE ) - { - SMESH_TNodeXYZ n1 = pyram->GetNode( iBase ); - SMESH_TNodeXYZ n2 = pyram->GetNode( ( iBase+1 ) % 4 ); - gp_Pnt middle = 0.5 * ( n1 + n2 ); - gp_Vec move( middle, medium ); - double moveSize2 = move.SquareMagnitude(); - if ( moveSize2 > maxMoveSize2 ) - maxMove = move, maxMoveSize2 = moveSize2; - } - } - } - - // move the apex - if ( maxMoveSize2 > 1e-20 ) - { - apex += maxMove.XYZ(); - GetMeshDS()->MoveNode( apex._node, apex.X(), apex.Y(), apex.Z()); - - // move medium nodes neighboring the apex to the middle - const int base2MediumShift_2 = 9; - for ( unsigned i = 0; i < pyramids.size(); ++i ) - for ( int iBase = 0; iBase < apexIndex; ++iBase ) - { - SMESH_TNodeXYZ base = pyramids[i]->GetNode( iBase ); - const SMDS_MeshNode* medium = pyramids[i]->GetNode( iBase + base2MediumShift_2 ); - gp_XYZ middle = 0.5 * ( apex + base ); - GetMeshDS()->MoveNode( medium, middle.X(), middle.Y(), middle.Z()); - } } } + + // Issue 0020982 + // Move the apex of pyramid together with the most curved link. + // TIDSortedNodeSet::iterator apexIt = apexOfPyramid.begin(); + // for ( ; apexIt != apexOfPyramid.end(); ++apexIt ) + // { + // SMESH_TNodeXYZ apex = *apexIt; + + // gp_Vec maxMove( 0,0,0 ); + // double maxMoveSize2 = 0; + + // // shift of node index to get medium nodes between the base nodes + // const int base2MediumShift = 5; + + // // find maximal movement of medium node + // SMDS_ElemIteratorPtr volIt = apex._node->GetInverseElementIterator( SMDSAbs_Volume ); + // vector< const SMDS_MeshElement* > pyramids; + // while ( volIt->more() ) + // { + // const SMDS_MeshElement* pyram = volIt->next(); + // if ( pyram->GetEntityType() != SMDSEntity_Quad_Pyramid ) continue; + // pyramids.push_back( pyram ); + + // for ( int iBase = 0; iBase < apexIndex; ++iBase ) + // { + // SMESH_TNodeXYZ medium = pyram->GetNode( iBase + base2MediumShift ); + // if ( medium._node->GetPosition()->GetTypeOfPosition() != SMDS_TOP_3DSPACE ) + // { + // SMESH_TNodeXYZ n1 = pyram->GetNode( iBase ); + // SMESH_TNodeXYZ n2 = pyram->GetNode( ( iBase+1 ) % 4 ); + // gp_Pnt middle = 0.5 * ( n1 + n2 ); + // gp_Vec move( middle, medium ); + // double moveSize2 = move.SquareMagnitude(); + // if ( moveSize2 > maxMoveSize2 ) + // maxMove = move, maxMoveSize2 = moveSize2; + // } + // } + // } + + // // move the apex + // if ( maxMoveSize2 > 1e-20 ) + // { + // apex += maxMove.XYZ(); + // GetMeshDS()->MoveNode( apex._node, apex.X(), apex.Y(), apex.Z()); + + // // move medium nodes neighboring the apex to the middle + // const int base2MediumShift_2 = 9; + // for ( unsigned i = 0; i < pyramids.size(); ++i ) + // for ( int iBase = 0; iBase < apexIndex; ++iBase ) + // { + // SMESH_TNodeXYZ base = pyramids[i]->GetNode( iBase ); + // const SMDS_MeshNode* medium = pyramids[i]->GetNode( iBase + base2MediumShift_2 ); + // gp_XYZ middle = 0.5 * ( apex + base ); + // GetMeshDS()->MoveNode( medium, middle.X(), middle.Y(), middle.Z()); + // } + // } + // } }