Salome HOME
#18782 EDF 20946 - Free nodes with biquadratic
[modules/smesh.git] / src / Controls / SMESH_Controls.cxx
index 06c40dae2ad9c8b21d5ec8098df9528899d990c8..9a3f828fbedbe34c0cde2cd530adbcd71e1b52b8 100644 (file)
@@ -142,13 +142,13 @@ namespace {
     //  Case 1          Case 2
     //  |     |      |        |      |
     //  |     |      |        |      |
-    //  +-----+------+  +-----+------+ 
+    //  +-----+------+  +-----+------+
     //  |            |  |            |
     //  |            |  |            |
     // result should be 2 in both cases
     //
     int aResult0 = 0, aResult1 = 0;
-     // last node, it is a medium one in a quadratic edge
+    // last node, it is a medium one in a quadratic edge
     const SMDS_MeshNode* aLastNode = anEdge->GetNode( anEdge->NbNodes() - 1 );
     const SMDS_MeshNode*    aNode0 = anEdge->GetNode( 0 );
     const SMDS_MeshNode*    aNode1 = anEdge->GetNode( 1 );
@@ -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,24 +884,24 @@ 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
     // Si - areas of the triangles
     const double alpha = sqrt( 1 / 32. );
     double L = Max( aLen[ 0 ],
-                 Max( aLen[ 1 ],
-                   Max( aLen[ 2 ],
-                     Max( aLen[ 3 ],
-                       Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
-    double C1 = sqrt( aLen[0] * aLen[0] +
-                        aLen[1] * aLen[1] +
-                        aLen[2] * aLen[2] +
-                        aLen[3] * aLen[3] ) / 4. );
+                    Max( aLen[ 1 ],
+                         Max( aLen[ 2 ],
+                              Max( aLen[ 3 ],
+                                   Max( aDia[ 0 ], aDia[ 1 ] ) ) ) ) );
+    double C1 = sqrt( aLen[0] * aLen[0] +
+                      aLen[1] * aLen[1] +
+                      aLen[2] * aLen[2] +
+                      aLen[3] * aLen[3] );
     double C2 = Min( anArea[ 0 ],
-                  Min( anArea[ 1 ],
-                    Min( anArea[ 2 ], anArea[ 3 ] ) ) );
+                     Min( anArea[ 1 ],
+                          Min( anArea[ 2 ], anArea[ 3 ] ) ) );
     if ( C2 <= theEps )
       return theInf;
     return alpha * L * C1 / C2;
@@ -1091,7 +1081,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
   case 5:{
     {
       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 3 ),P( 5 )};
-      aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+      aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4]));
     }
     {
       gp_XYZ aXYZ[4] = {P( 1 ),P( 3 ),P( 4 ),P( 5 )};
@@ -1110,7 +1100,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
   case 6:{
     {
       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 6 )};
-      aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+      aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4]));
     }
     {
       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 4 ),P( 3 )};
@@ -1137,7 +1127,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
   case 8:{
     {
       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 3 )};
-      aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4])),aQuality);
+      aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[4]));
     }
     {
       gp_XYZ aXYZ[4] = {P( 1 ),P( 2 ),P( 5 ),P( 4 )};
@@ -1272,7 +1262,7 @@ double AspectRatio3D::GetValue( const TSequenceOfXYZ& P )
   case 12:
     {
       gp_XYZ aXYZ[8] = {P( 1 ),P( 2 ),P( 4 ),P( 5 ),P( 7 ),P( 8 ),P( 10 ),P( 11 )};
-      aQuality = std::max(GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8])),aQuality);
+      aQuality = GetValue(TSequenceOfXYZ(&aXYZ[0],&aXYZ[8]));
     }
     {
       gp_XYZ aXYZ[8] = {P( 2 ),P( 3 ),P( 5 ),P( 6 ),P( 8 ),P( 9 ),P( 11 ),P( 12 )};
@@ -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
@@ -2743,14 +2745,14 @@ 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();
   TItrMapOfVolume volEnd = mapOfVol.end();
   for ( ; volItr != volEnd; ++volItr )
     if ( (*volItr).second >= nbNode )
-       nbVol++;
+      nbVol++;
   // face is not free if number of volumes constructed on their nodes more than one
   return (nbVol < 2);
 }
@@ -3159,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 ));
   }
 }
@@ -4260,6 +4262,7 @@ private:
   bool isOutOfFace  (const gp_Pnt& p);
   bool isOutOfEdge  (const gp_Pnt& p);
   bool isOutOfVertex(const gp_Pnt& p);
+  bool isOutOfNone  (const gp_Pnt& p) { return true; }
   bool isBox        (const TopoDS_Shape& s);
 
   bool (Classifier::*          myIsOutFun)(const gp_Pnt& p);
@@ -4581,7 +4584,7 @@ bool ElementsOnShape::IsSatisfy (const SMDS_MeshNode* node,
         {
           isNodeOut = false;
           if ( okShape )
-            *okShape = myWorkClassifiers[i]->Shape();
+            *okShape = myClassifiers[i].Shape();
           break;
         }
     }
@@ -4619,17 +4622,27 @@ void ElementsOnShape::Classifier::Init( const TopoDS_Shape& theShape,
   {
     Standard_Real u1,u2,v1,v2;
     Handle(Geom_Surface) surf = BRep_Tool::Surface( TopoDS::Face( theShape ));
-    surf->Bounds( u1,u2,v1,v2 );
-    myProjFace.Init(surf, u1,u2, v1,v2, myTol );
-    myIsOutFun = & ElementsOnShape::Classifier::isOutOfFace;
+    if ( surf.IsNull() )
+      myIsOutFun = & ElementsOnShape::Classifier::isOutOfNone;
+    else
+    {
+      surf->Bounds( u1,u2,v1,v2 );
+      myProjFace.Init(surf, u1,u2, v1,v2, myTol );
+      myIsOutFun = & ElementsOnShape::Classifier::isOutOfFace;
+    }
     break;
   }
   case TopAbs_EDGE:
   {
     Standard_Real u1, u2;
     Handle(Geom_Curve) curve = BRep_Tool::Curve( TopoDS::Edge( theShape ), u1, u2);
-    myProjEdge.Init(curve, u1, u2);
-    myIsOutFun = & ElementsOnShape::Classifier::isOutOfEdge;
+    if ( curve.IsNull() )
+      myIsOutFun = & ElementsOnShape::Classifier::isOutOfNone;
+    else
+    {
+      myProjEdge.Init(curve, u1, u2);
+      myIsOutFun = & ElementsOnShape::Classifier::isOutOfEdge;
+    }
     break;
   }
   case TopAbs_VERTEX: