Salome HOME
#18611 EDF 20669 - Quality of quadrangles
[modules/smesh.git] / src / Controls / SMESH_Controls.cxx
index ecd12a0b5d493ef48d1f6c28a1cc4a912620bc73..c5f8fbf67a3d4e10cbe82042b433ed54803751a6 100644 (file)
@@ -142,7 +142,7 @@ namespace {
     //  Case 1          Case 2
     //  |     |      |        |      |
     //  |     |      |        |      |
-    //  +-----+------+  +-----+------+ 
+    //  +-----+------+  +-----+------+
     //  |            |  |            |
     //  |            |  |            |
     // result should be 2 in both cases
@@ -491,7 +491,7 @@ double MaxElementLength2D::GetValue( const TSequenceOfXYZ& P )
   //     }
   // }
   // { // polygons
-    
+
   // }
 
   if( myPrecision >= 0 )
@@ -765,19 +765,9 @@ double AspectRatio::GetValue( long theId )
 {
   double aVal = 0;
   myCurrElement = myMesh->FindElement( theId );
-  if ( myCurrElement && myCurrElement->GetVtkType() == VTK_QUAD )
-  {
-    // issue 21723
-    vtkUnstructuredGrid* grid = const_cast<SMDS_Mesh*>( myMesh )->GetGrid();
-    if ( vtkCell* avtkCell = grid->GetCell( myCurrElement->GetVtkID() ))
-      aVal = Round( vtkMeshQuality::QuadAspectRatio( avtkCell ));
-  }
-  else
-  {
-    TSequenceOfXYZ P;
-    if ( GetPoints( myCurrElement, P ))
-      aVal = Round( GetValue( P ));
-  }
+  TSequenceOfXYZ P;
+  if ( GetPoints( myCurrElement, P ))
+    aVal = Round( GetValue( P ));
   return aVal;
 }
 
@@ -850,7 +840,7 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P )
     //
     // alpha = sqrt( 1/32 )
     // L = max( L1, L2, L3, L4, D1, D2 )
-    // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
+    // C1 = sqrt( L1^2 + L1^2 + L1^2 + L1^2 )
     // C2 = min( S1, S2, S3, S4 )
     // Li - lengths of the edges
     // Di - lengths of the diagonals
@@ -861,10 +851,10 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P )
                          Max( aLen[ 2 ],
                               Max( aLen[ 3 ],
                                    Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
-    double C1 = sqrt( aLen[0] * aLen[0] +
-                        aLen[1] * aLen[1] +
-                        aLen[2] * aLen[2] +
-                        aLen[3] * aLen[3] ) / 4. );
+    double C1 = sqrt( aLen[0] * aLen[0] +
+                      aLen[1] * aLen[1] +
+                      aLen[2] * aLen[2] +
+                      aLen[3] * aLen[3] );
     double C2 = Min( anArea[ 0 ],
                      Min( anArea[ 1 ],
                           Min( anArea[ 2 ], anArea[ 3 ] ) ) );
@@ -894,7 +884,7 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P )
     //
     // alpha = sqrt( 1/32 )
     // L = max( L1, L2, L3, L4, D1, D2 )
-    // C1 = sqrt( ( L1^2 + L1^2 + L1^2 + L1^2 ) / 4 )
+    // C1 = sqrt( L1^2 + L1^2 + L1^2 + L1^2 )
     // C2 = min( S1, S2, S3, S4 )
     // Li - lengths of the edges
     // Di - lengths of the diagonals
@@ -905,10 +895,10 @@ double AspectRatio::GetValue( const TSequenceOfXYZ& P )
                    Max( aLen[ 2 ],
                      Max( aLen[ 3 ],
                        Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
-    double C1 = sqrt( aLen[0] * aLen[0] +
-                        aLen[1] * aLen[1] +
-                        aLen[2] * aLen[2] +
-                        aLen[3] * aLen[3] ) / 4. );
+    double C1 = sqrt( aLen[0] * aLen[0] +
+                      aLen[1] * aLen[1] +
+                      aLen[2] * aLen[2] +
+                      aLen[3] * aLen[3] );
     double C2 = Min( anArea[ 0 ],
                   Min( anArea[ 1 ],
                     Min( anArea[ 2 ], anArea[ 3 ] ) ) );
@@ -2286,7 +2276,19 @@ bool BadOrientedVolume::IsSatisfy( long theId )
     return false;
 
   SMDS_VolumeTool vTool( myMesh->FindElement( theId ));
-  return !vTool.IsForward();
+
+  bool isOk = true;
+  if ( vTool.IsPoly() )
+  {
+    isOk = true;
+    for ( int i = 0; i < vTool.NbFaces() && isOk; ++i )
+      isOk = vTool.IsFaceExternal( i );
+  }
+  else
+  {
+    isOk = vTool.IsForward();
+  }
+  return !isOk;
 }
 
 SMDSAbs_ElementType BadOrientedVolume::GetType() const
@@ -2436,6 +2438,15 @@ SMDSAbs_ElementType CoincidentNodes::GetType() const
   return SMDSAbs_Node;
 }
 
+void CoincidentNodes::SetTolerance( const double theToler )
+{
+  if ( myToler != theToler )
+  {
+    SetMesh(0);
+    myToler = theToler;
+  }
+}
+
 void CoincidentNodes::SetMesh( const SMDS_Mesh* theMesh )
 {
   myMeshModifTracer.SetMesh( theMesh );
@@ -2565,18 +2576,14 @@ void FreeEdges::SetMesh( const SMDS_Mesh* theMesh )
 
 bool FreeEdges::IsFreeEdge( const SMDS_MeshNode** theNodes, const int theFaceId  )
 {
-  TColStd_MapOfInteger aMap;
-  for ( int i = 0; i < 2; i++ )
+  SMDS_ElemIteratorPtr anElemIter = theNodes[ 0 ]->GetInverseElementIterator(SMDSAbs_Face);
+  while( anElemIter->more() )
   {
-    SMDS_ElemIteratorPtr anElemIter = theNodes[ i ]->GetInverseElementIterator(SMDSAbs_Face);
-    while( anElemIter->more() )
+    if ( const SMDS_MeshElement* anElem = anElemIter->next())
     {
-      if ( const SMDS_MeshElement* anElem = anElemIter->next())
-      {
-        const int anId = anElem->GetID();
-        if ( anId != theFaceId && !aMap.Add( anId ))
-          return false;
-      }
+      const int anId = anElem->GetID();
+      if ( anId != theFaceId && anElem->GetNodeIndex( theNodes[1] ) >= 0 )
+        return false;
     }
   }
   return true;
@@ -2738,7 +2745,7 @@ bool FreeFaces::IsSatisfy( long theId )
       SMDS_MeshElement* aVol = (SMDS_MeshElement*)volItr->next();
       TItrMapOfVolume    itr = mapOfVol.insert( std::make_pair( aVol, 0 )).first;
       (*itr).second++;
-    } 
+    }
   }
   int nbVol = 0;
   TItrMapOfVolume volItr = mapOfVol.begin();
@@ -3154,7 +3161,7 @@ namespace
     double dot = v1 * v2; // cos * |v1| * |v2|
     double l1  = v1.SquareMagnitude();
     double l2  = v2.SquareMagnitude();
-    return (( dot * cos >= 0 ) && 
+    return (( dot * cos >= 0 ) &&
             ( dot * dot ) / l1 / l2 >= ( cos * cos ));
   }
 }
@@ -4276,6 +4283,8 @@ struct ElementsOnShape::OctreeClassifier : public SMESH_Octree
                     std::vector< ElementsOnShape::Classifier >&       cls );
   void GetClassifiersAtPoint( const gp_XYZ& p,
                               std::vector< ElementsOnShape::Classifier* >& classifiers );
+  size_t GetSize();
+
 protected:
   OctreeClassifier() {}
   SMESH_Octree* newChild() const { return new OctreeClassifier; }
@@ -4301,6 +4310,21 @@ ElementsOnShape::~ElementsOnShape()
 
 Predicate* ElementsOnShape::clone() const
 {
+  size_t size = sizeof( *this );
+  if ( myOctree )
+    size += myOctree->GetSize();
+  if ( !myClassifiers.empty() )
+    size += sizeof( myClassifiers[0] ) * myClassifiers.size();
+  if ( !myWorkClassifiers.empty() )
+    size += sizeof( myWorkClassifiers[0] ) * myWorkClassifiers.size();
+  if ( size > 1e+9 ) // 1G
+  {
+#ifdef _DEBUG_
+    std::cout << "Avoid ElementsOnShape::clone(), too large: " << size << " bytes " << std::endl;
+#endif
+    return 0;
+  }
+
   ElementsOnShape* cln = new ElementsOnShape();
   cln->SetAllNodes ( myAllNodesFlag );
   cln->SetTolerance( myToler );
@@ -4457,6 +4481,8 @@ bool ElementsOnShape::IsSatisfy (const SMDS_MeshElement* elem)
     for ( size_t i = 0; i < myClassifiers.size(); ++i )
       myWorkClassifiers[ i ] = & myClassifiers[ i ];
     myOctree = new OctreeClassifier( myWorkClassifiers );
+
+    SMESHUtils::FreeVector( myWorkClassifiers );
   }
 
   for ( int i = 0, nb = elem->NbNodes(); i < nb  && (isSatisfy == myAllNodesFlag); ++i )
@@ -4627,7 +4653,15 @@ void ElementsOnShape::Classifier::Init( const TopoDS_Shape& theShape,
     else
     {
       Bnd_Box box;
-      BRepBndLib::Add( myShape, box );
+      if ( myShape.ShapeType() == TopAbs_FACE )
+      {
+        BRepAdaptor_Surface SA( TopoDS::Face( myShape ), /*useBoundaries=*/false );
+        if ( SA.GetType() == GeomAbs_BSplineSurface )
+          BRepBndLib::AddOptimal( myShape, box,
+                                  /*useTriangulation=*/true, /*useShapeTolerance=*/true );
+      }
+      if ( box.IsVoid() )
+        BRepBndLib::Add( myShape, box );
       myBox.Clear();
       myBox.Add( box.CornerMin() );
       myBox.Add( box.CornerMax() );
@@ -4647,19 +4681,19 @@ ElementsOnShape::Classifier::~Classifier()
   delete mySolidClfr; mySolidClfr = 0;
 }
 
-bool ElementsOnShape::Classifier::isOutOfSolid (const gp_Pnt& p)
+bool ElementsOnShape::Classifier::isOutOfSolid( const gp_Pnt& p )
 {
   if ( isOutOfBox( p )) return true;
   mySolidClfr->Perform( p, myTol );
   return ( mySolidClfr->State() != TopAbs_IN && mySolidClfr->State() != TopAbs_ON );
 }
 
-bool ElementsOnShape::Classifier::isOutOfBox (const gp_Pnt& p)
+bool ElementsOnShape::Classifier::isOutOfBox( const gp_Pnt& p )
 {
   return myBox.IsOut( p.XYZ() );
 }
 
-bool ElementsOnShape::Classifier::isOutOfFace  (const gp_Pnt& p)
+bool ElementsOnShape::Classifier::isOutOfFace( const gp_Pnt& p )
 {
   if ( isOutOfBox( p )) return true;
   myProjFace.Perform( p );
@@ -4676,19 +4710,19 @@ bool ElementsOnShape::Classifier::isOutOfFace  (const gp_Pnt& p)
   return true;
 }
 
-bool ElementsOnShape::Classifier::isOutOfEdge  (const gp_Pnt& p)
+bool ElementsOnShape::Classifier::isOutOfEdge( const gp_Pnt& p )
 {
   if ( isOutOfBox( p )) return true;
   myProjEdge.Perform( p );
   return ! ( myProjEdge.NbPoints() > 0 && myProjEdge.LowerDistance() <= myTol );
 }
 
-bool ElementsOnShape::Classifier::isOutOfVertex(const gp_Pnt& p)
+bool ElementsOnShape::Classifier::isOutOfVertex( const gp_Pnt& p )
 {
   return ( myVertexXYZ.Distance( p ) > myTol );
 }
 
-bool ElementsOnShape::Classifier::isBox (const TopoDS_Shape& theShape)
+bool ElementsOnShape::Classifier::isBox(const TopoDS_Shape& theShape )
 {
   TopTools_IndexedMapOfShape vMap;
   TopExp::MapShapes( theShape, TopAbs_VERTEX, vMap );
@@ -4770,6 +4804,19 @@ OctreeClassifier::GetClassifiersAtPoint( const gp_XYZ& point,
   }
 }
 
+size_t ElementsOnShape::OctreeClassifier::GetSize()
+{
+  size_t res = sizeof( *this );
+  if ( !myClassifiers.empty() )
+    res += sizeof( myClassifiers[0] ) * myClassifiers.size();
+
+  if ( !isLeaf() )
+    for (int i = 0; i < nbChildren(); i++)
+      res += ((OctreeClassifier*) myChildren[i])->GetSize();
+
+  return res;
+}
+
 void ElementsOnShape::OctreeClassifier::buildChildrenData()
 {
   // distribute myClassifiers among myChildren
@@ -4844,8 +4891,13 @@ BelongToGeom::BelongToGeom()
 
 Predicate* BelongToGeom::clone() const
 {
-  BelongToGeom* cln = new BelongToGeom( *this );
-  cln->myElementsOnShapePtr.reset( static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ));
+  BelongToGeom* cln = 0;
+  if ( myElementsOnShapePtr )
+    if ( ElementsOnShape* eos = static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ))
+    {
+      cln = new BelongToGeom( *this );
+      cln->myElementsOnShapePtr.reset( eos );
+    }
   return cln;
 }
 
@@ -4856,6 +4908,8 @@ void BelongToGeom::SetMesh( const SMDS_Mesh* theMesh )
     myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
     init();
   }
+  if ( myElementsOnShapePtr )
+    myElementsOnShapePtr->SetMesh( myMeshDS );
 }
 
 void BelongToGeom::SetGeom( const TopoDS_Shape& theShape )
@@ -5015,8 +5069,13 @@ LyingOnGeom::LyingOnGeom()
 
 Predicate* LyingOnGeom::clone() const
 {
-  LyingOnGeom* cln = new LyingOnGeom( *this );
-  cln->myElementsOnShapePtr.reset( static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ));
+  LyingOnGeom* cln = 0;
+  if ( myElementsOnShapePtr )
+    if ( ElementsOnShape* eos = static_cast<ElementsOnShape*>( myElementsOnShapePtr->clone() ))
+    {
+      cln = new LyingOnGeom( *this );
+      cln->myElementsOnShapePtr.reset( eos );
+    }
   return cln;
 }
 
@@ -5027,6 +5086,8 @@ void LyingOnGeom::SetMesh( const SMDS_Mesh* theMesh )
     myMeshDS = dynamic_cast<const SMESHDS_Mesh*>(theMesh);
     init();
   }
+  if ( myElementsOnShapePtr )
+    myElementsOnShapePtr->SetMesh( myMeshDS );
 }
 
 void LyingOnGeom::SetGeom( const TopoDS_Shape& theShape )