-// Copyright (C) 2003 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
-// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
+// Copyright (C) 2007-2013 CEA/DEN, EDF R&D, OPEN CASCADE
//
-// This library is free software; you can redistribute it and/or
-// modify it under the terms of the GNU Lesser General Public
-// License as published by the Free Software Foundation; either
-// version 2.1 of the License.
+// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
+// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
//
-// This library is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-// Lesser General Public License for more details.
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License.
//
-// You should have received a copy of the GNU Lesser General Public
-// License along with this library; if not, write to the Free Software
-// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// Lesser General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
//
-// See http://www.opencascade.org/SALOME/ or email : webmaster.salome@opencascade.org
#include "SMESH_ControlsDef.hxx"
-#include <set>
+#include "SMDS_BallElement.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_Mesh.hxx"
+#include "SMESH_OctreeNode.hxx"
+#include "SMESH_MeshAlgos.hxx"
+#include <BRepAdaptor_Surface.hxx>
+#include <BRepClass_FaceClassifier.hxx>
#include <BRep_Tool.hxx>
-#include <gp_Ax3.hxx>
-#include <gp_Cylinder.hxx>
-#include <gp_Dir.hxx>
-#include <gp_Pnt.hxx>
-#include <gp_Pln.hxx>
-#include <gp_Vec.hxx>
-#include <gp_XYZ.hxx>
-#include <Geom_Plane.hxx>
#include <Geom_CylindricalSurface.hxx>
+#include <Geom_Plane.hxx>
+#include <Geom_Surface.hxx>
#include <Precision.hxx>
-#include <TColgp_Array1OfXYZ.hxx>
+#include <TColStd_MapIteratorOfMapOfInteger.hxx>
#include <TColStd_MapOfInteger.hxx>
#include <TColStd_SequenceOfAsciiString.hxx>
-#include <TColStd_MapIteratorOfMapOfInteger.hxx>
+#include <TColgp_Array1OfXYZ.hxx>
#include <TopAbs.hxx>
#include <TopoDS.hxx>
+#include <TopoDS_Edge.hxx>
#include <TopoDS_Face.hxx>
+#include <TopoDS_Iterator.hxx>
#include <TopoDS_Shape.hxx>
+#include <TopoDS_Vertex.hxx>
+#include <gp_Ax3.hxx>
+#include <gp_Cylinder.hxx>
+#include <gp_Dir.hxx>
+#include <gp_Pln.hxx>
+#include <gp_Pnt.hxx>
+#include <gp_Vec.hxx>
+#include <gp_XYZ.hxx>
-#include "SMDS_Mesh.hxx"
-#include "SMDS_Iterator.hxx"
-#include "SMDS_MeshElement.hxx"
-#include "SMDS_MeshNode.hxx"
-#include "SMDS_VolumeTool.hxx"
-#include "SMDS_QuadraticFaceOfNodes.hxx"
-#include "SMDS_QuadraticEdge.hxx"
+#include <vtkMeshQuality.h>
+
+#include <set>
+#include <limits>
+#include <Basics_Utils.hxx>
/*
AUXILIARY METHODS
*/
-namespace{
+namespace {
+
+ const double theEps = 1e-100;
+ const double theInf = 1e+100;
+
+ inline gp_XYZ gpXYZ(const SMDS_MeshNode* aNode )
+ {
+ return gp_XYZ(aNode->X(), aNode->Y(), aNode->Z() );
+ }
+
inline double getAngle( const gp_XYZ& P1, const gp_XYZ& P2, const gp_XYZ& P3 )
{
gp_Vec v1( P1 - P2 ), v2( P3 - P2 );
}
}
}
- int aResult = max ( aResult0, aResult1 );
+ int aResult = std::max ( aResult0, aResult1 );
// TColStd_MapOfInteger aMap;
// SMDS_ElemIteratorPtr anIter = anEdge->nodesIterator();
// if ( anIter != 0 ) {
// while( anIter->more() ) {
-// const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
-// if ( aNode == 0 )
-// return 0;
-// SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
-// while( anElemIter->more() ) {
-// const SMDS_MeshElement* anElem = anElemIter->next();
-// if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
-// int anId = anElem->GetID();
-
-// if ( anIter->more() ) // i.e. first node
-// aMap.Add( anId );
-// else if ( aMap.Contains( anId ) )
-// aResult++;
-// }
-// }
+// const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
+// if ( aNode == 0 )
+// return 0;
+// SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
+// while( anElemIter->more() ) {
+// const SMDS_MeshElement* anElem = anElemIter->next();
+// if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
+// int anId = anElem->GetID();
+
+// if ( anIter->more() ) // i.e. first node
+// aMap.Add( anId );
+// else if ( aMap.Contains( anId ) )
+// aResult++;
+// }
+// }
// }
// }
return aResult;
}
+ gp_XYZ getNormale( const SMDS_MeshFace* theFace, bool* ok=0 )
+ {
+ int aNbNode = theFace->NbNodes();
+
+ gp_XYZ q1 = gpXYZ( theFace->GetNode(1)) - gpXYZ( theFace->GetNode(0));
+ gp_XYZ q2 = gpXYZ( theFace->GetNode(2)) - gpXYZ( theFace->GetNode(0));
+ gp_XYZ n = q1 ^ q2;
+ if ( aNbNode > 3 ) {
+ gp_XYZ q3 = gpXYZ( theFace->GetNode(3)) - gpXYZ( theFace->GetNode(0));
+ n += q2 ^ q3;
+ }
+ double len = n.Modulus();
+ bool zeroLen = ( len <= numeric_limits<double>::min());
+ if ( !zeroLen )
+ n /= len;
+
+ if (ok) *ok = !zeroLen;
+
+ return n;
+ }
}
using namespace SMESH::Controls;
/*
- FUNCTORS
-*/
+ * FUNCTORS
+ */
+//================================================================================
/*
Class : NumericalFunctor
Description : Base class for numerical functors
*/
+//================================================================================
+
NumericalFunctor::NumericalFunctor():
myMesh(NULL)
{
if ( myMesh == 0 )
return false;
- return GetPoints( myMesh->FindElement( theId ), theRes );
+ const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
+ if ( !anElem || anElem->GetType() != this->GetType() )
+ return false;
+
+ return GetPoints( anElem, theRes );
}
bool NumericalFunctor::GetPoints(const SMDS_MeshElement* anElem,
{
theRes.clear();
- if ( anElem == 0)
+ if ( anElem == 0 )
return false;
theRes.reserve( anElem->NbNodes() );
if ( anElem->IsQuadratic() ) {
switch ( anElem->GetType() ) {
case SMDSAbs_Edge:
- anIter = static_cast<const SMDS_QuadraticEdge*>
+ anIter = dynamic_cast<const SMDS_VtkEdge*>
(anElem)->interlacedNodesElemIterator();
break;
case SMDSAbs_Face:
- anIter = static_cast<const SMDS_QuadraticFaceOfNodes*>
+ anIter = dynamic_cast<const SMDS_VtkFace*>
(anElem)->interlacedNodesElemIterator();
break;
default:
void NumericalFunctor::SetPrecision( const long thePrecision )
{
myPrecision = thePrecision;
+ myPrecisionValue = pow( 10., (double)( myPrecision ) );
}
double NumericalFunctor::GetValue( long theId )
{
+ double aVal = 0;
+
+ myCurrElement = myMesh->FindElement( theId );
+
TSequenceOfXYZ P;
if ( GetPoints( theId, P ))
+ aVal = Round( GetValue( P ));
+
+ return aVal;
+}
+
+double NumericalFunctor::Round( const double & aVal )
+{
+ return ( myPrecision >= 0 ) ? floor( aVal * myPrecisionValue + 0.5 ) / myPrecisionValue : aVal;
+}
+
+//================================================================================
+/*!
+ * \brief Return histogram of functor values
+ * \param nbIntervals - number of intervals
+ * \param nbEvents - number of mesh elements having values within i-th interval
+ * \param funValues - boundaries of intervals
+ * \param elements - elements to check vulue of; empty list means "of all"
+ * \param minmax - boundaries of diapason of values to divide into intervals
+ */
+//================================================================================
+
+void NumericalFunctor::GetHistogram(int nbIntervals,
+ std::vector<int>& nbEvents,
+ std::vector<double>& funValues,
+ const vector<int>& elements,
+ const double* minmax,
+ const bool isLogarithmic)
+{
+ if ( nbIntervals < 1 ||
+ !myMesh ||
+ !myMesh->GetMeshInfo().NbElements( GetType() ))
+ return;
+ nbEvents.resize( nbIntervals, 0 );
+ funValues.resize( nbIntervals+1 );
+
+ // get all values sorted
+ std::multiset< double > values;
+ if ( elements.empty() )
{
- double aVal = GetValue( P );
- if ( myPrecision >= 0 )
+ SMDS_ElemIteratorPtr elemIt = myMesh->elementsIterator(GetType());
+ while ( elemIt->more() )
+ values.insert( GetValue( elemIt->next()->GetID() ));
+ }
+ else
+ {
+ vector<int>::const_iterator id = elements.begin();
+ for ( ; id != elements.end(); ++id )
+ values.insert( GetValue( *id ));
+ }
+
+ if ( minmax )
+ {
+ funValues[0] = minmax[0];
+ funValues[nbIntervals] = minmax[1];
+ }
+ else
+ {
+ funValues[0] = *values.begin();
+ funValues[nbIntervals] = *values.rbegin();
+ }
+ // case nbIntervals == 1
+ if ( nbIntervals == 1 )
+ {
+ nbEvents[0] = values.size();
+ return;
+ }
+ // case of 1 value
+ if (funValues.front() == funValues.back())
+ {
+ nbEvents.resize( 1 );
+ nbEvents[0] = values.size();
+ funValues[1] = funValues.back();
+ funValues.resize( 2 );
+ }
+ // generic case
+ std::multiset< double >::iterator min = values.begin(), max;
+ for ( int i = 0; i < nbIntervals; ++i )
+ {
+ // find end value of i-th interval
+ double r = (i+1) / double(nbIntervals);
+ if (isLogarithmic && funValues.front() > 1e-07 && funValues.back() > 1e-07) {
+ double logmin = log10(funValues.front());
+ double lval = logmin + r * (log10(funValues.back()) - logmin);
+ funValues[i+1] = pow(10.0, lval);
+ }
+ else {
+ funValues[i+1] = funValues.front() * (1-r) + funValues.back() * r;
+ }
+
+ // count values in the i-th interval if there are any
+ if ( min != values.end() && *min <= funValues[i+1] )
{
- double prec = pow( 10., (double)( myPrecision ) );
- aVal = floor( aVal * prec + 0.5 ) / prec;
+ // find the first value out of the interval
+ max = values.upper_bound( funValues[i+1] ); // max is greater than funValues[i+1], or end()
+ nbEvents[i] = std::distance( min, max );
+ min = max;
}
- return aVal;
}
-
- return 0.;
+ // add values larger than minmax[1]
+ nbEvents.back() += std::distance( min, values.end() );
}
//=======================================================================
-//function : GetValue
-//purpose :
-//=======================================================================
+/*
+ Class : Volume
+ Description : Functor calculating volume of a 3D element
+*/
+//================================================================================
double Volume::GetValue( long theElementId )
{
return 0;
}
-//=======================================================================
-//function : GetBadRate
-//purpose : meaningless as it is not quality control functor
-//=======================================================================
-
double Volume::GetBadRate( double Value, int /*nbNodes*/ ) const
{
return Value;
}
+SMDSAbs_ElementType Volume::GetType() const
+{
+ return SMDSAbs_Volume;
+}
+
//=======================================================================
-//function : GetType
-//purpose :
+/*
+ Class : MaxElementLength2D
+ Description : Functor calculating maximum length of 2D element
+*/
+//================================================================================
+
+double MaxElementLength2D::GetValue( const TSequenceOfXYZ& P )
+{
+ if(P.size() == 0)
+ return 0.;
+ double aVal = 0;
+ int len = P.size();
+ 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 = Max(L1,Max(L2,L3));
+ }
+ else 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 ));
+ double D1 = getDistance(P( 1 ),P( 3 ));
+ double D2 = getDistance(P( 2 ),P( 4 ));
+ aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
+ }
+ else 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 = Max(L1,Max(L2,L3));
+ }
+ else if( len == 8 || len == 9 ) { // 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 ));
+ double D1 = getDistance(P( 1 ),P( 5 ));
+ double D2 = getDistance(P( 3 ),P( 7 ));
+ aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(D1,D2));
+ }
+
+ if( myPrecision >= 0 )
+ {
+ double prec = pow( 10., (double)myPrecision );
+ aVal = floor( aVal * prec + 0.5 ) / prec;
+ }
+ return aVal;
+}
+
+double MaxElementLength2D::GetValue( long theElementId )
+{
+ TSequenceOfXYZ P;
+ return GetPoints( theElementId, P ) ? GetValue(P) : 0.0;
+}
+
+double MaxElementLength2D::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+ return Value;
+}
+
+SMDSAbs_ElementType MaxElementLength2D::GetType() const
+{
+ return SMDSAbs_Face;
+}
+
//=======================================================================
+/*
+ Class : MaxElementLength3D
+ Description : Functor calculating maximum length of 3D element
+*/
+//================================================================================
-SMDSAbs_ElementType Volume::GetType() const
+double MaxElementLength3D::GetValue( long theElementId )
{
- return SMDSAbs_Volume;
+ TSequenceOfXYZ P;
+ if( GetPoints( theElementId, P ) ) {
+ double aVal = 0;
+ const SMDS_MeshElement* aElem = myMesh->FindElement( theElementId );
+ SMDSAbs_ElementType aType = aElem->GetType();
+ int len = P.size();
+ switch( aType ) {
+ case SMDSAbs_Volume:
+ if( len == 4 ) { // tetras
+ 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ break;
+ }
+ else if( len == 5 ) { // pyramids
+ 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(L7,L8));
+ break;
+ }
+ else if( len == 6 ) { // pentas
+ 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(Max(L7,L8),L9));
+ break;
+ }
+ else if( len == 8 ) { // hexas
+ 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 ));
+ double D1 = getDistance(P( 1 ),P( 7 ));
+ double D2 = getDistance(P( 2 ),P( 8 ));
+ double D3 = getDistance(P( 3 ),P( 5 ));
+ double D4 = getDistance(P( 4 ),P( 6 ));
+ aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
+ aVal = Max(aVal,Max(L11,L12));
+ aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
+ break;
+ }
+ else if( len == 12 ) { // hexagonal prism
+ for ( int i1 = 1; i1 < 12; ++i1 )
+ for ( int i2 = i1+1; i1 <= 12; ++i1 )
+ aVal = Max( aVal, getDistance(P( i1 ),P( i2 )));
+ break;
+ }
+ else if( len == 10 ) { // quadratic tetras
+ 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ break;
+ }
+ else if( len == 13 ) { // quadratic pyramids
+ 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(L7,L8));
+ break;
+ }
+ else if( len == 15 ) { // 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 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(Max(L7,L8),L9));
+ break;
+ }
+ else if( len == 20 || len == 27 ) { // quadratic hexas
+ 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 ));
+ double D1 = getDistance(P( 1 ),P( 7 ));
+ double D2 = getDistance(P( 2 ),P( 8 ));
+ double D3 = getDistance(P( 3 ),P( 5 ));
+ double D4 = getDistance(P( 4 ),P( 6 ));
+ aVal = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
+ aVal = Max(aVal,Max(L11,L12));
+ aVal = Max(aVal,Max(Max(D1,D2),Max(D3,D4)));
+ break;
+ }
+ else if( len > 1 && aElem->IsPoly() ) { // polys
+ // get the maximum distance between all pairs of nodes
+ for( int i = 1; i <= len; i++ ) {
+ for( int j = 1; j <= len; j++ ) {
+ if( j > i ) { // optimization of the loop
+ double D = getDistance( P(i), P(j) );
+ aVal = Max( aVal, D );
+ }
+ }
+ }
+ }
+ }
+
+ if( myPrecision >= 0 )
+ {
+ double prec = pow( 10., (double)myPrecision );
+ aVal = floor( aVal * prec + 0.5 ) / prec;
+ }
+ return aVal;
+ }
+ return 0.;
+}
+
+double MaxElementLength3D::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+ return Value;
}
+SMDSAbs_ElementType MaxElementLength3D::GetType() const
+{
+ return SMDSAbs_Volume;
+}
+//=======================================================================
/*
Class : MinimumAngle
Description : Functor for calculation of minimum angle
*/
+//================================================================================
double MinimumAngle::GetValue( const TSequenceOfXYZ& P )
{
aMin = Min(aMin,A0);
}
- return aMin * 180.0 / PI;
+ return aMin * 180.0 / M_PI;
}
double MinimumAngle::GetBadRate( double Value, int nbNodes ) const
}
+//================================================================================
/*
Class : AspectRatio
Description : Functor for calculating aspect ratio
*/
+//================================================================================
+
+double AspectRatio::GetValue( long theId )
+{
+ 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 ));
+ }
+ return aVal;
+}
+
double AspectRatio::GetValue( const TSequenceOfXYZ& P )
{
// According to "Mesh quality control" by Nadir Bouhamau referring to
if ( nbNodes < 3 )
return 0;
- // Compute lengths of the sides
-
- 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 ) );
-
// Compute aspect ratio
- if ( nbNodes == 3 )
- {
+ 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 ) );
// 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 ) );
- if ( anArea <= Precision::Confusion() )
- return 0.;
-
+ if ( anArea <= theEps )
+ return theInf;
return alfa * maxLen * half_perimeter / anArea;
}
- else
- {
- // return aspect ratio of the worst triange which can be built
+ 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) );
+ 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);
+ 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);
+ 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);
+ 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) );
+ anArea[3] = getArea( P(2), P(3), P(4) );
+ // Q = alpha * L * C1 / C2, where
+ //
+ // alpha = sqrt( 1/32 )
+ // L = max( L1, L2, L3, L4, D1, D2 )
+ // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
+ // 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. );
+ double C2 = Min( anArea[ 0 ],
+ 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);
+ 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);
+ 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
- TSequenceOfXYZ triaPnts(3);
- // triangle on nodes 1 3 2
- triaPnts(1) = P(1);
- triaPnts(2) = P(3);
- triaPnts(3) = P(2);
- double ar = GetValue( triaPnts );
- // triangle on nodes 1 3 4
- triaPnts(3) = P(4);
- ar = Max ( ar, GetValue( triaPnts ));
- // triangle on nodes 1 2 4
- triaPnts(2) = P(2);
- ar = Max ( ar, GetValue( triaPnts ));
- // triangle on nodes 3 2 4
- triaPnts(1) = P(3);
- ar = Max ( ar, GetValue( triaPnts ));
-
- return ar;
+ std::vector< 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) );
+ anArea[3] = getArea( P(3), P(5), P(7) );
+ // Q = alpha * L * C1 / C2, where
+ //
+ // alpha = sqrt( 1/32 )
+ // L = max( L1, L2, L3, L4, D1, D2 )
+ // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
+ // 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. );
+ double C2 = Min( anArea[ 0 ],
+ Min( anArea[ 1 ],
+ Min( anArea[ 2 ], anArea[ 3 ] ) ) );
+ if ( C2 <= theEps )
+ return theInf;
+ return alpha * L * C1 / C2;
}
+ return 0;
}
double AspectRatio::GetBadRate( double Value, int /*nbNodes*/ ) const
{
// the aspect ratio is in the range [1.0,infinity]
+ // < 1.0 = very bad, zero area
// 1.0 = good
// infinity = bad
- return Value / 1000.;
+ return ( Value < 0.9 ) ? 1000 : Value / 1000.;
}
SMDSAbs_ElementType AspectRatio::GetType() const
}
+//================================================================================
/*
Class : AspectRatio3D
Description : Functor for calculating aspect ratio
*/
+//================================================================================
+
namespace{
inline double getHalfPerimeter(double theTria[3]){
inline double getArea(double theHalfPerim, double theTria[3]){
return sqrt(theHalfPerim*
- (theHalfPerim-theTria[0])*
- (theHalfPerim-theTria[1])*
- (theHalfPerim-theTria[2]));
+ (theHalfPerim-theTria[0])*
+ (theHalfPerim-theTria[1])*
+ (theHalfPerim-theTria[2]));
}
inline double getVolume(double theLen[6]){
inline double getMaxHeight(double theLen[6])
{
- double aHeight = max(theLen[0],theLen[1]);
- aHeight = max(aHeight,theLen[2]);
- aHeight = max(aHeight,theLen[3]);
- aHeight = max(aHeight,theLen[4]);
- aHeight = max(aHeight,theLen[5]);
+ double aHeight = std::max(theLen[0],theLen[1]);
+ aHeight = std::max(aHeight,theLen[2]);
+ aHeight = std::max(aHeight,theLen[3]);
+ aHeight = std::max(aHeight,theLen[4]);
+ aHeight = std::max(aHeight,theLen[5]);
return aHeight;
}
}
+double AspectRatio3D::GetValue( long theId )
+{
+ double aVal = 0;
+ myCurrElement = myMesh->FindElement( theId );
+ if ( myCurrElement && myCurrElement->GetVtkType() == VTK_TETRA )
+ {
+ // 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() ))
+ aVal = Round( vtkMeshQuality::TetAspectRatio( avtkCell ));
+ }
+ else
+ {
+ TSequenceOfXYZ P;
+ if ( GetPoints( myCurrElement, P ))
+ aVal = Round( GetValue( P ));
+ }
+ return aVal;
+}
+
double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
{
double aQuality = 0.0;
+ if(myCurrElement->IsPoly()) return aQuality;
+
int nbNodes = P.size();
- switch(nbNodes){
+
+ 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 return aQuality;
+ }
+
+ switch(nbNodes) {
case 4:{
double aLen[6] = {
getDistance(P( 1 ),P( 2 )), // a
//double aVolume = getVolume(aLen);
double aHeight = getMaxHeight(aLen);
static double aCoeff = sqrt(2.0)/12.0;
- aQuality = aCoeff*aHeight*aSumArea/aVolume;
+ if ( aVolume > DBL_MIN )
+ aQuality = aCoeff*aHeight*aSumArea/aVolume;
break;
}
case 5:{
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 5 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 4 ),P( 5 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
break;
}
case 6:{
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 6 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 6 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 2 ),P( 5 ),P( 4 ),P( 3 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
break;
}
case 8:{
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 7 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 8 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 3 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 4 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 7 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 6 ),P( 8 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 3 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 4 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 7 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 2 ),P( 6 ),P( 5 ),P( 8 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 1 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 2 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 5 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 8 ),P( 6 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 1 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 2 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 5 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 7 ),P( 6 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 1 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 5 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 6 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 4 ),P( 8 ),P( 7 ),P( 2 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 4 ),P( 5 ),P( 8 ),P( 2 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 4 ),P( 5 ),P( 3 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 3 ),P( 6 ),P( 7 ),P( 1 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 2 ),P( 3 ),P( 6 ),P( 4 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 5 ),P( 6 ),P( 8 ),P( 3 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 7 ),P( 8 ),P( 6 ),P( 1 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 7 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
{
gp_XYZ aXYZ[4] = {P( 3 ),P( 4 ),P( 2 ),P( 5 )};
- aQuality = max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
}
break;
}
- }
+ 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);
+ }
+ {
+ gp_XYZ aXYZ[8] = {P( 2 ),P( 3 ),P( 5 ),P( 6 ),P( 8 ),P( 9 ),P( 11 ),P( 12 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
+ }
+ {
+ gp_XYZ aXYZ[8] = {P( 3 ),P( 4 ),P( 6 ),P( 1 ),P( 9 ),P( 10 ),P( 12 ),P( 7 )};
+ aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
+ }
+ break;
+ } // switch(nbNodes)
+
if ( nbNodes > 4 ) {
// avaluate aspect ratio of quadranle faces
AspectRatio aspect2D;
const int* pInd = SMDS_VolumeTool::GetFaceNodesIndices( type, i, true );
for ( int p = 0; p < 4; ++p ) // loop on nodes of a quadranle face
points( p + 1 ) = P( pInd[ p ] + 1 );
- aQuality = max( aQuality, aspect2D.GetValue( points ));
+ aQuality = std::max( aQuality, aspect2D.GetValue( points ));
}
}
return aQuality;
}
+//================================================================================
/*
Class : Warping
Description : Functor for calculating warping
*/
+//================================================================================
+
double Warping::GetValue( const TSequenceOfXYZ& P )
{
if ( P.size() != 4 )
return 0;
- gp_XYZ G = ( P( 1 ) + P( 2 ) + P( 3 ) + P( 4 ) ) / 4;
+ 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 );
- return Max( Max( A1, A2 ), Max( A3, A4 ) );
+ 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::ComputeA( const gp_XYZ& thePnt1,
double aLen1 = gp_Pnt( thePnt1 ).Distance( gp_Pnt( thePnt2 ) );
double aLen2 = gp_Pnt( thePnt2 ).Distance( gp_Pnt( thePnt3 ) );
double L = Min( aLen1, aLen2 ) * 0.5;
- if ( L < Precision::Confusion())
- return 0.;
+ if ( L < theEps )
+ return theInf;
gp_XYZ GI = ( thePnt2 + thePnt1 ) / 2. - theG;
gp_XYZ GJ = ( thePnt3 + thePnt2 ) / 2. - theG;
gp_XYZ N = GI.Crossed( GJ );
if ( N.Modulus() < gp::Resolution() )
- return PI / 2;
+ return M_PI / 2;
N.Normalize();
double H = ( thePnt2 - theG ).Dot( N );
- return asin( fabs( H / L ) ) * 180 / PI;
+ return asin( fabs( H / L ) ) * 180. / M_PI;
}
double Warping::GetBadRate( double Value, int /*nbNodes*/ ) const
}
+//================================================================================
/*
Class : Taper
Description : Functor for calculating taper
*/
+//================================================================================
+
double Taper::GetValue( const TSequenceOfXYZ& P )
{
if ( P.size() != 4 )
- return 0;
+ return 0.;
// Compute taper
- double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ) / 2;
- double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ) / 2;
- double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ) / 2;
- double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2;
+ double J1 = getArea( P( 4 ), P( 1 ), P( 2 ) ) / 2.;
+ double J2 = getArea( P( 3 ), P( 1 ), P( 2 ) ) / 2.;
+ double J3 = getArea( P( 2 ), P( 3 ), P( 4 ) ) / 2.;
+ double J4 = getArea( P( 3 ), P( 4 ), P( 1 ) ) / 2.;
double JA = 0.25 * ( J1 + J2 + J3 + J4 );
- if ( JA <= Precision::Confusion() )
- return 0.;
+ if ( JA <= theEps )
+ return theInf;
double T1 = fabs( ( J1 - JA ) / JA );
double T2 = fabs( ( J2 - JA ) / JA );
double T3 = fabs( ( J3 - JA ) / JA );
double T4 = fabs( ( J4 - JA ) / JA );
- return Max( Max( T1, T2 ), Max( T3, T4 ) );
+ double val = Max( Max( T1, T2 ), Max( T3, T4 ) );
+
+ const double eps = 0.01;
+
+ return val < eps ? 0. : val;
}
double Taper::GetBadRate( double Value, int /*nbNodes*/ ) const
return SMDSAbs_Face;
}
-
+//================================================================================
/*
Class : Skew
Description : Functor for calculating skew in degrees
*/
+//================================================================================
+
static inline double skewAngle( const gp_XYZ& p1, const gp_XYZ& p2, const gp_XYZ& p3 )
{
- gp_XYZ p12 = ( p2 + p1 ) / 2;
- gp_XYZ p23 = ( p3 + p2 ) / 2;
- gp_XYZ p31 = ( p3 + p1 ) / 2;
+ gp_XYZ p12 = ( p2 + p1 ) / 2.;
+ gp_XYZ p23 = ( p3 + p2 ) / 2.;
+ gp_XYZ p31 = ( p3 + p1 ) / 2.;
gp_Vec v1( p31 - p2 ), v2( p12 - p23 );
- return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0 : v1.Angle( v2 );
+ return v1.Magnitude() < gp::Resolution() || v2.Magnitude() < gp::Resolution() ? 0. : v1.Angle( v2 );
}
double Skew::GetValue( const TSequenceOfXYZ& P )
{
if ( P.size() != 3 && P.size() != 4 )
- return 0;
+ return 0.;
// Compute skew
- static double PI2 = PI / 2;
+ const double PI2 = M_PI / 2.;
if ( P.size() == 3 )
{
double A0 = fabs( PI2 - skewAngle( P( 3 ), P( 1 ), P( 2 ) ) );
double A1 = fabs( PI2 - skewAngle( P( 1 ), P( 2 ), P( 3 ) ) );
double A2 = fabs( PI2 - skewAngle( P( 2 ), P( 3 ), P( 1 ) ) );
- return Max( A0, Max( A1, A2 ) ) * 180 / PI;
+ return Max( A0, Max( A1, A2 ) ) * 180. / M_PI;
}
else
{
- gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2;
- gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2;
- gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2;
- gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2;
+ gp_XYZ p12 = ( P( 1 ) + P( 2 ) ) / 2.;
+ gp_XYZ p23 = ( P( 2 ) + P( 3 ) ) / 2.;
+ gp_XYZ p34 = ( P( 3 ) + P( 4 ) ) / 2.;
+ gp_XYZ p41 = ( P( 4 ) + P( 1 ) ) / 2.;
gp_Vec v1( p34 - p12 ), v2( p23 - p41 );
double A = v1.Magnitude() <= gp::Resolution() || v2.Magnitude() <= gp::Resolution()
- ? 0 : fabs( PI2 - v1.Angle( v2 ) );
+ ? 0. : fabs( PI2 - v1.Angle( v2 ) );
- return A * 180 / PI;
+ double val = A * 180. / M_PI;
+
+ const double eps = 0.1; // val is in degrees
+
+ return val < eps ? 0. : val;
}
}
}
+//================================================================================
/*
Class : Area
Description : Functor for calculating area
*/
+//================================================================================
+
double Area::GetValue( const TSequenceOfXYZ& P )
{
- gp_Vec aVec1( P(2) - P(1) );
- gp_Vec aVec2( P(3) - P(1) );
- gp_Vec SumVec = aVec1 ^ aVec2;
- for (int i=4; i<=P.size(); i++) {
- gp_Vec aVec1( P(i-1) - P(1) );
- gp_Vec aVec2( P(i) - P(1) );
- gp_Vec tmp = aVec1 ^ aVec2;
- SumVec.Add(tmp);
+ double val = 0.0;
+ if ( P.size() > 2 ) {
+ gp_Vec aVec1( P(2) - P(1) );
+ gp_Vec aVec2( P(3) - P(1) );
+ gp_Vec SumVec = aVec1 ^ aVec2;
+ for (int i=4; i<=P.size(); i++) {
+ gp_Vec aVec1( P(i-1) - P(1) );
+ gp_Vec aVec2( P(i) - P(1) );
+ gp_Vec tmp = aVec1 ^ aVec2;
+ SumVec.Add(tmp);
+ }
+ val = SumVec.Magnitude() * 0.5;
}
- return SumVec.Magnitude() * 0.5;
+ return val;
}
double Area::GetBadRate( double Value, int /*nbNodes*/ ) const
return SMDSAbs_Face;
}
-
+//================================================================================
/*
Class : Length
- Description : Functor for calculating length off edge
+ Description : Functor for calculating length of edge
*/
+//================================================================================
+
double Length::GetValue( const TSequenceOfXYZ& P )
{
switch ( P.size() ) {
return SMDSAbs_Edge;
}
+//================================================================================
/*
Class : Length2D
Description : Functor for calculating length of edge
*/
+//================================================================================
double Length2D::GetValue( long theElementId)
{
case SMDSAbs_Node:
case SMDSAbs_Edge:
if (len == 2){
- aVal = getDistance( P( 1 ), P( 2 ) );
+ aVal = getDistance( P( 1 ), P( 2 ) );
break;
}
else if (len == 3){ // quadratic edge
- aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
+ aVal = getDistance(P( 1 ),P( 3 )) + getDistance(P( 3 ),P( 2 ));
break;
}
case SMDSAbs_Face:
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 = Max(L1,Max(L2,L3));
- break;
+ double L1 = getDistance(P( 1 ),P( 2 ));
+ double L2 = getDistance(P( 2 ),P( 3 ));
+ double L3 = getDistance(P( 3 ),P( 1 ));
+ aVal = Max(L1,Max(L2,L3));
+ break;
}
else 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 = Max(Max(L1,L2),Max(L3,L4));
- break;
+ 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 = Max(Max(L1,L2),Max(L3,L4));
+ break;
}
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 = Max(L1,Max(L2,L3));
+ 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 = Max(L1,Max(L2,L3));
//cout<<"L1="<<L1<<" L2="<<L2<<"L3="<<L3<<" aVal="<<aVal<<endl;
- break;
+ break;
}
else 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 = Max(Max(L1,L2),Max(L3,L4));
- break;
+ 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 = Max(Max(L1,L2),Max(L3,L4));
+ break;
}
case SMDSAbs_Volume:
if (len == 4){ // tetraidrs
- 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
- break;
+ 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ break;
}
else if (len == 5){ // piramids
- 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( 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
- aVal = Max(aVal,Max(L7,L8));
- break;
+ 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(L7,L8));
+ break;
}
else 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
- aVal = Max(aVal,Max(Max(L7,L8),L9));
- break;
+ 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(Max(L7,L8),L9));
+ break;
}
else if (len == 8){ // hexaider
- 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
- aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
- aVal = Max(aVal,Max(L11,L12));
- break;
+ 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
+ aVal = Max(aVal,Max(L11,L12));
+ break;
}
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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
- break;
+ 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ break;
}
else 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( 1 ));
- 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
- aVal = Max(aVal,Max(L7,L8));
- break;
+ 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(L7,L8));
+ break;
}
else 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
- aVal = Max(aVal,Max(Max(L7,L8),L9));
- break;
+ 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(Max(L7,L8),L9));
+ break;
}
else 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
- aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
- aVal = Max(aVal,Max(L11,L12));
- break;
+ 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 = Max(Max(Max(L1,L2),Max(L3,L4)),Max(L5,L6));
+ aVal = Max(aVal,Max(Max(L7,L8),Max(L9,L10)));
+ aVal = Max(aVal,Max(L11,L12));
+ break;
}
const SMDS_MeshFace* anElem = anIter->next();
if(anElem->IsQuadratic()) {
- const SMDS_QuadraticFaceOfNodes* F =
- static_cast<const SMDS_QuadraticFaceOfNodes*>(anElem);
+ const SMDS_VtkFace* F =
+ dynamic_cast<const SMDS_VtkFace*>(anElem);
// use special nodes iterator
- SMDS_NodeIteratorPtr anIter = F->interlacedNodesIterator();
+ SMDS_ElemIteratorPtr anIter = F->interlacedNodesElemIterator();
long aNodeId[4];
gp_Pnt P[4];
}
}
+//================================================================================
/*
Class : MultiConnection
Description : Functor for calculating number of faces conneted to the edge
*/
+//================================================================================
+
double MultiConnection::GetValue( const TSequenceOfXYZ& P )
{
return 0;
return SMDSAbs_Edge;
}
+//================================================================================
/*
Class : MultiConnection2D
Description : Functor for calculating number of faces conneted to the edge
*/
+//================================================================================
+
double MultiConnection2D::GetValue( const TSequenceOfXYZ& P )
{
return 0;
double MultiConnection2D::GetValue( long theElementId )
{
- TSequenceOfXYZ P;
int aResult = 0;
- if (GetPoints(theElementId,P)){
- const SMDS_MeshElement* anFaceElem = myMesh->FindElement( theElementId );
- SMDSAbs_ElementType aType = anFaceElem->GetType();
-
- int len = P.size();
-
- TColStd_MapOfInteger aMap;
- int aResult = 0;
+ const SMDS_MeshElement* aFaceElem = myMesh->FindElement(theElementId);
+ SMDSAbs_ElementType aType = aFaceElem->GetType();
- switch (aType){
- case SMDSAbs_All:
- case SMDSAbs_Node:
- case SMDSAbs_Edge:
- case SMDSAbs_Face:
- if (len == 3){ // triangles
- int Nb[3] = {0,0,0};
-
- int i=0;
- SMDS_ElemIteratorPtr anIter = anFaceElem->nodesIterator();
- if ( anIter != 0 ) {
- while( anIter->more() ) {
- const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
- if ( aNode == 0 ){
- break;
- }
- SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
- while( anElemIter->more() ) {
- const SMDS_MeshElement* anElem = anElemIter->next();
- if ( anElem != 0 && anElem->GetType() != SMDSAbs_Edge ) {
- int anId = anElem->GetID();
-
- if ( anIter->more() ) // i.e. first node
- aMap.Add( anId );
- else if ( aMap.Contains( anId ) ){
- Nb[i]++;
- }
- }
- else if ( anElem != 0 && anElem->GetType() == SMDSAbs_Edge ) i++;
- }
- }
- }
-
- aResult = Max(Max(Nb[0],Nb[1]),Nb[2]);
+ switch (aType) {
+ case SMDSAbs_Face:
+ {
+ int i = 0, len = aFaceElem->NbNodes();
+ SMDS_ElemIteratorPtr anIter = aFaceElem->nodesIterator();
+ if (!anIter) break;
+
+ const SMDS_MeshNode *aNode, *aNode0;
+ TColStd_MapOfInteger aMap, aMapPrev;
+
+ for (i = 0; i <= len; i++) {
+ aMapPrev = aMap;
+ aMap.Clear();
+
+ int aNb = 0;
+ if (anIter->more()) {
+ aNode = (SMDS_MeshNode*)anIter->next();
+ } else {
+ if (i == len)
+ aNode = aNode0;
+ else
+ break;
+ }
+ if (!aNode) break;
+ if (i == 0) aNode0 = aNode;
+
+ SMDS_ElemIteratorPtr anElemIter = aNode->GetInverseElementIterator();
+ while (anElemIter->more()) {
+ const SMDS_MeshElement* anElem = anElemIter->next();
+ if (anElem != 0 && anElem->GetType() == SMDSAbs_Face) {
+ int anId = anElem->GetID();
+
+ aMap.Add(anId);
+ if (aMapPrev.Contains(anId)) {
+ aNb++;
+ }
+ }
+ }
+ aResult = Max(aResult, aNb);
}
- break;
- case SMDSAbs_Volume:
- default: aResult=0;
}
-
+ break;
+ default:
+ aResult = 0;
}
- return aResult;//getNbMultiConnection( myMesh, theId );
+
+ return aResult;
}
double MultiConnection2D::GetBadRate( double Value, int /*nbNodes*/ ) const
const SMDS_MeshFace* anElem = anIter->next();
SMDS_ElemIteratorPtr aNodesIter;
if ( anElem->IsQuadratic() )
- aNodesIter = static_cast<const SMDS_QuadraticFaceOfNodes*>
+ aNodesIter = dynamic_cast<const SMDS_VtkFace*>
(anElem)->interlacedNodesElemIterator();
else
aNodesIter = anElem->nodesIterator();
Value aValue(aNodeId[1],aNodeId[2]);
MValues::iterator aItr = theValues.find(aValue);
if (aItr != theValues.end()){
- aItr->second += 1;
- //aNbConnects = nb;
+ aItr->second += 1;
+ //aNbConnects = nb;
}
else {
- theValues[aValue] = 1;
- //aNbConnects = 1;
+ theValues[aValue] = 1;
+ //aNbConnects = 1;
}
//cout << "NodeIds: "<<aNodeId[1]<<","<<aNodeId[2]<<" nbconn="<<aNbConnects<<endl;
aNodeId[1] = aNodeId[2];
}
+//================================================================================
+/*
+ Class : BallDiameter
+ Description : Functor returning diameter of a ball element
+*/
+//================================================================================
+
+double BallDiameter::GetValue( long theId )
+{
+ double diameter = 0;
+
+ if ( const SMDS_BallElement* ball =
+ dynamic_cast<const SMDS_BallElement*>( myMesh->FindElement( theId )))
+ {
+ diameter = ball->GetDiameter();
+ }
+ return diameter;
+}
+
+double BallDiameter::GetBadRate( double Value, int /*nbNodes*/ ) const
+{
+ // meaningless as it is not a quality control functor
+ return Value;
+}
+
+SMDSAbs_ElementType BallDiameter::GetType() const
+{
+ return SMDSAbs_Ball;
+}
+
+
/*
PREDICATES
*/
+//================================================================================
/*
Class : BadOrientedVolume
Description : Predicate bad oriented volumes
*/
+//================================================================================
BadOrientedVolume::BadOrientedVolume()
{
return SMDSAbs_Volume;
}
+/*
+ Class : BareBorderVolume
+*/
+bool BareBorderVolume::IsSatisfy(long theElementId )
+{
+ SMDS_VolumeTool myTool;
+ if ( myTool.Set( myMesh->FindElement(theElementId)))
+ {
+ for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
+ if ( myTool.IsFreeFace( iF ))
+ {
+ const SMDS_MeshNode** n = myTool.GetFaceNodes(iF);
+ vector< const SMDS_MeshNode*> nodes( n, n+myTool.NbFaceNodes(iF));
+ if ( !myMesh->FindElement( nodes, SMDSAbs_Face, /*Nomedium=*/false))
+ return true;
+ }
+ }
+ return false;
+}
+
+//================================================================================
+/*
+ Class : BareBorderFace
+*/
+//================================================================================
+
+bool BareBorderFace::IsSatisfy(long theElementId )
+{
+ bool ok = false;
+ if ( const SMDS_MeshElement* face = myMesh->FindElement(theElementId))
+ {
+ if ( face->GetType() == SMDSAbs_Face )
+ {
+ int nbN = face->NbCornerNodes();
+ for ( int i = 0; i < nbN && !ok; ++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 )
+ {
+ const int iQuad = face->IsQuadratic();
+ myLinkNodes.resize( 2 + iQuad);
+ myLinkNodes[0] = n1;
+ myLinkNodes[1] = n2;
+ if ( iQuad )
+ myLinkNodes[2] = face->GetNode( i+nbN );
+ ok = !myMesh->FindElement( myLinkNodes, SMDSAbs_Edge, /*noMedium=*/false);
+ }
+ }
+ }
+ }
+ return ok;
+}
+
+//================================================================================
+/*
+ Class : OverConstrainedVolume
+*/
+//================================================================================
+
+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)))
+ {
+ int nbSharedFaces = 0;
+ for ( int iF = 0; iF < myTool.NbFaces(); ++iF )
+ if ( !myTool.IsFreeFace( iF ) && ++nbSharedFaces > 1 )
+ break;
+ return ( nbSharedFaces == 1 );
+ }
+ return false;
+}
+//================================================================================
+/*
+ Class : OverConstrainedFace
+*/
+//================================================================================
+
+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 )
+ {
+ int nbSharedBorders = 0;
+ int nbN = face->NbCornerNodes();
+ 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;
+ }
+ return ( nbSharedBorders == 1 );
+ }
+ return false;
+}
+
+//================================================================================
+/*
+ Class : CoincidentNodes
+ Description : Predicate of Coincident nodes
+*/
+//================================================================================
+
+CoincidentNodes::CoincidentNodes()
+{
+ myToler = 1e-5;
+}
+
+bool CoincidentNodes::IsSatisfy( long theElementId )
+{
+ return myCoincidentIDs.Contains( theElementId );
+}
+
+SMDSAbs_ElementType CoincidentNodes::GetType() const
+{
+ return SMDSAbs_Node;
+}
+
+void CoincidentNodes::SetMesh( const SMDS_Mesh* theMesh )
+{
+ myMeshModifTracer.SetMesh( theMesh );
+ if ( myMeshModifTracer.IsMeshModified() )
+ {
+ TIDSortedNodeSet nodesToCheck;
+ SMDS_NodeIteratorPtr nIt = theMesh->nodesIterator(/*idInceasingOrder=*/true);
+ while ( nIt->more() )
+ nodesToCheck.insert( nodesToCheck.end(), nIt->next() );
+
+ list< list< const SMDS_MeshNode*> > nodeGroups;
+ SMESH_OctreeNode::FindCoincidentNodes ( nodesToCheck, &nodeGroups, myToler );
+
+ myCoincidentIDs.Clear();
+ list< list< const SMDS_MeshNode*> >::iterator groupIt = nodeGroups.begin();
+ for ( ; groupIt != nodeGroups.end(); ++groupIt )
+ {
+ list< const SMDS_MeshNode*>& coincNodes = *groupIt;
+ list< const SMDS_MeshNode*>::iterator n = coincNodes.begin();
+ for ( ; n != coincNodes.end(); ++n )
+ myCoincidentIDs.Add( (*n)->GetID() );
+ }
+ }
+}
+
+//================================================================================
+/*
+ Class : CoincidentElements
+ Description : Predicate of Coincident Elements
+ Note : This class is suitable only for visualization of Coincident Elements
+*/
+//================================================================================
+
+CoincidentElements::CoincidentElements()
+{
+ myMesh = 0;
+}
+
+void CoincidentElements::SetMesh( const SMDS_Mesh* theMesh )
+{
+ myMesh = theMesh;
+}
+
+bool CoincidentElements::IsSatisfy( long theElementId )
+{
+ if ( !myMesh ) return false;
+
+ if ( const SMDS_MeshElement* e = myMesh->FindElement( theElementId ))
+ {
+ if ( e->GetType() != GetType() ) return false;
+ set< const SMDS_MeshNode* > elemNodes( e->begin_nodes(), e->end_nodes() );
+ const int nbNodes = e->NbNodes();
+ SMDS_ElemIteratorPtr invIt = (*elemNodes.begin())->GetInverseElementIterator( GetType() );
+ while ( invIt->more() )
+ {
+ const SMDS_MeshElement* e2 = invIt->next();
+ if ( e2 == e || e2->NbNodes() != nbNodes ) continue;
+
+ bool sameNodes = true;
+ for ( size_t i = 0; i < elemNodes.size() && sameNodes; ++i )
+ sameNodes = ( elemNodes.count( e2->GetNode( i )));
+ if ( sameNodes )
+ return true;
+ }
+ }
+ return false;
+}
+
+SMDSAbs_ElementType CoincidentElements1D::GetType() const
+{
+ return SMDSAbs_Edge;
+}
+SMDSAbs_ElementType CoincidentElements2D::GetType() const
+{
+ return SMDSAbs_Face;
+}
+SMDSAbs_ElementType CoincidentElements3D::GetType() const
+{
+ return SMDSAbs_Volume;
+}
+
+
+//================================================================================
/*
Class : FreeBorders
Description : Predicate for free borders
*/
+//================================================================================
FreeBorders::FreeBorders()
{
}
+//================================================================================
/*
Class : FreeEdges
Description : Predicate for free Edges
*/
+//================================================================================
+
FreeEdges::FreeEdges()
{
myMesh = 0;
TColStd_MapOfInteger aMap;
for ( int i = 0; i < 2; i++ )
{
- SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator();
+ SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator(SMDSAbs_Face);
while( anElemIter->more() )
{
- const SMDS_MeshElement* anElem = anElemIter->next();
- if ( anElem != 0 && anElem->GetType() == SMDSAbs_Face )
+ if ( const SMDS_MeshElement* anElem = anElemIter->next())
{
- int anId = anElem->GetID();
-
- if ( i == 0 )
- aMap.Add( anId );
- else if ( aMap.Contains( anId ) && anId != theFaceId )
+ const int anId = anElem->GetID();
+ if ( anId != theFaceId && !aMap.Add( anId ))
return false;
}
}
SMDS_ElemIteratorPtr anIter;
if ( aFace->IsQuadratic() ) {
- anIter = static_cast<const SMDS_QuadraticFaceOfNodes*>
+ anIter = dynamic_cast<const SMDS_VtkFace*>
(aFace)->interlacedNodesElemIterator();
}
else {
anIter = aFace->nodesIterator();
}
- if ( anIter != 0 )
+ if ( !anIter )
return false;
int i = 0, nbNodes = aFace->NbNodes();
- vector <const SMDS_MeshNode*> aNodes( nbNodes+1 );
+ std::vector <const SMDS_MeshNode*> aNodes( nbNodes+1 );
while( anIter->more() )
{
const SMDS_MeshNode* aNode = (SMDS_MeshNode*)anIter->next();
return false;
}
-SMDSAbs_ElementType FreeEdges::GetType() const
-{
- return SMDSAbs_Face;
-}
+SMDSAbs_ElementType FreeEdges::GetType() const
+{
+ return SMDSAbs_Face;
+}
+
+FreeEdges::Border::Border(long theElemId, long thePntId1, long thePntId2):
+ myElemId(theElemId)
+{
+ myPntId[0] = thePntId1; myPntId[1] = thePntId2;
+ if(thePntId1 > thePntId2){
+ myPntId[1] = thePntId1; myPntId[0] = thePntId2;
+ }
+}
+
+bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
+ if(myPntId[0] < x.myPntId[0]) return true;
+ if(myPntId[0] == x.myPntId[0])
+ if(myPntId[1] < x.myPntId[1]) return true;
+ return false;
+}
+
+inline void UpdateBorders(const FreeEdges::Border& theBorder,
+ FreeEdges::TBorders& theRegistry,
+ FreeEdges::TBorders& theContainer)
+{
+ if(theRegistry.find(theBorder) == theRegistry.end()){
+ theRegistry.insert(theBorder);
+ theContainer.insert(theBorder);
+ }else{
+ theContainer.erase(theBorder);
+ }
+}
+
+void FreeEdges::GetBoreders(TBorders& theBorders)
+{
+ TBorders aRegistry;
+ SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
+ for(; 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();
+ long aNodeId[2];
+ 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);
+ }
+ Border aBorder(anElemId,aNodeId[0],aNodeId[1]);
+ UpdateBorders(aBorder,aRegistry,theBorders);
+ }
+}
+
+//================================================================================
+/*
+ Class : FreeNodes
+ Description : Predicate for free nodes
+*/
+//================================================================================
+
+FreeNodes::FreeNodes()
+{
+ myMesh = 0;
+}
+
+void FreeNodes::SetMesh( const SMDS_Mesh* theMesh )
+{
+ myMesh = theMesh;
+}
+
+bool FreeNodes::IsSatisfy( long theNodeId )
+{
+ const SMDS_MeshNode* aNode = myMesh->FindNode( theNodeId );
+ if (!aNode)
+ return false;
+
+ return (aNode->NbInverseElements() < 1);
+}
+
+SMDSAbs_ElementType FreeNodes::GetType() const
+{
+ return SMDSAbs_Node;
+}
+
+
+//================================================================================
+/*
+ Class : FreeFaces
+ Description : Predicate for free faces
+*/
+//================================================================================
+
+FreeFaces::FreeFaces()
+{
+ myMesh = 0;
+}
+
+void FreeFaces::SetMesh( const SMDS_Mesh* theMesh )
+{
+ myMesh = theMesh;
+}
+
+bool FreeFaces::IsSatisfy( long theId )
+{
+ if (!myMesh) return false;
+ // check that faces nodes refers to less than two common volumes
+ const SMDS_MeshElement* aFace = myMesh->FindElement( theId );
+ if ( !aFace || aFace->GetType() != SMDSAbs_Face )
+ return false;
+
+ int nbNode = aFace->NbNodes();
+
+ // collect volumes check that number of volumss with count equal nbNode not less than 2
+ typedef map< SMDS_MeshElement*, int > TMapOfVolume; // map of volume counters
+ typedef map< SMDS_MeshElement*, int >::iterator TItrMapOfVolume; // iterator
+ TMapOfVolume mapOfVol;
+
+ SMDS_ElemIteratorPtr nodeItr = aFace->nodesIterator();
+ while ( nodeItr->more() ) {
+ const SMDS_MeshNode* aNode = static_cast<const SMDS_MeshNode*>(nodeItr->next());
+ if ( !aNode ) continue;
+ SMDS_ElemIteratorPtr volItr = aNode->GetInverseElementIterator(SMDSAbs_Volume);
+ while ( volItr->more() ) {
+ SMDS_MeshElement* aVol = (SMDS_MeshElement*)volItr->next();
+ TItrMapOfVolume itr = mapOfVol.insert(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++;
+ // face is not free if number of volumes constructed on thier nodes more than one
+ return (nbVol < 2);
+}
+
+SMDSAbs_ElementType FreeFaces::GetType() const
+{
+ return SMDSAbs_Face;
+}
+
+//================================================================================
+/*
+ Class : LinearOrQuadratic
+ Description : Predicate to verify whether a mesh element is linear
+*/
+//================================================================================
+
+LinearOrQuadratic::LinearOrQuadratic()
+{
+ myMesh = 0;
+}
+
+void LinearOrQuadratic::SetMesh( const SMDS_Mesh* theMesh )
+{
+ myMesh = theMesh;
+}
+
+bool LinearOrQuadratic::IsSatisfy( long theId )
+{
+ if (!myMesh) return false;
+ const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
+ if ( !anElem || (myType != SMDSAbs_All && anElem->GetType() != myType) )
+ return false;
+ return (!anElem->IsQuadratic());
+}
+
+void LinearOrQuadratic::SetType( SMDSAbs_ElementType theType )
+{
+ myType = theType;
+}
+
+SMDSAbs_ElementType LinearOrQuadratic::GetType() const
+{
+ return myType;
+}
+
+//================================================================================
+/*
+ Class : GroupColor
+ Description : Functor for check color of group to whic mesh element belongs to
+*/
+//================================================================================
+
+GroupColor::GroupColor()
+{
+}
+
+bool GroupColor::IsSatisfy( long theId )
+{
+ return (myIDs.find( theId ) != myIDs.end());
+}
+
+void GroupColor::SetType( SMDSAbs_ElementType theType )
+{
+ myType = theType;
+}
+
+SMDSAbs_ElementType GroupColor::GetType() const
+{
+ return myType;
+}
+
+static bool isEqual( const Quantity_Color& theColor1,
+ const Quantity_Color& theColor2 )
+{
+ // tolerance to compare colors
+ const double tol = 5*1e-3;
+ return ( fabs( theColor1.Red() - theColor2.Red() ) < tol &&
+ fabs( theColor1.Green() - theColor2.Green() ) < tol &&
+ fabs( theColor1.Blue() - theColor2.Blue() ) < tol );
+}
+
+
+void GroupColor::SetMesh( const SMDS_Mesh* theMesh )
+{
+ myIDs.clear();
+
+ const SMESHDS_Mesh* aMesh = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
+ if ( !aMesh )
+ return;
+
+ int nbGrp = aMesh->GetNbGroups();
+ if ( !nbGrp )
+ return;
+
+ // iterates on groups and find necessary elements ids
+ const std::set<SMESHDS_GroupBase*>& aGroups = aMesh->GetGroups();
+ set<SMESHDS_GroupBase*>::const_iterator GrIt = aGroups.begin();
+ for (; GrIt != aGroups.end(); GrIt++) {
+ SMESHDS_GroupBase* aGrp = (*GrIt);
+ if ( !aGrp )
+ continue;
+ // check type and color of group
+ if ( !isEqual( myColor, aGrp->GetColor() ) )
+ continue;
+ if ( myType != SMDSAbs_All && myType != (SMDSAbs_ElementType)aGrp->GetType() )
+ continue;
+
+ 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++)
+ myIDs.insert( aGrp->GetID(i+1) );
+ }
+ }
+}
+
+void GroupColor::SetColorStr( const TCollection_AsciiString& theStr )
+{
+ Kernel_Utils::Localizer loc;
+ TCollection_AsciiString aStr = theStr;
+ aStr.RemoveAll( ' ' );
+ aStr.RemoveAll( '\t' );
+ for ( int aPos = aStr.Search( ";;" ); aPos != -1; aPos = aStr.Search( ";;" ) )
+ aStr.Remove( aPos, 2 );
+ Standard_Real clr[3];
+ clr[0] = clr[1] = clr[2] = 0.;
+ for ( int i = 0; i < 3; i++ ) {
+ TCollection_AsciiString tmpStr = aStr.Token( ";", i+1 );
+ if ( !tmpStr.IsEmpty() && tmpStr.IsRealValue() )
+ clr[i] = tmpStr.RealValue();
+ }
+ myColor = Quantity_Color( clr[0], clr[1], clr[2], Quantity_TOC_RGB );
+}
+
+//=======================================================================
+// name : GetRangeStr
+// Purpose : Get range as a string.
+// Example: "1,2,3,50-60,63,67,70-"
+//=======================================================================
+
+void GroupColor::GetColorStr( TCollection_AsciiString& theResStr ) const
+{
+ theResStr.Clear();
+ theResStr += TCollection_AsciiString( myColor.Red() );
+ theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Green() );
+ theResStr += TCollection_AsciiString( ";" ) + TCollection_AsciiString( myColor.Blue() );
+}
+
+//================================================================================
+/*
+ Class : ElemGeomType
+ Description : Predicate to check element geometry type
+*/
+//================================================================================
+
+ElemGeomType::ElemGeomType()
+{
+ myMesh = 0;
+ myType = SMDSAbs_All;
+ myGeomType = SMDSGeom_TRIANGLE;
+}
+
+void ElemGeomType::SetMesh( const SMDS_Mesh* theMesh )
+{
+ myMesh = theMesh;
+}
+
+bool ElemGeomType::IsSatisfy( long theId )
+{
+ if (!myMesh) return false;
+ const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
+ if ( !anElem )
+ return false;
+ const SMDSAbs_ElementType anElemType = anElem->GetType();
+ if ( myType != SMDSAbs_All && anElemType != myType )
+ return false;
+ bool isOk = ( anElem->GetGeomType() == myGeomType );
+ return isOk;
+}
+
+void ElemGeomType::SetType( SMDSAbs_ElementType theType )
+{
+ myType = theType;
+}
+
+SMDSAbs_ElementType ElemGeomType::GetType() const
+{
+ return myType;
+}
+
+void ElemGeomType::SetGeomType( SMDSAbs_GeometryType theType )
+{
+ myGeomType = theType;
+}
+
+SMDSAbs_GeometryType ElemGeomType::GetGeomType() const
+{
+ return myGeomType;
+}
+
+//================================================================================
+/*
+ Class : ElemEntityType
+ Description : Predicate to check element entity type
+*/
+//================================================================================
+
+ElemEntityType::ElemEntityType():
+ myMesh( 0 ),
+ myType( SMDSAbs_All ),
+ myEntityType( SMDSEntity_0D )
+{
+}
+
+void ElemEntityType::SetMesh( const SMDS_Mesh* theMesh )
+{
+ myMesh = theMesh;
+}
+
+bool ElemEntityType::IsSatisfy( long theId )
+{
+ if ( !myMesh ) return false;
+ if ( myType == SMDSAbs_Node )
+ return myMesh->FindNode( theId );
+ const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
+ return ( anElem &&
+ myEntityType == anElem->GetEntityType() );
+}
+
+void ElemEntityType::SetType( SMDSAbs_ElementType theType )
+{
+ myType = theType;
+}
+
+SMDSAbs_ElementType ElemEntityType::GetType() const
+{
+ return myType;
+}
+
+void ElemEntityType::SetElemEntityType( SMDSAbs_EntityType theEntityType )
+{
+ myEntityType = theEntityType;
+}
+
+SMDSAbs_EntityType ElemEntityType::GetElemEntityType() const
+{
+ return myEntityType;
+}
+
+//================================================================================
+/*!
+ * \brief Class ConnectedElements
+ */
+//================================================================================
+
+ConnectedElements::ConnectedElements():
+ myNodeID(0), myType( SMDSAbs_All ), myOkIDsReady( false ) {}
+
+SMDSAbs_ElementType ConnectedElements::GetType() const
+{ return myType; }
+
+int ConnectedElements::GetNode() const
+{ return myXYZ.empty() ? myNodeID : 0; } // myNodeID can be found by myXYZ
+
+std::vector<double> ConnectedElements::GetPoint() const
+{ return myXYZ; }
+
+void ConnectedElements::clearOkIDs()
+{ myOkIDsReady = false; myOkIDs.clear(); }
+
+void ConnectedElements::SetType( SMDSAbs_ElementType theType )
+{
+ if ( myType != theType || myMeshModifTracer.IsMeshModified() )
+ clearOkIDs();
+ myType = theType;
+}
+
+void ConnectedElements::SetMesh( const SMDS_Mesh* theMesh )
+{
+ myMeshModifTracer.SetMesh( theMesh );
+ if ( myMeshModifTracer.IsMeshModified() )
+ {
+ clearOkIDs();
+ if ( !myXYZ.empty() )
+ SetPoint( myXYZ[0], myXYZ[1], myXYZ[2] ); // find a node near myXYZ it in a new mesh
+ }
+}
+
+void ConnectedElements::SetNode( int nodeID )
+{
+ myNodeID = nodeID;
+ myXYZ.clear();
+
+ bool isSameDomain = false;
+ if ( myOkIDsReady && myMeshModifTracer.GetMesh() && !myMeshModifTracer.IsMeshModified() )
+ if ( const SMDS_MeshNode* n = myMeshModifTracer.GetMesh()->FindNode( myNodeID ))
+ {
+ SMDS_ElemIteratorPtr eIt = n->GetInverseElementIterator( myType );
+ while ( !isSameDomain && eIt->more() )
+ isSameDomain = IsSatisfy( eIt->next()->GetID() );
+ }
+ if ( !isSameDomain )
+ clearOkIDs();
+}
+
+void ConnectedElements::SetPoint( double x, double y, double z )
+{
+ myXYZ.resize(3);
+ myXYZ[0] = x;
+ myXYZ[1] = y;
+ myXYZ[2] = z;
+ myNodeID = 0;
+
+ bool isSameDomain = false;
+
+ // find myNodeID by myXYZ if possible
+ if ( myMeshModifTracer.GetMesh() )
+ {
+ auto_ptr<SMESH_ElementSearcher> searcher
+ ( SMESH_MeshAlgos::GetElementSearcher( (SMDS_Mesh&) *myMeshModifTracer.GetMesh() ));
+
+ vector< const SMDS_MeshElement* > foundElems;
+ searcher->FindElementsByPoint( gp_Pnt(x,y,z), SMDSAbs_All, foundElems );
+
+ if ( !foundElems.empty() )
+ {
+ myNodeID = foundElems[0]->GetNode(0)->GetID();
+ if ( myOkIDsReady && !myMeshModifTracer.IsMeshModified() )
+ isSameDomain = IsSatisfy( foundElems[0]->GetID() );
+ }
+ }
+ if ( !isSameDomain )
+ clearOkIDs();
+}
+
+bool ConnectedElements::IsSatisfy( long theElementId )
+{
+ // Here we do NOT check if the mesh has changed, we do it in Set...() only!!!
+
+ if ( !myOkIDsReady )
+ {
+ if ( !myMeshModifTracer.GetMesh() )
+ return false;
+ const SMDS_MeshNode* node0 = myMeshModifTracer.GetMesh()->FindNode( myNodeID );
+ if ( !node0 )
+ return false;
+
+ list< const SMDS_MeshNode* > nodeQueue( 1, node0 );
+ std::set< int > checkedNodeIDs;
+ // algo:
+ // foreach node in nodeQueue:
+ // foreach element sharing a node:
+ // add ID of an element of myType to myOkIDs;
+ // push all element nodes absent from checkedNodeIDs to nodeQueue;
+ while ( !nodeQueue.empty() )
+ {
+ const SMDS_MeshNode* node = nodeQueue.front();
+ nodeQueue.pop_front();
+
+ // loop on elements sharing the node
+ SMDS_ElemIteratorPtr eIt = node->GetInverseElementIterator();
+ while ( eIt->more() )
+ {
+ // keep elements of myType
+ const SMDS_MeshElement* element = eIt->next();
+ if ( element->GetType() == myType )
+ myOkIDs.insert( myOkIDs.end(), element->GetID() );
+
+ // enqueue nodes of the element
+ SMDS_ElemIteratorPtr nIt = element->nodesIterator();
+ while ( nIt->more() )
+ {
+ const SMDS_MeshNode* n = static_cast< const SMDS_MeshNode* >( nIt->next() );
+ if ( checkedNodeIDs.insert( n->GetID() ).second )
+ nodeQueue.push_back( n );
+ }
+ }
+ }
+ if ( myType == SMDSAbs_Node )
+ std::swap( myOkIDs, checkedNodeIDs );
+
+ size_t totalNbElems = myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType );
+ if ( myOkIDs.size() == totalNbElems )
+ myOkIDs.clear();
-FreeEdges::Border::Border(long theElemId, long thePntId1, long thePntId2):
- myElemId(theElemId)
-{
- myPntId[0] = thePntId1; myPntId[1] = thePntId2;
- if(thePntId1 > thePntId2){
- myPntId[1] = thePntId1; myPntId[0] = thePntId2;
+ myOkIDsReady = true;
}
-}
-bool FreeEdges::Border::operator<(const FreeEdges::Border& x) const{
- if(myPntId[0] < x.myPntId[0]) return true;
- if(myPntId[0] == x.myPntId[0])
- if(myPntId[1] < x.myPntId[1]) return true;
- return false;
+ return myOkIDs.empty() ? true : myOkIDs.count( theElementId );
}
-inline void UpdateBorders(const FreeEdges::Border& theBorder,
- FreeEdges::TBorders& theRegistry,
- FreeEdges::TBorders& theContainer)
+//================================================================================
+/*!
+ * \brief Class CoplanarFaces
+ */
+//================================================================================
+
+CoplanarFaces::CoplanarFaces()
+ : myFaceID(0), myToler(0)
{
- if(theRegistry.find(theBorder) == theRegistry.end()){
- theRegistry.insert(theBorder);
- theContainer.insert(theBorder);
- }else{
- theContainer.erase(theBorder);
- }
}
-
-void FreeEdges::GetBoreders(TBorders& theBorders)
+void CoplanarFaces::SetMesh( const SMDS_Mesh* theMesh )
{
- TBorders aRegistry;
- SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
- for(; anIter->more(); ){
- const SMDS_MeshFace* anElem = anIter->next();
- long anElemId = anElem->GetID();
- SMDS_ElemIteratorPtr aNodesIter;
- if ( anElem->IsQuadratic() )
- aNodesIter = static_cast<const SMDS_QuadraticFaceOfNodes*>(anElem)->
- interlacedNodesElemIterator();
- else
- aNodesIter = anElem->nodesIterator();
- long aNodeId[2];
- 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;
- //std::cout<<aBorder.myPntId[0]<<"; "<<aBorder.myPntId[1]<<"; "<<aBorder.myElemId<<endl;
- UpdateBorders(aBorder,aRegistry,theBorders);
+ myMeshModifTracer.SetMesh( theMesh );
+ if ( myMeshModifTracer.IsMeshModified() )
+ {
+ // Build a set of coplanar face ids
+
+ myCoplanarIDs.clear();
+
+ if ( !myMeshModifTracer.GetMesh() || !myFaceID || !myToler )
+ return;
+
+ const SMDS_MeshElement* face = myMeshModifTracer.GetMesh()->FindElement( myFaceID );
+ if ( !face || face->GetType() != SMDSAbs_Face )
+ return;
+
+ bool normOK;
+ gp_Vec myNorm = getNormale( static_cast<const SMDS_MeshFace*>(face), &normOK );
+ if (!normOK)
+ return;
+
+ const double radianTol = myToler * M_PI / 180.;
+ std::set< SMESH_TLink > checkedLinks;
+
+ std::list< pair< const SMDS_MeshElement*, gp_Vec > > faceQueue;
+ faceQueue.push_back( make_pair( face, myNorm ));
+ while ( !faceQueue.empty() )
+ {
+ face = faceQueue.front().first;
+ myNorm = faceQueue.front().second;
+ faceQueue.pop_front();
+
+ for ( int i = 0, nbN = face->NbCornerNodes(); i < nbN; ++i )
+ {
+ const SMDS_MeshNode* n1 = face->GetNode( i );
+ const SMDS_MeshNode* n2 = face->GetNode(( i+1 )%nbN);
+ if ( !checkedLinks.insert( SMESH_TLink( n1, n2 )).second )
+ continue;
+ SMDS_ElemIteratorPtr fIt = n1->GetInverseElementIterator(SMDSAbs_Face);
+ while ( fIt->more() )
+ {
+ const SMDS_MeshElement* f = fIt->next();
+ if ( f->GetNodeIndex( n2 ) > -1 )
+ {
+ gp_Vec norm = getNormale( static_cast<const SMDS_MeshFace*>(f), &normOK );
+ if (!normOK || myNorm.Angle( norm ) <= radianTol)
+ {
+ myCoplanarIDs.insert( f->GetID() );
+ faceQueue.push_back( make_pair( f, norm ));
+ }
+ }
+ }
+ }
}
- Border aBorder(anElemId,aNodeId[0],aNodeId[1]);
- //std::cout<<aBorder.myPntId[0]<<"; "<<aBorder.myPntId[1]<<"; "<<aBorder.myElemId<<endl;
- UpdateBorders(aBorder,aRegistry,theBorders);
}
- //std::cout<<"theBorders.size() = "<<theBorders.size()<<endl;
+}
+bool CoplanarFaces::IsSatisfy( long theElementId )
+{
+ return myCoplanarIDs.count( theElementId );
}
/*
- Class : RangeOfIds
- Description : Predicate for Range of Ids.
- Range may be specified with two ways.
- 1. Using AddToRange method
- 2. With SetRangeStr method. Parameter of this method is a string
- like as "1,2,3,50-60,63,67,70-"
+ *Class : RangeOfIds
+ *Description : Predicate for Range of Ids.
+ * Range may be specified with two ways.
+ * 1. Using AddToRange method
+ * 2. With SetRangeStr method. Parameter of this method is a string
+ * like as "1,2,3,50-60,63,67,70-"
*/
//=======================================================================
while ( aMinStr.Search( "-" ) != -1 ) aMinStr.RemoveAll( '-' );
while ( aMaxStr.Search( "-" ) != -1 ) aMaxStr.RemoveAll( '-' );
- if ( !aMinStr.IsEmpty() && !aMinStr.IsIntegerValue() ||
- !aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue() )
+ if ( (!aMinStr.IsEmpty() && !aMinStr.IsIntegerValue()) ||
+ (!aMaxStr.IsEmpty() && !aMaxStr.IsIntegerValue()) )
return false;
myMin.Append( aMinStr.IsEmpty() ? IntegerFirst() : aMinStr.IntegerValue() );
else
{
const SMDS_MeshElement* anElem = myMesh->FindElement( theId );
- if ( anElem == 0 || myType != anElem->GetType() && myType != SMDSAbs_All )
+ if ( anElem == 0 || (myType != anElem->GetType() && myType != SMDSAbs_All ))
return false;
}
return
myPredicate1 &&
myPredicate2 &&
- myPredicate1->IsSatisfy( theId ) ||
- myPredicate2->IsSatisfy( theId );
+ (myPredicate1->IsSatisfy( theId ) ||
+ myPredicate2->IsSatisfy( theId ));
}
myPredicate = thePredicate;
}
-template<class TElement, class TIterator, class TPredicate>
-inline void FillSequence(const TIterator& theIterator,
- TPredicate& thePredicate,
- Filter::TIdSequence& theSequence)
-{
- if ( theIterator ) {
- while( theIterator->more() ) {
- TElement anElem = theIterator->next();
- long anId = anElem->GetID();
- if ( thePredicate->IsSatisfy( anId ) )
- theSequence.push_back( anId );
- }
- }
-}
-
-void
-Filter::
-GetElementsId( const SMDS_Mesh* theMesh,
- PredicatePtr thePredicate,
- TIdSequence& theSequence )
+void Filter::GetElementsId( const SMDS_Mesh* theMesh,
+ PredicatePtr thePredicate,
+ TIdSequence& theSequence )
{
theSequence.clear();
thePredicate->SetMesh( theMesh );
- SMDSAbs_ElementType aType = thePredicate->GetType();
- switch(aType){
- case SMDSAbs_Node:
- FillSequence<const SMDS_MeshNode*>(theMesh->nodesIterator(),thePredicate,theSequence);
- break;
- case SMDSAbs_Edge:
- FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),thePredicate,theSequence);
- break;
- case SMDSAbs_Face:
- FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),thePredicate,theSequence);
- break;
- case SMDSAbs_Volume:
- FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),thePredicate,theSequence);
- break;
- case SMDSAbs_All:
- FillSequence<const SMDS_MeshElement*>(theMesh->edgesIterator(),thePredicate,theSequence);
- FillSequence<const SMDS_MeshElement*>(theMesh->facesIterator(),thePredicate,theSequence);
- FillSequence<const SMDS_MeshElement*>(theMesh->volumesIterator(),thePredicate,theSequence);
- break;
+ 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 );
+ }
}
}
-void
-Filter::GetElementsId( const SMDS_Mesh* theMesh,
- Filter::TIdSequence& theSequence )
+void Filter::GetElementsId( const SMDS_Mesh* theMesh,
+ Filter::TIdSequence& theSequence )
{
GetElementsId(theMesh,myPredicate,theSequence);
}
}
}
-static gp_XYZ getNormale( const SMDS_MeshFace* theFace )
-{
- gp_XYZ n;
- int aNbNode = theFace->NbNodes();
- TColgp_Array1OfXYZ anArrOfXYZ(1,4);
- SMDS_ElemIteratorPtr aNodeItr = theFace->nodesIterator();
- int i = 1;
- for ( ; aNodeItr->more() && i <= 4; i++ ) {
- SMDS_MeshNode* aNode = (SMDS_MeshNode*)aNodeItr->next();
- anArrOfXYZ.SetValue(i, gp_XYZ( aNode->X(), aNode->Y(), aNode->Z() ) );
- }
-
- gp_XYZ q1 = anArrOfXYZ.Value(2) - anArrOfXYZ.Value(1);
- gp_XYZ q2 = anArrOfXYZ.Value(3) - anArrOfXYZ.Value(1);
- n = q1 ^ q2;
- if ( aNbNode > 3 ) {
- gp_XYZ q3 = anArrOfXYZ.Value(4) - anArrOfXYZ.Value(1);
- n += q2 ^ q3;
- }
- double len = n.Modulus();
- if ( len > 0 )
- n /= len;
-
- return n;
-}
-
bool ManifoldPart::findConnected
( const ManifoldPart::TDataMapFacePtrInt& theAllFacePtrInt,
SMDS_MeshFace* theStartFace,
void ManifoldPart::getFacesByLink( const ManifoldPart::Link& theLink,
ManifoldPart::TVectorOfFacePtr& theFaces ) const
{
- SMDS_Mesh::SetOfFaces aSetOfFaces;
+ 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.Add( aFace );
+ aSetOfFaces.insert( aFace );
}
// take all faces that shared second node
anItr = theLink.myNode2->facesIterator();
for ( ; anItr->more(); )
{
SMDS_MeshFace* aFace = (SMDS_MeshFace*)anItr->next();
- if ( aSetOfFaces.Contains( aFace ) )
+ if ( aSetOfFaces.count( aFace ) )
theFaces.push_back( aFace );
}
}
ElementsOnSurface::ElementsOnSurface()
{
- myMesh = 0;
myIds.Clear();
myType = SMDSAbs_All;
mySurf.Nullify();
myToler = Precision::Confusion();
+ myUseBoundaries = false;
}
ElementsOnSurface::~ElementsOnSurface()
{
- myMesh = 0;
}
void ElementsOnSurface::SetMesh( const SMDS_Mesh* theMesh )
{
- if ( myMesh == theMesh )
- return;
- myMesh = theMesh;
- myIds.Clear();
- process();
+ myMeshModifTracer.SetMesh( theMesh );
+ if ( myMeshModifTracer.IsMeshModified())
+ process();
}
bool ElementsOnSurface::IsSatisfy( long theElementId )
{ return myType; }
void ElementsOnSurface::SetTolerance( const double theToler )
-{ myToler = theToler; }
+{
+ if ( myToler != theToler )
+ myIds.Clear();
+ myToler = theToler;
+}
double ElementsOnSurface::GetTolerance() const
+{ return myToler; }
+
+void ElementsOnSurface::SetUseBoundaries( bool theUse )
{
- return myToler;
+ if ( myUseBoundaries != theUse ) {
+ myUseBoundaries = theUse;
+ SetSurface( mySurf, myType );
+ }
}
void ElementsOnSurface::SetSurface( const TopoDS_Shape& theShape,
const SMDSAbs_ElementType theType )
{
+ myIds.Clear();
myType = theType;
mySurf.Nullify();
if ( theShape.IsNull() || theShape.ShapeType() != TopAbs_FACE )
- {
- mySurf.Nullify();
return;
- }
- TopoDS_Face aFace = TopoDS::Face( theShape );
- mySurf = BRep_Tool::Surface( aFace );
+ mySurf = TopoDS::Face( theShape );
+ BRepAdaptor_Surface SA( mySurf, myUseBoundaries );
+ Standard_Real
+ u1 = SA.FirstUParameter(),
+ u2 = SA.LastUParameter(),
+ v1 = SA.FirstVParameter(),
+ v2 = SA.LastVParameter();
+ Handle(Geom_Surface) surf = BRep_Tool::Surface( mySurf );
+ myProjector.Init( surf, u1,u2, v1,v2 );
+ process();
}
void ElementsOnSurface::process()
if ( mySurf.IsNull() )
return;
- if ( myMesh == 0 )
+ if ( !myMeshModifTracer.GetMesh() )
return;
- if ( myType == SMDSAbs_Face || myType == SMDSAbs_All )
- {
- SMDS_FaceIteratorPtr anIter = myMesh->facesIterator();
- for(; anIter->more(); )
- process( anIter->next() );
- }
-
- if ( myType == SMDSAbs_Edge || myType == SMDSAbs_All )
- {
- SMDS_EdgeIteratorPtr anIter = myMesh->edgesIterator();
- for(; anIter->more(); )
- process( anIter->next() );
- }
+ myIds.ReSize( myMeshModifTracer.GetMesh()->GetMeshInfo().NbElements( myType ));
- if ( myType == SMDSAbs_Node )
- {
- SMDS_NodeIteratorPtr anIter = myMesh->nodesIterator();
- for(; anIter->more(); )
- process( anIter->next() );
- }
+ SMDS_ElemIteratorPtr anIter = myMeshModifTracer.GetMesh()->elementsIterator( myType );
+ for(; anIter->more(); )
+ process( anIter->next() );
}
void ElementsOnSurface::process( const SMDS_MeshElement* theElemPtr )
myIds.Add( theElemPtr->GetID() );
}
-bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode ) const
+bool ElementsOnSurface::isOnSurface( const SMDS_MeshNode* theNode )
{
if ( mySurf.IsNull() )
return false;
gp_Pnt aPnt( theNode->X(), theNode->Y(), theNode->Z() );
- double aToler2 = myToler * myToler;
- if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane)))
- {
- gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln();
- if ( aPln.SquareDistance( aPnt ) > aToler2 )
- return false;
+ // double aToler2 = myToler * myToler;
+// if ( mySurf->IsKind(STANDARD_TYPE(Geom_Plane)))
+// {
+// gp_Pln aPln = Handle(Geom_Plane)::DownCast(mySurf)->Pln();
+// if ( aPln.SquareDistance( aPnt ) > aToler2 )
+// return false;
+// }
+// else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)))
+// {
+// gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder();
+// double aRad = aCyl.Radius();
+// gp_Ax3 anAxis = aCyl.Position();
+// gp_XYZ aLoc = aCyl.Location().XYZ();
+// double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc );
+// double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc );
+// if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 )
+// return false;
+// }
+// else
+// return false;
+ myProjector.Perform( aPnt );
+ bool isOn = ( myProjector.IsDone() && myProjector.LowerDistance() <= myToler );
+
+ return isOn;
+}
+
+
+/*
+ ElementsOnShape
+*/
+
+ElementsOnShape::ElementsOnShape()
+ : //myMesh(0),
+ myType(SMDSAbs_All),
+ myToler(Precision::Confusion()),
+ myAllNodesFlag(false)
+{
+}
+
+ElementsOnShape::~ElementsOnShape()
+{
+ clearClassifiers();
+}
+
+SMDSAbs_ElementType ElementsOnShape::GetType() const
+{
+ return myType;
+}
+
+void ElementsOnShape::SetTolerance (const double theToler)
+{
+ if (myToler != theToler) {
+ myToler = theToler;
+ SetShape(myShape, myType);
}
- else if ( mySurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)))
+}
+
+double ElementsOnShape::GetTolerance() const
+{
+ return myToler;
+}
+
+void ElementsOnShape::SetAllNodes (bool theAllNodes)
+{
+ myAllNodesFlag = theAllNodes;
+}
+
+void ElementsOnShape::SetMesh (const SMDS_Mesh* theMesh)
+{
+ myMesh = theMesh;
+}
+
+void ElementsOnShape::SetShape (const TopoDS_Shape& theShape,
+ const SMDSAbs_ElementType theType)
+{
+ myType = theType;
+ myShape = theShape;
+ if ( myShape.IsNull() ) return;
+
+ TopTools_IndexedMapOfShape shapesMap;
+ TopAbs_ShapeEnum shapeTypes[4] = { TopAbs_SOLID, TopAbs_FACE, TopAbs_EDGE, TopAbs_VERTEX };
+ TopExp_Explorer sub;
+ for ( int i = 0; i < 4; ++i )
{
- gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(mySurf)->Cylinder();
- double aRad = aCyl.Radius();
- gp_Ax3 anAxis = aCyl.Position();
- gp_XYZ aLoc = aCyl.Location().XYZ();
- double aXDist = anAxis.XDirection().XYZ() * ( aPnt.XYZ() - aLoc );
- double aYDist = anAxis.YDirection().XYZ() * ( aPnt.XYZ() - aLoc );
- if ( fabs(aXDist*aXDist + aYDist*aYDist - aRad*aRad) > aToler2 )
- return false;
+ if ( shapesMap.IsEmpty() )
+ for ( sub.Init( myShape, shapeTypes[i] ); sub.More(); sub.Next() )
+ shapesMap.Add( sub.Current() );
+ if ( i > 0 )
+ for ( sub.Init( myShape, shapeTypes[i], shapeTypes[i-1] ); sub.More(); sub.Next() )
+ shapesMap.Add( sub.Current() );
}
- else
+
+ clearClassifiers();
+ myClassifiers.resize( shapesMap.Extent() );
+ for ( int i = 0; i < shapesMap.Extent(); ++i )
+ myClassifiers[ i ] = new TClassifier( shapesMap( i+1 ), myToler );
+}
+
+void ElementsOnShape::clearClassifiers()
+{
+ for ( size_t i = 0; i < myClassifiers.size(); ++i )
+ delete myClassifiers[ i ];
+ myClassifiers.clear();
+}
+
+bool ElementsOnShape::IsSatisfy (long elemId)
+{
+ const SMDS_MeshElement* elem =
+ ( myType == SMDSAbs_Node ? myMesh->FindNode( elemId ) : myMesh->FindElement( elemId ));
+ if ( !elem || myClassifiers.empty() )
return false;
+ for ( size_t i = 0; i < myClassifiers.size(); ++i )
+ {
+ SMDS_ElemIteratorPtr aNodeItr = elem->nodesIterator();
+ bool isSatisfy = myAllNodesFlag;
+
+ gp_XYZ centerXYZ (0, 0, 0);
+
+ while (aNodeItr->more() && (isSatisfy == myAllNodesFlag))
+ {
+ SMESH_TNodeXYZ aPnt ( aNodeItr->next() );
+ centerXYZ += aPnt;
+ isSatisfy = ! myClassifiers[i]->IsOut( aPnt );
+ }
+
+ // Check the center point for volumes MantisBug 0020168
+ if (isSatisfy &&
+ myAllNodesFlag &&
+ myClassifiers[i]->ShapeType() == TopAbs_SOLID)
+ {
+ centerXYZ /= elem->NbNodes();
+ isSatisfy = ! myClassifiers[i]->IsOut( centerXYZ );
+ }
+ if ( isSatisfy )
+ return true;
+ }
+
+ return false;
+}
+
+TopAbs_ShapeEnum ElementsOnShape::TClassifier::ShapeType() const
+{
+ return myShape.ShapeType();
+}
+
+bool ElementsOnShape::TClassifier::IsOut(const gp_Pnt& p)
+{
+ return (this->*myIsOutFun)( p );
+}
+
+void ElementsOnShape::TClassifier::Init (const TopoDS_Shape& theShape, double theTol)
+{
+ myShape = theShape;
+ myTol = theTol;
+ switch ( myShape.ShapeType() )
+ {
+ case TopAbs_SOLID: {
+ mySolidClfr.Load(theShape);
+ myIsOutFun = & ElementsOnShape::TClassifier::isOutOfSolid;
+ break;
+ }
+ case TopAbs_FACE: {
+ 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::TClassifier::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::TClassifier::isOutOfEdge;
+ break;
+ }
+ case TopAbs_VERTEX:{
+ myVertexXYZ = BRep_Tool::Pnt( TopoDS::Vertex( theShape ) );
+ myIsOutFun = & ElementsOnShape::TClassifier::isOutOfVertex;
+ break;
+ }
+ default:
+ throw SALOME_Exception("Programmer error in usage of ElementsOnShape::TClassifier");
+ }
+}
+
+bool ElementsOnShape::TClassifier::isOutOfSolid (const gp_Pnt& p)
+{
+ mySolidClfr.Perform( p, myTol );
+ return ( mySolidClfr.State() != TopAbs_IN && mySolidClfr.State() != TopAbs_ON );
+}
+
+bool ElementsOnShape::TClassifier::isOutOfFace (const gp_Pnt& p)
+{
+ myProjFace.Perform( p );
+ if ( myProjFace.IsDone() && myProjFace.LowerDistance() <= myTol )
+ {
+ // check relatively to the face
+ Quantity_Parameter u, v;
+ myProjFace.LowerDistanceParameters(u, v);
+ gp_Pnt2d aProjPnt (u, v);
+ BRepClass_FaceClassifier aClsf ( TopoDS::Face( myShape ), aProjPnt, myTol );
+ if ( aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON )
+ return false;
+ }
return true;
}
+
+bool ElementsOnShape::TClassifier::isOutOfEdge (const gp_Pnt& p)
+{
+ myProjEdge.Perform( p );
+ return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol );
+}
+
+bool ElementsOnShape::TClassifier::isOutOfVertex(const gp_Pnt& p)
+{
+ return ( myVertexXYZ.Distance( p ) > myTol );
+}
+
+
+TSequenceOfXYZ::TSequenceOfXYZ()
+{}
+
+TSequenceOfXYZ::TSequenceOfXYZ(size_type n) : myArray(n)
+{}
+
+TSequenceOfXYZ::TSequenceOfXYZ(size_type n, const gp_XYZ& t) : myArray(n,t)
+{}
+
+TSequenceOfXYZ::TSequenceOfXYZ(const TSequenceOfXYZ& theSequenceOfXYZ) : myArray(theSequenceOfXYZ.myArray)
+{}
+
+template <class InputIterator>
+TSequenceOfXYZ::TSequenceOfXYZ(InputIterator theBegin, InputIterator theEnd): myArray(theBegin,theEnd)
+{}
+
+TSequenceOfXYZ::~TSequenceOfXYZ()
+{}
+
+TSequenceOfXYZ& TSequenceOfXYZ::operator=(const TSequenceOfXYZ& theSequenceOfXYZ)
+{
+ myArray = theSequenceOfXYZ.myArray;
+ return *this;
+}
+
+gp_XYZ& TSequenceOfXYZ::operator()(size_type n)
+{
+ return myArray[n-1];
+}
+
+const gp_XYZ& TSequenceOfXYZ::operator()(size_type n) const
+{
+ return myArray[n-1];
+}
+
+void TSequenceOfXYZ::clear()
+{
+ myArray.clear();
+}
+
+void TSequenceOfXYZ::reserve(size_type n)
+{
+ myArray.reserve(n);
+}
+
+void TSequenceOfXYZ::push_back(const gp_XYZ& v)
+{
+ myArray.push_back(v);
+}
+
+TSequenceOfXYZ::size_type TSequenceOfXYZ::size() const
+{
+ return myArray.size();
+}
+
+TMeshModifTracer::TMeshModifTracer():
+ myMeshModifTime(0), myMesh(0)
+{
+}
+void TMeshModifTracer::SetMesh( const SMDS_Mesh* theMesh )
+{
+ if ( theMesh != myMesh )
+ myMeshModifTime = 0;
+ myMesh = theMesh;
+}
+bool TMeshModifTracer::IsMeshModified()
+{
+ bool modified = false;
+ if ( myMesh )
+ {
+ modified = ( myMeshModifTime != myMesh->GetMTime() );
+ myMeshModifTime = myMesh->GetMTime();
+ }
+ return modified;
+}
