diff --git a/doc/salome/gui/SMESH/input/extrusion_along_path.doc b/doc/salome/gui/SMESH/input/extrusion_along_path.doc index 03717c77f..2253d08d8 100644 --- a/doc/salome/gui/SMESH/input/extrusion_along_path.doc +++ b/doc/salome/gui/SMESH/input/extrusion_along_path.doc @@ -153,7 +153,7 @@ The following dialog will appear: Linear variation of the angles option allows defining the angle of gradual rotation for the whole path. At each step the elements will -be rotated by angle / nb. of steps. +be rotated by ( angle / nb. of steps ). diff --git a/src/SMDS/SMDS_Mesh.cxx b/src/SMDS/SMDS_Mesh.cxx index f37887482..86524eb8d 100644 --- a/src/SMDS/SMDS_Mesh.cxx +++ b/src/SMDS/SMDS_Mesh.cxx @@ -3228,7 +3228,6 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement * elem, list& removedNodes, bool removenodes) { - //MESSAGE("SMDS_Mesh::RemoveElement " << elem->getVtkId() << " " << removenodes); // get finite elements built on elem set * s1; if ( (elem->GetType() == SMDSAbs_0DElement) @@ -3273,19 +3272,16 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement * elem, // Remove element from of its nodes SMDS_ElemIteratorPtr itn = (*it)->nodesIterator(); while (itn->more()) - { - SMDS_MeshNode * n = static_cast (const_cast (itn->next())); - n->RemoveInverseElement((*it)); - } + { + SMDS_MeshNode * n = static_cast (const_cast (itn->next())); + n->RemoveInverseElement((*it)); + } int IdToRemove = (*it)->GetID(); int vtkid = (*it)->getVtkId(); - //MESSAGE("elem Id to remove " << IdToRemove << " vtkid " << vtkid << - // " vtktype " << (*it)->GetVtkType() << " type " << (*it)->GetType()); switch ((*it)->GetType()) { case SMDSAbs_Node: - MYASSERT("Internal Error: This should not happen") - ; + MYASSERT("Internal Error: This should not happen"); break; case SMDSAbs_0DElement: if (IdToRemove >= 0) @@ -3305,8 +3301,10 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement * elem, } removedElems.push_back((*it)); myElementIDFactory->ReleaseID(IdToRemove, vtkid); - if (const SMDS_VtkEdge* vtkElem = dynamic_cast(*it)) + if (const SMDS_VtkEdge* vtkElem = dynamic_cast(*it)) { myEdgePool->destroy((SMDS_VtkEdge*) vtkElem); + ((SMDS_MeshElement*) *it)->init( -1, -1, -1 ); // avoid reuse + } else delete (*it); break; @@ -3318,8 +3316,10 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement * elem, } removedElems.push_back((*it)); myElementIDFactory->ReleaseID(IdToRemove, vtkid); - if (const SMDS_VtkFace* vtkElem = dynamic_cast(*it)) + if (const SMDS_VtkFace* vtkElem = dynamic_cast(*it)) { myFacePool->destroy((SMDS_VtkFace*) vtkElem); + ((SMDS_MeshElement*) *it)->init( -1, -1, -1 ); // avoid reuse + } else delete (*it); break; @@ -3331,8 +3331,10 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement * elem, } removedElems.push_back((*it)); myElementIDFactory->ReleaseID(IdToRemove, vtkid); - if (const SMDS_VtkVolume* vtkElem = dynamic_cast(*it)) + if (const SMDS_VtkVolume* vtkElem = dynamic_cast(*it)) { myVolumePool->destroy((SMDS_VtkVolume*) vtkElem); + ((SMDS_MeshElement*) *it)->init( -1, -1, -1 ); // avoid reuse + } else delete (*it); break; @@ -3344,18 +3346,19 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement * elem, } removedElems.push_back((*it)); myElementIDFactory->ReleaseID(IdToRemove, vtkid); - if (const SMDS_BallElement* vtkElem = dynamic_cast(*it)) + if (const SMDS_BallElement* vtkElem = dynamic_cast(*it)) { myBallPool->destroy(const_cast( vtkElem )); + ((SMDS_MeshElement*) *it)->init( -1, -1, -1 ); // avoid reuse + } else delete (*it); break; - case SMDSAbs_All: + case SMDSAbs_All: // avoid compilation warning case SMDSAbs_NbElementTypes: break; } if (vtkid >= 0) { - //MESSAGE("VTK_EMPTY_CELL in " << vtkid); this->myGrid->GetCellTypesArray()->SetValue(vtkid, VTK_EMPTY_CELL); } it++; @@ -3368,7 +3371,6 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement * elem, while (it != s2->end()) { int IdToRemove = (*it)->GetID(); - //MESSAGE( "SMDS: RM node " << IdToRemove); if (IdToRemove >= 0) { myNodes[IdToRemove] = 0; @@ -3430,11 +3432,12 @@ void SMDS_Mesh::RemoveFreeElement(const SMDS_MeshElement * elem) } // in meshes without descendants elements are always free - switch (aType) { + switch (aType) { case SMDSAbs_0DElement: myCells[elemId] = 0; myInfo.remove(elem); delete elem; + elem = 0; break; case SMDSAbs_Edge: myCells[elemId] = 0; @@ -3463,6 +3466,9 @@ void SMDS_Mesh::RemoveFreeElement(const SMDS_MeshElement * elem) this->myGrid->GetCellTypesArray()->SetValue(vtkId, VTK_EMPTY_CELL); // --- to do: keep vtkid in a list of reusable cells + + if ( elem ) + ((SMDS_MeshElement*) elem)->init( -1, -1, -1 ); // avoid reuse } } diff --git a/src/SMDS/SMDS_MeshNode.cxx b/src/SMDS/SMDS_MeshNode.cxx index 3524480cc..7b0f6a97c 100644 --- a/src/SMDS/SMDS_MeshNode.cxx +++ b/src/SMDS/SMDS_MeshNode.cxx @@ -182,17 +182,14 @@ public: MESSAGE("SMDS_MeshNode_MyInvIterator problem Null element"); throw SALOME_Exception("SMDS_MeshNode_MyInvIterator problem Null element"); } - //MESSAGE("vtkId " << vtkId << " smdsId " << smdsId << " " << elem->GetType()); iter++; return elem; } }; -SMDS_ElemIteratorPtr SMDS_MeshNode:: -GetInverseElementIterator(SMDSAbs_ElementType type) const +SMDS_ElemIteratorPtr SMDS_MeshNode::GetInverseElementIterator(SMDSAbs_ElementType type) const { vtkCellLinks::Link l = SMDS_Mesh::_meshList[myMeshId]->getGrid()->GetCellLinks()->GetLink(myVtkID); - //MESSAGE("myID " << myID << " ncells " << l.ncells); return SMDS_ElemIteratorPtr(new SMDS_MeshNode_MyInvIterator(SMDS_Mesh::_meshList[myMeshId], l.cells, l.ncells, type)); } diff --git a/src/SMESHGUI/SMESHGUI_RemoveElementsDlg.cxx b/src/SMESHGUI/SMESHGUI_RemoveElementsDlg.cxx index 1ebfd0622..e60f386b6 100644 --- a/src/SMESHGUI/SMESHGUI_RemoveElementsDlg.cxx +++ b/src/SMESHGUI/SMESHGUI_RemoveElementsDlg.cxx @@ -42,6 +42,7 @@ #include #include #include +#include #include #include @@ -225,7 +226,10 @@ void SMESHGUI_RemoveElementsDlg::ClickOnApply() if (mySMESHGUI->isActiveStudyLocked()) return; - if (myNbOkElements) { + if (myNbOkElements) + { + SUIT_OverrideCursor wc; + QStringList aListId = myEditCurrentArgument->text().split(" ", QString::SkipEmptyParts); SMESH::long_array_var anArrayOfIdeces = new SMESH::long_array; anArrayOfIdeces->length(aListId.count()); @@ -233,7 +237,8 @@ void SMESHGUI_RemoveElementsDlg::ClickOnApply() anArrayOfIdeces[i] = aListId[ i ].toInt(); bool aResult = false; - try { + try + { SMESH::SMESH_MeshEditor_var aMeshEditor = myMesh->GetMeshEditor(); aResult = aMeshEditor->RemoveElements(anArrayOfIdeces.in()); diff --git a/src/SMESHUtils/SMESH_MeshAlgos.cxx b/src/SMESHUtils/SMESH_MeshAlgos.cxx index c42f02023..a897decd6 100644 --- a/src/SMESHUtils/SMESH_MeshAlgos.cxx +++ b/src/SMESHUtils/SMESH_MeshAlgos.cxx @@ -458,6 +458,10 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher void GetElementsNearLine( const gp_Ax1& line, SMDSAbs_ElementType type, vector< const SMDS_MeshElement* >& foundElems); + void GetElementsInSphere( const gp_XYZ& center, + const double radius, + SMDSAbs_ElementType type, + vector< const SMDS_MeshElement* >& foundElems); double getTolerance(); bool getIntersParamOnLine(const gp_Lin& line, const SMDS_MeshElement* face, const double tolerance, double & param); @@ -466,6 +470,7 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher { return _outerFaces.empty() || _outerFaces.count(face); } + struct TInters //!< data of intersection of the line and the mesh face (used in GetPointState()) { const SMDS_MeshElement* _face; @@ -646,6 +651,7 @@ void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerF set< const SMDS_MeshElement*, TIDCompare >::const_iterator face = faces.begin(); for ( ; face != faces.end(); ++face ) { + if ( *face == outerFace ) continue; if ( !SMESH_MeshAlgos::FaceNormal( *face, fNorm, /*normalized=*/false )) continue; gp_Vec dirInF = gp_Vec( fNorm ) ^ n1n2; @@ -660,8 +666,8 @@ void SMESH_ElementSearcherImpl::findOuterBoundary(const SMDS_MeshElement* outerF // store the found outer face and add its links to continue seaching from if ( outerFace2 ) { - _outerFaces.insert( outerFace ); - int nbNodes = outerFace2->NbNodes()/( outerFace2->IsQuadratic() ? 2 : 1 ); + _outerFaces.insert( outerFace2 ); + int nbNodes = outerFace2->NbCornerNodes(); for ( int i = 0; i < nbNodes; ++i ) { SMESH_TLink link2( outerFace2->GetNode(i), outerFace2->GetNode((i+1)%nbNodes)); @@ -1078,6 +1084,27 @@ void SMESH_ElementSearcherImpl::GetElementsNearLine( const gp_Ax1& foundElems.assign( suspectFaces.begin(), suspectFaces.end()); } +//======================================================================= +/* + * Return elements whose bounding box intersects a sphere + */ +//======================================================================= + +void SMESH_ElementSearcherImpl::GetElementsInSphere( const gp_XYZ& center, + const double radius, + SMDSAbs_ElementType type, + vector< const SMDS_MeshElement* >& foundElems) +{ + if ( !_ebbTree || _elementType != type ) + { + if ( _ebbTree ) delete _ebbTree; + _ebbTree = new ElementBndBoxTree( *_mesh, _elementType = type, _meshPartIt ); + } + TIDSortedElemSet suspectFaces; // elements possibly intersecting the line + _ebbTree->getElementsInSphere( center, radius, suspectFaces ); + foundElems.assign( suspectFaces.begin(), suspectFaces.end() ); +} + //======================================================================= /*! * \brief Return true if the point is IN or ON of the element diff --git a/src/SMESHUtils/SMESH_MeshAlgos.hxx b/src/SMESHUtils/SMESH_MeshAlgos.hxx index 9b860a6ab..f8ea69a6b 100644 --- a/src/SMESHUtils/SMESH_MeshAlgos.hxx +++ b/src/SMESHUtils/SMESH_MeshAlgos.hxx @@ -88,6 +88,13 @@ struct SMESHUtils_EXPORT SMESH_ElementSearcher virtual void GetElementsNearLine( const gp_Ax1& line, SMDSAbs_ElementType type, std::vector< const SMDS_MeshElement* >& foundElems) = 0; + /*! + * \brief Return elements whose bounding box intersects a sphere + */ + virtual void GetElementsInSphere( const gp_XYZ& center, + const double radius, + SMDSAbs_ElementType type, + std::vector< const SMDS_MeshElement* >& foundElems) = 0; /*! * \brief Find out if the given point is out of closed 2D mesh. */ diff --git a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx index f4c519ae7..a286e0e48 100644 --- a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx +++ b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx @@ -24,6 +24,7 @@ #include "StdMeshers_QuadToTriaAdaptor.hxx" +#include "SMDS_IteratorOnIterators.hxx" #include "SMDS_SetIterator.hxx" #include "SMESHDS_GroupBase.hxx" #include "SMESH_Algo.hxx" @@ -115,20 +116,20 @@ namespace gp_Vec nJ = baseVec.Crossed( baJ ); // Check angle between normals - double angle = nI.Angle( nJ ); + double angle = nI.Angle( nJ ); bool tooClose = ( angle < 15. * M_PI / 180. ); // Check if pyramids collide if ( !tooClose && ( baI * baJ > 0 ) && ( nI * nJ > 0 )) { // find out if nI points outside of PrmI or inside - int dInd = baseNodesIndI[1] - baseNodesIndI[0]; + int dInd = baseNodesIndI[1] - baseNodesIndI[0]; bool isOutI = ( abs(dInd)==1 ) ? dInd < 0 : dInd > 0; - // find out sign of projection of nJ to baI + // find out sign of projection of baI to nJ double proj = baI * nJ; - tooClose = isOutI ? proj > 0 : proj < 0; + tooClose = ( isOutI ? proj > 0 : proj < 0 ); } // Check if PrmI and PrmJ are in same domain @@ -170,8 +171,9 @@ namespace continue; // f is a base quadrangle // check projections of face direction (baOFN) to triange normals (nI and nJ) - gp_Vec baOFN( base1, SMESH_TNodeXYZ( otherFaceNode )); - if ( nI * baOFN > 0 && nJ * baOFN > 0 ) + gp_Vec baOFN( base2, SMESH_TNodeXYZ( otherFaceNode )); + if ( nI * baOFN > 0 && nJ * baOFN > 0 && + baI* baOFN > 0 && baJ* baOFN > 0 ) // issue 0023212 { tooClose = false; // f is between pyramids break; @@ -253,7 +255,6 @@ namespace } } } - } //================================================================================ @@ -266,6 +267,8 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr const SMDS_MeshElement* PrmJ, set & nodesToMove) { + // cout << endl << "Merge " << PrmI->GetID() << " " << PrmJ->GetID() << " " + // << PrmI->GetNode(4) << PrmJ->GetNode(4) << endl; const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove //int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume ); SMESH_TNodeXYZ Pj( Nrem ); @@ -288,7 +291,7 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr vector< const SMDS_MeshElement* > inverseElems // copy inverse elements to avoid iteration on changing container ( TStdElemIterator( CommonNode->GetInverseElementIterator(SMDSAbs_Face)), itEnd); - for ( unsigned i = 0; i < inverseElems.size(); ++i ) + for ( size_t i = 0; i < inverseElems.size(); ++i ) { const SMDS_MeshElement* FI = inverseElems[i]; const SMDS_MeshElement* FJEqual = 0; @@ -309,11 +312,12 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr } // set the common apex node to pyramids and triangles merged with J + vector< const SMDS_MeshNode* > nodes; inverseElems.assign( TStdElemIterator( Nrem->GetInverseElementIterator()), itEnd ); - for ( unsigned i = 0; i < inverseElems.size(); ++i ) + for ( size_t i = 0; i < inverseElems.size(); ++i ) { const SMDS_MeshElement* elem = inverseElems[i]; - vector< const SMDS_MeshNode* > nodes( elem->begin_nodes(), elem->end_nodes() ); + nodes.assign( elem->begin_nodes(), elem->end_nodes() ); nodes[ elem->GetType() == SMDSAbs_Volume ? PYRAM_APEX : TRIA_APEX ] = CommonNode; GetMeshDS()->ChangeElementNodes( elem, &nodes[0], nodes.size()); } @@ -330,33 +334,34 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr //================================================================================ void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement* PrmI, - set& nodesToMove) + set& nodesToMove, + const bool isRecursion) { TIDSortedElemSet adjacentPyrams; bool mergedPyrams = false; - for(int k=0; k<4; k++) // loop on 4 base nodes of PrmI + for ( int k=0; k<4; k++ ) // loop on 4 base nodes of PrmI { - const SMDS_MeshNode* n = PrmI->GetNode(k); + const SMDS_MeshNode* n = PrmI->GetNode(k); SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); while ( vIt->more() ) { const SMDS_MeshElement* PrmJ = vIt->next(); - if ( PrmJ->NbCornerNodes() != 5 || !adjacentPyrams.insert( PrmJ ).second ) + if ( PrmJ == PrmI || PrmJ->NbCornerNodes() != 5 || !adjacentPyrams.insert( PrmJ ).second ) continue; - if ( PrmI != PrmJ && TooCloseAdjacent( PrmI, PrmJ, GetMesh()->HasShapeToMesh() )) + if ( TooCloseAdjacent( PrmI, PrmJ, GetMesh()->HasShapeToMesh() )) { MergePiramids( PrmI, PrmJ, nodesToMove ); mergedPyrams = true; // container of inverse elements can change - vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); + // vIt = n->GetInverseElementIterator( SMDSAbs_Volume ); -- iterator re-implemented } } } - if ( mergedPyrams ) + if ( mergedPyrams && !isRecursion ) { TIDSortedElemSet::iterator prm; for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++prm) - MergeAdjacent( *prm, nodesToMove ); + MergeAdjacent( *prm, nodesToMove, true ); } } @@ -414,79 +419,58 @@ static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2, //======================================================================= //function : HasIntersection3 -//purpose : Auxilare for HasIntersection() -// find intersection point between triangle (P1,P2,P3) -// and segment [PC,P] +//purpose : Find intersection point between a triangle (P1,P2,P3) +// and a segment [PC,P] //======================================================================= static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint, const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3) { - //cout<<"HasIntersection3"< preci ) || - ( (PC.Y()-PIn.Y())*(P.Y()-PIn.Y()) > preci ) || - ( (PC.Z()-PIn.Z())*(P.Z()-PIn.Z()) > preci ); - if(IsExternal) { - return false; - } - // check if this point is internal for triangle (P1,P2,P3) - gp_Vec V1(PIn,P1); - gp_Vec V2(PIn,P2); - gp_Vec V3(PIn,P3); - if( V1.Magnitude() -EPSILON && det < EPSILON) + return false; + + /* calculate U parameter and test bounds */ + double u = ( tvec * pvec ) / det; + //if (u < 0.0 || u > 1.0) + if (u < -EPSILON || u > 1.0 + EPSILON) + return false; + + /* prepare to test V parameter */ + gp_XYZ qvec = tvec ^ edge1; + + /* calculate V parameter and test bounds */ + double v = (dir * qvec) / det; + //if ( v < 0.0 || u + v > 1.0 ) + if ( v < -EPSILON || u + v > 1.0 + EPSILON) + return false; + + /* calculate t, ray intersects triangle */ + double t = (edge2 * qvec) / det; + + Pint = orig + dir * t; + + return ( t > 0. && t < segLen ); } //======================================================================= @@ -521,54 +505,99 @@ static bool HasIntersection(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint, //================================================================================ /*! - * \brief Checks if a line segment (P,PC) intersects any mesh face. - * \param P - first segment end - * \param PC - second segment end (it is a gravity center of quadrangle) - * \param Pint - (out) intersection point + * \brief Return allowed height of a pyramid + * \param Papex - optimal pyramid apex + * \param PC - gravity center of a quadrangle + * \param PN - four nodes of the quadrangle * \param aMesh - mesh - * \param aShape - shape to check faces on - * \param NotCheckedFace - mesh face not to check - * \retval bool - true if there is an intersection + * \param NotCheckedFace - the quadrangle face + * \retval double - pyramid height */ //================================================================================ -bool StdMeshers_QuadToTriaAdaptor::CheckIntersection (const gp_Pnt& P, - const gp_Pnt& PC, - gp_Pnt& Pint, - SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape, - const SMDS_MeshElement* NotCheckedFace) +void StdMeshers_QuadToTriaAdaptor::LimitHeight (gp_Pnt& Papex, + const gp_Pnt& PC, + const TColgp_Array1OfPnt& PN, + const vector& FNodes, + SMESH_Mesh& aMesh, + const SMDS_MeshElement* NotCheckedFace, + const bool UseApexRay) { if ( !myElemSearcher ) myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *aMesh.GetMeshDS() ); SMESH_ElementSearcher* searcher = const_cast(myElemSearcher); - bool res = false; - double dist = RealLast(); // find intersection closest to PC - gp_Pnt Pres; - - gp_Ax1 line( P, gp_Vec(P,PC)); - vector< const SMDS_MeshElement* > suspectElems; - searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems); + // Find intersection of faces with (P,PC) segment elongated 3 times + double height = Papex.Distance( PC ); + gp_Ax1 line( PC, gp_Vec( PC, Papex )); + gp_Pnt Pint, Ptest; + vector< const SMDS_MeshElement* > suspectFaces; TColgp_SequenceOfPnt aContour; - for ( size_t iF = 0; iF < suspectElems.size(); ++iF ) + + if ( UseApexRay ) { - const SMDS_MeshElement* face = suspectElems[iF]; - if ( face == NotCheckedFace ) continue; - aContour.Clear(); - for ( int i = 0; i < face->NbCornerNodes(); ++i ) - aContour.Append( SMESH_TNodeXYZ( face->GetNode(i) )); - if ( HasIntersection(P, PC, Pres, aContour)) { - res = true; - double tmp = PC.Distance(Pres); - if ( tmp < dist ) { - Pint = Pres; - dist = tmp; + // find intersection closest to PC + Ptest = PC.XYZ() + line.Direction().XYZ() * height * 3; + + searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectFaces ); + for ( size_t iF = 0; iF < suspectFaces.size(); ++iF ) + { + const SMDS_MeshElement* face = suspectFaces[iF]; + if ( face == NotCheckedFace ) continue; + + aContour.Clear(); + for ( int i = 0, nb = face->NbCornerNodes(); i < nb; ++i ) + aContour.Append( SMESH_TNodeXYZ( face->GetNode(i) )); + + if ( HasIntersection( Ptest, PC, Pint, aContour )) + { + double dInt = PC.Distance( Pint ); + height = Min( height, dInt / 3. ); } } } - return res; + + // Find faces intersecting triangular facets of the pyramid (issue 23212) + + gp_XYZ center = PC.XYZ() + line.Direction().XYZ() * height * 0.5; + double diameter = Max( PN(1).Distance(PN(3)), PN(2).Distance(PN(4))); + suspectFaces.clear(); + searcher->GetElementsInSphere( center, diameter * 0.6, SMDSAbs_Face, suspectFaces); + + const double upShift = 1.5; + Ptest = PC.XYZ() + line.Direction().XYZ() * height * upShift; // tmp apex + + for ( size_t iF = 0; iF < suspectFaces.size(); ++iF ) + { + const SMDS_MeshElement* face = suspectFaces[iF]; + if ( face == NotCheckedFace ) continue; + if ( face->GetNodeIndex( FNodes[0] ) >= 0 || + face->GetNodeIndex( FNodes[1] ) >= 0 || + face->GetNodeIndex( FNodes[2] ) >= 0 || + face->GetNodeIndex( FNodes[3] ) >= 0 ) + continue; // neighbor face of the quadrangle + + // limit height using points of intersection of face links with pyramid facets + int nbN = face->NbCornerNodes(); + gp_Pnt P1 = SMESH_TNodeXYZ( face->GetNode( nbN-1 )); // 1st link end + for ( int i = 0; i < nbN; ++i ) + { + gp_Pnt P2 = SMESH_TNodeXYZ( face->GetNode(i) ); // 2nd link end + + for ( int iN = 1; iN <= 4; ++iN ) // loop on pyramid facets + { + if ( HasIntersection3( P1, P2, Pint, PN(iN), PN(iN+1), Ptest )) + { + height = Min( height, gp_Vec( PC, Pint ) * line.Direction() ); + //Ptest = PC.XYZ() + line.Direction().XYZ() * height * upShift; // new tmp apex + } + } + P1 = P2; + } + } + + Papex = PC.XYZ() + line.Direction().XYZ() * height; } //================================================================================ @@ -720,25 +749,36 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, vector myPyramids; + const SMESHDS_SubMesh * aSubMeshDSFace; SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); SMESH_MesherHelper helper(aMesh); helper.IsQuadraticSubMesh(aShape); helper.SetElementsOnShape( true ); if ( myElemSearcher ) delete myElemSearcher; + vector< SMDS_ElemIteratorPtr > itVec; if ( aProxyMesh ) - myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, aProxyMesh->GetFaces(aShape)); + { + itVec.push_back( aProxyMesh->GetFaces( aShape )); + } else - myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS ); + { + for ( TopExp_Explorer exp(aShape,TopAbs_FACE); exp.More(); exp.Next() ) + if (( aSubMeshDSFace = aProxyMesh->GetSubMesh( exp.Current() ))) + itVec.push_back( aSubMeshDSFace->GetElements() ); + } + typedef + SMDS_IteratorOnIterators< const SMDS_MeshElement*, vector< SMDS_ElemIteratorPtr > > TIter; + SMDS_ElemIteratorPtr faceIt( new TIter( itVec )); + myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, faceIt ); - const SMESHDS_SubMesh * aSubMeshDSFace; TColgp_Array1OfPnt PN(1,5); TColgp_Array1OfVec VN(1,4); vector FNodes(5); gp_Pnt PC; gp_Vec VNorm; - for (TopExp_Explorer exp(aShape,TopAbs_FACE);exp.More();exp.Next()) + for ( TopExp_Explorer exp(aShape,TopAbs_FACE); exp.More(); exp.Next() ) { const TopoDS_Shape& aShapeFace = exp.Current(); if ( aProxyMesh ) @@ -809,17 +849,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh, } else { // check possible intersection with other faces - gp_Pnt Pint; - gp_Vec VB(PC,PCbest); - gp_Pnt PCbestTmp = PC.XYZ() + VB.XYZ() * 3.0; - bool hasInters = CheckIntersection(PCbestTmp, PC, Pint, aMesh, aShape, face); - if ( hasInters ) { - double dist = PC.Distance(Pint)/3.; - if ( dist < height ) { - gp_Dir aDir( VB ); - PCbest = PC.XYZ() + aDir.XYZ() * dist; - } - } + LimitHeight( PCbest, PC, PN, FNodes, aMesh, face, /*UseApexRay=*/true ); } // create node for PCbest SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() ); @@ -971,11 +1001,11 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) gp_Pnt Pres1,Pres2; gp_Ax1 line( PC, VNorm ); - vector< const SMDS_MeshElement* > suspectElems; - searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems); + vector< const SMDS_MeshElement* > suspectFaces; + searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectFaces); - for ( size_t iF = 0; iF < suspectElems.size(); ++iF ) { - const SMDS_MeshElement* F = suspectElems[iF]; + for ( size_t iF = 0; iF < suspectFaces.size(); ++iF ) { + const SMDS_MeshElement* F = suspectFaces[iF]; if ( F == face ) continue; aContour.Clear(); for ( int i = 0; i < 4; ++i ) @@ -1026,7 +1056,9 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) continue; } + // ----------------------------------- // Case of non-degenerated quadrangle + // ----------------------------------- // Find pyramid peak @@ -1059,12 +1091,12 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) gp_Pnt intPnt[2]; gp_Ax1 line( PC, tmpDir ); - vector< const SMDS_MeshElement* > suspectElems; - searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems); + vector< const SMDS_MeshElement* > suspectFaces; + searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectFaces); - for ( size_t iF = 0; iF < suspectElems.size(); ++iF ) + for ( size_t iF = 0; iF < suspectFaces.size(); ++iF ) { - const SMDS_MeshElement* F = suspectElems[iF]; + const SMDS_MeshElement* F = suspectFaces[iF]; if ( F == face ) continue; aContour.Clear(); int nbN = F->NbCornerNodes(); @@ -1109,13 +1141,15 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) { if( !intersected[isRev] ) continue; double pyramidH = Min( height, PC.Distance(intPnt[isRev])/3.); - PCbest = PC.XYZ() + tmpDir.XYZ() * (isRev ? -pyramidH : pyramidH); + gp_Pnt Papex = PC.XYZ() + tmpDir.XYZ() * (isRev ? -pyramidH : pyramidH); - // create node for PCbest - SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() ); + LimitHeight( Papex, PC, PN, FNodes, aMesh, face, /*UseApexRay=*/false ); + + // create node for Papex + SMDS_MeshNode* NewNode = helper.AddNode( Papex.X(), Papex.Y(), Papex.Z() ); // add triangles to result map - for(i=0; i<4; i++) { + for ( i = 0; i < 4; i++) { SMDS_MeshFace* NewFace; if(isRev) NewFace = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] ); @@ -1146,16 +1180,19 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh) bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh, const vector& myPyramids) { - if(myPyramids.empty()) + if ( myPyramids.empty() ) return true; SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); size_t i, j, k; - int myShapeID = myPyramids[0]->GetNode(4)->getshapeId(); - - if ( myElemSearcher ) delete myElemSearcher; - myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS ); - SMESH_ElementSearcher* searcher = const_cast(myElemSearcher); + //int myShapeID = myPyramids[0]->GetNode(4)->getshapeId(); + { + SMDS_ElemIteratorPtr + pyramIt( new SMDS_ElementVectorIterator( myPyramids.begin(), myPyramids.end() )); + if ( myElemSearcher ) delete myElemSearcher; + myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, pyramIt ); + } + SMESH_ElementSearcher* searcher = const_cast( myElemSearcher ); set nodesToMove; @@ -1167,17 +1204,19 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& MergeAdjacent( PrmI, nodesToMove ); } - // iterate on all pyramids + // iterate on all new pyramids + vector< const SMDS_MeshElement* > suspectPyrams; for ( i = 0; i < myPyramids.size(); ++i ) { - const SMDS_MeshElement* PrmI = myPyramids[i]; + const SMDS_MeshElement* PrmI = myPyramids[i]; + const SMDS_MeshNode* apexI = PrmI->GetNode( PYRAM_APEX ); // compare PrmI with all the rest pyramids // collect adjacent pyramids and nodes coordinates of PrmI set checkedPyrams; - vector PsI(5); - for(k=0; k<5; k++) // loop on 4 base nodes of PrmI + gp_Pnt PsI[5]; + for ( k = 0; k < 5; k++ ) { const SMDS_MeshNode* n = PrmI->GetNode(k); PsI[k] = SMESH_TNodeXYZ( n ); @@ -1190,70 +1229,77 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& } } + // get pyramids to check + gp_XYZ PC = ( PsI[0].XYZ() + PsI[1].XYZ() + PsI[2].XYZ() + PsI[3].XYZ() ) / 4.; + gp_XYZ ray = PsI[4].XYZ() - PC; + gp_XYZ center = PC + 0.5 * ray; + double diameter = Max( PsI[0].Distance(PsI[2]), PsI[1].Distance(PsI[3])); + suspectPyrams.clear(); + searcher->GetElementsInSphere( center, diameter * 0.6, SMDSAbs_Volume, suspectPyrams); + // check intersection with distant pyramids - for(k=0; k<4; k++) // loop on 4 base nodes of PrmI + for ( j = 0; j < suspectPyrams.size(); ++j ) { - gp_Vec Vtmp(PsI[k],PsI[4]); - gp_Ax1 line( PsI[k], Vtmp ); - vector< const SMDS_MeshElement* > suspectPyrams; - searcher->GetElementsNearLine( line, SMDSAbs_Volume, suspectPyrams); + const SMDS_MeshElement* PrmJ = suspectPyrams[j]; + if ( PrmJ == PrmI ) + continue; + if ( apexI == PrmJ->GetNode( PYRAM_APEX )) + continue; // pyramids PrmI and PrmJ already merged + if ( !checkedPyrams.insert( PrmJ ).second ) + continue; // already checked - for ( j = 0; j < suspectPyrams.size(); ++j ) + gp_Pnt PsJ[5]; + for ( k = 0; k < 5; k++ ) + PsJ[k] = SMESH_TNodeXYZ( PrmJ->GetNode(k) ); + + if ( ray * ( PsJ[4].XYZ() - PC ) < 0. ) + continue; // PrmJ is below PrmI + + for ( k = 0; k < 4; k++ ) // loop on 4 base nodes of PrmI { - const SMDS_MeshElement* PrmJ = suspectPyrams[j]; - if ( PrmJ == PrmI || PrmJ->NbCornerNodes() != 5 ) - continue; - if ( myShapeID != PrmJ->GetNode(4)->getshapeId()) - continue; // pyramid from other SOLID - if ( PrmI->GetNode(4) == PrmJ->GetNode(4) ) - continue; // pyramids PrmI and PrmJ already merged - if ( !checkedPyrams.insert( PrmJ ).second ) - continue; // already checked - - TXyzIterator xyzIt( PrmJ->nodesIterator() ); - vector PsJ( xyzIt, TXyzIterator() ); - gp_Pnt Pint; bool hasInt=false; - for(k=0; k<4 && !hasInt; k++) { - gp_Vec Vtmp(PsI[k],PsI[4]); + for ( k = 0; k < 4 && !hasInt; k++ ) + { + gp_Vec Vtmp( PsI[k], PsI[ PYRAM_APEX ]); gp_Pnt Pshift = PsI[k].XYZ() + Vtmp.XYZ() * 0.01; // base node moved a bit to apex hasInt = - ( HasIntersection3( Pshift, PsI[4], Pint, PsJ[0], PsJ[1], PsJ[4]) || - HasIntersection3( Pshift, PsI[4], Pint, PsJ[1], PsJ[2], PsJ[4]) || - HasIntersection3( Pshift, PsI[4], Pint, PsJ[2], PsJ[3], PsJ[4]) || - HasIntersection3( Pshift, PsI[4], Pint, PsJ[3], PsJ[0], PsJ[4]) ); + ( HasIntersection3( Pshift, PsI[4], Pint, PsJ[0], PsJ[1], PsJ[PYRAM_APEX]) || + HasIntersection3( Pshift, PsI[4], Pint, PsJ[1], PsJ[2], PsJ[PYRAM_APEX]) || + HasIntersection3( Pshift, PsI[4], Pint, PsJ[2], PsJ[3], PsJ[PYRAM_APEX]) || + HasIntersection3( Pshift, PsI[4], Pint, PsJ[3], PsJ[0], PsJ[PYRAM_APEX]) ); } - for(k=0; k<4 && !hasInt; k++) { - gp_Vec Vtmp(PsJ[k],PsJ[4]); + for ( k = 0; k < 4 && !hasInt; k++ ) + { + gp_Vec Vtmp( PsJ[k], PsJ[ PYRAM_APEX ]); gp_Pnt Pshift = PsJ[k].XYZ() + Vtmp.XYZ() * 0.01; hasInt = - ( HasIntersection3( Pshift, PsJ[4], Pint, PsI[0], PsI[1], PsI[4]) || - HasIntersection3( Pshift, PsJ[4], Pint, PsI[1], PsI[2], PsI[4]) || - HasIntersection3( Pshift, PsJ[4], Pint, PsI[2], PsI[3], PsI[4]) || - HasIntersection3( Pshift, PsJ[4], Pint, PsI[3], PsI[0], PsI[4]) ); + ( HasIntersection3( Pshift, PsJ[4], Pint, PsI[0], PsI[1], PsI[PYRAM_APEX]) || + HasIntersection3( Pshift, PsJ[4], Pint, PsI[1], PsI[2], PsI[PYRAM_APEX]) || + HasIntersection3( Pshift, PsJ[4], Pint, PsI[2], PsI[3], PsI[PYRAM_APEX]) || + HasIntersection3( Pshift, PsJ[4], Pint, PsI[3], PsI[0], PsI[PYRAM_APEX]) ); } if ( hasInt ) { // count common nodes of base faces of two pyramids int nbc = 0; - for (k=0; k<4; k++) + for ( k = 0; k < 4; k++ ) nbc += int ( PrmI->GetNodeIndex( PrmJ->GetNode(k) ) >= 0 ); if ( nbc == 4 ) continue; // pyrams have a common base face - if(nbc>0) + if ( nbc > 0 ) { // Merge the two pyramids and others already merged with them MergePiramids( PrmI, PrmJ, nodesToMove ); } - else { // nbc==0 - + else // nbc==0 + { // decrease height of pyramids gp_XYZ PCi(0,0,0), PCj(0,0,0); - for(k=0; k<4; k++) { + for ( k = 0; k < 4; k++ ) { PCi += PsI[k].XYZ(); PCj += PsJ[k].XYZ(); } @@ -1272,9 +1318,9 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& VN1.Scale(coef1); VN2.Scale(coef2); - SMDS_MeshNode* aNode1 = const_cast(PrmI->GetNode(4)); + SMDS_MeshNode* aNode1 = const_cast( apexI ); aNode1->setXYZ( PCi.X()+VN1.X(), PCi.Y()+VN1.Y(), PCi.Z()+VN1.Z() ); - SMDS_MeshNode* aNode2 = const_cast(PrmJ->GetNode(4)); + SMDS_MeshNode* aNode2 = const_cast(PrmJ->GetNode( PYRAM_APEX )); aNode2->setXYZ( PCj.X()+VN2.X(), PCj.Y()+VN2.Y(), PCj.Z()+VN2.Z() ); nodesToMove.insert( aNode1 ); nodesToMove.insert( aNode2 ); diff --git a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.hxx b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.hxx index 44470b8e9..3c2807a14 100644 --- a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.hxx +++ b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.hxx @@ -71,10 +71,13 @@ protected: gp_Pnt& PC, gp_Vec& VNorm, const SMDS_MeshElement** volumes=0); - bool CheckIntersection(const gp_Pnt& P, const gp_Pnt& PC, - gp_Pnt& Pint, SMESH_Mesh& aMesh, - const TopoDS_Shape& aShape, - const SMDS_MeshElement* NotCheckedFace); + void LimitHeight (gp_Pnt& Papex, + const gp_Pnt& PC, + const TColgp_Array1OfPnt& PN, + const std::vector& FNodes, + SMESH_Mesh& aMesh, + const SMDS_MeshElement* NotCheckedFace, + const bool UseApexRay); bool Compute2ndPart(SMESH_Mesh& aMesh, const std::vector& pyramids); @@ -85,7 +88,8 @@ protected: std::set & nodesToMove); void MergeAdjacent(const SMDS_MeshElement* PrmI, - std::set& nodesToMove); + std::set& nodesToMove, + const bool isRecursion = false); TopoDS_Shape myShape;