-// Copyright (C) 2007-2016 CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2024 CEA, EDF, 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_ControlsDef.hxx"
#include "SMDS_BallElement.hxx"
+#include "SMDS_FacePosition.hxx"
#include "SMDS_Iterator.hxx"
#include "SMDS_Mesh.hxx"
#include "SMDS_MeshElement.hxx"
#include "SMDS_MeshNode.hxx"
-#include "SMDS_QuadraticEdge.hxx"
-#include "SMDS_QuadraticFaceOfNodes.hxx"
#include "SMDS_VolumeTool.hxx"
#include "SMESHDS_GroupBase.hxx"
#include "SMESHDS_GroupOnFilter.hxx"
#include "SMESHDS_Mesh.hxx"
#include "SMESH_MeshAlgos.hxx"
#include "SMESH_OctreeNode.hxx"
+#include "SMESH_Comment.hxx"
+#include <GEOMUtils.hxx>
#include <Basics_Utils.hxx>
#include <BRepAdaptor_Surface.hxx>
#include <BRepBndLib.hxx>
#include <BRepBuilderAPI_Copy.hxx>
+#include <BRepClass3d_SolidClassifier.hxx>
#include <BRepClass_FaceClassifier.hxx>
#include <BRep_Tool.hxx>
+#include <GeomLib_IsPlanarSurface.hxx>
#include <Geom_CylindricalSurface.hxx>
#include <Geom_Plane.hxx>
#include <Geom_Surface.hxx>
#include <NCollection_Map.hxx>
#include <Precision.hxx>
+#include <ShapeAnalysis_Surface.hxx>
#include <TColStd_MapIteratorOfMapOfInteger.hxx>
#include <TColStd_MapOfInteger.hxx>
#include <TColStd_SequenceOfAsciiString.hxx>
v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
}
+ inline double getCos2( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
+ {
+ gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
+ double dot = v1 * v2, len1 = v1.SquareMagnitude(), len2 = v2.SquareMagnitude();
+
+ return ( dot < 0 || len1 < gp::Resolution() || len2 < gp::Resolution() ? -1 :
+ dot * dot / len1 / len2 );
+ }
+
inline double getArea( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
{
gp_Vec aVec1( P2 - P1 );
return aDist;
}
- int getNbMultiConnection( const SMDS_Mesh* theMesh, const int theId )
+ int getNbMultiConnection( const SMDS_Mesh* theMesh, const smIdType theId )
{
if ( theMesh == 0 )
return 0;
// 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 );
myMesh = theMesh;
}
-bool NumericalFunctor::GetPoints(const int theId,
+bool NumericalFunctor::GetPoints(const smIdType theId,
TSequenceOfXYZ& theRes ) const
{
theRes.clear();
return false;
const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
- if ( !anElem || anElem->GetType() != this->GetType() )
+ if ( !IsApplicable( anElem ))
return false;
return GetPoints( anElem, theRes );
theRes.setElement( anElem );
// Get nodes of the element
- SMDS_ElemIteratorPtr anIter;
-
- if ( anElem->IsQuadratic() ) {
- switch ( anElem->GetType() ) {
- case SMDSAbs_Edge:
- anIter = dynamic_cast<const SMDS_VtkEdge*>
- (anElem)->interlacedNodesElemIterator();
- break;
- case SMDSAbs_Face:
- anIter = dynamic_cast<const SMDS_VtkFace*>
- (anElem)->interlacedNodesElemIterator();
- break;
- default:
- anIter = anElem->nodesIterator();
- }
- }
- else {
- anIter = anElem->nodesIterator();
- }
-
+ SMDS_NodeIteratorPtr anIter= anElem->interlacedNodesIterator();
if ( anIter ) {
- double xyz[3];
+ SMESH_NodeXYZ p;
while( anIter->more() ) {
- if ( const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>( anIter->next() ))
- {
- aNode->GetXYZ( xyz );
- theRes.push_back( gp_XYZ( xyz[0], xyz[1], xyz[2] ));
- }
+ if ( p.Set( anIter->next() ))
+ theRes.push_back( p );
}
}
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
*/
//================================================================================
-void NumericalFunctor::GetHistogram(int nbIntervals,
- std::vector<int>& nbEvents,
- std::vector<double>& funValues,
- const std::vector<int>& elements,
- const double* minmax,
- const bool isLogarithmic)
+void NumericalFunctor::GetHistogram(int nbIntervals,
+ std::vector<int>& nbEvents,
+ std::vector<double>& funValues,
+ const std::vector<smIdType>& elements,
+ const double* minmax,
+ const bool isLogarithmic)
{
if ( nbIntervals < 1 ||
!myMesh ||
}
else
{
- std::vector<int>::const_iterator id = elements.begin();
+ std::vector<smIdType>::const_iterator id = elements.begin();
for ( ; id != elements.end(); ++id )
values.insert( GetValue( *id ));
}
// }
// }
// { // polygons
-
+
// }
if( myPrecision >= 0 )
aVal = Max(aVal,Max(L7,L8));
break;
}
- case SMDSEntity_Quad_Penta: { // quadratic pentas
+ case SMDSEntity_Quad_Penta:
+ case SMDSEntity_BiQuad_Penta: { // quadratic pentas
double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
{
- double aMin;
-
- if (P.size() <3)
+ if ( P.size() < 3 )
return 0.;
- aMin = getAngle(P( P.size() ), P( 1 ), P( 2 ));
- aMin = Min(aMin,getAngle(P( P.size()-1 ), P( P.size() ), P( 1 )));
+ double aMaxCos2;
+
+ aMaxCos2 = getCos2( P( P.size() ), P( 1 ), P( 2 ));
+ aMaxCos2 = Max( aMaxCos2, getCos2( P( P.size()-1 ), P( P.size() ), P( 1 )));
for ( size_t i = 2; i < P.size(); i++ )
{
- double A0 = getAngle( P( i-1 ), P( i ), P( i+1 ) );
- aMin = Min(aMin,A0);
+ double A0 = getCos2( P( i-1 ), P( i ), P( i+1 ) );
+ aMaxCos2 = Max( aMaxCos2, A0 );
}
+ if ( aMaxCos2 < 0 )
+ return 0; // all nodes coincide
- return aMin * 180.0 / M_PI;
+ double cos = sqrt( aMaxCos2 );
+ if ( cos >= 1 ) return 0;
+ return acos( cos ) * 180.0 / M_PI;
}
double MinimumAngle::GetBadRate( double Value, int nbNodes ) const
{
double aVal = 0;
myCurrElement = myMesh->FindElement( theId );
- if ( myCurrElement && myCurrElement->GetVtkType() == VTK_QUAD )
- {
- // issue 21723
- vtkUnstructuredGrid* grid = SMDS_Mesh::_meshList[myCurrElement->getMeshId()]->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;
}
if ( nbNodes == 3 ) {
// Compute lengths of the sides
- std::vector< double > aLen (nbNodes);
- for ( int i = 0; i < nbNodes - 1; i++ )
- aLen[ i ] = getDistance( P( i + 1 ), P( i + 2 ) );
- aLen[ nbNodes - 1 ] = getDistance( P( 1 ), P( nbNodes ) );
+ double aLen1 = getDistance( P( 1 ), P( 2 ));
+ double aLen2 = getDistance( P( 2 ), P( 3 ));
+ double aLen3 = getDistance( P( 3 ), P( 1 ));
// Q = alfa * h * p / S, where
//
// alfa = sqrt( 3 ) / 6
// h - length of the longest edge
// p - half perimeter
// S - triangle surface
- const double alfa = sqrt( 3. ) / 6.;
- double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
- double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.;
- double anArea = getArea( P( 1 ), P( 2 ), P( 3 ) );
+ const double alfa = sqrt( 3. ) / 6.;
+ double maxLen = Max( aLen1, Max( aLen2, aLen3 ));
+ double half_perimeter = ( aLen1 + aLen2 + aLen3 ) / 2.;
+ double anArea = getArea( P( 1 ), P( 2 ), P( 3 ));
if ( anArea <= theEps )
return theInf;
return alfa * maxLen * half_perimeter / anArea;
}
else if ( nbNodes == 6 ) { // quadratic triangles
// Compute lengths of the sides
- std::vector< double > aLen (3);
- aLen[0] = getDistance( P(1), P(3) );
- aLen[1] = getDistance( P(3), P(5) );
- aLen[2] = getDistance( P(5), P(1) );
- // Q = alfa * h * p / S, where
- //
- // alfa = sqrt( 3 ) / 6
- // h - length of the longest edge
- // p - half perimeter
- // S - triangle surface
- const double alfa = sqrt( 3. ) / 6.;
- double maxLen = Max( aLen[ 0 ], Max( aLen[ 1 ], aLen[ 2 ] ) );
- double half_perimeter = ( aLen[0] + aLen[1] + aLen[2] ) / 2.;
- double anArea = getArea( P(1), P(3), P(5) );
+ double aLen1 = getDistance( P( 1 ), P( 3 ));
+ double aLen2 = getDistance( P( 3 ), P( 5 ));
+ double aLen3 = getDistance( P( 5 ), P( 1 ));
+ // algo same as for the linear triangle
+ const double alfa = sqrt( 3. ) / 6.;
+ double maxLen = Max( aLen1, Max( aLen2, aLen3 ));
+ double half_perimeter = ( aLen1 + aLen2 + aLen3 ) / 2.;
+ double anArea = getArea( P( 1 ), P( 3 ), P( 5 ));
if ( anArea <= theEps )
return theInf;
return alfa * maxLen * half_perimeter / anArea;
}
else if( nbNodes == 4 ) { // quadrangle
// Compute lengths of the sides
- std::vector< double > aLen (4);
+ double aLen[4];
aLen[0] = getDistance( P(1), P(2) );
aLen[1] = getDistance( P(2), P(3) );
aLen[2] = getDistance( P(3), P(4) );
aLen[3] = getDistance( P(4), P(1) );
// Compute lengths of the diagonals
- std::vector< double > aDia (2);
+ double aDia[2];
aDia[0] = getDistance( P(1), P(3) );
aDia[1] = getDistance( P(2), P(4) );
// Compute areas of all triangles which can be built
// taking three nodes of the quadrangle
- std::vector< double > anArea (4);
+ double anArea[4];
anArea[0] = getArea( P(1), P(2), P(3) );
anArea[1] = getArea( P(1), P(2), P(4) );
anArea[2] = getArea( P(1), P(3), P(4) );
//
// 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;
}
else if( nbNodes == 8 || nbNodes == 9 ) { // nbNodes==8 - quadratic quadrangle
// Compute lengths of the sides
- std::vector< double > aLen (4);
+ double aLen[4];
aLen[0] = getDistance( P(1), P(3) );
aLen[1] = getDistance( P(3), P(5) );
aLen[2] = getDistance( P(5), P(7) );
aLen[3] = getDistance( P(7), P(1) );
// Compute lengths of the diagonals
- std::vector< double > aDia (2);
+ double aDia[2];
aDia[0] = getDistance( P(1), P(5) );
aDia[1] = getDistance( P(3), P(7) );
// Compute areas of all triangles which can be built
// taking three nodes of the quadrangle
- std::vector< double > anArea (4);
+ double anArea[4];
anArea[0] = getArea( P(1), P(3), P(5) );
anArea[1] = getArea( P(1), P(3), P(7) );
anArea[2] = getArea( P(1), P(5), P(7) );
//
// 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 )
// Action from CoTech | ACTION 31.3:
// EURIWARE BO: Homogenize the formulas used to calculate the Controls in SMESH to fit with
// those of ParaView. The library used by ParaView for those calculations can be reused in SMESH.
- vtkUnstructuredGrid* grid = SMDS_Mesh::_meshList[myCurrElement->getMeshId()]->getGrid();
- if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->getVtkId() ))
+ vtkUnstructuredGrid* grid = const_cast<SMDS_Mesh*>( myMesh )->GetGrid();
+ if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->GetVtkID() ))
aVal = Round( vtkMeshQuality::TetAspectRatio( avtkCell ));
}
else
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( myCurrElement->IsQuadratic() ) {
+ 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 )};
} // switch(nbNodes)
if ( nbNodes > 4 ) {
- // avaluate aspect ratio of quadranle faces
+ // evaluate aspect ratio of quadrangle faces
AspectRatio aspect2D;
SMDS_VolumeTool::VolumeType type = SMDS_VolumeTool::GetType( nbNodes );
int nbFaces = SMDS_VolumeTool::NbFaces( type );
if ( SMDS_VolumeTool::NbFaceNodes( type, i ) != 4 )
continue;
const int* pInd = SMDS_VolumeTool::GetFaceNodesIndices( type, i, true );
- for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadranle face
+ for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadrangle face
points( p + 1 ) = P( pInd[ p ] + 1 );
aQuality = std::max( aQuality, aspect2D.GetValue( points ));
}
*/
//================================================================================
-double Warping::GetValue( const TSequenceOfXYZ& P )
+bool Warping::IsApplicable( const SMDS_MeshElement* element ) const
{
- if ( P.size() != 4 )
- return 0;
-
- gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4.;
-
- double A1 = ComputeA( P( 1 ), P( 2 ), P( 3 ), G );
- double A2 = ComputeA( P( 2 ), P( 3 ), P( 4 ), G );
- double A3 = ComputeA( P( 3 ), P( 4 ), P( 1 ), G );
- double A4 = ComputeA( P( 4 ), P( 1 ), P( 2 ), G );
-
- double val = Max( Max( A1, A2 ), Max( A3, A4 ) );
-
- const double eps = 0.1; // val is in degrees
+ return NumericalFunctor::IsApplicable( element ) && element->NbNodes() == 4;
+}
- return val < eps ? 0. : val;
+double Warping::GetValue( const TSequenceOfXYZ& P )
+{
+ return ComputeValue(P);
}
double Warping::ComputeA( const gp_XYZ& thePnt1,
return asin( fabs( H / L ) ) * 180. / M_PI;
}
+double Warping::ComputeValue(const TSequenceOfXYZ& thePoints) const
+{
+ if (thePoints.size() != 4)
+ return 0;
+
+ gp_XYZ G = (thePoints(1) + thePoints(2) + thePoints(3) + thePoints(4)) / 4.;
+
+ double A1 = ComputeA(thePoints(1), thePoints(2), thePoints(3), G);
+ double A2 = ComputeA(thePoints(2), thePoints(3), thePoints(4), G);
+ double A3 = ComputeA(thePoints(3), thePoints(4), thePoints(1), G);
+ double A4 = ComputeA(thePoints(4), thePoints(1), thePoints(2), G);
+
+ double val = Max(Max(A1, A2), Max(A3, A4));
+
+ const double eps = 0.1; // val is in degrees
+
+ return val < eps ? 0. : val;
+}
+
double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const
{
// the warp is in the range [0.0,PI/2]
}
+//================================================================================
+/*
+ Class : Warping3D
+ Description : Functor for calculating warping
+*/
+//================================================================================
+
+bool Warping3D::IsApplicable(const SMDS_MeshElement* element) const
+{
+ return NumericalFunctor::IsApplicable(element);//&& element->NbNodes() == 4;
+}
+
+double Warping3D::GetValue(long theId)
+{
+ double aVal = 0;
+ myCurrElement = myMesh->FindElement(theId);
+ if (myCurrElement)
+ {
+ WValues aValues;
+ ProcessVolumeELement(aValues);
+ for (const auto& aValue: aValues)
+ {
+ aVal = Max(aVal, aValue.myWarp);
+ }
+ }
+ return aVal;
+}
+
+double Warping3D::GetValue(const TSequenceOfXYZ& P)
+{
+ return ComputeValue(P);
+}
+
+SMDSAbs_ElementType Warping3D::GetType() const
+{
+ return SMDSAbs_Volume;
+}
+
+bool Warping3D::Value::operator<(const Warping3D::Value& x) const
+{
+ if (myPntIds.size() != x.myPntIds.size())
+ return myPntIds.size() < x.myPntIds.size();
+
+ for (int anInd = 0; anInd < myPntIds.size(); ++anInd)
+ if (myPntIds[anInd] != x.myPntIds[anInd])
+ return myPntIds[anInd] != x.myPntIds[anInd];
+
+ return false;
+}
+
+// Compute value on each face of volume
+void Warping3D::ProcessVolumeELement(WValues& theValues)
+{
+ SMDS_VolumeTool aVTool(myCurrElement);
+ double aCoord[3];
+ for (int aFaceID = 0; aFaceID < aVTool.NbFaces(); ++aFaceID)
+ {
+ TSequenceOfXYZ aPoints;
+ std::set<const SMDS_MeshNode*> aNodes;
+ std::vector<long> aNodeIds;
+ const SMDS_MeshNode** aNodesPtr = aVTool.GetFaceNodes(aFaceID);
+
+ if (aNodesPtr)
+ {
+ for (int i = 0; i < aVTool.NbFaceNodes(aFaceID); ++i)
+ {
+ aNodesPtr[i]->GetXYZ(aCoord);
+ aPoints.push_back(gp_XYZ{ aCoord[0], aCoord[1], aCoord[2] });
+ aNodeIds.push_back(aNodesPtr[i]->GetID());
+ }
+ double aWarp = GetValue(aPoints);
+ Value aVal{ aWarp, aNodeIds };
+
+ theValues.push_back(aVal);
+ }
+ }
+}
+
+void Warping3D::GetValues(WValues& theValues)
+{
+ for (SMDS_VolumeIteratorPtr anIter = myMesh->volumesIterator(); anIter->more(); )
+ {
+ myCurrElement = anIter->next();
+ ProcessVolumeELement(theValues);
+ }
+}
+
//================================================================================
/*
Class : Taper
*/
//================================================================================
+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
*/
//================================================================================
-double Length2D::GetValue( long theElementId )
+Length2D::Length2D( SMDSAbs_ElementType type ):
+ myType ( type )
{
- TSequenceOfXYZ P;
+}
- if ( GetPoints( theElementId, P ))
- {
- double aVal = 0;
- int len = P.size();
- SMDSAbs_EntityType aType = P.getElementEntity();
+bool Length2D::IsApplicable( const SMDS_MeshElement* element ) const
+{
+ return ( NumericalFunctor::IsApplicable( element ) &&
+ element->GetEntityType() != SMDSEntity_Polyhedra );
+}
- switch (aType) {
- case SMDSEntity_Edge:
- if (len == 2)
- aVal = getDistance( P( 1 ), P( 2 ) );
- break;
- case SMDSEntity_Quad_Edge:
- if (len == 3) // quadratic edge
- aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
- break;
- case SMDSEntity_Triangle:
- if (len == 3){ // triangles
- double L1 = getDistance(P( 1 ),P( 2 ));
- double L2 = getDistance(P( 2 ),P( 3 ));
- double L3 = getDistance(P( 3 ),P( 1 ));
- aVal = Min(L1,Min(L2,L3));
- }
- break;
- case SMDSEntity_Quadrangle:
- if (len == 4){ // quadrangles
- double L1 = getDistance(P( 1 ),P( 2 ));
- double L2 = getDistance(P( 2 ),P( 3 ));
- double L3 = getDistance(P( 3 ),P( 4 ));
- double L4 = getDistance(P( 4 ),P( 1 ));
- aVal = Min(Min(L1,L2),Min(L3,L4));
- }
- break;
- case SMDSEntity_Quad_Triangle:
- case SMDSEntity_BiQuad_Triangle:
- if (len >= 6){ // quadratic triangles
- double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
- double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
- double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
- aVal = Min(L1,Min(L2,L3));
- }
- break;
- case SMDSEntity_Quad_Quadrangle:
- case SMDSEntity_BiQuad_Quadrangle:
- if (len >= 8){ // quadratic quadrangles
- double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
- double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
- double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
- double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
- aVal = Min(Min(L1,L2),Min(L3,L4));
- }
- break;
- case SMDSEntity_Tetra:
- if (len == 4){ // tetrahedra
- double L1 = getDistance(P( 1 ),P( 2 ));
- double L2 = getDistance(P( 2 ),P( 3 ));
- double L3 = getDistance(P( 3 ),P( 1 ));
- double L4 = getDistance(P( 1 ),P( 4 ));
- double L5 = getDistance(P( 2 ),P( 4 ));
- double L6 = getDistance(P( 3 ),P( 4 ));
- aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
- }
- break;
- case SMDSEntity_Pyramid:
- if (len == 5){ // piramids
- double L1 = getDistance(P( 1 ),P( 2 ));
- double L2 = getDistance(P( 2 ),P( 3 ));
- double L3 = getDistance(P( 3 ),P( 4 ));
- double L4 = getDistance(P( 4 ),P( 1 ));
- double L5 = getDistance(P( 1 ),P( 5 ));
- double L6 = getDistance(P( 2 ),P( 5 ));
- double L7 = getDistance(P( 3 ),P( 5 ));
- double L8 = getDistance(P( 4 ),P( 5 ));
-
- aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
- aVal = Min(aVal,Min(L7,L8));
- }
- break;
- case SMDSEntity_Penta:
- if (len == 6) { // pentaidres
- double L1 = getDistance(P( 1 ),P( 2 ));
- double L2 = getDistance(P( 2 ),P( 3 ));
- double L3 = getDistance(P( 3 ),P( 1 ));
- double L4 = getDistance(P( 4 ),P( 5 ));
- double L5 = getDistance(P( 5 ),P( 6 ));
- double L6 = getDistance(P( 6 ),P( 4 ));
- double L7 = getDistance(P( 1 ),P( 4 ));
- double L8 = getDistance(P( 2 ),P( 5 ));
- double L9 = getDistance(P( 3 ),P( 6 ));
-
- aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
- aVal = Min(aVal,Min(Min(L7,L8),L9));
- }
- break;
- case SMDSEntity_Hexa:
- if (len == 8){ // hexahedron
- double L1 = getDistance(P( 1 ),P( 2 ));
- double L2 = getDistance(P( 2 ),P( 3 ));
- double L3 = getDistance(P( 3 ),P( 4 ));
- double L4 = getDistance(P( 4 ),P( 1 ));
- double L5 = getDistance(P( 5 ),P( 6 ));
- double L6 = getDistance(P( 6 ),P( 7 ));
- double L7 = getDistance(P( 7 ),P( 8 ));
- double L8 = getDistance(P( 8 ),P( 5 ));
- double L9 = getDistance(P( 1 ),P( 5 ));
- double L10= getDistance(P( 2 ),P( 6 ));
- double L11= getDistance(P( 3 ),P( 7 ));
- double L12= getDistance(P( 4 ),P( 8 ));
-
- aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
- aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10)));
- aVal = Min(aVal,Min(L11,L12));
- }
- break;
- case SMDSEntity_Quad_Tetra:
- if (len == 10){ // quadratic tetraidrs
- double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
- double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
- double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
- double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
- double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
- double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
- aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
- }
- break;
- case SMDSEntity_Quad_Pyramid:
- if (len == 13){ // quadratic piramids
- double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
- double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
- double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
- double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
- double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
- double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
- double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
- double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
- aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
- aVal = Min(aVal,Min(L7,L8));
- }
- break;
- case SMDSEntity_Quad_Penta:
- if (len == 15){ // quadratic pentaidres
- double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
- double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
- double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
- double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
- double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
- double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
- double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
- double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
- double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
- aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
- aVal = Min(aVal,Min(Min(L7,L8),L9));
- }
- break;
- case SMDSEntity_Quad_Hexa:
- case SMDSEntity_TriQuad_Hexa:
- if (len >= 20) { // quadratic hexaider
- double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
- double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
- double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
- double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
- double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
- double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
- double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
- double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
- double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
- double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
- double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
- double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
- aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
- aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10)));
- aVal = Min(aVal,Min(L11,L12));
- }
- break;
- case SMDSEntity_Polygon:
- if ( len > 1 ) {
- aVal = getDistance( P(1), P( P.size() ));
- for ( size_t i = 1; i < P.size(); ++i )
- aVal = Min( aVal, getDistance( P( i ), P( i+1 )));
- }
- break;
- case SMDSEntity_Quad_Polygon:
- if ( len > 2 ) {
- aVal = getDistance( P(1), P( P.size() )) + getDistance( P(P.size()), P( P.size()-1 ));
- for ( size_t i = 1; i < P.size()-1; i += 2 )
- aVal = Min( aVal, getDistance( P( i ), P( i+1 )) + getDistance( P( i+1 ), P( i+2 )));
- }
- break;
- case SMDSEntity_Hexagonal_Prism:
- if (len == 12) { // hexagonal prism
- double L1 = getDistance(P( 1 ),P( 2 ));
- double L2 = getDistance(P( 2 ),P( 3 ));
- double L3 = getDistance(P( 3 ),P( 4 ));
- double L4 = getDistance(P( 4 ),P( 5 ));
- double L5 = getDistance(P( 5 ),P( 6 ));
- double L6 = getDistance(P( 6 ),P( 1 ));
-
- double L7 = getDistance(P( 7 ), P( 8 ));
- double L8 = getDistance(P( 8 ), P( 9 ));
- double L9 = getDistance(P( 9 ), P( 10 ));
- double L10= getDistance(P( 10 ),P( 11 ));
- double L11= getDistance(P( 11 ),P( 12 ));
- double L12= getDistance(P( 12 ),P( 7 ));
-
- double L13 = getDistance(P( 1 ),P( 7 ));
- double L14 = getDistance(P( 2 ),P( 8 ));
- double L15 = getDistance(P( 3 ),P( 9 ));
- double L16 = getDistance(P( 4 ),P( 10 ));
- double L17 = getDistance(P( 5 ),P( 11 ));
- double L18 = getDistance(P( 6 ),P( 12 ));
- aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
- aVal = Min(aVal, Min(Min(Min(L7,L8),Min(L9,L10)),Min(L11,L12)));
- aVal = Min(aVal, Min(Min(Min(L13,L14),Min(L15,L16)),Min(L17,L18)));
- }
- break;
- case SMDSEntity_Polyhedra:
- {
+double Length2D::GetValue( const TSequenceOfXYZ& P )
+{
+ double aVal = 0;
+ int len = P.size();
+ SMDSAbs_EntityType aType = P.getElementEntity();
+
+ switch (aType) {
+ case SMDSEntity_Edge:
+ if (len == 2)
+ aVal = getDistance( P( 1 ), P( 2 ) );
+ break;
+ case SMDSEntity_Quad_Edge:
+ if (len == 3) // quadratic edge
+ aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
+ break;
+ case SMDSEntity_Triangle:
+ if (len == 3){ // triangles
+ double L1 = getDistance(P( 1 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 1 ));
+ aVal = Min(L1,Min(L2,L3));
+ }
+ break;
+ case SMDSEntity_Quadrangle:
+ if (len == 4){ // quadrangles
+ double L1 = getDistance(P( 1 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 4 ));
+ double L4 = getDistance(P( 4 ),P( 1 ));
+ aVal = Min(Min(L1,L2),Min(L3,L4));
+ }
+ break;
+ case SMDSEntity_Quad_Triangle:
+ case SMDSEntity_BiQuad_Triangle:
+ if (len >= 6){ // quadratic triangles
+ double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
+ double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
+ double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 1 ));
+ aVal = Min(L1,Min(L2,L3));
+ }
+ break;
+ case SMDSEntity_Quad_Quadrangle:
+ case SMDSEntity_BiQuad_Quadrangle:
+ if (len >= 8){ // quadratic quadrangles
+ double L1 = getDistance(P( 1 ),P( 2 )) + getDistance(P( 2 ),P( 3 ));
+ double L2 = getDistance(P( 3 ),P( 4 )) + getDistance(P( 4 ),P( 5 ));
+ double L3 = getDistance(P( 5 ),P( 6 )) + getDistance(P( 6 ),P( 7 ));
+ double L4 = getDistance(P( 7 ),P( 8 )) + getDistance(P( 8 ),P( 1 ));
+ aVal = Min(Min(L1,L2),Min(L3,L4));
+ }
+ break;
+ case SMDSEntity_Tetra:
+ if (len == 4){ // tetrahedra
+ double L1 = getDistance(P( 1 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 1 ));
+ double L4 = getDistance(P( 1 ),P( 4 ));
+ double L5 = getDistance(P( 2 ),P( 4 ));
+ double L6 = getDistance(P( 3 ),P( 4 ));
+ aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+ }
+ break;
+ case SMDSEntity_Pyramid:
+ if (len == 5){ // pyramid
+ double L1 = getDistance(P( 1 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 4 ));
+ double L4 = getDistance(P( 4 ),P( 1 ));
+ double L5 = getDistance(P( 1 ),P( 5 ));
+ double L6 = getDistance(P( 2 ),P( 5 ));
+ double L7 = getDistance(P( 3 ),P( 5 ));
+ double L8 = getDistance(P( 4 ),P( 5 ));
+
+ aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+ aVal = Min(aVal,Min(L7,L8));
+ }
+ break;
+ case SMDSEntity_Penta:
+ if (len == 6) { // pentahedron
+ double L1 = getDistance(P( 1 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 1 ));
+ double L4 = getDistance(P( 4 ),P( 5 ));
+ double L5 = getDistance(P( 5 ),P( 6 ));
+ double L6 = getDistance(P( 6 ),P( 4 ));
+ double L7 = getDistance(P( 1 ),P( 4 ));
+ double L8 = getDistance(P( 2 ),P( 5 ));
+ double L9 = getDistance(P( 3 ),P( 6 ));
+
+ aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+ aVal = Min(aVal,Min(Min(L7,L8),L9));
+ }
+ break;
+ case SMDSEntity_Hexa:
+ if (len == 8){ // hexahedron
+ double L1 = getDistance(P( 1 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 4 ));
+ double L4 = getDistance(P( 4 ),P( 1 ));
+ double L5 = getDistance(P( 5 ),P( 6 ));
+ double L6 = getDistance(P( 6 ),P( 7 ));
+ double L7 = getDistance(P( 7 ),P( 8 ));
+ double L8 = getDistance(P( 8 ),P( 5 ));
+ double L9 = getDistance(P( 1 ),P( 5 ));
+ double L10= getDistance(P( 2 ),P( 6 ));
+ double L11= getDistance(P( 3 ),P( 7 ));
+ double L12= getDistance(P( 4 ),P( 8 ));
+
+ aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+ aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10)));
+ aVal = Min(aVal,Min(L11,L12));
+ }
+ break;
+ case SMDSEntity_Quad_Tetra:
+ if (len == 10){ // quadratic tetrahedron
+ double L1 = getDistance(P( 1 ),P( 5 )) + getDistance(P( 5 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 6 )) + getDistance(P( 6 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 7 )) + getDistance(P( 7 ),P( 1 ));
+ double L4 = getDistance(P( 1 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
+ double L5 = getDistance(P( 2 ),P( 9 )) + getDistance(P( 9 ),P( 4 ));
+ double L6 = getDistance(P( 3 ),P( 10 )) + getDistance(P( 10 ),P( 4 ));
+ aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+ }
+ break;
+ case SMDSEntity_Quad_Pyramid:
+ if (len == 13){ // quadratic pyramid
+ double L1 = getDistance(P( 1 ),P( 6 )) + getDistance(P( 6 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 7 )) + getDistance(P( 7 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 8 )) + getDistance(P( 8 ),P( 4 ));
+ double L4 = getDistance(P( 4 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
+ double L5 = getDistance(P( 1 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
+ double L6 = getDistance(P( 2 ),P( 11 )) + getDistance(P( 11 ),P( 5 ));
+ double L7 = getDistance(P( 3 ),P( 12 )) + getDistance(P( 12 ),P( 5 ));
+ double L8 = getDistance(P( 4 ),P( 13 )) + getDistance(P( 13 ),P( 5 ));
+ aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+ aVal = Min(aVal,Min(L7,L8));
+ }
+ break;
+ case SMDSEntity_Quad_Penta:
+ case SMDSEntity_BiQuad_Penta:
+ if (len >= 15){ // quadratic pentahedron
+ double L1 = getDistance(P( 1 ),P( 7 )) + getDistance(P( 7 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 8 )) + getDistance(P( 8 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 9 )) + getDistance(P( 9 ),P( 1 ));
+ double L4 = getDistance(P( 4 ),P( 10 )) + getDistance(P( 10 ),P( 5 ));
+ double L5 = getDistance(P( 5 ),P( 11 )) + getDistance(P( 11 ),P( 6 ));
+ double L6 = getDistance(P( 6 ),P( 12 )) + getDistance(P( 12 ),P( 4 ));
+ double L7 = getDistance(P( 1 ),P( 13 )) + getDistance(P( 13 ),P( 4 ));
+ double L8 = getDistance(P( 2 ),P( 14 )) + getDistance(P( 14 ),P( 5 ));
+ double L9 = getDistance(P( 3 ),P( 15 )) + getDistance(P( 15 ),P( 6 ));
+ aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+ aVal = Min(aVal,Min(Min(L7,L8),L9));
+ }
+ break;
+ case SMDSEntity_Quad_Hexa:
+ case SMDSEntity_TriQuad_Hexa:
+ if (len >= 20) { // quadratic hexahedron
+ double L1 = getDistance(P( 1 ),P( 9 )) + getDistance(P( 9 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 10 )) + getDistance(P( 10 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 11 )) + getDistance(P( 11 ),P( 4 ));
+ double L4 = getDistance(P( 4 ),P( 12 )) + getDistance(P( 12 ),P( 1 ));
+ double L5 = getDistance(P( 5 ),P( 13 )) + getDistance(P( 13 ),P( 6 ));
+ double L6 = getDistance(P( 6 ),P( 14 )) + getDistance(P( 14 ),P( 7 ));
+ double L7 = getDistance(P( 7 ),P( 15 )) + getDistance(P( 15 ),P( 8 ));
+ double L8 = getDistance(P( 8 ),P( 16 )) + getDistance(P( 16 ),P( 5 ));
+ double L9 = getDistance(P( 1 ),P( 17 )) + getDistance(P( 17 ),P( 5 ));
+ double L10= getDistance(P( 2 ),P( 18 )) + getDistance(P( 18 ),P( 6 ));
+ double L11= getDistance(P( 3 ),P( 19 )) + getDistance(P( 19 ),P( 7 ));
+ double L12= getDistance(P( 4 ),P( 20 )) + getDistance(P( 20 ),P( 8 ));
+ aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+ aVal = Min(aVal,Min(Min(L7,L8),Min(L9,L10)));
+ aVal = Min(aVal,Min(L11,L12));
}
break;
- default:
- return 0;
+ case SMDSEntity_Polygon:
+ if ( len > 1 ) {
+ aVal = getDistance( P(1), P( P.size() ));
+ for ( size_t i = 1; i < P.size(); ++i )
+ aVal = Min( aVal, getDistance( P( i ), P( i+1 )));
}
-
- if (aVal < 0 ) {
- return 0.;
+ break;
+ case SMDSEntity_Quad_Polygon:
+ if ( len > 2 ) {
+ aVal = getDistance( P(1), P( P.size() )) + getDistance( P(P.size()), P( P.size()-1 ));
+ for ( size_t i = 1; i < P.size()-1; i += 2 )
+ aVal = Min( aVal, getDistance( P( i ), P( i+1 )) + getDistance( P( i+1 ), P( i+2 )));
}
-
- if ( myPrecision >= 0 )
- {
- double prec = pow( 10., (double)( myPrecision ) );
- aVal = floor( aVal * prec + 0.5 ) / prec;
+ break;
+ case SMDSEntity_Hexagonal_Prism:
+ if (len == 12) { // hexagonal prism
+ double L1 = getDistance(P( 1 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 4 ));
+ double L4 = getDistance(P( 4 ),P( 5 ));
+ double L5 = getDistance(P( 5 ),P( 6 ));
+ double L6 = getDistance(P( 6 ),P( 1 ));
+
+ double L7 = getDistance(P( 7 ), P( 8 ));
+ double L8 = getDistance(P( 8 ), P( 9 ));
+ double L9 = getDistance(P( 9 ), P( 10 ));
+ double L10= getDistance(P( 10 ),P( 11 ));
+ double L11= getDistance(P( 11 ),P( 12 ));
+ double L12= getDistance(P( 12 ),P( 7 ));
+
+ double L13 = getDistance(P( 1 ),P( 7 ));
+ double L14 = getDistance(P( 2 ),P( 8 ));
+ double L15 = getDistance(P( 3 ),P( 9 ));
+ double L16 = getDistance(P( 4 ),P( 10 ));
+ double L17 = getDistance(P( 5 ),P( 11 ));
+ double L18 = getDistance(P( 6 ),P( 12 ));
+ aVal = Min(Min(Min(L1,L2),Min(L3,L4)),Min(L5,L6));
+ aVal = Min(aVal, Min(Min(Min(L7,L8),Min(L9,L10)),Min(L11,L12)));
+ aVal = Min(aVal, Min(Min(Min(L13,L14),Min(L15,L16)),Min(L17,L18)));
}
+ break;
+ case SMDSEntity_Polyhedra:
+ {
+ }
+ break;
+ default:
+ return 0;
+ }
- return aVal;
+ if (aVal < 0 ) {
+ return 0.;
+ }
+ if ( myPrecision >= 0 )
+ {
+ double prec = pow( 10., (double)( myPrecision ) );
+ aVal = floor( aVal * prec + 0.5 ) / prec;
}
- return 0.;
+
+ return aVal;
}
double Length2D::GetBadRate( double Value, int /*nbNodes*/ ) const
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;
- SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
- for(; anIter->more(); ){
- const SMDS_MeshFace* anElem = anIter->next();
+ if ( myType == SMDSAbs_Face )
+ {
+ for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); )
+ {
+ const SMDS_MeshFace* anElem = anIter->next();
+ if ( anElem->IsQuadratic() )
+ {
+ // use special nodes iterator
+ SMDS_NodeIteratorPtr anIter = anElem->interlacedNodesIterator();
+ smIdType aNodeId[4] = { 0,0,0,0 };
+ gp_Pnt P[4];
- if(anElem->IsQuadratic()) {
- const SMDS_VtkFace* F =
- dynamic_cast<const SMDS_VtkFace*>(anElem);
- // use special nodes iterator
- SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
- long aNodeId[4] = { 0,0,0,0 };
- gp_Pnt P[4];
-
- double aLength = 0;
- const SMDS_MeshElement* aNode;
- if(anIter->more()){
- aNode = anIter->next();
- const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
- P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
- aNodeId[0] = aNodeId[1] = aNode->GetID();
- aLength = 0;
- }
- for(; anIter->more(); ){
- const SMDS_MeshNode* N1 = static_cast<const SMDS_MeshNode*> (anIter->next());
- P[2] = gp_Pnt(N1->X(),N1->Y(),N1->Z());
- aNodeId[2] = N1->GetID();
- aLength = P[1].Distance(P[2]);
- if(!anIter->more()) break;
- const SMDS_MeshNode* N2 = static_cast<const SMDS_MeshNode*> (anIter->next());
- P[3] = gp_Pnt(N2->X(),N2->Y(),N2->Z());
- aNodeId[3] = N2->GetID();
- aLength += P[2].Distance(P[3]);
+ 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_ElemIteratorPtr aNodesIter = anElem->nodesIterator();
- long aNodeId[2] = {0,0};
- gp_Pnt P[3];
-
- double aLength;
- const SMDS_MeshElement* aNode;
- if(aNodesIter->more()){
- aNode = aNodesIter->next();
- const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
- P[0] = P[1] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
- aNodeId[0] = aNodeId[1] = aNode->GetID();
- aLength = 0;
- }
- for(; aNodesIter->more(); ){
- aNode = aNodesIter->next();
- const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode;
- long anId = aNode->GetID();
-
- P[2] = gp_Pnt(aNodes->X(),aNodes->Y(),aNodes->Z());
-
- aLength = P[1].Distance(P[2]);
-
- Value aValue(aLength,aNodeId[1],anId);
- aNodeId[1] = anId;
- P[1] = P[2];
+ else {
+ SMDS_NodeIteratorPtr aNodesIter = anElem->nodeIterator();
+ smIdType 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();
+ smIdType anId = aNode->GetID();
+
+ 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[0].Distance(P[1]);
+
+ Value aValue(aLength,aNodeId[0],aNodeId[1]);
theValues.insert(aValue);
}
+ }
+ }
+ else
+ {
+ // not implemented
+ }
+}
+
+//================================================================================
+/*
+ Class : Deflection2D
+ Description : computes distance between a face center and an underlying surface
+*/
+//================================================================================
+
+double Deflection2D::GetValue( const TSequenceOfXYZ& P )
+{
+ if ( myMesh && P.getElement() )
+ {
+ // get underlying surface
+ if ( myShapeIndex != P.getElement()->getshapeId() )
+ {
+ mySurface.Nullify();
+ myShapeIndex = P.getElement()->getshapeId();
+ const TopoDS_Shape& S =
+ static_cast< const SMESHDS_Mesh* >( myMesh )->IndexToShape( myShapeIndex );
+ if ( !S.IsNull() && S.ShapeType() == TopAbs_FACE )
+ {
+ mySurface = new ShapeAnalysis_Surface( BRep_Tool::Surface( TopoDS::Face( S )));
+
+ GeomLib_IsPlanarSurface isPlaneCheck( mySurface->Surface() );
+ if ( isPlaneCheck.IsPlanar() )
+ myPlane.reset( new gp_Pln( isPlaneCheck.Plan() ));
+ else
+ myPlane.reset();
+ }
+ }
+ // project gravity center to the surface
+ if ( !mySurface.IsNull() )
+ {
+ gp_XYZ gc(0,0,0);
+ gp_XY uv(0,0);
+ int nbUV = 0;
+ for ( size_t i = 0; i < P.size(); ++i )
+ {
+ gc += P(i+1);
- aLength = P[0].Distance(P[1]);
+ if ( SMDS_FacePositionPtr fPos = P.getElement()->GetNode( i )->GetPosition() )
+ {
+ uv.ChangeCoord(1) += fPos->GetUParameter();
+ uv.ChangeCoord(2) += fPos->GetVParameter();
+ ++nbUV;
+ }
+ }
+ gc /= P.size();
+ if ( nbUV ) uv /= nbUV;
- Value aValue(aLength,aNodeId[0],aNodeId[1]);
- theValues.insert(aValue);
+ double maxLen = MaxElementLength2D().GetValue( P );
+ double tol = 1e-3 * maxLen;
+ double dist;
+ if ( myPlane )
+ {
+ dist = myPlane->Distance( gc );
+ if ( dist < tol )
+ dist = 0;
+ }
+ else
+ {
+ if ( uv.X() != 0 && uv.Y() != 0 ) // faster way
+ mySurface->NextValueOfUV( uv, gc, tol, 0.5 * maxLen );
+ else
+ mySurface->ValueOfUV( gc, tol );
+ dist = mySurface->Gap();
+ }
+ return Round( dist );
}
}
+ return 0;
+}
+
+void Deflection2D::SetMesh( const SMDS_Mesh* theMesh )
+{
+ NumericalFunctor::SetMesh( dynamic_cast<const SMESHDS_Mesh* >( theMesh ));
+ myShapeIndex = -100;
+ myPlane.reset();
+}
+
+SMDSAbs_ElementType Deflection2D::GetType() const
+{
+ return SMDSAbs_Face;
+}
+
+double Deflection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+ // meaningless as it is not quality control functor
+ return Value;
}
//================================================================================
*/
//================================================================================
-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;
}
if (!anIter) break;
const SMDS_MeshNode *aNode, *aNode0 = 0;
- TColStd_MapOfInteger aMap, aMapPrev;
+ NCollection_Map< smIdType, smIdHasher > aMap, aMapPrev;
for (i = 0; i <= len; i++) {
aMapPrev = aMap;
while (anElemIter->more()) {
const SMDS_MeshElement* anElem = anElemIter->next();
if (anElem != 0 && anElem->GetType() == SMDSAbs_Face) {
- int anId = anElem->GetID();
+ smIdType anId = anElem->GetID();
aMap.Add(anId);
if (aMapPrev.Contains(anId)) {
void MultiConnection2D::GetValues(MValues& theValues)
{
if ( !myMesh ) return;
- SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
- for(; anIter->more(); ){
- const SMDS_MeshFace* anElem = anIter->next();
- SMDS_ElemIteratorPtr aNodesIter;
- if ( anElem->IsQuadratic() )
- aNodesIter = dynamic_cast<const SMDS_VtkFace*>
- (anElem)->interlacedNodesElemIterator();
- else
- aNodesIter = anElem->nodesIterator();
- long aNodeId[3] = {0,0,0};
+ for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); )
+ {
+ const SMDS_MeshFace* anElem = anIter->next();
+ SMDS_NodeIteratorPtr aNodesIter = anElem->interlacedNodesIterator();
- //int aNbConnects=0;
- const SMDS_MeshNode* aNode0;
- const SMDS_MeshNode* aNode1;
+ const SMDS_MeshNode* aNode1 = anElem->GetNode( anElem->NbNodes() - 1 );
const SMDS_MeshNode* aNode2;
- if(aNodesIter->more()){
- aNode0 = (SMDS_MeshNode*) aNodesIter->next();
- aNode1 = aNode0;
- const SMDS_MeshNode* aNodes = (SMDS_MeshNode*) aNode1;
- aNodeId[0] = aNodeId[1] = aNodes->GetID();
- }
- for(; aNodesIter->more(); ) {
- aNode2 = (SMDS_MeshNode*) aNodesIter->next();
- long anId = aNode2->GetID();
- aNodeId[2] = anId;
-
- Value aValue(aNodeId[1],aNodeId[2]);
- MValues::iterator aItr = theValues.find(aValue);
- if (aItr != theValues.end()){
- aItr->second += 1;
- //aNbConnects = nb;
- }
- else {
- theValues[aValue] = 1;
- //aNbConnects = 1;
- }
- //cout << "NodeIds: "<<aNodeId[1]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
- aNodeId[1] = aNodeId[2];
+ for ( ; aNodesIter->more(); )
+ {
+ aNode2 = aNodesIter->next();
+
+ Value aValue ( aNode1->GetID(), aNode2->GetID() );
+ MValues::iterator aItr = theValues.insert( std::make_pair( aValue, 0 )).first;
+ aItr->second++;
aNode1 = aNode2;
}
- Value aValue(aNodeId[0],aNodeId[2]);
- MValues::iterator aItr = theValues.find(aValue);
- if (aItr != theValues.end()) {
- aItr->second += 1;
- //aNbConnects = nb;
- }
- else {
- theValues[aValue] = 1;
- //aNbConnects = 1;
- }
- //cout << "NodeIds: "<<aNodeId[0]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
}
-
+ return;
}
//================================================================================
double diameter = 0;
if ( const SMDS_BallElement* ball =
- dynamic_cast<const SMDS_BallElement*>( myMesh->FindElement( theId )))
+ myMesh->DownCast< SMDS_BallElement >( myMesh->FindElement( theId )))
{
diameter = ball->GetDiameter();
}
return SMDSAbs_Node;
}
+//================================================================================
+/*
+ Class : ScaledJacobian
+ Description : Functor returning the ScaledJacobian for volumetric elements
+*/
+//================================================================================
+
+double ScaledJacobian::GetValue( long theElementId )
+{
+ if ( theElementId && myMesh ) {
+ SMDS_VolumeTool aVolumeTool;
+ if ( aVolumeTool.Set( myMesh->FindElement( theElementId )))
+ return aVolumeTool.GetScaledJacobian();
+ }
+ return 0;
+
+ /*
+ //VTK version not used because lack of implementation for HEXAGONAL_PRISM.
+ //Several mesh quality measures implemented in vtkMeshQuality can be accessed left here as reference
+ double aVal = 0;
+ myCurrElement = myMesh->FindElement( theElementId );
+ if ( myCurrElement )
+ {
+ VTKCellType cellType = myCurrElement->GetVtkType();
+ vtkUnstructuredGrid* grid = const_cast<SMDS_Mesh*>( myMesh )->GetGrid();
+ vtkCell* avtkCell = grid->GetCell( myCurrElement->GetVtkID() );
+ switch ( cellType )
+ {
+ case VTK_QUADRATIC_TETRA:
+ case VTK_TETRA:
+ aVal = Round( vtkMeshQuality::TetScaledJacobian( avtkCell ));
+ break;
+ case VTK_QUADRATIC_HEXAHEDRON:
+ case VTK_HEXAHEDRON:
+ aVal = Round( vtkMeshQuality::HexScaledJacobian( avtkCell ));
+ break;
+ case VTK_QUADRATIC_WEDGE:
+ case VTK_WEDGE: //Pentahedron
+ aVal = Round( vtkMeshQuality::WedgeScaledJacobian( avtkCell ));
+ break;
+ case VTK_QUADRATIC_PYRAMID:
+ case VTK_PYRAMID:
+ aVal = Round( vtkMeshQuality::PyramidScaledJacobian( avtkCell ));
+ break;
+ case VTK_HEXAGONAL_PRISM:
+ case VTK_POLYHEDRON:
+ default:
+ break;
+ }
+ }
+ return aVal;
+ */
+}
+
+double ScaledJacobian::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+ return Value;
+}
+
+SMDSAbs_ElementType ScaledJacobian::GetType() const
+{
+ return SMDSAbs_Volume;
+}
+
/*
PREDICATES
*/
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
bool OverConstrainedVolume::IsSatisfy(long theElementId )
{
- // An element is over-constrained if it has N-1 free borders where
- // N is the number of edges/faces for a 2D/3D element.
- SMDS_VolumeTool myTool;
- if ( myTool.Set( myMesh->FindElement(theElementId)))
+ // An element is over-constrained if all its nodes are on the boundary.
+ // A node is on the boundary if it is connected to one or more faces.
+ SMDS_VolumeTool myTool;
+ if (myTool.Set(myMesh->FindElement(theElementId)))
{
- int nbSharedFaces = 0;
- for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
- if ( !myTool.IsFreeFace( iF ) && ++nbSharedFaces > 1 )
- break;
- return ( nbSharedFaces == 1 );
+ auto nodes = myTool.GetNodes();
+
+ for (int i = 0; i < myTool.NbNodes(); ++i)
+ {
+ auto node = nodes[i];
+ if (node->NbInverseElements(SMDSAbs_Face) == 0)
+ {
+ return false;
+ }
+ }
+ return true;
}
return false;
}
bool OverConstrainedFace::IsSatisfy(long theElementId )
{
- // An element is over-constrained if it has N-1 free borders where
- // N is the number of edges/faces for a 2D/3D element.
- if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
- if ( face->GetType() == SMDSAbs_Face )
+ // An element is over-constrained if all its nodes are on the boundary.
+ // A node is on the boundary if it is connected to one or more faces.
+ if (const SMDS_MeshElement *face = myMesh->FindElement(theElementId))
+ if (face->GetType() == SMDSAbs_Face)
{
- int nbSharedBorders = 0;
int nbN = face->NbCornerNodes();
- for ( int i = 0; i < nbN; ++i )
+ for (int i = 0; i < nbN; ++i)
{
- // check if a link is shared by another face
- const SMDS_MeshNode* n1 = face->GetNode( i );
- const SMDS_MeshNode* n2 = face->GetNode( (i+1)%nbN );
- SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator( SMDSAbs_Face );
- bool isShared = false;
- while ( !isShared && fIt->more() )
- {
- const SMDS_MeshElement* f = fIt->next();
- isShared = ( f != face && f->GetNodeIndex(n2) != -1 );
- }
- if ( isShared && ++nbSharedBorders > 1 )
- break;
+ const SMDS_MeshNode *n1 = face->GetNode(i);
+ if (n1->NbInverseElements(SMDSAbs_Edge) == 0)
+ return false;
}
- return ( nbSharedBorders == 1 );
+ return true;
}
return false;
}
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 );
if ( myMeshModifTracer.IsMeshModified() )
{
TIDSortedNodeSet nodesToCheck;
- SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator(/*idInceasingOrder=*/true);
+ SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator();
while ( nIt->more() )
nodesToCheck.insert( nodesToCheck.end(), nIt->next() );
myMesh = theMesh;
}
-bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId )
+bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const smIdType 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 smIdType anId = anElem->GetID();
+ if ( anId != theFaceId && anElem->GetNodeIndex( theNodes[1] ) >= 0 )
+ return false;
}
}
return true;
void FreeEdges::GetBoreders(TBorders& theBorders)
{
TBorders aRegistry;
- SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
- for(; anIter->more(); ){
+ for ( SMDS_FaceIteratorPtr anIter = myMesh->facesIterator(); anIter->more(); )
+ {
const SMDS_MeshFace* anElem = anIter->next();
long anElemId = anElem->GetID();
- SMDS_ElemIteratorPtr aNodesIter;
- if ( anElem->IsQuadratic() )
- aNodesIter = static_cast<const SMDS_VtkFace*>(anElem)->
- interlacedNodesElemIterator();
- else
- aNodesIter = anElem->nodesIterator();
+ SMDS_NodeIteratorPtr aNodesIter = anElem->interlacedNodesIterator();
+ if ( !aNodesIter->more() ) continue;
long aNodeId[2] = {0,0};
- const SMDS_MeshElement* aNode;
- if(aNodesIter->more()){
- aNode = aNodesIter->next();
- aNodeId[0] = aNodeId[1] = aNode->GetID();
- }
- for(; aNodesIter->more(); ){
- aNode = aNodesIter->next();
- long anId = aNode->GetID();
- Border aBorder(anElemId,aNodeId[1],anId);
- aNodeId[1] = anId;
- UpdateBorders(aBorder,aRegistry,theBorders);
+ aNodeId[0] = anElem->GetNode( anElem->NbNodes()-1 )->GetID();
+ for ( ; aNodesIter->more(); )
+ {
+ aNodeId[1] = aNodesIter->next()->GetID();
+ Border aBorder( anElemId, aNodeId[0], aNodeId[1] );
+ UpdateBorders( aBorder, aRegistry, theBorders );
+ aNodeId[0] = aNodeId[1];
}
- Border aBorder(anElemId,aNodeId[0],aNodeId[1]);
- UpdateBorders(aBorder,aRegistry,theBorders);
}
}
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);
}
SMDSAbs_ElementType aGrpElType = (SMDSAbs_ElementType)aGrp->GetType();
if ( myType == aGrpElType || (myType == SMDSAbs_All && aGrpElType != SMDSAbs_Node) ) {
// add elements IDS into control
- int aSize = aGrp->Extent();
- for (int i = 0; i < aSize; i++)
+ smIdType aSize = aGrp->Extent();
+ for (smIdType i = 0; i < aSize; i++)
myIDs.insert( aGrp->GetID(i+1) );
}
}
SMDSAbs_ElementType ConnectedElements::GetType() const
{ return myType; }
-int ConnectedElements::GetNode() const
+smIdType ConnectedElements::GetNode() const
{ return myXYZ.empty() ? myNodeID : 0; } // myNodeID can be found by myXYZ
std::vector<double> ConnectedElements::GetPoint() const
}
}
-void ConnectedElements::SetNode( int nodeID )
+void ConnectedElements::SetNode( smIdType nodeID )
{
myNodeID = nodeID;
myXYZ.clear();
return false;
std::list< const SMDS_MeshNode* > nodeQueue( 1, node0 );
- std::set< int > checkedNodeIDs;
+ std::set< smIdType > checkedNodeIDs;
// algo:
// foreach node in nodeQueue:
// foreach element sharing a node:
{
// keep elements of myType
const SMDS_MeshElement* element = eIt->next();
- if ( element->GetType() == myType )
+ if ( myType == SMDSAbs_All || element->GetType() == myType )
myOkIDs.insert( myOkIDs.end(), element->GetID() );
// enqueue nodes of the element
while ( nIt->more() )
{
const SMDS_MeshNode* n = static_cast< const SMDS_MeshNode* >( nIt->next() );
- if ( checkedNodeIDs.insert( n->GetID() ).second )
+ if ( checkedNodeIDs.insert( n->GetID()).second )
nodeQueue.push_back( n );
}
}
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 ));
}
}
{
theResStr.Clear();
- TColStd_SequenceOfInteger anIntSeq;
- TColStd_SequenceOfAsciiString aStrSeq;
+ TIDsSeq anIntSeq;
+ NCollection_Sequence< std::string > aStrSeq;
- TColStd_MapIteratorOfMapOfInteger anIter( myIds );
+ TIDsMap::Iterator anIter( myIds );
for ( ; anIter.More(); anIter.Next() )
{
- int anId = anIter.Key();
- TCollection_AsciiString aStr( anId );
+ smIdType anId = anIter.Key();
+ SMESH_Comment aStr( anId );
anIntSeq.Append( anId );
aStrSeq.Append( aStr );
}
- for ( int i = 1, n = myMin.Length(); i <= n; i++ )
+ for ( smIdType i = 1, n = myMin.size(); i <= n; i++ )
{
- int aMinId = myMin( i );
- int aMaxId = myMax( i );
+ smIdType aMinId = myMin[i];
+ smIdType aMaxId = myMax[i];
- TCollection_AsciiString aStr;
+ SMESH_Comment aStr;
if ( aMinId != IntegerFirst() )
- aStr += aMinId;
+ aStr << aMinId;
- aStr += "-";
+ aStr << "-";
- if ( aMaxId != IntegerLast() )
- aStr += aMaxId;
+ if ( aMaxId != std::numeric_limits<smIdType>::max() )
+ aStr << aMaxId;
// find position of the string in result sequence and insert string in it
if ( anIntSeq.Length() == 0 )
{
anIntSeq.Append( aMinId );
- aStrSeq.Append( aStr );
+ aStrSeq.Append( (const char*)aStr );
}
else
{
if ( aMinId < anIntSeq.First() )
{
anIntSeq.Prepend( aMinId );
- aStrSeq.Prepend( aStr );
+ aStrSeq.Prepend( (const char*)aStr );
}
else if ( aMinId > anIntSeq.Last() )
{
anIntSeq.Append( aMinId );
- aStrSeq.Append( aStr );
+ aStrSeq.Append( (const char*)aStr );
}
else
for ( int j = 1, k = anIntSeq.Length(); j <= k; j++ )
if ( aMinId < anIntSeq( j ) )
{
anIntSeq.InsertBefore( j, aMinId );
- aStrSeq.InsertBefore( j, aStr );
+ aStrSeq.InsertBefore( j, (const char*)aStr );
break;
}
}
if ( aStrSeq.Length() == 0 )
return;
-
- theResStr = aStrSeq( 1 );
+ std::string aResStr;
+ aResStr = aStrSeq( 1 );
for ( int j = 2, k = aStrSeq.Length(); j <= k; j++ )
{
- theResStr += ",";
- theResStr += aStrSeq( j );
+ aResStr += ",";
+ aResStr += aStrSeq( j );
}
+ theResStr = aResStr.c_str();
}
//=======================================================================
//=======================================================================
bool RangeOfIds::SetRangeStr( const TCollection_AsciiString& theStr )
{
- myMin.Clear();
- myMax.Clear();
+ myMin.clear();
+ myMax.clear();
myIds.Clear();
TCollection_AsciiString aStr = theStr;
(!aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue()) )
return false;
- myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
- myMax.Append( aMaxStr.IsEmpty() ? IntegerLast() : aMaxStr.IntegerValue() );
+ myMin.push_back( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
+ myMax.push_back( aMaxStr.IsEmpty() ? IntegerLast() : aMaxStr.IntegerValue() );
}
}
if ( myIds.Contains( theId ) )
return true;
- for ( int i = 1, n = myMin.Length(); i <= n; i++ )
- if ( theId >= myMin( i ) && theId <= myMax( i ) )
+ for ( size_t i = 0; i < myMin.size(); i++ )
+ if ( theId >= myMin[i] && theId <= myMax[i] )
return true;
return false;
myPredicate = thePredicate;
}
-void Filter::GetElementsId( const SMDS_Mesh* theMesh,
- PredicatePtr thePredicate,
- TIdSequence& theSequence )
+void Filter::GetElementsId( const SMDS_Mesh* theMesh,
+ PredicatePtr thePredicate,
+ TIdSequence& theSequence,
+ SMDS_ElemIteratorPtr theElements )
{
theSequence.clear();
thePredicate->SetMesh( theMesh );
- SMDS_ElemIteratorPtr elemIt = theMesh->elementsIterator( thePredicate->GetType() );
- if ( elemIt ) {
- while ( elemIt->more() ) {
- const SMDS_MeshElement* anElem = elemIt->next();
- long anId = anElem->GetID();
- if ( thePredicate->IsSatisfy( anId ) )
- theSequence.push_back( anId );
+ if ( !theElements )
+ theElements = theMesh->elementsIterator( thePredicate->GetType() );
+
+ if ( theElements ) {
+ while ( theElements->more() ) {
+ const SMDS_MeshElement* anElem = theElements->next();
+ if ( thePredicate->GetType() == SMDSAbs_All ||
+ thePredicate->GetType() == anElem->GetType() )
+ {
+ long anId = anElem->GetID();
+ if ( thePredicate->IsSatisfy( anId ) )
+ theSequence.push_back( anId );
+ }
}
}
}
void Filter::GetElementsId( const SMDS_Mesh* theMesh,
- Filter::TIdSequence& theSequence )
+ Filter::TIdSequence& theSequence,
+ SMDS_ElemIteratorPtr theElements )
{
- GetElementsId(theMesh,myPredicate,theSequence);
+ GetElementsId(theMesh,myPredicate,theSequence,theElements);
}
/*
// the map of non manifold links and bad geometry
TMapOfLink aMapOfNonManifold;
- TColStd_MapOfInteger aMapOfTreated;
+ TIDsMap aMapOfTreated;
// begin cycle on faces from start index and run on vector till the end
// and from begin to start index to cover whole vector
// as result next time when fi will be equal to aStartIndx
SMDS_MeshFace* aFacePtr = myAllFacePtr[ fi ];
- if ( aMapOfTreated.Contains( aFacePtr->GetID() ) )
+ if ( aMapOfTreated.Contains( aFacePtr->GetID()) )
continue;
aMapOfTreated.Add( aFacePtr->GetID() );
- TColStd_MapOfInteger aResFaces;
+ TIDsMap aResFaces;
if ( !findConnected( myAllFacePtrIntDMap, aFacePtr,
aMapOfNonManifold, aResFaces ) )
continue;
- TColStd_MapIteratorOfMapOfInteger anItr( aResFaces );
+ TIDsMap::Iterator anItr( aResFaces );
for ( ; anItr.More(); anItr.Next() )
{
- int aFaceId = anItr.Key();
+ smIdType aFaceId = anItr.Key();
aMapOfTreated.Add( aFaceId );
myMapIds.Add( aFaceId );
}
( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
SMDS_MeshFace* theStartFace,
ManifoldPart::TMapOfLink& theNonManifold,
- TColStd_MapOfInteger& theResFaces )
+ TIDsMap& theResFaces )
{
theResFaces.Clear();
if ( !theAllFacePtrInt.size() )
SMDS_MeshFace* aNextFace = *pFace;
if ( aPrevFace == aNextFace )
continue;
- int anNextFaceID = aNextFace->GetID();
+ smIdType anNextFaceID = aNextFace->GetID();
if ( myIsOnlyManifold && theResFaces.Contains( anNextFaceID ) )
// should not be with non manifold restriction. probably bad topology
continue;
void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
ManifoldPart::TVectorOfFacePtr& theFaces ) const
{
- std::set<SMDS_MeshCell *> aSetOfFaces;
+
// take all faces that shared first node
- SMDS_ElemIteratorPtr anItr = theLink.myNode1->facesIterator();
- for ( ; anItr->more(); )
- {
- SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
- if ( !aFace )
- continue;
- aSetOfFaces.insert( aFace );
- }
+ SMDS_ElemIteratorPtr anItr = theLink.myNode1->GetInverseElementIterator( SMDSAbs_Face );
+ SMDS_StdIterator< const SMDS_MeshElement*, SMDS_ElemIteratorPtr > faces( anItr ), facesEnd;
+ std::set<const SMDS_MeshElement *> aSetOfFaces( faces, facesEnd );
+
// take all faces that shared second node
- anItr = theLink.myNode2->facesIterator();
+ anItr = theLink.myNode2->GetInverseElementIterator( SMDSAbs_Face );
// find the common part of two sets
for ( ; anItr->more(); )
{
- SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
- if ( aSetOfFaces.count( aFace ) )
- theFaces.push_back( aFace );
+ const SMDS_MeshElement* aFace = anItr->next();
+ if ( aSetOfFaces.count( aFace ))
+ theFaces.push_back( (SMDS_MeshFace*) aFace );
}
}
void ElementsOnSurface::SetTolerance( const double theToler )
{
if ( myToler != theToler )
- myIds.Clear();
- myToler = theToler;
+ {
+ myToler = theToler;
+ process();
+ }
}
double ElementsOnSurface::GetTolerance() const
if ( !myMeshModifTracer.GetMesh() )
return;
- myIds.ReSize( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType ));
+ int nbElems = FromSmIdType<int>( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType ));
+ if ( nbElems > 0 )
+ myIds.ReSize( nbElems );
SMDS_ElemIteratorPtr anIter = myMeshModifTracer.GetMesh()->elementsIterator( myType );
for(; anIter->more(); )
struct ElementsOnShape::Classifier
{
- Classifier() { mySolidClfr = 0; myFlags = 0; }
+ Classifier(): mySolidClfr(0), myProjFace(0), myProjEdge(0), myFlags(0) { myU = myV = 1e100; }
~Classifier();
void Init(const TopoDS_Shape& s, double tol, const Bnd_B3d* box = 0 );
bool IsOut(const gp_Pnt& p) { return SetChecked( true ), (this->*myIsOutFun)( p ); }
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 ); }
void SetFlag ( int flag ) { myFlags |= flag; }
void UnsetFlag( int flag ) { myFlags &= ~flag; }
+ void GetParams( double & u, double & v ) const { u = myU; v = myV; }
private:
bool isOutOfSolid (const gp_Pnt& p);
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
+ BRepClass3d_SolidClassifier* mySolidClfr;
Bnd_B3d myBox;
- GeomAPI_ProjectPointOnSurf myProjFace;
- GeomAPI_ProjectPointOnCurve myProjEdge;
+ GeomAPI_ProjectPointOnSurf* myProjFace;
+ GeomAPI_ProjectPointOnCurve* myProjEdge;
gp_Pnt myVertexXYZ;
TopoDS_Shape myShape;
double myTol;
+ double myU, myV; // result of isOutOfFace() and isOutOfEdge()
int myFlags;
};
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
+ {
+
+ if (SALOME::VerbosityActivated())
+ std::cout << "Avoid ElementsOnShape::clone(), too large: " << size << " bytes " << std::endl;
+
+ return 0;
+ }
+
ElementsOnShape* cln = new ElementsOnShape();
cln->SetAllNodes ( myAllNodesFlag );
cln->SetTolerance( myToler );
void ElementsOnShape::SetTolerance (const double theToler)
{
- if (myToler != theToler) {
+ if (myToler != theToler)
+ {
myToler = theToler;
- SetShape(myShape, myType);
+ TopoDS_Shape s = myShape;
+ myShape.Nullify();
+ SetShape( s, myType );
}
}
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;
}
+//================================================================================
+/*!
+ * \brief Check and optionally return a satisfying shape
+ */
+//================================================================================
+
bool ElementsOnShape::IsSatisfy (const SMDS_MeshNode* node,
TopoDS_Shape* okShape)
{
isNodeOut = false;
if ( okShape )
*okShape = myWorkClassifiers[i]->Shape();
+ myWorkClassifiers[i]->GetParams( myU, myV );
break;
}
}
{
isNodeOut = false;
if ( okShape )
- *okShape = myWorkClassifiers[i]->Shape();
+ *okShape = myClassifiers[i].Shape();
+ myClassifiers[i].GetParams( myU, myV );
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 = new GeomAPI_ProjectPointOnSurf;
+ 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 = new GeomAPI_ProjectPointOnCurve;
+ 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() );
ElementsOnShape::Classifier::~Classifier()
{
delete mySolidClfr; mySolidClfr = 0;
+ delete myProjFace; myProjFace = 0;
+ delete myProjEdge; myProjEdge = 0;
+}
+
+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)
+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 )
{
- myProjFace.Perform( p );
- if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol )
+ if ( isOutOfBox( p )) return true;
+ myProjFace->Perform( p );
+ if ( myProjFace->IsDone() && myProjFace->LowerDistance() <= myTol )
{
// check relatively to the face
- Standard_Real u, v;
- myProjFace.LowerDistanceParameters(u, v);
- gp_Pnt2d aProjPnt (u, v);
+ myProjFace->LowerDistanceParameters( myU, myV );
+ gp_Pnt2d aProjPnt( myU, myV );
BRepClass_FaceClassifier aClsf ( TopoDS::Face( myShape ), aProjPnt, myTol );
if ( aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON )
return false;
return true;
}
-bool ElementsOnShape::Classifier::isOutOfEdge (const gp_Pnt& p)
+bool ElementsOnShape::Classifier::isOutOfEdge( const gp_Pnt& p )
{
- myProjEdge.Perform( p );
- return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol );
+ if ( isOutOfBox( p )) return true;
+ myProjEdge->Perform( p );
+ bool isOn = ( myProjEdge->NbPoints() > 0 && myProjEdge->LowerDistance() <= myTol );
+ if ( isOn )
+ myU = myProjEdge->LowerDistanceParameter();
+ return !isOn;
}
-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 )
{
if ( const SMDS_MeshElement* anElem = myMeshDS->FindElement( theId ))
{
- if ( anElem->GetType() == myType )
+ if ( myType == SMDSAbs_All || anElem->GetType() == myType )
{
if ( anElem->getshapeId() < 1 )
return myElementsOnShapePtr->IsSatisfy(theId);
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 )
if ( mySubShapesIDs.Contains( elem->getshapeId() ))
return true;
- if ( elem->GetType() != SMDSAbs_Node && elem->GetType() == myType )
+ if (( elem->GetType() != SMDSAbs_Node ) &&
+ ( myType == SMDSAbs_All || elem->GetType() == myType ))
{
SMDS_ElemIteratorPtr nodeItr = elem->nodesIterator();
while ( nodeItr->more() )
