0021825: Error in the example of "Projection Algorithms" in the user's guide

[modules/smesh.git] / src / StdMeshers / StdMeshers_QuadToTriaAdaptor.cxx

diff --git a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx
index db32a5f1a5f92274ef99ed98d721e89459bd9711..5764cc6eac0520dda5a7f2a24ae226c59a627b0d 100644 (file)
--- a/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx
+++ b/src/StdMeshers/StdMeshers_QuadToTriaAdaptor.cxx
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2011  CEA/DEN, EDF R&D, OPEN CASCADE
+// Copyright (C) 2007-2012  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
@@ -21,7 +21,7 @@
 // Module    : SMESH
 // Created   : Wen May 07 16:37:07 2008
 // Author    : Sergey KUUL (skl)
-//
+
 #include "StdMeshers_QuadToTriaAdaptor.hxx"
 
 #include "SMDS_SetIterator.hxx"
@@ -40,6 +40,7 @@
 #include <gp_Pln.hxx>
 
 #include <numeric>
+#include <limits>
 
 using namespace std;
 
@@ -110,7 +111,7 @@ namespace
 
     // Check angle between normals
     double angle = nI.Angle( nJ );
-    bool tooClose = ( angle < 15 * PI180 );
+    bool tooClose = ( angle < 15. * M_PI / 180. );
 
     // Check if pyramids collide
     if ( !tooClose && baI * baJ > 0 )
@@ -378,37 +379,41 @@ StdMeshers_QuadToTriaAdaptor::~StdMeshers_QuadToTriaAdaptor()
   myElemSearcher=0;
 }
 
-
 //=======================================================================
 //function : FindBestPoint
 //purpose  : Return a point P laying on the line (PC,V) so that triangle
 //           (P, P1, P2) to be equilateral as much as possible
 //           V - normal to (P1,P2,PC)
 //=======================================================================
+
 static gp_Pnt FindBestPoint(const gp_Pnt& P1, const gp_Pnt& P2,
                             const gp_Pnt& PC, const gp_Vec& V)
 {
-  double a = P1.Distance(P2);
-  double b = P1.Distance(PC);
-  double c = P2.Distance(PC);
+  gp_Pnt Pbest = PC;
+  const double a = P1.Distance(P2);
+  const double b = P1.Distance(PC);
+  const double c = P2.Distance(PC);
   if( a < (b+c)/2 )
-    return PC;
+    return Pbest;
   else {
     // find shift along V in order a to became equal to (b+c)/2
-    double shift = sqrt( a*a + (b*b-c*c)*(b*b-c*c)/16/a/a - (b*b+c*c)/2 );
-    gp_Dir aDir(V);
-    gp_Pnt Pbest = PC.XYZ() + aDir.XYZ() * shift;
-    return Pbest;
+    const double Vsize = V.Magnitude();
+    if ( fabs( Vsize ) > std::numeric_limits<double>::min() )
+    {
+      const double shift = sqrt( a*a + (b*b-c*c)*(b*b-c*c)/16/a/a - (b*b+c*c)/2 );
+      Pbest.ChangeCoord() += shift * V.XYZ() / Vsize;
+    }
   }
+  return Pbest;
 }
 
-
 //=======================================================================
 //function : HasIntersection3
 //purpose  : Auxilare for HasIntersection()
 //           find intersection point between triangle (P1,P2,P3)
 //           and segment [PC,P]
 //=======================================================================
+
 static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
                              const gp_Pnt& P1, const gp_Pnt& P2, const gp_Pnt& P3)
 {
@@ -479,7 +484,6 @@ static bool HasIntersection3(const gp_Pnt& P, const gp_Pnt& PC, gp_Pnt& Pint,
   return false;
 }
 
-
 //=======================================================================
 //function : HasIntersection
 //purpose  : Auxilare for CheckIntersection()
@@ -701,9 +705,9 @@ int StdMeshers_QuadToTriaAdaptor::Preparation(const SMDS_MeshElement*       face
 //purpose  : 
 //=======================================================================
 
-bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&           aMesh,
-                                           const TopoDS_Shape&   aShape,
-                                           SMESH_ProxyMesh* aProxyMesh)
+bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&         aMesh,
+                                           const TopoDS_Shape& aShape,
+                                           SMESH_ProxyMesh*    aProxyMesh)
 {
   SMESH_ProxyMesh::setMesh( aMesh );
 
@@ -711,6 +715,8 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh&           aMesh,
        aShape.ShapeType() != TopAbs_SHELL )
     return false;
 
+  myShape = aShape;
+
   vector<const SMDS_MeshElement*> myPyramids;
 
   SMESHDS_Mesh * meshDS = aMesh.GetMeshDS();
@@ -1010,10 +1016,12 @@ bool StdMeshers_QuadToTriaAdaptor::Compute(SMESH_Mesh& aMesh)
     PCbest /= 4;
 
     double height = PC.Distance(PCbest); // pyramid height to precise
-    if(height<1.e-6) {
+    if ( height < 1.e-6 ) {
       // create new PCbest using a bit shift along VNorm
       PCbest = PC.XYZ() + VNorm.XYZ() * 0.001;
       height = PC.Distance(PCbest);
+      if ( height < std::numeric_limits<double>::min() )
+        return false; // batterfly element
     }
 
     // Restrict pyramid height by intersection with other faces
@@ -1214,8 +1222,8 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&
             gp_Vec VI2(PCj,Pint);
             double ang1 = fabs(VN1.Angle(VI1));
             double ang2 = fabs(VN2.Angle(VI2));
-            double coef1 = 0.5 - (( ang1<PI/3 ) ? cos(ang1)*0.25 : 0 );
-            double coef2 = 0.5 - (( ang2<PI/3 ) ? cos(ang2)*0.25 : 0 ); // cos(ang2) ?
+            double coef1 = 0.5 - (( ang1 < M_PI/3. ) ? cos(ang1)*0.25 : 0 );
+            double coef2 = 0.5 - (( ang2 < M_PI/3. ) ? cos(ang2)*0.25 : 0 ); // cos(ang2) ?
 //             double coef2 = 0.5;
 //             if(ang2<PI/3)
 //               coef2 -= cos(ang1)*0.25;
@@ -1279,4 +1287,3 @@ bool StdMeshers_QuadToTriaAdaptor::Compute2ndPart(SMESH_Mesh&
 
   return true;
 }
-