ForkMaster/smesh
1
0
Fork 0
mirror of https://git.salome-platform.org/gitpub/modules/smesh.git synced 2025-01-26 23:10:34 +05:00
Code Issues Packages Projects Releases Wiki Activity

0023212: EDF 12054 - Problem with a pyramidal layer

Browse Source
This commit is contained in:
eap 2016-01-20 16:02:17 +03:00
parent d1ed915c68
commit 3c33b14157
8 changed files with 316 additions and 224 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
doc/salome/gui/SMESH/input/extrusion_along_path.doc
View File

@ -153,7 +153,7 @@ The following dialog will appear:
<b>Linear variation of the angles</b> option allows defining the angle <b>Linear variation of the angles</b> option allows defining the angle
of gradual rotation for the whole path. At each step the elements will of gradual rotation for the whole path. At each step the elements will
be rotated by <code>angle / nb. of steps</code>. be rotated by <code>( angle / nb. of steps )</code>.
</li> </li>
</ul> </ul>
</li> </li>

40
src/SMDS/SMDS_Mesh.cxx
View File

@ -3228,7 +3228,6 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement * elem,
list<const SMDS_MeshElement *>& removedNodes, list<const SMDS_MeshElement *>& removedNodes,
bool removenodes) bool removenodes)
{ {
//MESSAGE("SMDS_Mesh::RemoveElement " << elem->getVtkId() << " " << removenodes);
// get finite elements built on elem // get finite elements built on elem
set<const SMDS_MeshElement*> * s1; set<const SMDS_MeshElement*> * s1;
if ( (elem->GetType() == SMDSAbs_0DElement) if ( (elem->GetType() == SMDSAbs_0DElement)
@ -3273,19 +3272,16 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement * elem,
// Remove element from <InverseElements> of its nodes // Remove element from <InverseElements> of its nodes
SMDS_ElemIteratorPtr itn = (*it)->nodesIterator(); SMDS_ElemIteratorPtr itn = (*it)->nodesIterator();
while (itn->more()) while (itn->more())
{ {
SMDS_MeshNode * n = static_cast<SMDS_MeshNode *> (const_cast<SMDS_MeshElement *> (itn->next())); SMDS_MeshNode * n = static_cast<SMDS_MeshNode *> (const_cast<SMDS_MeshElement *> (itn->next()));
n->RemoveInverseElement((*it)); n->RemoveInverseElement((*it));
} }
int IdToRemove = (*it)->GetID(); int IdToRemove = (*it)->GetID();
int vtkid = (*it)->getVtkId(); int vtkid = (*it)->getVtkId();
//MESSAGE("elem Id to remove " << IdToRemove << " vtkid " << vtkid <<
// " vtktype " << (*it)->GetVtkType() << " type " << (*it)->GetType());
switch ((*it)->GetType()) switch ((*it)->GetType())
{ {
case SMDSAbs_Node: case SMDSAbs_Node:
MYASSERT("Internal Error: This should not happen") MYASSERT("Internal Error: This should not happen");
;
break; break;
case SMDSAbs_0DElement: case SMDSAbs_0DElement:
if (IdToRemove >= 0) if (IdToRemove >= 0)
@ -3305,8 +3301,10 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement * elem,
} }
removedElems.push_back((*it)); removedElems.push_back((*it));
myElementIDFactory->ReleaseID(IdToRemove, vtkid); myElementIDFactory->ReleaseID(IdToRemove, vtkid);
if (const SMDS_VtkEdge* vtkElem = dynamic_cast<const SMDS_VtkEdge*>(*it)) if (const SMDS_VtkEdge* vtkElem = dynamic_cast<const SMDS_VtkEdge*>(*it)) {
myEdgePool->destroy((SMDS_VtkEdge*) vtkElem); myEdgePool->destroy((SMDS_VtkEdge*) vtkElem);
((SMDS_MeshElement*) *it)->init( -1, -1, -1 ); // avoid reuse
}
else else
delete (*it); delete (*it);
break; break;
@ -3318,8 +3316,10 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement * elem,
} }
removedElems.push_back((*it)); removedElems.push_back((*it));
myElementIDFactory->ReleaseID(IdToRemove, vtkid); myElementIDFactory->ReleaseID(IdToRemove, vtkid);
if (const SMDS_VtkFace* vtkElem = dynamic_cast<const SMDS_VtkFace*>(*it)) if (const SMDS_VtkFace* vtkElem = dynamic_cast<const SMDS_VtkFace*>(*it)) {
myFacePool->destroy((SMDS_VtkFace*) vtkElem); myFacePool->destroy((SMDS_VtkFace*) vtkElem);
((SMDS_MeshElement*) *it)->init( -1, -1, -1 ); // avoid reuse
}
else else
delete (*it); delete (*it);
break; break;
@ -3331,8 +3331,10 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement * elem,
} }
removedElems.push_back((*it)); removedElems.push_back((*it));
myElementIDFactory->ReleaseID(IdToRemove, vtkid); myElementIDFactory->ReleaseID(IdToRemove, vtkid);
if (const SMDS_VtkVolume* vtkElem = dynamic_cast<const SMDS_VtkVolume*>(*it)) if (const SMDS_VtkVolume* vtkElem = dynamic_cast<const SMDS_VtkVolume*>(*it)) {
myVolumePool->destroy((SMDS_VtkVolume*) vtkElem); myVolumePool->destroy((SMDS_VtkVolume*) vtkElem);
((SMDS_MeshElement*) *it)->init( -1, -1, -1 ); // avoid reuse
}
else else
delete (*it); delete (*it);
break; break;
@ -3344,18 +3346,19 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement * elem,
} }
removedElems.push_back((*it)); removedElems.push_back((*it));
myElementIDFactory->ReleaseID(IdToRemove, vtkid); myElementIDFactory->ReleaseID(IdToRemove, vtkid);
if (const SMDS_BallElement* vtkElem = dynamic_cast<const SMDS_BallElement*>(*it)) if (const SMDS_BallElement* vtkElem = dynamic_cast<const SMDS_BallElement*>(*it)) {
myBallPool->destroy(const_cast<SMDS_BallElement*>( vtkElem )); myBallPool->destroy(const_cast<SMDS_BallElement*>( vtkElem ));
((SMDS_MeshElement*) *it)->init( -1, -1, -1 ); // avoid reuse
}
else else
delete (*it); delete (*it);
break; break;
case SMDSAbs_All: case SMDSAbs_All: // avoid compilation warning
case SMDSAbs_NbElementTypes: break; case SMDSAbs_NbElementTypes: break;
} }
if (vtkid >= 0) if (vtkid >= 0)
{ {
//MESSAGE("VTK_EMPTY_CELL in " << vtkid);
this->myGrid->GetCellTypesArray()->SetValue(vtkid, VTK_EMPTY_CELL); this->myGrid->GetCellTypesArray()->SetValue(vtkid, VTK_EMPTY_CELL);
} }
it++; it++;
@ -3368,7 +3371,6 @@ void SMDS_Mesh::RemoveElement(const SMDS_MeshElement * elem,
while (it != s2->end()) while (it != s2->end())
{ {
int IdToRemove = (*it)->GetID(); int IdToRemove = (*it)->GetID();
//MESSAGE( "SMDS: RM node " << IdToRemove);
if (IdToRemove >= 0) if (IdToRemove >= 0)
{ {
myNodes[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 // in meshes without descendants elements are always free
switch (aType) { switch (aType) {
case SMDSAbs_0DElement: case SMDSAbs_0DElement:
myCells[elemId] = 0; myCells[elemId] = 0;
myInfo.remove(elem); myInfo.remove(elem);
delete elem; delete elem;
elem = 0;
break; break;
case SMDSAbs_Edge: case SMDSAbs_Edge:
myCells[elemId] = 0; myCells[elemId] = 0;
@ -3463,6 +3466,9 @@ void SMDS_Mesh::RemoveFreeElement(const SMDS_MeshElement * elem)
this->myGrid->GetCellTypesArray()->SetValue(vtkId, VTK_EMPTY_CELL); this->myGrid->GetCellTypesArray()->SetValue(vtkId, VTK_EMPTY_CELL);
// --- to do: keep vtkid in a list of reusable cells // --- to do: keep vtkid in a list of reusable cells
if ( elem )
((SMDS_MeshElement*) elem)->init( -1, -1, -1 ); // avoid reuse
} }
} }

5
src/SMDS/SMDS_MeshNode.cxx
View File

@ -182,17 +182,14 @@ public:
MESSAGE("SMDS_MeshNode_MyInvIterator problem Null element"); MESSAGE("SMDS_MeshNode_MyInvIterator problem Null element");
throw SALOME_Exception("SMDS_MeshNode_MyInvIterator problem Null element"); throw SALOME_Exception("SMDS_MeshNode_MyInvIterator problem Null element");
} }
//MESSAGE("vtkId " << vtkId << " smdsId " << smdsId << " " << elem->GetType());
iter++; iter++;
return elem; return elem;
} }
}; };
SMDS_ElemIteratorPtr SMDS_MeshNode:: SMDS_ElemIteratorPtr SMDS_MeshNode::GetInverseElementIterator(SMDSAbs_ElementType type) const
GetInverseElementIterator(SMDSAbs_ElementType type) const
{ {
vtkCellLinks::Link l = SMDS_Mesh::_meshList[myMeshId]->getGrid()->GetCellLinks()->GetLink(myVtkID); 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)); return SMDS_ElemIteratorPtr(new SMDS_MeshNode_MyInvIterator(SMDS_Mesh::_meshList[myMeshId], l.cells, l.ncells, type));
} }

9
src/SMESHGUI/SMESHGUI_RemoveElementsDlg.cxx
View File

@ -42,6 +42,7 @@
#include <SUIT_Desktop.h> #include <SUIT_Desktop.h>
#include <SUIT_Session.h> #include <SUIT_Session.h>
#include <SUIT_MessageBox.h> #include <SUIT_MessageBox.h>
#include <SUIT_OverrideCursor.h>
#include <LightApp_Application.h> #include <LightApp_Application.h>
#include <LightApp_SelectionMgr.h> #include <LightApp_SelectionMgr.h>
@ -225,7 +226,10 @@ void SMESHGUI_RemoveElementsDlg::ClickOnApply()
if (mySMESHGUI->isActiveStudyLocked()) if (mySMESHGUI->isActiveStudyLocked())
return; return;
if (myNbOkElements) { if (myNbOkElements)
{
SUIT_OverrideCursor wc;
QStringList aListId = myEditCurrentArgument->text().split(" ", QString::SkipEmptyParts); QStringList aListId = myEditCurrentArgument->text().split(" ", QString::SkipEmptyParts);
SMESH::long_array_var anArrayOfIdeces = new SMESH::long_array; SMESH::long_array_var anArrayOfIdeces = new SMESH::long_array;
anArrayOfIdeces->length(aListId.count()); anArrayOfIdeces->length(aListId.count());
@ -233,7 +237,8 @@ void SMESHGUI_RemoveElementsDlg::ClickOnApply()
anArrayOfIdeces[i] = aListId[ i ].toInt(); anArrayOfIdeces[i] = aListId[ i ].toInt();
bool aResult = false; bool aResult = false;
try { try
{
SMESH::SMESH_MeshEditor_var aMeshEditor = myMesh->GetMeshEditor(); SMESH::SMESH_MeshEditor_var aMeshEditor = myMesh->GetMeshEditor();
aResult = aMeshEditor->RemoveElements(anArrayOfIdeces.in()); aResult = aMeshEditor->RemoveElements(anArrayOfIdeces.in());

31
src/SMESHUtils/SMESH_MeshAlgos.cxx
View File

@ -458,6 +458,10 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
void GetElementsNearLine( const gp_Ax1& line, void GetElementsNearLine( const gp_Ax1& line,
SMDSAbs_ElementType type, SMDSAbs_ElementType type,
vector< const SMDS_MeshElement* >& foundElems); vector< const SMDS_MeshElement* >& foundElems);
void GetElementsInSphere( const gp_XYZ& center,
const double radius,
SMDSAbs_ElementType type,
vector< const SMDS_MeshElement* >& foundElems);
double getTolerance(); double getTolerance();
bool getIntersParamOnLine(const gp_Lin& line, const SMDS_MeshElement* face, bool getIntersParamOnLine(const gp_Lin& line, const SMDS_MeshElement* face,
const double tolerance, double & param); const double tolerance, double & param);
@ -466,6 +470,7 @@ struct SMESH_ElementSearcherImpl: public SMESH_ElementSearcher
{ {
return _outerFaces.empty() || _outerFaces.count(face); return _outerFaces.empty() || _outerFaces.count(face);
} }
struct TInters //!< data of intersection of the line and the mesh face (used in GetPointState()) struct TInters //!< data of intersection of the line and the mesh face (used in GetPointState())
{ {
const SMDS_MeshElement* _face; 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(); set< const SMDS_MeshElement*, TIDCompare >::const_iterator face = faces.begin();
for ( ; face != faces.end(); ++face ) for ( ; face != faces.end(); ++face )
{ {
if ( *face == outerFace ) continue;
if ( !SMESH_MeshAlgos::FaceNormal( *face, fNorm, /*normalized=*/false )) if ( !SMESH_MeshAlgos::FaceNormal( *face, fNorm, /*normalized=*/false ))
continue; continue;
gp_Vec dirInF = gp_Vec( fNorm ) ^ n1n2; 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 // store the found outer face and add its links to continue seaching from
if ( outerFace2 ) if ( outerFace2 )
{ {
_outerFaces.insert( outerFace ); _outerFaces.insert( outerFace2 );
int nbNodes = outerFace2->NbNodes()/( outerFace2->IsQuadratic() ? 2 : 1 ); int nbNodes = outerFace2->NbCornerNodes();
for ( int i = 0; i < nbNodes; ++i ) for ( int i = 0; i < nbNodes; ++i )
{ {
SMESH_TLink link2( outerFace2->GetNode(i), outerFace2->GetNode((i+1)%nbNodes)); 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()); 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 * \brief Return true if the point is IN or ON of the element

7
src/SMESHUtils/SMESH_MeshAlgos.hxx
View File

@ -88,6 +88,13 @@ struct SMESHUtils_EXPORT SMESH_ElementSearcher
virtual void GetElementsNearLine( const gp_Ax1& line, virtual void GetElementsNearLine( const gp_Ax1& line,
SMDSAbs_ElementType type, SMDSAbs_ElementType type,
std::vector< const SMDS_MeshElement* >& foundElems) = 0; 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. * \brief Find out if the given point is out of closed 2D mesh.
*/ */

432
src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx
View File

@ -24,6 +24,7 @@
#include "StdMeshers_QuadToTriaAdaptor.hxx" #include "StdMeshers_QuadToTriaAdaptor.hxx"
#include "SMDS_IteratorOnIterators.hxx"
#include "SMDS_SetIterator.hxx" #include "SMDS_SetIterator.hxx"
#include "SMESHDS_GroupBase.hxx" #include "SMESHDS_GroupBase.hxx"
#include "SMESH_Algo.hxx" #include "SMESH_Algo.hxx"
@ -115,20 +116,20 @@ namespace
gp_Vec nJ = baseVec.Crossed( baJ ); gp_Vec nJ = baseVec.Crossed( baJ );
// Check angle between normals // Check angle between normals
double angle = nI.Angle( nJ ); double angle = nI.Angle( nJ );
bool tooClose = ( angle < 15. * M_PI / 180. ); bool tooClose = ( angle < 15. * M_PI / 180. );
// Check if pyramids collide // Check if pyramids collide
if ( !tooClose && ( baI * baJ > 0 ) && ( nI * nJ > 0 )) if ( !tooClose && ( baI * baJ > 0 ) && ( nI * nJ > 0 ))
{ {
// find out if nI points outside of PrmI or inside // 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; 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; 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 // Check if PrmI and PrmJ are in same domain
@ -170,8 +171,9 @@ namespace
continue; // f is a base quadrangle continue; // f is a base quadrangle
// check projections of face direction (baOFN) to triange normals (nI and nJ) // check projections of face direction (baOFN) to triange normals (nI and nJ)
gp_Vec baOFN( base1, SMESH_TNodeXYZ( otherFaceNode )); gp_Vec baOFN( base2, SMESH_TNodeXYZ( otherFaceNode ));
if ( nI * baOFN > 0 && nJ * baOFN > 0 ) if ( nI * baOFN > 0 && nJ * baOFN > 0 &&
baI* baOFN > 0 && baJ* baOFN > 0 ) // issue 0023212
{ {
tooClose = false; // f is between pyramids tooClose = false; // f is between pyramids
break; break;
@ -253,7 +255,6 @@ namespace
} }
} }
} }
} }
//================================================================================ //================================================================================
@ -266,6 +267,8 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr
const SMDS_MeshElement* PrmJ, const SMDS_MeshElement* PrmJ,
set<const SMDS_MeshNode*> & nodesToMove) set<const SMDS_MeshNode*> & 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 const SMDS_MeshNode* Nrem = PrmJ->GetNode(4); // node to remove
//int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume ); //int nbJ = Nrem->NbInverseElements( SMDSAbs_Volume );
SMESH_TNodeXYZ Pj( Nrem ); SMESH_TNodeXYZ Pj( Nrem );
@ -288,7 +291,7 @@ void StdMeshers_QuadToTriaAdaptor::MergePiramids( const SMDS_MeshElement* Pr
vector< const SMDS_MeshElement* > inverseElems vector< const SMDS_MeshElement* > inverseElems
// copy inverse elements to avoid iteration on changing container // copy inverse elements to avoid iteration on changing container
( TStdElemIterator( CommonNode->GetInverseElementIterator(SMDSAbs_Face)), itEnd); ( 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* FI = inverseElems[i];
const SMDS_MeshElement* FJEqual = 0; 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 // set the common apex node to pyramids and triangles merged with J
vector< const SMDS_MeshNode* > nodes;
inverseElems.assign( TStdElemIterator( Nrem->GetInverseElementIterator()), itEnd ); 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]; 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; nodes[ elem->GetType() == SMDSAbs_Volume ? PYRAM_APEX : TRIA_APEX ] = CommonNode;
GetMeshDS()->ChangeElementNodes( elem, &nodes[0], nodes.size()); 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, void StdMeshers_QuadToTriaAdaptor::MergeAdjacent(const SMDS_MeshElement* PrmI,
set<const SMDS_MeshNode*>& nodesToMove) set<const SMDS_MeshNode*>& nodesToMove,
const bool isRecursion)
{ {
TIDSortedElemSet adjacentPyrams; TIDSortedElemSet adjacentPyrams;
bool mergedPyrams = false; 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 ); SMDS_ElemIteratorPtr vIt = n->GetInverseElementIterator( SMDSAbs_Volume );
while ( vIt->more() ) while ( vIt->more() )
{ {
const SMDS_MeshElement* PrmJ = vIt->next(); 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; continue;
if ( PrmI != PrmJ && TooCloseAdjacent( PrmI, PrmJ, GetMesh()->HasShapeToMesh() )) if ( TooCloseAdjacent( PrmI, PrmJ, GetMesh()->HasShapeToMesh() ))
{ {
MergePiramids( PrmI, PrmJ, nodesToMove ); MergePiramids( PrmI, PrmJ, nodesToMove );
mergedPyrams = true; mergedPyrams = true;
// container of inverse elements can change // 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; TIDSortedElemSet::iterator prm;
for (prm = adjacentPyrams.begin(); prm != adjacentPyrams.end(); ++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 //function : HasIntersection3
//purpose : Auxilare for HasIntersection() //purpose : Find intersection point between a triangle (P1,P2,P3)
// find intersection point between triangle (P1,P2,P3) // and a segment [PC,P]
// and segment [PC,P]
//======================================================================= //=======================================================================
static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint, 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) const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3)
{ {
//cout<<"HasIntersection3"<<endl; const double EPSILON = 1e-6;
//cout<<" PC("<<PC.X()<<","<<PC.Y()<<","<<PC.Z()<<")"<<endl; double segLen = P.Distance( PC );
//cout<<" P("<<P.X()<<","<<P.Y()<<","<<P.Z()<<")"<<endl;
//cout<<" P1("<<P1.X()<<","<<P1.Y()<<","<<P1.Z()<<")"<<endl;
//cout<<" P2("<<P2.X()<<","<<P2.Y()<<","<<P2.Z()<<")"<<endl;
//cout<<" P3("<<P3.X()<<","<<P3.Y()<<","<<P3.Z()<<")"<<endl;
gp_Vec VP1(P1,P2);
gp_Vec VP2(P1,P3);
IntAna_Quadric IAQ(gp_Pln(P1,VP1.Crossed(VP2)));
IntAna_IntConicQuad IAICQ(gp_Lin(PC,gp_Dir(gp_Vec(PC,P))),IAQ);
if(IAICQ.IsDone()) {
if( IAICQ.IsInQuadric() )
return false;
if( IAICQ.NbPoints() == 1 ) {
gp_Pnt PIn = IAICQ.Point(1);
const double preci = 1.e-10 * P.Distance(PC);
// check if this point is internal for segment [PC,P]
bool IsExternal =
( (PC.X()-PIn.X())*(P.X()-PIn.X()) > 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()<preci ||
V2.Magnitude()<preci ||
V3.Magnitude()<preci ) {
Pint = PIn;
return true;
}
const double angularTol = 1e-6;
gp_Vec VC1 = V1.Crossed(V2);
gp_Vec VC2 = V2.Crossed(V3);
gp_Vec VC3 = V3.Crossed(V1);
if(VC1.Magnitude()<gp::Resolution()) {
if(VC2.IsOpposite(VC3,angularTol)) {
return false;
}
}
else if(VC2.Magnitude()<gp::Resolution()) {
if(VC1.IsOpposite(VC3,angularTol)) {
return false;
}
}
else if(VC3.Magnitude()<gp::Resolution()) {
if(VC1.IsOpposite(VC2,angularTol)) {
return false;
}
}
else {
if( VC1.IsOpposite(VC2,angularTol) || VC1.IsOpposite(VC3,angularTol) ||
VC2.IsOpposite(VC3,angularTol) ) {
return false;
}
}
Pint = PIn;
return true;
}
}
return false; gp_XYZ orig = PC.XYZ();
gp_XYZ dir = ( P.XYZ() - PC.XYZ() ) / segLen;
gp_XYZ vert0 = P1.XYZ();
gp_XYZ vert1 = P2.XYZ();
gp_XYZ vert2 = P3.XYZ();
/* calculate distance from vert0 to ray origin */
gp_XYZ tvec = orig - vert0;
gp_XYZ edge1 = vert1 - vert0;
gp_XYZ edge2 = vert2 - vert0;
/* begin calculating determinant - also used to calculate U parameter */
gp_XYZ pvec = dir ^ edge2;
/* if determinant is near zero, ray lies in plane of triangle */
double det = edge1 * pvec;
if (det > -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. * \brief Return allowed height of a pyramid
* \param P - first segment end * \param Papex - optimal pyramid apex
* \param PC - second segment end (it is a gravity center of quadrangle) * \param PC - gravity center of a quadrangle
* \param Pint - (out) intersection point * \param PN - four nodes of the quadrangle
* \param aMesh - mesh * \param aMesh - mesh
* \param aShape - shape to check faces on * \param NotCheckedFace - the quadrangle face
* \param NotCheckedFace - mesh face not to check * \retval double - pyramid height
* \retval bool - true if there is an intersection
*/ */
//================================================================================ //================================================================================
bool StdMeshers_QuadToTriaAdaptor::CheckIntersection (const gp_Pnt& P, void StdMeshers_QuadToTriaAdaptor::LimitHeight (gp_Pnt& Papex,
const gp_Pnt& PC, const gp_Pnt& PC,
gp_Pnt& Pint, const TColgp_Array1OfPnt& PN,
SMESH_Mesh& aMesh, const vector<const SMDS_MeshNode*>& FNodes,
const TopoDS_Shape& aShape, SMESH_Mesh& aMesh,
const SMDS_MeshElement* NotCheckedFace) const SMDS_MeshElement* NotCheckedFace,
const bool UseApexRay)
{ {
if ( !myElemSearcher ) if ( !myElemSearcher )
myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *aMesh.GetMeshDS() ); myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *aMesh.GetMeshDS() );
SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher); SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher);
bool res = false; // Find intersection of faces with (P,PC) segment elongated 3 times
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);
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; TColgp_SequenceOfPnt aContour;
for ( size_t iF = 0; iF < suspectElems.size(); ++iF )
if ( UseApexRay )
{ {
const SMDS_MeshElement* face = suspectElems[iF]; // find intersection closest to PC
if ( face == NotCheckedFace ) continue; Ptest = PC.XYZ() + line.Direction().XYZ() * height * 3;
aContour.Clear();
for ( int i = 0; i < face->NbCornerNodes(); ++i ) searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectFaces );
aContour.Append( SMESH_TNodeXYZ( face->GetNode(i) )); for ( size_t iF = 0; iF < suspectFaces.size(); ++iF )
if ( HasIntersection(P, PC, Pres, aContour)) { {
res = true; const SMDS_MeshElement* face = suspectFaces[iF];
double tmp = PC.Distance(Pres); if ( face == NotCheckedFace ) continue;
if ( tmp < dist ) {
Pint = Pres; aContour.Clear();
dist = tmp; 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<const SMDS_MeshElement*> myPyramids; vector<const SMDS_MeshElement*> myPyramids;
const SMESHDS_SubMesh * aSubMeshDSFace;
SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
SMESH_MesherHelper helper(aMesh); SMESH_MesherHelper helper(aMesh);
helper.IsQuadraticSubMesh(aShape); helper.IsQuadraticSubMesh(aShape);
helper.SetElementsOnShape( true ); helper.SetElementsOnShape( true );
if ( myElemSearcher ) delete myElemSearcher; if ( myElemSearcher ) delete myElemSearcher;
vector< SMDS_ElemIteratorPtr > itVec;
if ( aProxyMesh ) if ( aProxyMesh )
myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, aProxyMesh->GetFaces(aShape)); {
itVec.push_back( aProxyMesh->GetFaces( aShape ));
}
else 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_Array1OfPnt PN(1,5);
TColgp_Array1OfVec VN(1,4); TColgp_Array1OfVec VN(1,4);
vector<const SMDS_MeshNode*> FNodes(5); vector<const SMDS_MeshNode*> FNodes(5);
gp_Pnt PC; gp_Pnt PC;
gp_Vec VNorm; 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(); const TopoDS_Shape& aShapeFace = exp.Current();
if ( aProxyMesh ) if ( aProxyMesh )
@ -809,17 +849,7 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh,
} }
else { else {
// check possible intersection with other faces // check possible intersection with other faces
gp_Pnt Pint; LimitHeight( PCbest, PC, PN, FNodes, aMesh, face, /*UseApexRay=*/true );
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;
}
}
} }
// create node for PCbest // create node for PCbest
SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() ); 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_Pnt Pres1,Pres2;
gp_Ax1 line( PC, VNorm ); gp_Ax1 line( PC, VNorm );
vector< const SMDS_MeshElement* > suspectElems; vector< const SMDS_MeshElement* > suspectFaces;
searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems); 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; if ( F == face ) continue;
aContour.Clear(); aContour.Clear();
for ( int i = 0; i < 4; ++i ) for ( int i = 0; i < 4; ++i )
@ -1026,7 +1056,9 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
continue; continue;
} }
// -----------------------------------
// Case of non-degenerated quadrangle // Case of non-degenerated quadrangle
// -----------------------------------
// Find pyramid peak // Find pyramid peak
@ -1059,12 +1091,12 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
gp_Pnt intPnt[2]; gp_Pnt intPnt[2];
gp_Ax1 line( PC, tmpDir ); gp_Ax1 line( PC, tmpDir );
vector< const SMDS_MeshElement* > suspectElems; vector< const SMDS_MeshElement* > suspectFaces;
searcher->GetElementsNearLine( line, SMDSAbs_Face, suspectElems); 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; if ( F == face ) continue;
aContour.Clear(); aContour.Clear();
int nbN = F->NbCornerNodes(); int nbN = F->NbCornerNodes();
@ -1109,13 +1141,15 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
{ {
if( !intersected[isRev] ) continue; if( !intersected[isRev] ) continue;
double pyramidH = Min( height, PC.Distance(intPnt[isRev])/3.); 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 LimitHeight( Papex, PC, PN, FNodes, aMesh, face, /*UseApexRay=*/false );
SMDS_MeshNode* NewNode = helper.AddNode( PCbest.X(), PCbest.Y(), PCbest.Z() );
// create node for Papex
SMDS_MeshNode* NewNode = helper.AddNode( Papex.X(), Papex.Y(), Papex.Z() );
// add triangles to result map // add triangles to result map
for(i=0; i<4; i++) { for ( i = 0; i < 4; i++) {
SMDS_MeshFace* NewFace; SMDS_MeshFace* NewFace;
if(isRev) if(isRev)
NewFace = meshDS->AddFace( NewNode, FNodes[i], FNodes[i+1] ); 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, bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh& aMesh,
const vector<const SMDS_MeshElement*>& myPyramids) const vector<const SMDS_MeshElement*>& myPyramids)
{ {
if(myPyramids.empty()) if ( myPyramids.empty() )
return true; return true;
SMESHDS_Mesh * meshDS = aMesh.GetMeshDS(); SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
size_t i, j, k; size_t i, j, k;
int myShapeID = myPyramids[0]->GetNode(4)->getshapeId(); //int myShapeID = myPyramids[0]->GetNode(4)->getshapeId();
{
if ( myElemSearcher ) delete myElemSearcher; SMDS_ElemIteratorPtr
myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS ); pyramIt( new SMDS_ElementVectorIterator( myPyramids.begin(), myPyramids.end() ));
SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>(myElemSearcher); if ( myElemSearcher ) delete myElemSearcher;
myElemSearcher = SMESH_MeshAlgos::GetElementSearcher( *meshDS, pyramIt );
}
SMESH_ElementSearcher* searcher = const_cast<SMESH_ElementSearcher*>( myElemSearcher );
set<const SMDS_MeshNode*> nodesToMove; set<const SMDS_MeshNode*> nodesToMove;
@ -1167,17 +1204,19 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&
MergeAdjacent( PrmI, nodesToMove ); MergeAdjacent( PrmI, nodesToMove );
} }
// iterate on all pyramids // iterate on all new pyramids
vector< const SMDS_MeshElement* > suspectPyrams;
for ( i = 0; i < myPyramids.size(); ++i ) 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 // compare PrmI with all the rest pyramids
// collect adjacent pyramids and nodes coordinates of PrmI // collect adjacent pyramids and nodes coordinates of PrmI
set<const SMDS_MeshElement*> checkedPyrams; set<const SMDS_MeshElement*> checkedPyrams;
vector<gp_Pnt> PsI(5); gp_Pnt PsI[5];
for(k=0; k<5; k++) // loop on 4 base nodes of PrmI for ( k = 0; k < 5; k++ )
{ {
const SMDS_MeshNode* n = PrmI->GetNode(k); const SMDS_MeshNode* n = PrmI->GetNode(k);
PsI[k] = SMESH_TNodeXYZ( n ); 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 // 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]); const SMDS_MeshElement* PrmJ = suspectPyrams[j];
gp_Ax1 line( PsI[k], Vtmp ); if ( PrmJ == PrmI )
vector< const SMDS_MeshElement* > suspectPyrams; continue;
searcher->GetElementsNearLine( line, SMDSAbs_Volume, suspectPyrams); 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<gp_Pnt> PsJ( xyzIt, TXyzIterator() );
gp_Pnt Pint; gp_Pnt Pint;
bool hasInt=false; bool hasInt=false;
for(k=0; k<4 && !hasInt; k++) { for ( k = 0; k < 4 && !hasInt; k++ )
gp_Vec Vtmp(PsI[k],PsI[4]); {
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 gp_Pnt Pshift = PsI[k].XYZ() + Vtmp.XYZ() * 0.01; // base node moved a bit to apex
hasInt = hasInt =
( HasIntersection3( Pshift, PsI[4], Pint, PsJ[0], PsJ[1], PsJ[4]) || ( HasIntersection3( Pshift, PsI[4], Pint, PsJ[0], PsJ[1], PsJ[PYRAM_APEX]) ||
HasIntersection3( Pshift, PsI[4], Pint, PsJ[1], PsJ[2], PsJ[4]) || HasIntersection3( Pshift, PsI[4], Pint, PsJ[1], PsJ[2], PsJ[PYRAM_APEX]) ||
HasIntersection3( Pshift, PsI[4], Pint, PsJ[2], PsJ[3], PsJ[4]) || HasIntersection3( Pshift, PsI[4], Pint, PsJ[2], PsJ[3], PsJ[PYRAM_APEX]) ||
HasIntersection3( Pshift, PsI[4], Pint, PsJ[3], PsJ[0], PsJ[4]) ); HasIntersection3( Pshift, PsI[4], Pint, PsJ[3], PsJ[0], PsJ[PYRAM_APEX]) );
} }
for(k=0; k<4 && !hasInt; k++) { for ( k = 0; k < 4 && !hasInt; k++ )
gp_Vec Vtmp(PsJ[k],PsJ[4]); {
gp_Vec Vtmp( PsJ[k], PsJ[ PYRAM_APEX ]);
gp_Pnt Pshift = PsJ[k].XYZ() + Vtmp.XYZ() * 0.01; gp_Pnt Pshift = PsJ[k].XYZ() + Vtmp.XYZ() * 0.01;
hasInt = hasInt =
( HasIntersection3( Pshift, PsJ[4], Pint, PsI[0], PsI[1], PsI[4]) || ( HasIntersection3( Pshift, PsJ[4], Pint, PsI[0], PsI[1], PsI[PYRAM_APEX]) ||
HasIntersection3( Pshift, PsJ[4], Pint, PsI[1], PsI[2], PsI[4]) || HasIntersection3( Pshift, PsJ[4], Pint, PsI[1], PsI[2], PsI[PYRAM_APEX]) ||
HasIntersection3( Pshift, PsJ[4], Pint, PsI[2], PsI[3], PsI[4]) || HasIntersection3( Pshift, PsJ[4], Pint, PsI[2], PsI[3], PsI[PYRAM_APEX]) ||
HasIntersection3( Pshift, PsJ[4], Pint, PsI[3], PsI[0], PsI[4]) ); HasIntersection3( Pshift, PsJ[4], Pint, PsI[3], PsI[0], PsI[PYRAM_APEX]) );
} }
if ( hasInt ) if ( hasInt )
{ {
// count common nodes of base faces of two pyramids // count common nodes of base faces of two pyramids
int nbc = 0; int nbc = 0;
for (k=0; k<4; k++) for ( k = 0; k < 4; k++ )
nbc += int ( PrmI->GetNodeIndex( PrmJ->GetNode(k) ) >= 0 ); nbc += int ( PrmI->GetNodeIndex( PrmJ->GetNode(k) ) >= 0 );
if ( nbc == 4 ) if ( nbc == 4 )
continue; // pyrams have a common base face continue; // pyrams have a common base face
if(nbc>0) if ( nbc > 0 )
{ {
// Merge the two pyramids and others already merged with them // Merge the two pyramids and others already merged with them
MergePiramids( PrmI, PrmJ, nodesToMove ); MergePiramids( PrmI, PrmJ, nodesToMove );
} }
else { // nbc==0 else // nbc==0
{
// decrease height of pyramids // decrease height of pyramids
gp_XYZ PCi(0,0,0), PCj(0,0,0); 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(); PCi += PsI[k].XYZ();
PCj += PsJ[k].XYZ(); PCj += PsJ[k].XYZ();
} }
@ -1272,9 +1318,9 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&
VN1.Scale(coef1); VN1.Scale(coef1);
VN2.Scale(coef2); VN2.Scale(coef2);
SMDS_MeshNode* aNode1 = const_cast<SMDS_MeshNode*>(PrmI->GetNode(4)); SMDS_MeshNode* aNode1 = const_cast<SMDS_MeshNode*>( apexI );
aNode1->setXYZ( PCi.X()+VN1.X(), PCi.Y()+VN1.Y(), PCi.Z()+VN1.Z() ); aNode1->setXYZ( PCi.X()+VN1.X(), PCi.Y()+VN1.Y(), PCi.Z()+VN1.Z() );
SMDS_MeshNode* aNode2 = const_cast<SMDS_MeshNode*>(PrmJ->GetNode(4)); SMDS_MeshNode* aNode2 = const_cast<SMDS_MeshNode*>(PrmJ->GetNode( PYRAM_APEX ));
aNode2->setXYZ( PCj.X()+VN2.X(), PCj.Y()+VN2.Y(), PCj.Z()+VN2.Z() ); aNode2->setXYZ( PCj.X()+VN2.X(), PCj.Y()+VN2.Y(), PCj.Z()+VN2.Z() );
nodesToMove.insert( aNode1 ); nodesToMove.insert( aNode1 );
nodesToMove.insert( aNode2 ); nodesToMove.insert( aNode2 );

14
src/StdMeshers/StdMeshers_QuadToTriaAdaptor.hxx
View File

@ -71,10 +71,13 @@ protected:
gp_Pnt& PC, gp_Vec& VNorm, gp_Pnt& PC, gp_Vec& VNorm,
const SMDS_MeshElement** volumes=0); const SMDS_MeshElement** volumes=0);
bool CheckIntersection(const gp_Pnt& P, const gp_Pnt& PC, void LimitHeight (gp_Pnt& Papex,
gp_Pnt& Pint, SMESH_Mesh& aMesh, const gp_Pnt& PC,
const TopoDS_Shape& aShape, const TColgp_Array1OfPnt& PN,
const SMDS_MeshElement* NotCheckedFace); const std::vector<const SMDS_MeshNode*>& FNodes,
SMESH_Mesh& aMesh,
const SMDS_MeshElement* NotCheckedFace,
const bool UseApexRay);
bool Compute2ndPart(SMESH_Mesh& aMesh, bool Compute2ndPart(SMESH_Mesh& aMesh,
const std::vector<const SMDS_MeshElement*>& pyramids); const std::vector<const SMDS_MeshElement*>& pyramids);
@ -85,7 +88,8 @@ protected:
std::set<const SMDS_MeshNode*> & nodesToMove); std::set<const SMDS_MeshNode*> & nodesToMove);
void MergeAdjacent(const SMDS_MeshElement* PrmI, void MergeAdjacent(const SMDS_MeshElement* PrmI,
std::set<const SMDS_MeshNode*>& nodesToMove); std::set<const SMDS_MeshNode*>& nodesToMove,
const bool isRecursion = false);
TopoDS_Shape myShape; TopoDS_Shape myShape;