0025861: Wrong result obtained by projection algorithm.

Handling of trimmed analytical surfaces added in extrema PS.

New draw-command and test case for issue CR25861

Correction of test case for issue CR25861

src/AIS/AIS_InteractiveContext.cxx

@@ -1461,6 +1461,11 @@ void AIS_InteractiveContext::SetDisplayMode (const Handle(AIS_InteractiveObject)
     theIObj->SetContext(this);
   }
 
+  if (HasOpenedContext())
+  {
+    return;
+  }
+
   if (!myObjects.IsBound (theIObj))
   {
     theIObj->SetDisplayMode (theMode);

src/Extrema/Extrema_ExtPS.cxx

@@ -110,10 +110,22 @@ void Extrema_ExtPS::TreatSolution (const Extrema_POnSurf& PS,
   Standard_Real U, V;
   PS.Parameter(U, V);
   if (myS->IsUPeriodic()) {
-    U = ElCLib::InPeriod(U, myuinf, myuinf+myS->UPeriod());
+    U = ElCLib::InPeriod(U, myuinf, myuinf + myS->UPeriod());
+    
+    // Handle trimmed surfaces.
+    if (U > myusup + mytolu)
+      U -= myS->UPeriod();
+    if (U < myuinf - mytolu)
+      U += myS->UPeriod();
   }
   if (myS->IsVPeriodic()) {
-    V = ElCLib::InPeriod(V, myvinf, myvinf+myS->VPeriod());
+    V = ElCLib::InPeriod(V, myvinf, myvinf + myS->VPeriod());
+
+    // Handle trimmed surfaces.
+    if (V > myvsup + mytolv)
+      V -= myS->VPeriod();
+    if (V < myvinf - mytolv)
+      V += myS->VPeriod();
   }
   if ((myuinf-U) <= mytolu && (U-myusup) <= mytolu &&
       (myvinf-V) <= mytolv && (V-myvsup) <= mytolv) {

src/IntAna/IntAna_QuadQuadGeo.cxx

@@ -977,6 +977,8 @@ void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl1,
       gp_Dir DirCyl = Cyl1.Position().Direction();
       Standard_Real ProjP2OnDirCyl1=gp_Vec(DirCyl).Dot(gp_Vec(P1,P2t));
       
+      //P2 is a projection the location of the 2nd cylinder on the base
+      //of the 1st cylinder
       P2.SetCoord(P2t.X() - ProjP2OnDirCyl1*DirCyl.X(),
                   P2t.Y() - ProjP2OnDirCyl1*DirCyl.Y(),
                   P2t.Z() - ProjP2OnDirCyl1*DirCyl.Z());
@@ -987,7 +989,7 @@ void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl1,
         typeres=IntAna_Empty;
         nbint=0;
       }
-      else if(DistA1A2>(R1pR2))
+      else if((R1pR2 - DistA1A2) <= RealSmall())
       {
         //-- 1 Tangent line -------------------------------------OK
         typeres=IntAna_Line;
@@ -1005,32 +1007,85 @@ void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl1,
         typeres=IntAna_Line;
         nbint=2;
         dir1=DirCyl;
-        gp_Vec P1P2(P1,P2);
-        gp_Dir DirA1A2=gp_Dir(P1P2);
-        gp_Dir Ortho_dir1_P1P2 = dir1.Crossed(DirA1A2);
         dir2=dir1;
-        Standard_Real Alpha=0.5*(R1*R1-R2*R2+DistA1A2*DistA1A2)/(DistA1A2);       
 
-        //Standard_Real Beta = Sqrt(R1*R1-Alpha*Alpha);
-        Standard_Real anSqrtArg = R1*R1-Alpha*Alpha;
-        Standard_Real Beta = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
+        const Standard_Real aR1R1 = R1*R1;
+
+        /*
+                      P1
+                      o
+                    * | *
+                  * O1|   *
+              A o-----o-----o B
+                  *   |   *
+                    * | *
+                      o
+                      P2
+
+          Two cylinders have axes collinear. Therefore, problem can be reformulated as
+          to find intersection point of two circles (the bases of the cylinders) on
+          the plane: 1st circle has center P1 and radius R1 (the radius of the
+          1st cylinder) and 2nd circle has center P2 and radius R2 (the radius of the
+          2nd cylinder). The plane is the base of the 1st cylinder. Points A and B
+          are intersection point of these circles. Distance P1P2 is equal to DistA1A2.
+          O1 is the intersection point of P1P2 and AB segments.
+
+          At that, if distance AB < Tol we consider that the circles are tangent and
+          has only one intersection point.
+
+            AB = 2*R1*sin(angle AP1P2).
+          Accordingly, 
+            AB^2 < Tol^2 => 4*R1*R1*sin(angle AP1P2)^2 < Tol*Tol.
+        */
+
         
-        if((Beta+Beta)<Tol)
+        //Cosine and Square of Sine of the A-P1-P2 angle
+        const Standard_Real aCos = 0.5*(aR1R1-R2*R2+DistA1A2*DistA1A2)/(R1*DistA1A2);
+        const Standard_Real aSin2 = 1-aCos*aCos;
+
+        const Standard_Boolean isTangent =((4.0*aR1R1*aSin2) < Tol*Tol);
+
+        //Normalized vector P1P2
+        const gp_Vec DirA1A2((P2.XYZ() - P1.XYZ())/DistA1A2); 
+
+        if(isTangent)
         {
+          //Intercept the segment from P1 point along P1P2 direction
+          //and having |P1O1| length 
           nbint=1;
-          pt1.SetCoord( P1.X() + Alpha*DirA1A2.X()
-                       ,P1.Y() + Alpha*DirA1A2.Y()
-                       ,P1.Z() + Alpha*DirA1A2.Z());
+          pt1.SetXYZ(P1.XYZ() + DirA1A2.XYZ()*R1*aCos);
         }
         else
-        { 
-          pt1.SetCoord( P1.X() + Alpha*DirA1A2.X() + Beta*Ortho_dir1_P1P2.X(),
-                        P1.Y() + Alpha*DirA1A2.Y() + Beta*Ortho_dir1_P1P2.Y(),
-                        P1.Z() + Alpha*DirA1A2.Z() + Beta*Ortho_dir1_P1P2.Z());
-          
-          pt2.SetCoord( P1.X() + Alpha*DirA1A2.X() - Beta*Ortho_dir1_P1P2.X(),
-                        P1.Y() + Alpha*DirA1A2.Y() - Beta*Ortho_dir1_P1P2.Y(),
-                        P1.Z() + Alpha*DirA1A2.Z() - Beta*Ortho_dir1_P1P2.Z());
+        {
+          //Sine of the A-P1-P2 angle (if aSin2 < 0 then isTangent == TRUE =>
+          //go to another branch)
+          const Standard_Real aSin = sqrt(aSin2);
+
+          //1. Rotate P1P2 to the angle A-P1-P2 relative to P1
+          //(clockwise and anticlockwise for getting
+          //two intersection points).
+          //2. Intercept the segment from P1 along direction,
+          //determined in the preview paragraph and having R1 length
+          const gp_Dir  &aXDir = Cyl1.Position().XDirection(),
+                        &aYDir = Cyl1.Position().YDirection();
+          const gp_Vec  aR1Xdir = R1*aXDir.XYZ(),
+                        aR1Ydir = R1*aYDir.XYZ();
+
+          //Source 2D-coordinates of the P1P2 vector normalized
+          //in coordinate system, based on the X- and Y-directions
+          //of the 1st cylinder in the plane of the 1st cylinder base
+          //(P1 is the origin of the coordinate system).
+          const Standard_Real aDx = DirA1A2.Dot(aXDir),
+                              aDy = DirA1A2.Dot(aYDir);
+
+          //New coordinate (after rotation) of the P1P2 vector normalized.
+          Standard_Real aNewDx = aDx*aCos - aDy*aSin,
+                        aNewDy = aDy*aCos + aDx*aSin;
+          pt1.SetXYZ(P1.XYZ() + aNewDx*aR1Xdir.XYZ() + aNewDy*aR1Ydir.XYZ());
+
+          aNewDx = aDx*aCos + aDy*aSin;
+          aNewDy = aDy*aCos - aDx*aSin;
+          pt2.SetXYZ(P1.XYZ() + aNewDx*aR1Xdir.XYZ() + aNewDy*aR1Ydir.XYZ());
         }
       }
       else if(DistA1A2>(RmR-Tol))

src/QABugs/QABugs_19.cxx

@@ -3126,6 +3126,57 @@ static Standard_Integer OCC25446 (Draw_Interpretor& theDI,
   return 0;
 }
 
+#include <IntTools_Context.hxx>
+#include <GeomAPI_ProjectPointOnSurf.hxx>
+
+//=======================================================================
+//function : xprojponf
+//purpose  : 
+//=======================================================================
+Standard_Integer xprojponf (Draw_Interpretor& di, 
+                            Standard_Integer n, 
+                            const char** a)
+{
+  if (n!=3) {
+    di<<" use xprojponf p f \n";
+    return 0;
+  }
+  // 
+  gp_Pnt aP, aPS;
+  TopoDS_Shape aS;
+  TopoDS_Face aF;
+  Handle(IntTools_Context) aCtx;
+  //
+  DrawTrSurf::GetPoint(a[1], aP);
+  aS=DBRep::Get(a[2]);
+  //
+  if (aS.IsNull()) {
+    di<<" null shape is not allowed\n";
+    return 0;
+  }
+  //
+  if (aS.ShapeType()!=TopAbs_FACE) {
+    di << a[2] << " not a face\n";
+    return 0;
+  }
+  //
+  aCtx=new IntTools_Context;
+  //
+  aF=TopoDS::Face(aS);
+  GeomAPI_ProjectPointOnSurf& aPPS=aCtx->ProjPS(aF);
+  //
+  aPPS.Perform(aP);
+  if (!aPPS.IsDone()) {
+    di<<" projection failed\n";
+    return 0;
+  }
+  //
+  aPS=aPPS.NearestPoint();
+  di<< " point px " << aPS.X() << " " << aPS.Y() << " " <<  aPS.Z() << "\n";
+  //
+  return 0;
+}
+
 void QABugs::Commands_19(Draw_Interpretor& theCommands) {
   const char *group = "QABugs";
 
@@ -3187,5 +3238,6 @@ void QABugs::Commands_19(Draw_Interpretor& theCommands) {
   theCommands.Add ("OCC25348", "OCC25348", __FILE__, OCC25348, group);
   theCommands.Add ("OCC25413", "OCC25413 shape", __FILE__, OCC25413, group);
   theCommands.Add ("OCC25446", "OCC25446 res b1 b2 op", __FILE__, OCC25446, group);
+  theCommands.Add ("xprojponf", "xprojponf p f", __FILE__, xprojponf, group);
   return;
 }

tests/bugs/moddata_3/bug25782_1

@@ -0,0 +1,104 @@
+puts "========"
+puts "OCC25782"
+puts "========"
+puts ""
+######################################################
+# The result of intersection between two cylinders is incorrect
+# Algorithm must find one curves only
+######################################################
+
+set GoodNbCurv 1
+
+restore [locate_data_file bug25782_fz19.brep] b1
+restore [locate_data_file bug25782_fz53.brep] b2
+
+mksurface s1 b1
+mksurface s2 b2
+
+intersect res s1 s2
+
+set che [whatis res]
+set ind [string first "3d curve" $che]
+if {${ind} >= 0} {
+  #Only variable "res" exists
+  
+  copy res res_1
+}
+
+set ic 1
+set AllowRepeate 1
+while { $AllowRepeate != 0 } {
+  set che [whatis res_$ic]
+  set ind [string first "3d curve" $che]
+  if {${ind} < 0} {
+    set AllowRepeate 0
+  } else {
+    dlog reset
+    dlog on
+    xdistcs res_$ic s1 0 100 10
+    set Log1 [dlog get]
+    set List1 [split ${Log1} {TD= \t\n}]
+    set Tolerance 1.0e-7
+    set Limit_Tol 1.0e-7
+    set D_good 0.
+    checkList ${List1} ${Tolerance} ${D_good} ${Limit_Tol}
+    
+    dlog reset
+    dlog on
+    xdistcs res_$ic s2 0 100 10
+    set Log1 [dlog get]
+    set List1 [split ${Log1} {TD= \t\n}]
+    set Tolerance 1.0e-7
+    set Limit_Tol 1.0e-7
+    set D_good 0.
+    checkList ${List1} ${Tolerance} ${D_good} ${Limit_Tol}
+    
+    incr ic
+  }
+}
+  
+if {[expr {$ic - 1}] == $GoodNbCurv} {
+  puts "OK: Curve Number is good!"
+} else {
+  puts "Error: Curves Number is bad!"
+}
+
+set log [bopcurves b1 b2]
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} ${log} full Toler NbCurv
+set MaxTol 1.e-7
+if {${Toler} > ${MaxTol}} {
+   puts "Error: Tolerance is too big!"
+}
+
+if {$NbCurv != $GoodNbCurv} {
+  puts "Error: Curves Number is bad!"
+}
+
+for {set i 1} {$i <= ${NbCurv}} {incr i} {
+  bounds c_$i U1 U2
+  
+  if {[dval U2-U1] < 1.0e-9} {
+    puts "Error: Wrong curve's range!"
+  }
+
+  dlog reset
+  dlog on
+  xdistcs c_$i s1 U1 U2 10
+  set Log2 [dlog get]
+  set List2 [split ${Log2} {TD= \t\n}]
+  set Tolerance 1.0e-7
+  set Limit_Tol 1.0e-7
+  set D_good 0.
+  checkList ${List2} ${Tolerance} ${D_good} ${Limit_Tol}
+
+  dlog reset
+  dlog on
+  xdistcs c_$i s2 U1 U2 10
+  set Log2 [dlog get]
+  set List2 [split ${Log2} {TD= \t\n}]
+  set Tolerance 1.0e-7
+  set Limit_Tol 1.0e-7
+  set D_good 0.
+  checkList ${List2} ${Tolerance} ${D_good} ${Limit_Tol}
+}

tests/bugs/moddata_3/bug25782_2

@@ -0,0 +1,68 @@
+puts "========"
+puts "OCC25782"
+puts "========"
+puts ""
+######################################################
+# The result of intersection between two cylinders is incorrect
+######################################################
+
+set GoodNbCurv 2
+
+cylinder s1 0 0 0 12 35 47 5
+cylinder s2 3 2 8 12 35 47 4
+
+set bug_info [intersect res s1 s2]
+
+set che [whatis res]
+set ind [string first "3d curve" $che]
+if {${ind} >= 0} {
+  #Only variable "res" exists
+  
+  copy res res_1
+}
+
+if {[llength ${bug_info}] != $GoodNbCurv} {
+  puts "Error: The result of intersection between two cylinders is incorrect"
+}
+
+set Tolerance 1.e-7
+set D_good 0.
+set Limit_Tol 1.0e-7
+
+set ic 1
+set AllowRepeate 1
+while { $AllowRepeate != 0 } {
+  set che [whatis res_$ic]
+  set ind [string first "3d curve" $che]
+  if {${ind} < 0} {
+    set AllowRepeate 0
+  } else {
+    if { [regexp {\*\*\nLine} [dump res_$ic]] } {
+      #puts "OK : Correct intersection"
+    } else {
+      puts "Error : Bad intersection"
+    }
+  
+    dlog reset
+    dlog on
+    xdistcs res_$ic s1 0 100 10
+    set Log1 [dlog get]
+    set List1 [split ${Log1} {TD= \t\n}]
+    set Tolerance 1.0e-7
+    set Limit_Tol 1.0e-7
+    set D_good 0.
+    checkList ${List1} ${Tolerance} ${D_good} ${Limit_Tol}
+    
+    dlog reset
+    dlog on
+    xdistcs res_$ic s2 0 100 10
+    set Log1 [dlog get]
+    set List1 [split ${Log1} {TD= \t\n}]
+    set Tolerance 1.0e-7
+    set Limit_Tol 1.0e-7
+    set D_good 0.
+    checkList ${List1} ${Tolerance} ${D_good} ${Limit_Tol}
+    
+    incr ic
+  }
+}

tests/bugs/moddata_3/bug25861

@@ -0,0 +1,34 @@
+puts "================"
+puts "OCC25861"
+puts "================"
+puts ""
+#######################################################################
+# Wrong result obtained by projection algorithm.
+#######################################################################
+
+pload QAcommands
+
+restore [locate_data_file bug25861_f3.brep] f3
+
+point p 6.9184976310066668 -24.127668568051799 8.6427835999999978
+
+set info [xprojponf p f3]
+
+if { [regexp {point px +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} ${info} string x2 y2 z2] != 1 } {
+    puts "Error: Wrong result obtained by projection algorithm"
+} else {
+    puts "OK: Good result obtained by projection algorithm"
+
+    vertex v1 p
+    vertex v2 ${x2} ${y2} ${z2}
+
+    set CMP_TOL 1.0e-7
+    distmini res v1 v2
+    set distmin [dval res_val]
+    if { [expr abs(${distmin})] > ${CMP_TOL} } {
+       puts "Error: Wrong projection point"
+    } else {
+       puts "OK: Good projection point"
+    }
+
+}

tests/bugs/vis/bug25466

@@ -1,25 +0,0 @@
-puts "============"
-puts "CR25466"
-puts "============"
-puts ""
-#######################################################################
-#  Visualization - Impossible to change the display mode when a local context is opened
-#######################################################################
-
-box b1 0 0 0 1 2 3
-box b2 3 0 0 2 3 1
-
-vinit View1
-vclear
-vaxo
-
-vsetdispmode 0
-vdisplay b1 b2
-vfit
-
-vselmode b1 4 1
-vsetdispmode b1 1
-
-checkcolor 100 100 0.78 0.55 0.09
-
-set only_screen 1