-// Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2021 CEA/DEN, EDF R&D, OPEN CASCADE
//
// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
#include "SMESH_MeshAlgos.hxx"
#include "SMESH_OctreeNode.hxx"
+#include <GEOMUtils.hxx>
#include <Basics_Utils.hxx>
#include <BRepAdaptor_Surface.hxx>
// Case 1 Case 2
// | | | | |
// | | | | |
- // +-----+------+ +-----+------+
+ // +-----+------+ +-----+------+
// | | | |
// | | | |
// result should be 2 in both cases
//
int aResult0 = 0, aResult1 = 0;
- // last node, it is a medium one in a quadratic edge
+ // last node, it is a medium one in a quadratic edge
const SMDS_MeshNode* aLastNode = anEdge->GetNode( anEdge->NbNodes() - 1 );
const SMDS_MeshNode* aNode0 = anEdge->GetNode( 0 );
const SMDS_MeshNode* aNode1 = anEdge->GetNode( 1 );
return false;
const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
- if ( !anElem || anElem->GetType() != this->GetType() )
+ if ( !IsApplicable( anElem ))
return false;
return GetPoints( anElem, theRes );
return ( myPrecision >= 0 ) ? floor( aVal * myPrecisionValue + 0.5 ) / myPrecisionValue : aVal;
}
+//================================================================================
+/*!
+ * \brief Return true if a value can be computed for a given element.
+ * Some NumericalFunctor's are meaningful for elements of a certain
+ * geometry only.
+ */
+//================================================================================
+
+bool NumericalFunctor::IsApplicable( const SMDS_MeshElement* element ) const
+{
+ return element && element->GetType() == this->GetType();
+}
+
+bool NumericalFunctor::IsApplicable( long theElementId ) const
+{
+ return IsApplicable( myMesh->FindElement( theElementId ));
+}
+
//================================================================================
/*!
* \brief Return histogram of functor values
// }
// }
// { // polygons
-
+
// }
if( myPrecision >= 0 )
{
double aVal = 0;
myCurrElement = myMesh->FindElement( theId );
- if ( myCurrElement && myCurrElement->GetVtkType() == VTK_QUAD )
- {
- // issue 21723
- vtkUnstructuredGrid* grid = const_cast<SMDS_Mesh*>( myMesh )->GetGrid();
- if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->GetVtkID() ))
- aVal = Round( vtkMeshQuality::QuadAspectRatio( avtkCell ));
- }
- else
- {
- TSequenceOfXYZ P;
- if ( GetPoints( myCurrElement, P ))
- aVal = Round( GetValue( P ));
- }
+ TSequenceOfXYZ P;
+ if ( GetPoints( myCurrElement, P ))
+ aVal = Round( GetValue( P ));
return aVal;
}
//
// alpha = sqrt( 1/32 )
// L = max( L1, L2, L3, L4, D1, D2 )
- // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
+ // C1 = sqrt( L1^2 + L1^2 + L1^2 + L1^2 )
// C2 = min( S1, S2, S3, S4 )
// Li - lengths of the edges
// Di - lengths of the diagonals
Max( aLen[ 2 ],
Max( aLen[ 3 ],
Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
- double C1 = sqrt( ( aLen[0] * aLen[0] +
- aLen[1] * aLen[1] +
- aLen[2] * aLen[2] +
- aLen[3] * aLen[3] ) / 4. );
+ double C1 = sqrt( aLen[0] * aLen[0] +
+ aLen[1] * aLen[1] +
+ aLen[2] * aLen[2] +
+ aLen[3] * aLen[3] );
double C2 = Min( anArea[ 0 ],
Min( anArea[ 1 ],
Min( anArea[ 2 ], anArea[ 3 ] ) ) );
//
// alpha = sqrt( 1/32 )
// L = max( L1, L2, L3, L4, D1, D2 )
- // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
+ // C1 = sqrt( L1^2 + L1^2 + L1^2 + L1^2 )
// C2 = min( S1, S2, S3, S4 )
// Li - lengths of the edges
// Di - lengths of the diagonals
// Si - areas of the triangles
const double alpha = sqrt( 1 / 32. );
double L = Max( aLen[ 0 ],
- Max( aLen[ 1 ],
- Max( aLen[ 2 ],
- Max( aLen[ 3 ],
- Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
- double C1 = sqrt( ( aLen[0] * aLen[0] +
- aLen[1] * aLen[1] +
- aLen[2] * aLen[2] +
- aLen[3] * aLen[3] ) / 4. );
+ Max( aLen[ 1 ],
+ Max( aLen[ 2 ],
+ Max( aLen[ 3 ],
+ Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
+ double C1 = sqrt( aLen[0] * aLen[0] +
+ aLen[1] * aLen[1] +
+ aLen[2] * aLen[2] +
+ aLen[3] * aLen[3] );
double C2 = Min( anArea[ 0 ],
- Min( anArea[ 1 ],
- Min( anArea[ 2 ], anArea[ 3 ] ) ) );
+ Min( anArea[ 1 ],
+ Min( anArea[ 2 ], anArea[ 3 ] ) ) );
if ( C2 <= theEps )
return theInf;
return alpha * L * C1 / C2;
return 0;
}
+bool AspectRatio::IsApplicable( const SMDS_MeshElement* element ) const
+{
+ return ( NumericalFunctor::IsApplicable( element ) && !element->IsPoly() );
+}
+
double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const
{
// the aspect ratio is in the range [1.0,infinity]
return aHeight;
}
+ //================================================================================
+ /*!
+ * \brief Standard quality of a tetrahedron but not normalized
+ */
+ //================================================================================
+
+ double tetQualityByHomardMethod( const gp_XYZ & p1,
+ const gp_XYZ & p2,
+ const gp_XYZ & p3,
+ const gp_XYZ & p4 )
+ {
+ gp_XYZ edgeVec[6];
+ edgeVec[0] = ( p1 - p2 );
+ edgeVec[1] = ( p2 - p3 );
+ edgeVec[2] = ( p3 - p1 );
+ edgeVec[3] = ( p4 - p1 );
+ edgeVec[4] = ( p4 - p2 );
+ edgeVec[5] = ( p4 - p3 );
+
+ double maxEdgeLen2 = edgeVec[0].SquareModulus();
+ maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[1].SquareModulus() );
+ maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[2].SquareModulus() );
+ maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[3].SquareModulus() );
+ maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[4].SquareModulus() );
+ maxEdgeLen2 = Max( maxEdgeLen2, edgeVec[5].SquareModulus() );
+ double maxEdgeLen = Sqrt( maxEdgeLen2 );
+
+ gp_XYZ cross01 = edgeVec[0] ^ edgeVec[1];
+ double sumArea = ( cross01 ).Modulus(); // actually double area
+ sumArea += ( edgeVec[0] ^ edgeVec[3] ).Modulus();
+ sumArea += ( edgeVec[1] ^ edgeVec[4] ).Modulus();
+ sumArea += ( edgeVec[2] ^ edgeVec[5] ).Modulus();
+
+ double sixVolume = Abs( cross01 * edgeVec[4] ); // 6 * volume
+ double quality = maxEdgeLen * sumArea / sixVolume; // not normalized!!!
+ return quality;
+ }
+
+ //================================================================================
+ /*!
+ * \brief HOMARD method of hexahedron quality
+ * 1. Decompose the hexa into 24 tetra: each face is splitted into 4 triangles by
+ * adding the diagonals and every triangle is connected to the center of the hexa.
+ * 2. Compute the quality of every tetra with the same formula as for the standard quality,
+ * except that the factor for the normalization is not the same because the final goal
+ * is to have a quality equal to 1 for a perfect cube. So the formula is:
+ * qual = max(lengthes of 6 edges) * (sum of surfaces of 4 faces) / (7.6569*6*volume)
+ * 3. The quality of the hexa is the highest value of the qualities of the 24 tetra
+ */
+ //================================================================================
+
+ double hexQualityByHomardMethod( const TSequenceOfXYZ& P )
+ {
+ gp_XYZ quadCenter[6];
+ quadCenter[0] = ( P(1) + P(2) + P(3) + P(4) ) / 4.;
+ quadCenter[1] = ( P(5) + P(6) + P(7) + P(8) ) / 4.;
+ quadCenter[2] = ( P(1) + P(2) + P(6) + P(5) ) / 4.;
+ quadCenter[3] = ( P(2) + P(3) + P(7) + P(6) ) / 4.;
+ quadCenter[4] = ( P(3) + P(4) + P(8) + P(7) ) / 4.;
+ quadCenter[5] = ( P(1) + P(4) + P(8) + P(5) ) / 4.;
+
+ gp_XYZ hexCenter = ( P(1) + P(2) + P(3) + P(4) + P(5) + P(6) + P(7) + P(8) ) / 8.;
+
+ // quad 1 ( 1 2 3 4 )
+ double quality = tetQualityByHomardMethod( P(1), P(2), quadCenter[0], hexCenter );
+ quality = Max( quality, tetQualityByHomardMethod( P(2), P(3), quadCenter[0], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(3), P(4), quadCenter[0], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(4), P(1), quadCenter[0], hexCenter ));
+ // quad 2 ( 5 6 7 8 )
+ quality = Max( quality, tetQualityByHomardMethod( P(5), P(6), quadCenter[1], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(6), P(7), quadCenter[1], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(7), P(8), quadCenter[1], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(8), P(5), quadCenter[1], hexCenter ));
+ // quad 3 ( 1 2 6 5 )
+ quality = Max( quality, tetQualityByHomardMethod( P(1), P(2), quadCenter[2], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(2), P(6), quadCenter[2], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(6), P(5), quadCenter[2], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(5), P(1), quadCenter[2], hexCenter ));
+ // quad 4 ( 2 3 7 6 )
+ quality = Max( quality, tetQualityByHomardMethod( P(2), P(3), quadCenter[3], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(3), P(7), quadCenter[3], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(7), P(6), quadCenter[3], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(6), P(2), quadCenter[3], hexCenter ));
+ // quad 5 ( 3 4 8 7 )
+ quality = Max( quality, tetQualityByHomardMethod( P(3), P(4), quadCenter[4], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(4), P(8), quadCenter[4], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(8), P(7), quadCenter[4], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(7), P(3), quadCenter[4], hexCenter ));
+ // quad 6 ( 1 4 8 5 )
+ quality = Max( quality, tetQualityByHomardMethod( P(1), P(4), quadCenter[5], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(4), P(8), quadCenter[5], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(8), P(5), quadCenter[5], hexCenter ));
+ quality = Max( quality, tetQualityByHomardMethod( P(5), P(1), quadCenter[5], hexCenter ));
+
+ return quality / 7.65685424949;
+ }
}
double AspectRatio3D::GetValue( long theId )
return aVal;
}
+bool AspectRatio3D::IsApplicable( const SMDS_MeshElement* element ) const
+{
+ return ( NumericalFunctor::IsApplicable( element ) && !element->IsPoly() );
+}
+
double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
{
double aQuality = 0.0;
int nbNodes = P.size();
if( myCurrElement->IsQuadratic() ) {
- if(nbNodes==10) nbNodes=4; // quadratic tetrahedron
+ if (nbNodes==10) nbNodes=4; // quadratic tetrahedron
else if(nbNodes==13) nbNodes=5; // quadratic pyramid
else if(nbNodes==15) nbNodes=6; // quadratic pentahedron
else if(nbNodes==20) nbNodes=8; // quadratic hexahedron
- else if(nbNodes==27) nbNodes=8; // quadratic hexahedron
+ else if(nbNodes==27) nbNodes=8; // tri-quadratic hexahedron
else return aQuality;
}
case 5:{
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )};
- aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4]));
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )};
case 6:{
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )};
- aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4]));
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )};
break;
}
case 8:{
+
+ return hexQualityByHomardMethod( P ); // bos #23982
+
+
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
- aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4]));
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )};
case 12:
{
gp_XYZ aXYZ[8] = {P( 1 ),P( 2 ),P( 4 ),P( 5 ),P( 7 ),P( 8 ),P( 10 ),P( 11 )};
- aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
+ aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8]));
}
{
gp_XYZ aXYZ[8] = {P( 2 ),P( 3 ),P( 5 ),P( 6 ),P( 8 ),P( 9 ),P( 11 ),P( 12 )};
*/
//================================================================================
+bool Warping::IsApplicable( const SMDS_MeshElement* element ) const
+{
+ return NumericalFunctor::IsApplicable( element ) && element->NbNodes() == 4;
+}
+
double Warping::GetValue( const TSequenceOfXYZ& P )
{
if ( P.size() != 4 )
*/
//================================================================================
+bool Taper::IsApplicable( const SMDS_MeshElement* element ) const
+{
+ return ( NumericalFunctor::IsApplicable( element ) && element->NbNodes() == 4 );
+}
+
double Taper::GetValue( const TSequenceOfXYZ& P )
{
if ( P.size() != 4 )
return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0. : v1.Angle( v2 );
}
+bool Skew::IsApplicable( const SMDS_MeshElement* element ) const
+{
+ return ( NumericalFunctor::IsApplicable( element ) && element->NbNodes() <= 4 );
+}
+
double Skew::GetValue( const TSequenceOfXYZ& P )
{
if ( P.size() != 3 && P.size() != 4 )
return SMDSAbs_Edge;
}
+//================================================================================
+/*
+ Class : Length3D
+ Description : Functor for calculating minimal length of element edge
+*/
+//================================================================================
+
+Length3D::Length3D():
+ Length2D ( SMDSAbs_Volume )
+{
+}
+
//================================================================================
/*
Class : Length2D
- Description : Functor for calculating minimal length of edge
+ Description : Functor for calculating minimal length of element edge
*/
//================================================================================
+Length2D::Length2D( SMDSAbs_ElementType type ):
+ myType ( type )
+{
+}
+
+bool Length2D::IsApplicable( const SMDS_MeshElement* element ) const
+{
+ return ( NumericalFunctor::IsApplicable( element ) &&
+ element->GetEntityType() != SMDSEntity_Polyhedra );
+}
+
double Length2D::GetValue( const TSequenceOfXYZ& P )
{
double aVal = 0;
SMDSAbs_ElementType Length2D::GetType() const
{
- return SMDSAbs_Face;
+ return myType;
}
Length2D::Value::Value(double theLength,long thePntId1, long thePntId2):
void Length2D::GetValues(TValues& theValues)
{
- TValues aValues;
- for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); )
+ if ( myType == SMDSAbs_Face )
{
- const SMDS_MeshFace* anElem = anIter->next();
- if ( anElem->IsQuadratic() )
+ for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); )
{
- // use special nodes iterator
- SMDS_NodeIteratorPtr anIter = anElem->interlacedNodesIterator();
- long aNodeId[4] = { 0,0,0,0 };
- gp_Pnt P[4];
-
- double aLength = 0;
- if ( anIter->more() )
- {
- const SMDS_MeshNode* aNode = anIter->next();
- P[0] = P[1] = SMESH_NodeXYZ( aNode );
- aNodeId[0] = aNodeId[1] = aNode->GetID();
- aLength = 0;
- }
- for ( ; anIter->more(); )
+ const SMDS_MeshFace* anElem = anIter->next();
+ if ( anElem->IsQuadratic() )
{
- const SMDS_MeshNode* N1 = anIter->next();
- P[2] = SMESH_NodeXYZ( N1 );
- aNodeId[2] = N1->GetID();
- aLength = P[1].Distance(P[2]);
- if(!anIter->more()) break;
- const SMDS_MeshNode* N2 = anIter->next();
- P[3] = SMESH_NodeXYZ( N2 );
- aNodeId[3] = N2->GetID();
- aLength += P[2].Distance(P[3]);
+ // use special nodes iterator
+ SMDS_NodeIteratorPtr anIter = anElem->interlacedNodesIterator();
+ long aNodeId[4] = { 0,0,0,0 };
+ gp_Pnt P[4];
+
+ double aLength = 0;
+ if ( anIter->more() )
+ {
+ const SMDS_MeshNode* aNode = anIter->next();
+ P[0] = P[1] = SMESH_NodeXYZ( aNode );
+ aNodeId[0] = aNodeId[1] = aNode->GetID();
+ aLength = 0;
+ }
+ for ( ; anIter->more(); )
+ {
+ const SMDS_MeshNode* N1 = anIter->next();
+ P[2] = SMESH_NodeXYZ( N1 );
+ aNodeId[2] = N1->GetID();
+ aLength = P[1].Distance(P[2]);
+ if(!anIter->more()) break;
+ const SMDS_MeshNode* N2 = anIter->next();
+ P[3] = SMESH_NodeXYZ( N2 );
+ aNodeId[3] = N2->GetID();
+ aLength += P[2].Distance(P[3]);
+ Value aValue1(aLength,aNodeId[1],aNodeId[2]);
+ Value aValue2(aLength,aNodeId[2],aNodeId[3]);
+ P[1] = P[3];
+ aNodeId[1] = aNodeId[3];
+ theValues.insert(aValue1);
+ theValues.insert(aValue2);
+ }
+ aLength += P[2].Distance(P[0]);
Value aValue1(aLength,aNodeId[1],aNodeId[2]);
- Value aValue2(aLength,aNodeId[2],aNodeId[3]);
- P[1] = P[3];
- aNodeId[1] = aNodeId[3];
+ Value aValue2(aLength,aNodeId[2],aNodeId[0]);
theValues.insert(aValue1);
theValues.insert(aValue2);
}
- aLength += P[2].Distance(P[0]);
- Value aValue1(aLength,aNodeId[1],aNodeId[2]);
- Value aValue2(aLength,aNodeId[2],aNodeId[0]);
- theValues.insert(aValue1);
- theValues.insert(aValue2);
- }
- else {
- SMDS_NodeIteratorPtr aNodesIter = anElem->nodeIterator();
- long aNodeId[2] = {0,0};
- gp_Pnt P[3];
+ else {
+ SMDS_NodeIteratorPtr aNodesIter = anElem->nodeIterator();
+ long aNodeId[2] = {0,0};
+ gp_Pnt P[3];
+
+ double aLength;
+ const SMDS_MeshElement* aNode;
+ if ( aNodesIter->more())
+ {
+ aNode = aNodesIter->next();
+ P[0] = P[1] = SMESH_NodeXYZ( aNode );
+ aNodeId[0] = aNodeId[1] = aNode->GetID();
+ aLength = 0;
+ }
+ for( ; aNodesIter->more(); )
+ {
+ aNode = aNodesIter->next();
+ long anId = aNode->GetID();
- double aLength;
- const SMDS_MeshElement* aNode;
- if ( aNodesIter->more())
- {
- aNode = aNodesIter->next();
- P[0] = P[1] = SMESH_NodeXYZ( aNode );
- aNodeId[0] = aNodeId[1] = aNode->GetID();
- aLength = 0;
- }
- for( ; aNodesIter->more(); )
- {
- aNode = aNodesIter->next();
- long anId = aNode->GetID();
+ P[2] = SMESH_NodeXYZ( aNode );
- P[2] = SMESH_NodeXYZ( aNode );
+ aLength = P[1].Distance(P[2]);
+
+ Value aValue(aLength,aNodeId[1],anId);
+ aNodeId[1] = anId;
+ P[1] = P[2];
+ theValues.insert(aValue);
+ }
- aLength = P[1].Distance(P[2]);
+ aLength = P[0].Distance(P[1]);
- Value aValue(aLength,aNodeId[1],anId);
- aNodeId[1] = anId;
- P[1] = P[2];
+ Value aValue(aLength,aNodeId[0],aNodeId[1]);
theValues.insert(aValue);
}
-
- aLength = P[0].Distance(P[1]);
-
- Value aValue(aLength,aNodeId[0],aNodeId[1]);
- theValues.insert(aValue);
}
}
+ else
+ {
+ // not implemented
+ }
}
//================================================================================
*/
//================================================================================
-double MultiConnection::GetValue( const TSequenceOfXYZ& P )
+double MultiConnection::GetValue( const TSequenceOfXYZ& /*P*/ )
{
return 0;
}
*/
//================================================================================
-double MultiConnection2D::GetValue( const TSequenceOfXYZ& P )
+double MultiConnection2D::GetValue( const TSequenceOfXYZ& /*P*/ )
{
return 0;
}
return false;
SMDS_VolumeTool vTool( myMesh->FindElement( theId ));
- return !vTool.IsForward();
+
+ bool isOk = true;
+ if ( vTool.IsPoly() )
+ {
+ isOk = true;
+ for ( int i = 0; i < vTool.NbFaces() && isOk; ++i )
+ isOk = vTool.IsFaceExternal( i );
+ }
+ else
+ {
+ isOk = vTool.IsForward();
+ }
+ return !isOk;
}
SMDSAbs_ElementType BadOrientedVolume::GetType() const
return SMDSAbs_Node;
}
+void CoincidentNodes::SetTolerance( const double theToler )
+{
+ if ( myToler != theToler )
+ {
+ SetMesh(0);
+ myToler = theToler;
+ }
+}
+
void CoincidentNodes::SetMesh( const SMDS_Mesh* theMesh )
{
myMeshModifTracer.SetMesh( theMesh );
bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId )
{
- TColStd_MapOfInteger aMap;
- for ( int i = 0; i < 2; i++ )
+ SMDS_ElemIteratorPtr anElemIter = theNodes[ 0 ]->GetInverseElementIterator(SMDSAbs_Face);
+ while( anElemIter->more() )
{
- SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator(SMDSAbs_Face);
- while( anElemIter->more() )
+ if ( const SMDS_MeshElement* anElem = anElemIter->next())
{
- if ( const SMDS_MeshElement* anElem = anElemIter->next())
- {
- const int anId = anElem->GetID();
- if ( anId != theFaceId && !aMap.Add( anId ))
- return false;
- }
+ const int anId = anElem->GetID();
+ if ( anId != theFaceId && anElem->GetNodeIndex( theNodes[1] ) >= 0 )
+ return false;
}
}
return true;
SMDS_MeshElement* aVol = (SMDS_MeshElement*)volItr->next();
TItrMapOfVolume itr = mapOfVol.insert( std::make_pair( aVol, 0 )).first;
(*itr).second++;
- }
+ }
}
int nbVol = 0;
TItrMapOfVolume volItr = mapOfVol.begin();
TItrMapOfVolume volEnd = mapOfVol.end();
for ( ; volItr != volEnd; ++volItr )
if ( (*volItr).second >= nbNode )
- nbVol++;
+ nbVol++;
// face is not free if number of volumes constructed on their nodes more than one
return (nbVol < 2);
}
double dot = v1 * v2; // cos * |v1| * |v2|
double l1 = v1.SquareMagnitude();
double l2 = v2.SquareMagnitude();
- return (( dot * cos >= 0 ) &&
+ return (( dot * cos >= 0 ) &&
( dot * dot ) / l1 / l2 >= ( cos * cos ));
}
}
TopAbs_ShapeEnum ShapeType() const { return myShape.ShapeType(); }
const TopoDS_Shape& Shape() const { return myShape; }
const Bnd_B3d* GetBndBox() const { return & myBox; }
+ double Tolerance() const { return myTol; }
bool IsChecked() { return myFlags & theIsCheckedFlag; }
bool IsSetFlag( int flag ) const { return myFlags & flag; }
void SetChecked( bool is ) { is ? SetFlag( theIsCheckedFlag ) : UnsetFlag( theIsCheckedFlag ); }
bool isOutOfFace (const gp_Pnt& p);
bool isOutOfEdge (const gp_Pnt& p);
bool isOutOfVertex(const gp_Pnt& p);
+ bool isOutOfNone (const gp_Pnt& /*p*/) { return true; }
bool isBox (const TopoDS_Shape& s);
+ TopoDS_Shape prepareSolid( const TopoDS_Shape& theSolid );
+
bool (Classifier::* myIsOutFun)(const gp_Pnt& p);
BRepClass3d_SolidClassifier* mySolidClfr; // ptr because of a run-time forbidden copy-constructor
Bnd_B3d myBox;
std::vector< ElementsOnShape::Classifier >& cls );
void GetClassifiersAtPoint( const gp_XYZ& p,
std::vector< ElementsOnShape::Classifier* >& classifiers );
+ size_t GetSize();
+
protected:
OctreeClassifier() {}
SMESH_Octree* newChild() const { return new OctreeClassifier; }
Predicate* ElementsOnShape::clone() const
{
+ size_t size = sizeof( *this );
+ if ( myOctree )
+ size += myOctree->GetSize();
+ if ( !myClassifiers.empty() )
+ size += sizeof( myClassifiers[0] ) * myClassifiers.size();
+ if ( !myWorkClassifiers.empty() )
+ size += sizeof( myWorkClassifiers[0] ) * myWorkClassifiers.size();
+ if ( size > 1e+9 ) // 1G
+ {
+#ifdef _DEBUG_
+ std::cout << "Avoid ElementsOnShape::clone(), too large: " << size << " bytes " << std::endl;
+#endif
+ return 0;
+ }
+
ElementsOnShape* cln = new ElementsOnShape();
cln->SetAllNodes ( myAllNodesFlag );
cln->SetTolerance( myToler );
for ( size_t i = 0; i < myClassifiers.size(); ++i )
myWorkClassifiers[ i ] = & myClassifiers[ i ];
myOctree = new OctreeClassifier( myWorkClassifiers );
+
+ SMESHUtils::FreeVector( myWorkClassifiers );
}
- SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator();
- while (aNodeItr->more() && (isSatisfy == myAllNodesFlag))
+ for ( int i = 0, nb = elem->NbNodes(); i < nb && (isSatisfy == myAllNodesFlag); ++i )
{
- SMESH_TNodeXYZ aPnt( aNodeItr->next() );
+ SMESH_TNodeXYZ aPnt( elem->GetNode( i ));
centerXYZ += aPnt;
isNodeOut = true;
centerXYZ /= elem->NbNodes();
isSatisfy = false;
if ( myOctree )
+ {
+ myWorkClassifiers.clear();
+ myOctree->GetClassifiersAtPoint( centerXYZ, myWorkClassifiers );
for ( size_t i = 0; i < myWorkClassifiers.size() && !isSatisfy; ++i )
isSatisfy = ! myWorkClassifiers[i]->IsOut( centerXYZ );
+ }
else
+ {
for ( size_t i = 0; i < myClassifiers.size() && !isSatisfy; ++i )
isSatisfy = ! myClassifiers[i].IsOut( centerXYZ );
+ }
}
return isSatisfy;
{
isNodeOut = false;
if ( okShape )
- *okShape = myWorkClassifiers[i]->Shape();
+ *okShape = myClassifiers[i].Shape();
break;
}
}
}
else
{
- mySolidClfr = new BRepClass3d_SolidClassifier(theShape);
+ mySolidClfr = new BRepClass3d_SolidClassifier( prepareSolid( theShape ));
myIsOutFun = & ElementsOnShape::Classifier::isOutOfSolid;
}
break;
{
Standard_Real u1,u2,v1,v2;
Handle(Geom_Surface) surf = BRep_Tool::Surface( TopoDS::Face( theShape ));
- surf->Bounds( u1,u2,v1,v2 );
- myProjFace.Init(surf, u1,u2, v1,v2, myTol );
- myIsOutFun = & ElementsOnShape::Classifier::isOutOfFace;
+ if ( surf.IsNull() )
+ myIsOutFun = & ElementsOnShape::Classifier::isOutOfNone;
+ else
+ {
+ surf->Bounds( u1,u2,v1,v2 );
+ myProjFace.Init(surf, u1,u2, v1,v2, myTol );
+ myIsOutFun = & ElementsOnShape::Classifier::isOutOfFace;
+ }
break;
}
case TopAbs_EDGE:
{
Standard_Real u1, u2;
Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( theShape ), u1, u2);
- myProjEdge.Init(curve, u1, u2);
- myIsOutFun = & ElementsOnShape::Classifier::isOutOfEdge;
+ if ( curve.IsNull() )
+ myIsOutFun = & ElementsOnShape::Classifier::isOutOfNone;
+ else
+ {
+ myProjEdge.Init(curve, u1, u2);
+ myIsOutFun = & ElementsOnShape::Classifier::isOutOfEdge;
+ }
break;
}
case TopAbs_VERTEX:
else
{
Bnd_Box box;
- BRepBndLib::Add( myShape, box );
+ if ( myShape.ShapeType() == TopAbs_FACE )
+ {
+ BRepAdaptor_Surface SA( TopoDS::Face( myShape ), /*useBoundaries=*/false );
+ if ( SA.GetType() == GeomAbs_BSplineSurface )
+ BRepBndLib::AddOptimal( myShape, box,
+ /*useTriangulation=*/true, /*useShapeTolerance=*/true );
+ }
+ if ( box.IsVoid() )
+ BRepBndLib::Add( myShape, box );
myBox.Clear();
myBox.Add( box.CornerMin() );
myBox.Add( box.CornerMax() );
delete mySolidClfr; mySolidClfr = 0;
}
-bool ElementsOnShape::Classifier::isOutOfSolid (const gp_Pnt& p)
+TopoDS_Shape ElementsOnShape::Classifier::prepareSolid( const TopoDS_Shape& theSolid )
+{
+ // try to limit tolerance of theSolid down to myTol (issue #19026)
+
+ // check if tolerance of theSolid is more than myTol
+ bool tolIsOk = true; // max tolerance is at VERTEXes
+ for ( TopExp_Explorer exp( theSolid, TopAbs_VERTEX ); exp.More() && tolIsOk; exp.Next() )
+ tolIsOk = ( myTol >= BRep_Tool::Tolerance( TopoDS::Vertex( exp.Current() )));
+ if ( tolIsOk )
+ return theSolid;
+
+ // make a copy to prevent the original shape from changes
+ TopoDS_Shape resultShape = BRepBuilderAPI_Copy( theSolid );
+
+ if ( !GEOMUtils::FixShapeTolerance( resultShape, TopAbs_SHAPE, myTol ))
+ return theSolid;
+ return resultShape;
+}
+
+bool ElementsOnShape::Classifier::isOutOfSolid( const gp_Pnt& p )
{
if ( isOutOfBox( p )) return true;
mySolidClfr->Perform( p, myTol );
return ( mySolidClfr->State() != TopAbs_IN && mySolidClfr->State() != TopAbs_ON );
}
-bool ElementsOnShape::Classifier::isOutOfBox (const gp_Pnt& p)
+bool ElementsOnShape::Classifier::isOutOfBox( const gp_Pnt& p )
{
return myBox.IsOut( p.XYZ() );
}
-bool ElementsOnShape::Classifier::isOutOfFace (const gp_Pnt& p)
+bool ElementsOnShape::Classifier::isOutOfFace( const gp_Pnt& p )
{
if ( isOutOfBox( p )) return true;
myProjFace.Perform( p );
return true;
}
-bool ElementsOnShape::Classifier::isOutOfEdge (const gp_Pnt& p)
+bool ElementsOnShape::Classifier::isOutOfEdge( const gp_Pnt& p )
{
if ( isOutOfBox( p )) return true;
myProjEdge.Perform( p );
return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol );
}
-bool ElementsOnShape::Classifier::isOutOfVertex(const gp_Pnt& p)
+bool ElementsOnShape::Classifier::isOutOfVertex( const gp_Pnt& p )
{
return ( myVertexXYZ.Distance( p ) > myTol );
}
-bool ElementsOnShape::Classifier::isBox (const TopoDS_Shape& theShape)
+bool ElementsOnShape::Classifier::isBox(const TopoDS_Shape& theShape )
{
TopTools_IndexedMapOfShape vMap;
TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap );
}
}
+size_t ElementsOnShape::OctreeClassifier::GetSize()
+{
+ size_t res = sizeof( *this );
+ if ( !myClassifiers.empty() )
+ res += sizeof( myClassifiers[0] ) * myClassifiers.size();
+
+ if ( !isLeaf() )
+ for (int i = 0; i < nbChildren(); i++)
+ res += ((OctreeClassifier*) myChildren[i])->GetSize();
+
+ return res;
+}
+
void ElementsOnShape::OctreeClassifier::buildChildrenData()
{
// distribute myClassifiers among myChildren
for ( int i = 0; i < nbChildren(); i++ )
{
OctreeClassifier* child = static_cast<OctreeClassifier*>( myChildren[ i ]);
- child->myIsLeaf = ( child->myClassifiers.size() <= 5 );
+ child->myIsLeaf = ( child->myClassifiers.size() <= 5 ||
+ child->maxSize() < child->myClassifiers[0]->Tolerance() );
}
}
Predicate* BelongToGeom::clone() const
{
- BelongToGeom* cln = new BelongToGeom( *this );
- cln->myElementsOnShapePtr.reset( static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ));
+ BelongToGeom* cln = 0;
+ if ( myElementsOnShapePtr )
+ if ( ElementsOnShape* eos = static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ))
+ {
+ cln = new BelongToGeom( *this );
+ cln->myElementsOnShapePtr.reset( eos );
+ }
return cln;
}
myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
init();
}
+ if ( myElementsOnShapePtr )
+ myElementsOnShapePtr->SetMesh( myMeshDS );
}
void BelongToGeom::SetGeom( const TopoDS_Shape& theShape )
Predicate* LyingOnGeom::clone() const
{
- LyingOnGeom* cln = new LyingOnGeom( *this );
- cln->myElementsOnShapePtr.reset( static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ));
+ LyingOnGeom* cln = 0;
+ if ( myElementsOnShapePtr )
+ if ( ElementsOnShape* eos = static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ))
+ {
+ cln = new LyingOnGeom( *this );
+ cln->myElementsOnShapePtr.reset( eos );
+ }
return cln;
}
myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
init();
}
+ if ( myElementsOnShapePtr )
+ myElementsOnShapePtr->SetMesh( myMeshDS );
}
void LyingOnGeom::SetGeom( const TopoDS_Shape& theShape )
