Patch to OCCT 6.7.1 version

adm/templates/env.sh.in

@@ -60,12 +60,12 @@ fi
 BIN_PATH="${WOKSTATION}${ARCH}/${COMPILER}/bin${CASDEB}"
 LIBS_PATH="${WOKSTATION}${ARCH}/${COMPILER}/lib${CASDEB}"
 
-export PATH="${CASROOT}/${BIN_PATH}:${PATH}"
+export PATH="${CASROOT}/${BIN_PATH}:${THRDPARTY_PATH}:${PATH}"
 
 if [ "$LD_LIBRARY_PATH" != "" ]; then
-  export LD_LIBRARY_PATH="${CASROOT}/${LIBS_PATH}:${THRDPARTY_PATH}:${LD_LIBRARY_PATH}"
+  export LD_LIBRARY_PATH="${CASROOT}/${LIBS_PATH}:${LD_LIBRARY_PATH}"
 else
-  export LD_LIBRARY_PATH="${CASROOT}/${LIBS_PATH}:${THRDPARTY_PATH}"
+  export LD_LIBRARY_PATH="${CASROOT}/${LIBS_PATH}"
 fi
 
 if [ "$WOKSTATION" == "mac" ]; then

samples/mfc/standard/01_Geometry/src/GeometryDoc.cpp

@@ -1118,7 +1118,8 @@ static Standard_Boolean fixParam(Standard_Real& theParam)
 void CGeometryDoc::OnSimplify() 
 {
   CString initfile(((OCC_App*) AfxGetApp())->GetInitDataDir());
-  initfile += "\\..\\..\\Data\\";
+  //initfile += "\\..\\..\\Data\\";
+  initfile += "\\..\\Data\\";
   initfile += "shell1.brep";
 
   TCollection_AsciiString Path((Standard_CString)(LPCTSTR)initfile);

samples/qt/IESample/IESample-vc9.sln

@@ -1,6 +1,6 @@
 Microsoft Visual Studio Solution File, Format Version 10.00
 # Visual Studio 2008
-Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "IESample", "IESample.vcproj", "{AD4DADCB-F15E-37FD-B3AA-88504FAAF4FD}"
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "IESample", "IESample.vcproj", "{29A31181-FF46-3DB1-827A-05ADEAD293B5}"
 EndProject
 Global
 	GlobalSection(SolutionConfigurationPlatforms) = preSolution
@@ -8,10 +8,10 @@ Global
 		Release|Win32 = Release|Win32
 	EndGlobalSection
 	GlobalSection(ProjectConfigurationPlatforms) = postSolution
-		{AD4DADCB-F15E-37FD-B3AA-88504FAAF4FD}.Debug|Win32.ActiveCfg = Debug|Win32
-		{AD4DADCB-F15E-37FD-B3AA-88504FAAF4FD}.Debug|Win32.Build.0 = Debug|Win32
-		{AD4DADCB-F15E-37FD-B3AA-88504FAAF4FD}.Release|Win32.ActiveCfg = Release|Win32
-		{AD4DADCB-F15E-37FD-B3AA-88504FAAF4FD}.Release|Win32.Build.0 = Release|Win32
+		{29A31181-FF46-3DB1-827A-05ADEAD293B5}.Debug|Win32.ActiveCfg = Debug|Win32
+		{29A31181-FF46-3DB1-827A-05ADEAD293B5}.Debug|Win32.Build.0 = Debug|Win32
+		{29A31181-FF46-3DB1-827A-05ADEAD293B5}.Release|Win32.ActiveCfg = Release|Win32
+		{29A31181-FF46-3DB1-827A-05ADEAD293B5}.Release|Win32.Build.0 = Release|Win32
 	EndGlobalSection
 	GlobalSection(SolutionProperties) = preSolution
 		HideSolutionNode = FALSE

samples/qt/Tutorial/Tutorial-vc9.sln

@@ -1,6 +1,6 @@
 Microsoft Visual Studio Solution File, Format Version 10.00
 # Visual Studio 2008
-Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "Tutorial", "Tutorial.vcproj", "{28627D0B-F82A-39D5-A15D-DDAFA11059E7}"
+Project("{8BC9CEB8-8B4A-11D0-8D11-00A0C91BC942}") = "Tutorial", "Tutorial.vcproj", "{C9C6457A-5B88-3C7F-9964-9D8EE43E7414}"
 EndProject
 Global
 	GlobalSection(SolutionConfigurationPlatforms) = preSolution
@@ -8,10 +8,10 @@ Global
 		Release|Win32 = Release|Win32
 	EndGlobalSection
 	GlobalSection(ProjectConfigurationPlatforms) = postSolution
-		{28627D0B-F82A-39D5-A15D-DDAFA11059E7}.Debug|Win32.ActiveCfg = Debug|Win32
-		{28627D0B-F82A-39D5-A15D-DDAFA11059E7}.Debug|Win32.Build.0 = Debug|Win32
-		{28627D0B-F82A-39D5-A15D-DDAFA11059E7}.Release|Win32.ActiveCfg = Release|Win32
-		{28627D0B-F82A-39D5-A15D-DDAFA11059E7}.Release|Win32.Build.0 = Release|Win32
+		{C9C6457A-5B88-3C7F-9964-9D8EE43E7414}.Debug|Win32.ActiveCfg = Debug|Win32
+		{C9C6457A-5B88-3C7F-9964-9D8EE43E7414}.Debug|Win32.Build.0 = Debug|Win32
+		{C9C6457A-5B88-3C7F-9964-9D8EE43E7414}.Release|Win32.ActiveCfg = Release|Win32
+		{C9C6457A-5B88-3C7F-9964-9D8EE43E7414}.Release|Win32.Build.0 = Release|Win32
 	EndGlobalSection
 	GlobalSection(SolutionProperties) = preSolution
 		HideSolutionNode = FALSE

src/BOPAlgo/BOPAlgo.cdl

@@ -52,6 +52,7 @@ is
       IncompatibilityOfFace, 
       OperationAborted,
       GeomAbs_C0,
+      InvalidCurveOnSurface,
       NotValid
     end CheckStatus;
 

src/BOPAlgo/BOPAlgo_ArgumentAnalyzer.cdl

@@ -115,6 +115,13 @@ is
     ---Purpose: Returns (modifiable) mode that means
     --          checking of problem of continuity of the shape.
 
+    CurveOnSurfaceMode(me: in out) 
+    returns Boolean from Standard;
+    ---C++: return &
+    ---C++: inline
+    ---Purpose: Returns (modifiable) mode that means
+    --          checking of problem of invalid curve on surface.
+
     ---
     Perform(me: out);
     ---Purpose: performs analysis
@@ -154,6 +161,9 @@ is
     is protected;
 
     TestContinuity(me: out)
+    is protected; 
+     
+    TestCurveOnSurface(me: out) 
     is protected;
 
 --  TestMergeFace(me: out)
@@ -162,20 +172,20 @@ is
 
 fields
 
-    myShape1           : Shape     from TopoDS;
-    myShape2           : Shape     from TopoDS;
-    myStopOnFirst      : Boolean   from Standard;
-    myOperation        : Operation from BOPAlgo;
-    myArgumentTypeMode : Boolean   from Standard;
-    mySelfInterMode    : Boolean   from Standard;
-    mySmallEdgeMode    : Boolean   from Standard;
-    myRebuildFaceMode  : Boolean   from Standard;
-    myTangentMode      : Boolean   from Standard;
-    myMergeVertexMode  : Boolean   from Standard;
-    myMergeEdgeMode    : Boolean   from Standard;
-    myContinuityMode   : Boolean   from Standard;
-    myEmpty1,myEmpty2  : Boolean   from Standard; 
-    myResult      : ListOfCheckResult from BOPAlgo; 
-    
-    
+    myShape1             : Shape     from TopoDS;
+    myShape2             : Shape     from TopoDS;
+    myStopOnFirst        : Boolean   from Standard;
+    myOperation          : Operation from BOPAlgo;
+    myArgumentTypeMode   : Boolean   from Standard;
+    mySelfInterMode      : Boolean   from Standard;
+    mySmallEdgeMode      : Boolean   from Standard;
+    myRebuildFaceMode    : Boolean   from Standard;
+    myTangentMode        : Boolean   from Standard;
+    myMergeVertexMode    : Boolean   from Standard;
+    myMergeEdgeMode      : Boolean   from Standard;
+    myContinuityMode     : Boolean   from Standard; 
+    myCurveOnSurfaceMode : Boolean   from Standard; 
+    myEmpty1, myEmpty2   : Boolean   from Standard; 
+    myResult     : ListOfCheckResult from BOPAlgo; 
+        
 end ArgumentAnalyzer;

src/BOPAlgo/BOPAlgo_ArgumentAnalyzer.cxx

@@ -73,6 +73,7 @@ myTangentMode(Standard_False),
 myMergeVertexMode(Standard_False),
 myMergeEdgeMode(Standard_False),
 myContinuityMode(Standard_False),
+myCurveOnSurfaceMode(Standard_False),
 myEmpty1(Standard_False),
 myEmpty2(Standard_False)
 // myMergeFaceMode(Standard_False)
@@ -196,6 +197,11 @@ void BOPAlgo_ArgumentAnalyzer::Perform()
       if(!(!myResult.IsEmpty() && myStopOnFirst))
         TestContinuity();
     }
+
+    if(myCurveOnSurfaceMode) {
+      if(!(!myResult.IsEmpty() && myStopOnFirst))
+        TestCurveOnSurface();
+    }
   }
   catch(Standard_Failure) {
     BOPAlgo_CheckResult aResult;
@@ -904,6 +910,57 @@ void BOPAlgo_ArgumentAnalyzer::TestContinuity()
   }
 }
 
+// ================================================================================
+// function: TestCurveOnSurface
+// purpose:
+// ================================================================================
+void BOPAlgo_ArgumentAnalyzer::TestCurveOnSurface()
+{
+  Standard_Integer i;
+  Standard_Real aT, aD, aTolE;
+  TopExp_Explorer aExpF, aExpE;
+  //
+  for(i = 0; i < 2; i++) {
+    const TopoDS_Shape& aS = (i == 0) ? myShape1 : myShape2;
+    if(aS.IsNull()) {
+      continue;
+    }
+    //
+    aExpF.Init(aS, TopAbs_FACE);
+    for (; aExpF.More(); aExpF.Next()) {
+      const TopoDS_Face& aF = *(TopoDS_Face*)&aExpF.Current();
+      //
+      aExpE.Init(aF, TopAbs_EDGE);
+      for (; aExpE.More(); aExpE.Next()) {
+        const TopoDS_Edge& aE = *(TopoDS_Edge*)&aExpE.Current();
+        //
+        if (BOPTools_AlgoTools::ComputeTolerance(aF, aE, aD, aT)) {
+          aTolE = BRep_Tool::Tolerance(aE);
+          if (aD > aTolE) {
+            BOPAlgo_CheckResult aResult;
+            aResult.SetCheckStatus(BOPAlgo_InvalidCurveOnSurface);
+            if(i == 0) {
+              aResult.SetShape1(myShape1);
+              aResult.AddFaultyShape1(aE);
+              aResult.AddFaultyShape1(aF);
+              aResult.SetMaxDistance1(aD);
+              aResult.SetMaxParameter1(aT);
+            }
+            else {
+              aResult.SetShape2(myShape2);
+              aResult.AddFaultyShape2(aE);
+              aResult.AddFaultyShape2(aF);
+              aResult.SetMaxDistance2(aD);
+              aResult.SetMaxParameter2(aT);
+            }
+            myResult.Append(aResult);
+          }
+        }
+      }
+    }
+  }
+}
+
 // ================================================================================
 // function: TestMergeFace
 // purpose:

src/BOPAlgo/BOPAlgo_ArgumentAnalyzer.lxx

@@ -57,6 +57,10 @@ inline Standard_Boolean& BOPAlgo_ArgumentAnalyzer::ContinuityMode()
   return myContinuityMode;
 }
 
+inline Standard_Boolean& BOPAlgo_ArgumentAnalyzer::CurveOnSurfaceMode() 
+{
+  return myCurveOnSurfaceMode;
+}
 // inline Standard_Boolean& BOPAlgo_ArgumentAnalyzer::MergeFaceMode() 
 // {
 //   return myMergeFaceMode;

src/BOPAlgo/BOPAlgo_CheckResult.cdl

@@ -39,7 +39,7 @@ is
     ---Purpose: sets ancestor shape (tool) for faulty sub-shapes
     
     AddFaultyShape2(me: in out; TheShape: Shape from TopoDS);
-    ---Purpose: adds faulty sub-shapes from tool to a list
+    ---Purpose: adds faulty sub-shapes from tool to a list 
     
     GetShape1(me)
     returns Shape from TopoDS;
@@ -66,7 +66,39 @@ is
     
     GetCheckStatus(me)
     returns CheckStatus from BOPAlgo;
-    ---Purpose: gets status of faulty
+    ---Purpose: gets status of faulty 
+     
+    SetMaxDistance1(me:out; 
+        theDist : Real from Standard); 
+    ---Purpose: Sets max distance for the first shape
+ 
+    SetMaxDistance2(me:out; 
+        theDist : Real from Standard); 
+    ---Purpose: Sets max distance for the second shape
+ 
+    SetMaxParameter1(me:out; 
+        thePar : Real from Standard); 
+    ---Purpose: Sets the parameter for the first shape
+ 
+    SetMaxParameter2(me:out; 
+        thePar : Real from Standard); 
+    ---Purpose: Sets the parameter for the second shape 
+     
+    GetMaxDistance1(me) 
+    returns Real from Standard; 
+    ---Purpose: Returns the distance for the first shape
+ 
+    GetMaxDistance2(me) 
+    returns Real from Standard; 
+    ---Purpose: Returns the distance for the second shape
+ 
+    GetMaxParameter1(me) 
+    returns Real from Standard; 
+    ---Purpose: Returns the parameter for the fircst shape
+ 
+    GetMaxParameter2(me) 
+    returns Real from Standard; 
+    ---Purpose: Returns the parameter for the second shape
   
 fields
 
@@ -74,6 +106,10 @@ fields
     myShape2 : Shape from TopoDS;
     myStatus : CheckStatus from BOPAlgo;
     myFaulty1 : ListOfShape from BOPCol;
-    myFaulty2 : ListOfShape from BOPCol;
+    myFaulty2 : ListOfShape from BOPCol; 
+    myMaxDist1 : Real from Standard;
+    myMaxDist2 : Real from Standard;
+    myMaxPar1  : Real from Standard;
+    myMaxPar2  : Real from Standard;
 
 end CheckResult;

src/BOPAlgo/BOPAlgo_CheckResult.cxx

@@ -19,7 +19,13 @@
 // function:  BOPAlgo_CheckResult()
 // purpose: 
 //=======================================================================
-BOPAlgo_CheckResult::BOPAlgo_CheckResult() : myStatus(BOPAlgo_CheckUnknown)
+BOPAlgo_CheckResult::BOPAlgo_CheckResult() 
+: 
+  myStatus(BOPAlgo_CheckUnknown),
+  myMaxDist1(0.),
+  myMaxDist2(0.),
+  myMaxPar1(0.),
+  myMaxPar2(0.)
 {
 }
 
@@ -72,3 +78,43 @@ BOPAlgo_CheckStatus BOPAlgo_CheckResult::GetCheckStatus() const
 {
   return myStatus;
 }
+
+void BOPAlgo_CheckResult::SetMaxDistance1(const Standard_Real theDist)
+{
+  myMaxDist1 = theDist;
+}
+
+void BOPAlgo_CheckResult::SetMaxDistance2(const Standard_Real theDist)
+{
+  myMaxDist2 = theDist;
+}
+
+void BOPAlgo_CheckResult::SetMaxParameter1(const Standard_Real thePar)
+{
+  myMaxPar1 = thePar;
+}
+
+void BOPAlgo_CheckResult::SetMaxParameter2(const Standard_Real thePar)
+{
+  myMaxPar2 = thePar;
+}
+
+Standard_Real BOPAlgo_CheckResult::GetMaxDistance1() const
+{
+  return myMaxDist1;
+}
+
+Standard_Real BOPAlgo_CheckResult::GetMaxDistance2() const
+{
+  return myMaxDist2;
+}
+
+Standard_Real BOPAlgo_CheckResult::GetMaxParameter1() const
+{
+  return myMaxPar1;
+}
+
+Standard_Real BOPAlgo_CheckResult::GetMaxParameter2() const
+{
+  return myMaxPar2;
+}

src/BOPTest/BOPTest_CheckCommands.cxx

@@ -15,6 +15,11 @@
 #include <BOPTest.ixx>
 #include <TCollection_AsciiString.hxx>
 
+#include <TopExp_Explorer.hxx>
+
+#include <TopTools_MapOfShape.hxx>
+#include <TopTools_ShapeMapHasher.hxx>
+
 #include <gp_Pnt.hxx>
 
 #include <TopoDS_Shape.hxx>
@@ -50,11 +55,26 @@
 #include <BOPAlgo_ArgumentAnalyzer.hxx>
 #include <BOPAlgo_CheckResult.hxx>
 
+#include <BOPTools_AlgoTools.hxx>
+
 static 
   Standard_Integer bopcheck (Draw_Interpretor&, Standard_Integer, const char** );
 
 static 
   Standard_Integer bopargcheck (Draw_Interpretor&, Standard_Integer, const char** );
+
+static 
+  Standard_Integer xdistef(Draw_Interpretor&, Standard_Integer, const char** );
+
+static 
+  Standard_Integer checkcurveonsurf (Draw_Interpretor&, Standard_Integer, const char**);
+
+static 
+  Standard_Integer csdeviation(Draw_Interpretor&, Standard_Integer, const char** );
+
+static 
+  Standard_Integer checkcsdeviation(Draw_Interpretor&, Standard_Integer, const char** );
+
 //
 
 //=======================================================================
@@ -77,6 +97,20 @@ static
   theCommands.Add("bopargcheck" , 
                   "Use bopargcheck without parameters to get ",  
                   __FILE__, bopargcheck, g);
+  theCommands.Add ("xdistef" ,
+                   "Use xdistef edge face",
+                   __FILE__, xdistef, g);
+  theCommands.Add("checkcurveonsurf",
+                  "checkcurveonsurf shape",
+                  __FILE__, checkcurveonsurf, g);
+
+  theCommands.Add("csdeviation",
+                  "csdeviation edge face",
+                  __FILE__, csdeviation, g);
+
+  theCommands.Add("checkcsdeviation",
+                  "checkcsdeviation shape",
+                  __FILE__, checkcsdeviation, g);
 }
 
 //=======================================================================
@@ -293,7 +327,10 @@ static void MakeShapeForFullOutput
    const Standard_Integer          aIndex,
    const BOPCol_ListOfShape &    aList,
    Standard_Integer&               aCount,
-   Draw_Interpretor&               di)
+   Draw_Interpretor&               di,
+   Standard_Boolean                bCurveOnSurf = Standard_False,
+   Standard_Real                   aMaxDist = 0.,
+   Standard_Real                   aMaxParameter = 0.)
 {
   TCollection_AsciiString aNum(aIndex);
   TCollection_AsciiString aName = aBaseName + aNum;
@@ -309,7 +346,15 @@ static void MakeShapeForFullOutput
     BB.Add(cmp, aS);
     aCount++;
   }
-  di << "Made faulty shape: " << name << "\n";
+  di << "Made faulty shape: " << name;
+  //
+  if (bCurveOnSurf) {
+    di << " (MaxDist = " << aMaxDist 
+       << ", MaxPar = " << aMaxParameter << ")";
+  }
+  //
+  di << "\n";
+  //
   DBRep::Set(name, cmp);
 }
 
@@ -324,7 +369,7 @@ Standard_Integer bopargcheck
   if (n<2) {
     di << "\n";
     di << " Use >bopargcheck Shape1 [[Shape2] ";
-    di << "[-F/O/C/T/S/U] [/R|F|T|V|E|I|P]] [#BF]" << "\n" << "\n";
+    di << "[-F/O/C/T/S/U] [/R|F|T|V|E|I|P|C|S]] [#BF]" << "\n" << "\n";
     di << " -<Boolean Operation>" << "\n";
     di << " F (fuse)" << "\n";
     di << " O (common)" << "\n";
@@ -344,6 +389,7 @@ Standard_Integer bopargcheck
     di << " I (disable self-interference test)" << "\n";
     di << " P (disable shape type test)" << "\n";
     di << " C (disable test for shape continuity)" << "\n";
+    di << " S (disable curve on surface check)" << "\n";
     di << " For example: \"bopargcheck s1 s2 /RI\" disables ";
     di << "small edge detection and self-intersection detection" << "\n";
     di << " default - all options are enabled" << "\n" << "\n";
@@ -427,11 +473,12 @@ Standard_Integer bopargcheck
   aChecker.SetShape1(aS1);
 
   // set default options (always tested!) for single and couple shapes
-  aChecker.ArgumentTypeMode() = Standard_True;
-  aChecker.SelfInterMode()    = Standard_True;
-  aChecker.SmallEdgeMode()    = Standard_True;
-  aChecker.RebuildFaceMode()  = Standard_True;
-  aChecker.ContinuityMode()   = Standard_True;
+  aChecker.ArgumentTypeMode()   = Standard_True;
+  aChecker.SelfInterMode()      = Standard_True;
+  aChecker.SmallEdgeMode()      = Standard_True;
+  aChecker.RebuildFaceMode()    = Standard_True;
+  aChecker.ContinuityMode()     = Standard_True;
+  aChecker.CurveOnSurfaceMode() = Standard_True;
 
   // test & set options and operation for two shapes
   if(!aS22.IsNull()) {
@@ -499,6 +546,9 @@ Standard_Integer bopargcheck
       else if(a[indxOP][ind] == 'C' || a[indxOP][ind] == 'c') {
         aChecker.ContinuityMode() = Standard_False;
       }
+      else if(a[indxOP][ind] == 'S' || a[indxOP][ind] == 's') {
+        aChecker.CurveOnSurfaceMode() = Standard_False;
+      }
       else {
         di << "Error: invalid test option(s)!" << "\n";
         di << "Type bopargcheck without arguments for more information" << "\n";
@@ -549,6 +599,7 @@ Standard_Integer bopargcheck
       Standard_Integer S2_SelfIntAll = 0, S2_SmalEAll = 0, S2_BadFAll = 0, S2_BadVAll = 0, S2_BadEAll = 0;
       Standard_Integer S1_OpAb = 0, S2_OpAb = 0;
       Standard_Integer S1_C0 = 0, S2_C0 = 0, S1_C0All = 0, S2_C0All = 0;
+      Standard_Integer S1_COnS = 0, S2_COnS = 0, S1_COnSAll = 0, S2_COnSAll = 0;
       Standard_Boolean hasUnknown = Standard_False;
 
       TCollection_AsciiString aS1SIBaseName("s1si_");
@@ -557,12 +608,14 @@ Standard_Integer bopargcheck
       TCollection_AsciiString aS1BVBaseName("s1bv_");
       TCollection_AsciiString aS1BEBaseName("s1be_");
       TCollection_AsciiString aS1C0BaseName("s1C0_");
+      TCollection_AsciiString aS1COnSBaseName("s1COnS_");
       TCollection_AsciiString aS2SIBaseName("s2si_");
       TCollection_AsciiString aS2SEBaseName("s2se_");
       TCollection_AsciiString aS2BFBaseName("s2bf_");
       TCollection_AsciiString aS2BVBaseName("s2bv_");
       TCollection_AsciiString aS2BEBaseName("s2be_");
       TCollection_AsciiString aS2C0BaseName("s2C0_");
+      TCollection_AsciiString aS2COnSBaseName("s2COnS_");
 
       for(; anIt.More(); anIt.Next()) {
         const BOPAlgo_CheckResult& aResult = anIt.Value();
@@ -667,6 +720,27 @@ Standard_Integer bopargcheck
           }
         }
           break;
+        case BOPAlgo_InvalidCurveOnSurface: {
+          if(!aSS1.IsNull()) {
+            S1_COnS++;
+            if(isL1) {
+              Standard_Real aMaxDist = aResult.GetMaxDistance1();
+              Standard_Real aMaxParameter = aResult.GetMaxParameter1();
+              MakeShapeForFullOutput(aS1COnSBaseName, S1_COnS, aLS1, S1_COnSAll, di,
+                                     Standard_True, aMaxDist, aMaxParameter);
+            }
+          }
+          if(!aSS2.IsNull()) {
+            S2_COnS++;
+            if(isL2) {
+              Standard_Real aMaxDist = aResult.GetMaxDistance2();
+              Standard_Real aMaxParameter = aResult.GetMaxParameter2();
+              MakeShapeForFullOutput(aS2COnSBaseName, S2_COnS, aLS2, S2_COnSAll, di,
+                                     Standard_True, aMaxDist, aMaxParameter);
+            }
+          }
+        }
+          break;
         case BOPAlgo_OperationAborted: {
           if(!aSS1.IsNull()) S1_OpAb++;
           if(!aSS2.IsNull()) S2_OpAb++;
@@ -680,9 +754,9 @@ Standard_Integer bopargcheck
         } // switch
       }// faulties
 
-      Standard_Integer FS1 = S1_SelfInt + S1_SmalE + S1_BadF + S1_BadV + S1_BadE + S1_OpAb + S1_C0;
+      Standard_Integer FS1 = S1_SelfInt + S1_SmalE + S1_BadF + S1_BadV + S1_BadE + S1_OpAb + S1_C0 + S1_COnS;
       FS1 += (S1_BadType != 0) ? 1 : 0;
-      Standard_Integer FS2 = S2_SelfInt + S2_SmalE + S2_BadF + S2_BadV + S2_BadE + S2_OpAb + S2_C0;
+      Standard_Integer FS2 = S2_SelfInt + S2_SmalE + S2_BadF + S2_BadV + S2_BadE + S2_OpAb + S2_C0 + S2_COnS;
       FS2 += (S2_BadType != 0) ? 1 : 0;
       
       // output for first shape
@@ -761,6 +835,17 @@ Standard_Integer bopargcheck
           di << "  Cases(" << S1_C0 << ")  Total shapes(" << S1_C0All << ")" << "\n";
         else
           di << "\n";
+
+        Standard_CString CString17;
+        if (S1_COnS != 0)
+          CString17="YES";
+        else
+          CString17="NO";
+        di << "Invalid Curve on Surface        : " << CString17;
+        if(S1_COnS != 0)
+          di << "  Cases(" << S1_COnS << ")  Total shapes(" << S1_COnSAll << ")" << "\n";
+        else
+          di << "\n";
       }
 
       // output for second shape
@@ -842,6 +927,17 @@ Standard_Integer bopargcheck
         else
           di << "\n";
 
+        Standard_CString CString18;
+        if (S2_COnS != 0)
+          CString18="YES";
+        else
+          CString18="NO";
+        di << "Invalid Curve on Surface        : " << CString18;
+        if(S2_COnS != 0)
+          di << "  Cases(" << S2_COnS << ")  Total shapes(" << S2_COnSAll << ")" << "\n";
+        else
+          di << "\n";
+
         // warning
         di << "\n";
         if(hasUnknown)
@@ -852,3 +948,328 @@ Standard_Integer bopargcheck
 
   return 0;
 }
+
+//=======================================================================
+//function : xdistef
+//purpose  : 
+//=======================================================================
+Standard_Integer xdistef(Draw_Interpretor& di,
+                         Standard_Integer n,
+                         const char** a)
+{
+  if(n < 3) {
+    di << "Use efmaxdist edge face\n";
+    return 1;
+  }
+  //
+  const TopoDS_Shape aS1 = DBRep::Get(a[1]);
+  const TopoDS_Shape aS2 = DBRep::Get(a[2]);
+  //
+  if (aS1.IsNull() || aS2.IsNull()) {
+    di << "null shapes\n";
+    return 1;
+  }
+  //
+  if (aS1.ShapeType() != TopAbs_EDGE || 
+      aS2.ShapeType() != TopAbs_FACE) {
+    di << "type mismatch\n";
+    return 1;
+  }
+  //
+  Standard_Real aMaxDist = 0.0, aMaxPar = 0.0;
+  //
+  const TopoDS_Edge& anEdge = *(TopoDS_Edge*)&aS1;
+  const TopoDS_Face& aFace  = *(TopoDS_Face*)&aS2;
+  //
+  if(!BOPTools_AlgoTools::ComputeTolerance
+     (aFace, anEdge, aMaxDist, aMaxPar)) {
+    di << "Tolerance cannot be computed\n";
+    return 1;
+  }
+  //
+  di << "Max Distance = " << aMaxDist 
+     << "; Parameter on curve = " << aMaxPar << "\n";
+  //
+  return 0;
+}
+
+//=======================================================================
+//function : checkcurveonsurf
+//purpose  : 
+//=======================================================================
+Standard_Integer checkcurveonsurf(Draw_Interpretor& di,
+                                  Standard_Integer n, 
+                                  const char** a)
+{
+  if (n != 2) {
+    di << "use checkcurveonsurf shape\n";
+    return 1;
+  }
+  //
+  TopoDS_Shape aS =  DBRep::Get(a[1]);
+  if (aS.IsNull()) {
+    di << "null shape\n";
+    return 1;
+  }
+  //
+  Standard_Integer nE, nF, anECounter, aFCounter;
+  Standard_Real aT, aTolE, aDMax;
+  TopExp_Explorer aExpF, aExpE;
+  char buf[200], aFName[10], anEName[10];
+  NCollection_DataMap<TopoDS_Shape, Standard_Real, TopTools_ShapeMapHasher> aDMETol;
+  BOPCol_DataMapOfShapeInteger aMSI;
+  //
+  anECounter = 0;
+  aFCounter  = 0;
+  //
+  aExpF.Init(aS, TopAbs_FACE);
+  for (; aExpF.More(); aExpF.Next()) {
+    const TopoDS_Face& aF = *(TopoDS_Face*)&aExpF.Current();
+    //
+    aExpE.Init(aF, TopAbs_EDGE);
+    for (; aExpE.More(); aExpE.Next()) {
+      const TopoDS_Edge& aE = *(TopoDS_Edge*)&aExpE.Current();
+      //
+      if (!BOPTools_AlgoTools::ComputeTolerance(aF, aE, aDMax, aT)) {
+        continue;
+      }
+      //
+      aTolE = BRep_Tool::Tolerance(aE);
+      if (aDMax < aTolE) {
+        continue;
+      }
+      //
+      if (aDMETol.IsBound(aE)) {
+        Standard_Real& aD = aDMETol.ChangeFind(aE);
+        if (aDMax > aD) {
+          aD = aDMax;
+        }
+      }
+      else {
+        aDMETol.Bind(aE, aDMax);
+      }
+      //
+      if (anECounter == 0) {
+        di << "Invalid curves on surface:\n";
+      }
+      //
+      if (aMSI.IsBound(aE)) {
+        nE = aMSI.Find(aE);
+      }
+      else {
+        nE = anECounter;
+        aMSI.Bind(aE, nE);
+        ++anECounter;
+      }
+      //
+      if (aMSI.IsBound(aF)) {
+        nF = aMSI.Find(aF);
+      } else {
+        nF = aFCounter;
+        aMSI.Bind(aF, nF);
+        ++aFCounter;
+      }
+      //
+      Sprintf(anEName, "e_%d", nE);
+      Sprintf(aFName , "f_%d", nF);
+      Sprintf(buf, "edge %s on face %s (max dist: %3.16f, parameter on curve: %3.16f)\n",
+              anEName, aFName, aDMax, aT);
+      di << buf;
+      //
+      DBRep::Set(anEName, aE);
+      DBRep::Set(aFName , aF);
+    }
+  }
+  //
+  if (anECounter > 0) {
+    di << "\n\nSugestions to fix the shape:\n";
+    di << "explode " << a[1] << " e;\n";
+    //
+    TopTools_MapOfShape M;
+    aExpE.Init(aS, TopAbs_EDGE);
+    for (anECounter = 0; aExpE.More(); aExpE.Next()) {
+      const TopoDS_Shape& aE = aExpE.Current();
+      if (!M.Add(aE)) {
+        continue;
+      }
+      ++anECounter;
+      //
+      if (!aDMETol.IsBound(aE)) {
+        continue;
+      }
+      //
+      aTolE = aDMETol.Find(aE);
+      aTolE *= 1.001;
+      sprintf(buf, "settolerance %s_%d %3.16f;\n", a[1], anECounter, aTolE);
+      di << buf;
+    }
+  }
+  else {
+    di << "This shape seems to be OK.\n";
+  }
+  //
+  return 0;
+}
+
+//=======================================================================
+//function : csdeviation
+//purpose  : 
+//=======================================================================
+Standard_Integer csdeviation(Draw_Interpretor& di,
+                             Standard_Integer n,
+                             const char** a)
+{
+  if(n != 3) {
+    di << "Use csdeviation edge face\n";
+    return 1;
+  }
+  //
+  const TopoDS_Shape aS1 = DBRep::Get(a[1]);
+  const TopoDS_Shape aS2 = DBRep::Get(a[2]);
+  //
+  if (aS1.IsNull() || aS2.IsNull()) {
+    di << "null shapes\n";
+    return 1;
+  }
+  //
+  if (aS1.ShapeType() != TopAbs_EDGE || 
+      aS2.ShapeType() != TopAbs_FACE) {
+    di << "type mismatch\n";
+    return 1;
+  }
+  //
+  Standard_Real aMaxDist = 0., aMaxPar = 0.;
+  //
+  const TopoDS_Edge& anEdge = *(TopoDS_Edge*)&aS1;
+  const TopoDS_Face& aFace  = *(TopoDS_Face*)&aS2;
+  //
+  if(!BOPTools_AlgoTools::Deviation
+     (anEdge, aFace, aMaxDist, aMaxPar)) {
+    di << "Deviation cannot be computed\n";
+    return 1;
+  }
+  //
+  di << "Max Distance = " << aMaxDist 
+     << "; Parameter on curve = " << aMaxPar << "\n";
+  //
+  return 0;
+}
+
+//=======================================================================
+//function : checkcsdeviation
+//purpose  : 
+//=======================================================================
+Standard_Integer checkcsdeviation(Draw_Interpretor& di,
+                                  Standard_Integer n, 
+                                  const char** a)
+{
+  if (n != 2) {
+    di << "use checkcsdeviation shape\n";
+    return 1;
+  }
+  //
+  TopoDS_Shape aS =  DBRep::Get(a[1]);
+  if (aS.IsNull()) {
+    di << "null shape\n";
+    return 1;
+  }
+  //
+  Standard_Integer nE, nF, anECounter, aFCounter;
+  Standard_Real aT, aTolE, aDMax;
+  TopExp_Explorer aExpF, aExpE;
+  char buf[200], aFName[10], anEName[10];
+  NCollection_DataMap<TopoDS_Shape, Standard_Real, TopTools_ShapeMapHasher> aDMETol;
+  BOPCol_DataMapOfShapeInteger aMSI;
+  //
+  anECounter = 0;
+  aFCounter  = 0;
+  //
+  aExpF.Init(aS, TopAbs_FACE);
+  for (; aExpF.More(); aExpF.Next()) {
+    const TopoDS_Face& aF = *(TopoDS_Face*)&aExpF.Current();
+    //
+    aExpE.Init(aF, TopAbs_EDGE);
+    for (; aExpE.More(); aExpE.Next()) {
+      const TopoDS_Edge& aE = *(TopoDS_Edge*)&aExpE.Current();
+      //
+      if (!BOPTools_AlgoTools::Deviation(aE, aF, aDMax, aT)) {
+        continue;
+      }
+      //
+      aTolE = BRep_Tool::Tolerance(aE);
+      if (aDMax <= aTolE) {
+        continue;
+      }
+      //
+      if (aDMETol.IsBound(aE)) {
+        Standard_Real& aD = aDMETol.ChangeFind(aE);
+        if (aDMax > aD) {
+          aD = aDMax;
+        }
+      }
+      else {
+        aDMETol.Bind(aE, aDMax);
+      }
+      //
+      if (anECounter == 0) {
+        di << "Invalid curves on surface:\n";
+      }
+      //
+      if (aMSI.IsBound(aE)) {
+        nE = aMSI.Find(aE);
+      }
+      else {
+        nE = anECounter;
+        aMSI.Bind(aE, nE);
+        ++anECounter;
+      }
+      //
+      if (aMSI.IsBound(aF)) {
+        nF = aMSI.Find(aF);
+      } else {
+        nF = aFCounter;
+        aMSI.Bind(aF, nF);
+        ++aFCounter;
+      }
+      //
+      Sprintf(anEName, "e_%d", nE);
+      Sprintf(aFName , "f_%d", nF);
+      Sprintf(buf, "edge %s on face %s (max dist: %3.16f, parameter on curve: %3.16f)\n",
+              anEName, aFName, aDMax, aT);
+      di << buf;
+      //
+      DBRep::Set(anEName, aE);
+      DBRep::Set(aFName , aF);
+    }
+  }
+  //
+  if (anECounter > 0) {
+    di << "\n\nSugestions to fix the shape:\n";
+    di << "explode " << a[1] << " e;\n";
+    //
+    TopTools_MapOfShape M;
+    aExpE.Init(aS, TopAbs_EDGE);
+    for (anECounter = 0; aExpE.More(); aExpE.Next()) {
+      const TopoDS_Shape& aE = aExpE.Current();
+      if (!M.Add(aE)) {
+        continue;
+      }
+      ++anECounter;
+      //
+      if (!aDMETol.IsBound(aE)) {
+        continue;
+      }
+      //
+      aTolE = aDMETol.Find(aE);
+      aTolE *= 1.001;
+      sprintf(buf, "settolerance %s_%d %3.16f;\n", a[1], anECounter, aTolE);
+      di << buf;
+    }
+  }
+  else {
+    di << "This shape seems to be OK.\n";
+  }
+  //
+  return 0;
+}
+

src/BOPTools/BOPTools_AlgoTools.cdl

@@ -33,6 +33,9 @@ uses
     Wire from TopoDS, 
     Shell from TopoDS,  
     Solid from TopoDS, 
+    Curve from Geom,
+    Curve from Geom2d,
+    Surface from Geom,
     --  
     BaseAllocator from BOPCol,
     ListOfShape from BOPCol,  
@@ -458,4 +461,44 @@ is
     	returns Boolean from Standard;  
     	---Purpose: Returns true if the solid <theSolid> is inverted 	
  
+    ComputeTolerance(myclass;
+        theCurve3D : Curve from Geom;
+        theCurve2D : Curve from Geom2d;
+        theSurf    : Surface from Geom;
+        theMaxDist : out Real from Standard; 
+        theMaxPar  : out Real from Standard)
+    returns Boolean from Standard; 
+    ---Purpose:  
+    -- Computes the max distance between points  
+    -- taken from 3D and 2D curves by the same parameter
+
+    ComputeTolerance(myclass;
+        theFace : Face from TopoDS;
+        theEdge : Edge from TopoDS;
+        theMaxDist : out Real from Standard; 
+        theMaxPar  : out Real from Standard)
+    returns Boolean from Standard;
+    ---Purpose:  
+    -- Computes the neccessary value of the tolerance for the edge
+ 
+    Deviation(myclass;
+        theEdge  : Edge from TopoDS; 
+        theFace  : Face from TopoDS;
+        theDist  : out Real from Standard;
+        theParam : out Real from Standard)
+    returns Boolean from Standard;
+    ---Purpose:  
+    -- Computes the max distance between edge and face
+ 
+    Deviation(myclass;
+        theCurve   : Curve from Geom;
+        theSurface : Surface from Geom;  
+        theFirst   : Real from Standard;
+        theLast    : Real from Standard;
+        theDist    : out Real from Standard;
+        theParam   : out Real from Standard)
+    returns Boolean from Standard;
+    ---Purpose:  
+    -- Computes the max distance between curve and surface
+
 end AlgoTools;

src/BOPTools/BOPTools_AlgoTools2D.cxx

@@ -547,7 +547,7 @@ void BOPTools_AlgoTools2D::Make2D (const TopoDS_Edge& aE,
   }
   
   //
-  aToler=.5*BRep_Tool::Tolerance(aE);
+  aToler = BRep_Tool::Tolerance(aE);
   BOPTools_AlgoTools2D::MakePCurveOnFace(aF, C3D2, f3d, l3d, aC2D, aToler);
   //
   aFirst = f3d; 
@@ -616,9 +616,11 @@ void BOPTools_AlgoTools2D::Make2D (const TopoDS_Edge& aE,
     }
   }
   TolReached2d=aTolR;
-  
-  BOPTools_AlgoTools2D::AdjustPCurveOnFace (aF, aFirst, aLast, aC2D, aC2DA);
-  aC2D=aC2DA;
+  //
+  if (!aC2D.IsNull()) {
+    BOPTools_AlgoTools2D::AdjustPCurveOnFace (aF, aFirst, aLast, aC2D, aC2DA);
+    aC2D=aC2DA;
+  }
 }
 
 //=======================================================================

src/BOPTools/BOPTools_AlgoTools_1.cxx

@@ -111,6 +111,378 @@ static
 static
   void UpdateVertices(const TopoDS_Edge& aE);
 
+//=======================================================================
+//class   : BOPTools_CheckCurveOnSurface
+//purpose : it is used to check the curve on the surface
+//=======================================================================
+#include <math_GlobOptMin.hxx>
+#include <math_MultipleVarFunctionWithHessian.hxx>
+#include <math_Matrix.hxx>
+#include <Geom2d_TrimmedCurve.hxx>
+
+class BOPTools_CheckCurveOnSurface : 
+  public math_MultipleVarFunctionWithHessian
+{
+ public:
+  BOPTools_CheckCurveOnSurface(BOPTools_CheckCurveOnSurface&);
+  BOPTools_CheckCurveOnSurface(const Handle(Geom_Curve)& theC3D,
+                              const Handle(Geom2d_Curve)& theC2D,
+                              const Handle(Geom_Surface)& theSurf)
+    :
+      my3DCurve(theC3D),
+      my2DCurve(theC2D),
+      mySurf(theSurf)
+  {
+  }
+  //
+  virtual Standard_Integer NbVariables() const {
+    return 1;
+  }
+  //
+  virtual Standard_Boolean Value(const math_Vector& theX,
+                                 Standard_Real& theFVal)
+  {
+    Standard_Real u = theX(1);
+    Standard_Real aFirst = my3DCurve->FirstParameter();
+    Standard_Real aLast = my3DCurve->LastParameter();
+
+    if ((u < aFirst) || (u > aLast))
+    {
+      return Standard_False;
+    }
+
+    try {
+      const Standard_Real aPar = theX(1);
+      gp_Pnt aP1, aP2;
+      gp_Pnt2d aP2d;
+      my3DCurve->D0(aPar, aP1);
+      my2DCurve->D0(aPar, aP2d);
+      mySurf->D0(aP2d.X(), aP2d.Y(), aP2);
+      //
+      theFVal = -1.0*aP1.SquareDistance(aP2);
+    }
+    catch(...) {
+      return Standard_False;
+    }
+    //
+    return Standard_True;
+  }
+  //
+  virtual Standard_Integer GetStateNumber() {
+    return 0;
+  }
+  //
+  virtual Standard_Boolean Gradient(const math_Vector& theX,
+                                    math_Vector& theGrad) {
+
+  Standard_Real u = theX(1);
+  Standard_Real aRange = (my3DCurve->LastParameter() - my3DCurve->FirstParameter());
+  Standard_Real aFirst = my3DCurve->FirstParameter() - 0.1*aRange;
+  Standard_Real aLast = my3DCurve->LastParameter() + 0.1*aRange;
+
+  if ((u < aFirst) || (u > aLast))
+  {
+    return Standard_False;
+  }
+
+    try {
+      const Standard_Real aPar = theX(1);
+      
+      gp_Pnt aP1, aP2;
+      gp_Vec aDC3D, aDSU, aDSV;
+      gp_Pnt2d aP2d;
+      gp_Vec2d aDC2D;
+      
+      my3DCurve->D1(aPar, aP1, aDC3D);
+      my2DCurve->D1(aPar, aP2d, aDC2D);
+      mySurf->D1(aP2d.X(), aP2d.Y(), aP2, aDSU, aDSV);
+      
+      aP1.SetXYZ(aP1.XYZ() - aP2.XYZ());
+      aP2.SetXYZ(aDC3D.XYZ() - aDC2D.X()*aDSU.XYZ() - aDC2D.Y()*aDSV.XYZ());
+      
+      theGrad(1) = -2.0*aP1.XYZ().Dot(aP2.XYZ());
+    }
+    catch(...) {
+      return Standard_False;
+    }
+    
+    return Standard_True;
+  }
+  //   
+  virtual Standard_Boolean Values(const math_Vector& theX,
+                                  Standard_Real& theVal,
+                                  math_Vector& theGrad) {
+    if(!Value(theX, theVal))
+      return Standard_False;
+    
+    if(!Gradient(theX, theGrad))
+      return Standard_False;
+    
+    return Standard_True;
+  }
+  //
+  virtual Standard_Boolean Values(const math_Vector& theX,
+                                  Standard_Real& theVal,
+                                  math_Vector& theGrad,
+                                  math_Matrix& theHessian) {
+    if(!Value(theX, theVal))
+      return Standard_False;
+    
+    if(!Gradient(theX, theGrad))
+      return Standard_False;
+    
+    theHessian(1,1) = theGrad(1);
+    
+    return Standard_True;
+  }
+  //
+ private:
+  Handle(Geom_Curve) my3DCurve;
+  Handle(Geom2d_Curve) my2DCurve;
+  Handle(Geom_Surface) mySurf;
+};
+
+//=======================================================================
+//class   : BOPTools_CheckCSDistance
+//purpose : the class is used to measure the deviation between 
+//          edge and face
+//=======================================================================
+///--------------------------------------------------
+#include <Precision.hxx>
+#include <GeomAPI_ProjectPointOnSurf.hxx>
+
+#ifdef HAVE_TBB
+#include <tbb\parallel_reduce.h>
+#include <tbb\blocked_range.h>
+
+using namespace tbb;
+#endif
+
+class BOPTools_CheckCSDistance
+{
+  BOPTools_CheckCSDistance operator=(BOPTools_CheckCSDistance&);
+public:
+  BOPTools_CheckCSDistance( const Handle(Geom_Curve)& theC3D,
+                            const Handle(Geom_Surface)& theSurface,
+                            const Standard_Boolean theIsNotOptimize,
+                            const Standard_Real theFirst,
+                            const Standard_Real theLast)
+    :
+      myGoldRatio((sqrt(5.)-1.)/2.),
+      myCurve(theC3D), 
+      mySurface(theSurface),
+      myFirst(theFirst),
+      myLast(theLast),
+      myIsNotOptimize(theIsNotOptimize),
+      myMaxStep(1.e-4*(myLast - myFirst))
+  {
+    Standard_Real aUMin = 0.0, aUMax = 0.0, aVMin = 0.0, aVMax = 0.0;
+    theSurface->Bounds(aUMin, aUMax, aVMin, aVMax);
+    myProjPS.Init(theSurface, aUMin, aUMax, aVMin, aVMax);
+
+    myMaxParam = theFirst;
+    Value(myMaxParam, myMaxValue);
+  }
+
+  //
+  Standard_Boolean Value(const Standard_Real thePar,
+                         Standard_Real& theFVal)
+  {
+    if ((thePar < myFirst) || (thePar > myLast)) {
+      return Standard_False;
+    }
+    //
+    gp_Pnt aPOnC;
+    myCurve->D0(thePar, aPOnC);
+    myProjPS.Perform(aPOnC);
+    theFVal = myProjPS.IsDone() ? myProjPS.LowerDistance() : 0.;
+    //
+    return Standard_True;
+  }
+  //GoldRatio method
+  void MinimizeByGoldenRatio(const Standard_Real theFirst, 
+                             const Standard_Real theLast)
+  {
+    
+    Standard_Real aFirst, aLast, aMidPar1, aMidPar2, aMidVal1, aMidVal2;
+    //
+    aFirst = theFirst;
+    aLast  = theLast;
+    aMidPar1 = aLast - (aLast - aFirst) * myGoldRatio;
+    aMidPar2 = aFirst + (aLast - aFirst) * myGoldRatio;
+    Value(aMidPar1, aMidVal1);
+    Value(aMidPar2, aMidVal2);
+    //
+    Standard_Boolean bIsEqual = Standard_False;
+    //
+    do {
+      if (Abs(aMidVal1 - aMidVal2) <= Precision::Confusion()) {
+        bIsEqual = Standard_True;
+        break;
+      }
+      //
+      if (aMidVal1 > aMidVal2) {
+        aLast = aMidPar2;
+        aMidPar2 = aMidPar1;
+        aMidVal2 = aMidVal1;
+        aMidPar1 = aLast - (aLast - aFirst) * myGoldRatio;
+        Value(aMidPar1, aMidVal1);
+      }
+      else
+      {
+        aFirst = aMidPar1;
+        aMidPar1 = aMidPar2;
+        aMidVal1 = aMidVal2;
+        aMidPar2 = aFirst + (aLast - aFirst) * myGoldRatio;
+        Value(aMidPar2, aMidVal2);
+      }
+    }
+    while ((aMidPar2 - aMidPar1) > myMaxStep);
+
+    const Standard_Real aTMax = (aLast + aFirst)*0.5;
+
+    if(myIsNotOptimize)
+    {
+      if (bIsEqual && ((aMidPar2-aMidPar1) > myMaxStep))
+      {
+        MinimizeByGoldenRatio(aMidPar1, aMidPar2);
+        MinimizeByGoldenRatio(aFirst, aMidPar1);
+        MinimizeByGoldenRatio(aMidPar2, aLast);
+      }
+    }
+
+    //
+    Standard_Real aVal = 0.;
+    Value(aTMax, aVal);
+    //
+    if (aVal > myMaxValue) {
+      myMaxParam = aTMax;
+      myMaxValue = aVal;
+    }
+
+    //if(!myIsNotOptimize)
+    //{
+    //  if( ((aTMax - theFirst) > myMaxStep) &&
+    //      ((theLast - aTMax) > myMaxStep))
+    //  {
+    //    const Standard_Real aSt = 0.5*myMaxStep;
+    //    MinimizeByGoldenRatio(theFirst, aTMax-aSt);
+    //    MinimizeByGoldenRatio(aTMax+aSt, theLast);
+    //  }
+    //}
+  }
+
+  Handle(Geom_Curve) GetCurve() const 
+  {
+    return myCurve;
+  }
+
+  Handle(Geom_Surface) GetSurface() const 
+  {
+    return mySurface;
+  }
+
+  Standard_Real GetFirstParam() const
+  {
+    return myFirst;
+  }
+
+  Standard_Real GetLastParam() const
+  {
+    return myLast;
+  }
+
+  Standard_Boolean GetOptimizeFlag() const
+  {
+    return myIsNotOptimize;
+  }
+  //
+  Standard_Real ParameterOnCurve() const
+  {
+    return myMaxParam;
+  }
+  //
+  Standard_Real MaxDistance() const
+  {
+    return myMaxValue;
+  }
+  //
+
+
+private:
+  const Handle(Geom_Curve) myCurve;
+  const Handle(Geom_Surface) mySurface;
+
+  const Standard_Real myFirst;
+  const Standard_Real myLast;
+  const Standard_Real myMaxStep;
+  const Standard_Real myGoldRatio;
+
+  const Standard_Boolean myIsNotOptimize;
+  //
+  GeomAPI_ProjectPointOnSurf myProjPS;
+  //
+  Standard_Real myMaxValue;
+  Standard_Real myMaxParam;
+};
+
+#ifdef HAVE_TBB
+class BOPTools_CSDParallelComputing
+{
+public:
+  BOPTools_CSDParallelComputing(const BOPTools_CheckCSDistance& theCD,
+    const Standard_Real theDeltaParam): myBTCheckDist(theCD.GetCurve(), theCD.GetSurface(), theCD.GetOptimizeFlag(),theCD.GetFirstParam(), theCD.GetLastParam()), myDeltaParam(theDeltaParam){}
+//
+  BOPTools_CSDParallelComputing(BOPTools_CSDParallelComputing& theCD, split ) :
+  myBTCheckDist(theCD.myBTCheckDist.GetCurve(), theCD.myBTCheckDist.GetSurface(), theCD.myBTCheckDist.GetOptimizeFlag(), theCD.myBTCheckDist.GetFirstParam(), theCD.myBTCheckDist.GetLastParam()),
+      myMaxValue(theCD.myMaxValue),
+      myMaxParam(theCD.myMaxParam), myDeltaParam(theCD.myDeltaParam)  {}
+
+  void join(const BOPTools_CSDParallelComputing& theCD)
+  {
+    if(theCD.myMaxValue > myMaxValue)
+    {
+      myMaxValue = theCD.myMaxValue;
+      myMaxParam = theCD.myMaxParam;
+    }
+  }
+
+  void operator()(const blocked_range<Standard_Integer>& theRange)
+  {
+    for(Standard_Integer anInd = theRange.begin(); anInd != theRange.end(); anInd++)
+    {
+      const Standard_Real aStart = myBTCheckDist.GetFirstParam() + anInd*myDeltaParam;
+      myBTCheckDist.MinimizeByGoldenRatio(aStart, aStart+myDeltaParam);
+
+      if(myBTCheckDist.MaxDistance() > myMaxValue)
+      {
+        myMaxValue = myBTCheckDist.MaxDistance();
+        myMaxParam = myBTCheckDist.ParameterOnCurve();
+      }
+    }
+  }
+
+  Standard_Real ParameterOnCurve() const
+  {
+    return myMaxParam;
+  }
+  //
+  Standard_Real MaxDistance() const
+  {
+    return myMaxValue;
+  }
+
+
+private:
+  BOPTools_CheckCSDistance myBTCheckDist;
+
+  Standard_Real myMaxValue;
+  Standard_Real myMaxParam;
+
+  const Standard_Real myDeltaParam;
+};
+#endif
+
 //=======================================================================
 // Function : CorrectTolerances
 // purpose : 
@@ -812,3 +1184,350 @@ void CheckEdge (const TopoDS_Edge& Ed, const Standard_Real aMaxTol)
     }
   }
 }
+
+//=======================================================================
+// Function : MinComputing
+// purpose : 
+//=======================================================================
+static Standard_Boolean MinComputing( BOPTools_CheckCurveOnSurface& theFunction,
+                                      const Standard_Real theFirst,
+                                      const Standard_Real theLast,
+                                      const Standard_Real theEpsilon, //1.0e-3
+                                      Standard_Real & theBestValue,
+                                      Standard_Real & theBestParameter)
+{
+  //Standard_Real aPrevValue = theBestValue;
+  const Standard_Real aStepMin = 1.0e-2/**Min(1.0, theLast-theFirst)*/;
+  math_Vector aFirstV(1, 1), aLastV(1, 1), anOutputParam(1, 1);
+  aFirstV(1) = theFirst;
+  aLastV(1) = theLast;
+
+  math_GlobOptMin aFinder(&theFunction, aFirstV, aLastV);
+  aFinder.SetTol(aStepMin, theEpsilon);
+  aFinder.Perform();
+
+  const Standard_Integer aNbExtr = aFinder.NbExtrema();
+  for(Standard_Integer i = 1; i <= aNbExtr; i++)
+  {
+    Standard_Real aValue = 0.0;
+    aFinder.Points(i, anOutputParam);
+    theFunction.Value(anOutputParam, aValue);
+    
+    if(aValue < theBestValue)
+    {
+      theBestValue = aValue;
+      theBestParameter = anOutputParam(1);
+    }
+  }
+
+  return Standard_True;
+}
+
+
+//=======================================================================
+// Function : ComputeTolerance
+// purpose : 
+//=======================================================================
+Standard_Boolean BOPTools_AlgoTools::
+          ComputeTolerance( const Handle(Geom_Curve)& theCurve3D,
+                            const Handle(Geom2d_Curve)& theCurve2D,
+                            const Handle(Geom_Surface)& theSurf,
+                            Standard_Real& theMaxDist,
+                            Standard_Real& theMaxPar)
+{
+  if (theCurve3D.IsNull() ||
+      theCurve2D.IsNull() ||
+      theSurf.IsNull()) {
+    return Standard_False;
+  }
+
+  const Standard_Real anEpsilonRange = 1.0e-3, aMinDelta = 1.0e-5;
+
+  //
+  try {
+    Standard_Real aFirst = theCurve3D->FirstParameter(),
+                  aLast = theCurve3D->LastParameter();
+
+    BOPTools_CheckCurveOnSurface aFunc(theCurve3D, theCurve2D, theSurf);
+    //
+    math_Vector anOutputParam(1, 1);
+    anOutputParam(1) = theMaxPar = aFirst;
+    //
+    theMaxDist = 0.;
+    MinComputing(aFunc, aFirst, aLast, anEpsilonRange, theMaxDist, theMaxPar);
+
+    Standard_Integer aNbIteration = 100;
+    Standard_Boolean aStatus = Standard_True;
+    while((aNbIteration-- >= 0) && aStatus)
+    {
+      Standard_Real aValue = theMaxDist, aParam = theMaxPar;
+      Standard_Real aBP = theMaxPar - aMinDelta;
+      MinComputing(aFunc, aFirst, aBP, anEpsilonRange, theMaxDist, theMaxPar);
+
+      if(theMaxDist < aValue)
+      {
+        aLast = aBP;
+        aStatus = Standard_True;
+      }
+      else 
+      {
+        theMaxDist = aValue;
+        theMaxPar = aParam;
+        aStatus = Standard_False;
+      }
+
+      if(!aStatus)
+      {
+        aBP = theMaxPar + aMinDelta;
+        MinComputing(aFunc, aBP, aLast, 1.0e-3, theMaxDist, theMaxPar);
+
+        if(theMaxDist < aValue)
+        {
+          aFirst = aBP;
+          aStatus = Standard_True;
+        }
+        else
+        {
+          theMaxDist = aValue;
+          theMaxPar = aParam;
+          aStatus = Standard_False;
+        }
+      }
+    }
+
+    theMaxDist = sqrt(Abs(theMaxDist));
+  }
+  catch (...) {
+    return Standard_False;
+  }
+  //
+  return Standard_True;
+}
+
+//=======================================================================
+// Function : ComputeTolerance
+// purpose : 
+//=======================================================================
+Standard_Boolean BOPTools_AlgoTools::ComputeTolerance
+  (const TopoDS_Face& theFace,
+   const TopoDS_Edge& theEdge,
+   Standard_Real& theMaxDist,
+   Standard_Real& theParameter)
+{
+  Standard_Boolean bRet;
+  Standard_Real aT, aD, aFirst, aLast;
+  TopLoc_Location aLocC, aLocS;
+  //
+  theMaxDist = 0.;
+  theParameter = 0.;
+  bRet = Standard_False;
+  //
+  const Handle(BRep_TEdge)& aTE = *((Handle(BRep_TEdge)*)&theEdge.TShape());
+  //The edge is considered to be same range and not degenerated
+  if ((!aTE->SameRange() && aTE->SameParameter()) ||
+      aTE->Degenerated()) {
+    return bRet;
+  }
+  //
+  Handle(Geom_Curve) aC = Handle(Geom_Curve)::
+    DownCast(BRep_Tool::Curve(theEdge, aLocC, aFirst, aLast)->Copy());
+  aC = new Geom_TrimmedCurve(aC, aFirst, aLast);
+  aC->Transform(aLocC.Transformation());
+  //
+  const Handle(Geom_Surface)& aSurfF = BRep_Tool::Surface(theFace, aLocS);
+  const Handle(Geom_Surface)& aSurf = Handle(Geom_Surface)::
+    DownCast(aSurfF->Copy()->Transformed(aLocS.Transformation()));
+  //
+  Standard_Boolean isPCurveFound = Standard_False;
+  BRep_ListIteratorOfListOfCurveRepresentation itcr(aTE->Curves());
+  for (; itcr.More(); itcr.Next()) {
+    const Handle(BRep_CurveRepresentation)& cr = itcr.Value();
+    if (!(cr->IsCurveOnSurface(aSurfF, aLocS.Predivided(theEdge.Location())))) {
+      continue;
+    }
+    isPCurveFound = Standard_True;
+    //
+    Handle(Geom2d_Curve) aC2d = Handle(Geom2d_Curve)::
+      DownCast(cr->PCurve()->Copy());
+    aC2d = new Geom2d_TrimmedCurve(aC2d, aFirst, aLast);
+      //
+    if(BOPTools_AlgoTools::ComputeTolerance
+       (aC, aC2d, aSurf, aD, aT)) {
+      bRet = Standard_True;
+      if (aD > theMaxDist) {
+        theMaxDist = aD;
+        theParameter = aT;
+      }
+    }
+    //
+    if (cr->IsCurveOnClosedSurface()) {
+      Handle(Geom2d_Curve) aC2d = Handle(Geom2d_Curve)::
+        DownCast(cr->PCurve2()->Copy());
+      aC2d = new Geom2d_TrimmedCurve(aC2d, aFirst, aLast);
+      //
+      if(BOPTools_AlgoTools::ComputeTolerance
+         (aC, aC2d, aSurf, aD, aT)) {
+        bRet = Standard_True;
+        if (aD > theMaxDist) {
+          theMaxDist = aD;
+          theParameter = aT;
+        }
+      }
+    }
+  }
+  //
+  if (isPCurveFound) {
+    return bRet;
+  }
+  //
+  Handle(Geom_Plane) aPlane;
+  Handle(Standard_Type) dtyp = aSurf->DynamicType();
+  //
+  if (dtyp == STANDARD_TYPE(Geom_RectangularTrimmedSurface)) {
+    aPlane = Handle(Geom_Plane)::
+      DownCast(Handle(Geom_RectangularTrimmedSurface)::
+               DownCast(aSurf)->BasisSurface()->Copy());
+  }
+  else {
+    aPlane = Handle(Geom_Plane)::DownCast(aSurf->Copy());
+  }
+  //
+  if (aPlane.IsNull()) { // not a plane
+    return bRet;
+  }
+  //
+  aPlane = Handle(Geom_Plane)::DownCast(aPlane);//
+  //
+  Handle(GeomAdaptor_HSurface) GAHS = new GeomAdaptor_HSurface(aPlane);
+  Handle(Geom_Curve) ProjOnPlane = 
+    GeomProjLib::ProjectOnPlane (new Geom_TrimmedCurve(aC, aFirst, aLast), 
+                                 aPlane, aPlane->Position().Direction(), 
+                                 Standard_True);
+  Handle(GeomAdaptor_HCurve) aHCurve = new GeomAdaptor_HCurve(ProjOnPlane);
+  //
+  ProjLib_ProjectedCurve proj(GAHS,aHCurve);
+  Handle(Geom2d_Curve) aC2d = Geom2dAdaptor::MakeCurve(proj);
+  aC2d = new Geom2d_TrimmedCurve(aC2d, aFirst, aLast);
+  //
+  if(BOPTools_AlgoTools::ComputeTolerance
+     (aC, aC2d, aPlane, aD, aT)) {
+    bRet = Standard_True;
+    if (aD > theMaxDist) {
+      theMaxDist = aD;
+      theParameter = aT;
+    }
+  }
+  //
+  return bRet;
+}
+
+//=======================================================================
+// Function : Deviation
+// purpose : 
+//=======================================================================
+Standard_Boolean BOPTools_AlgoTools::Deviation
+  (const TopoDS_Edge& theEdge,
+   const TopoDS_Face& theFace,
+   Standard_Real& theMaxDist,
+   Standard_Real& theParam)
+{
+  Standard_Boolean bFound = Standard_False;
+  //
+  if (theEdge.IsNull() || theFace.IsNull()) {
+    return bFound;
+  }
+  //
+  if (BRep_Tool::Degenerated(theEdge) ||
+      !BRep_Tool::IsGeometric(theEdge)) {
+    return bFound;
+  }
+  //
+  Standard_Real aFirst, aLast;
+  const Handle(Geom_Curve)& aC = BRep_Tool::Curve(theEdge, aFirst, aLast);
+  const Handle(Geom_Surface)& aS = BRep_Tool::Surface(theFace);
+  //
+  bFound = BOPTools_AlgoTools::Deviation
+    (aC, aS, aFirst, aLast, theMaxDist, theParam);
+  //
+  return bFound;
+}
+
+//=======================================================================
+// Function : Deviation
+// purpose : 
+//=======================================================================
+
+
+
+Standard_Boolean BOPTools_AlgoTools::Deviation
+  (const Handle(Geom_Curve)& theCurve,
+   const Handle(Geom_Surface)& theSurface,
+   const Standard_Real theFirst,
+   const Standard_Real theLast,
+   Standard_Real& theMaxDist,
+   Standard_Real& theParam)
+{
+  if (theCurve.IsNull() || theSurface.IsNull()) {
+    return Standard_False;
+  }
+  //
+  const Standard_Integer aNbParts = 10;
+  const Standard_Real aDeltaPar = (theLast - theFirst) / aNbParts;
+  theMaxDist = 0.;
+  theParam = theFirst;
+
+  const Standard_Boolean isNotOptimize =  ((theCurve->Continuity() == GeomAbs_C0) || 
+                                          (theSurface->Continuity() == GeomAbs_C0));
+
+
+  BOPTools_CheckCSDistance aFunc(theCurve, theSurface, isNotOptimize, theFirst, theLast);
+
+
+  if(isNotOptimize)
+  {
+#ifdef HAVE_TBB
+    BOPTools_CSDParallelComputing aParall(aFunc, aDeltaPar);
+    parallel_reduce(blocked_range<Standard_Integer>(0,aNbParts-2), aParall);
+    const Standard_Real aDist = aParall.MaxDistance();
+
+    if (aDist > theMaxDist)
+    {
+      theMaxDist = aDist;
+      theParam = aParall.ParameterOnCurve();
+    }
+#else
+    Standard_Real aT1, aT2 = theFirst;
+    for (;;)
+    {
+      aT1 = aT2;
+      aT2 += aDeltaPar;
+
+      if (aT2 > theLast)
+      {
+        break;
+      }
+
+      aFunc.MinimizeByGoldenRatio(aT1, aT2);
+      
+      const Standard_Real aDist = aFunc.MaxDistance();
+      if (aDist > theMaxDist)
+      {
+        theMaxDist = aDist;
+        theParam = aFunc.ParameterOnCurve();
+      }
+    }
+#endif
+  }
+  else
+  {
+    aFunc.MinimizeByGoldenRatio(theFirst, theLast);
+    const Standard_Real aDist = aFunc.MaxDistance();
+    if (aDist > theMaxDist)
+    {
+      theMaxDist = aDist;
+      theParam = aFunc.ParameterOnCurve();
+    }
+  }
+
+  return Standard_True;
+}

src/BRepClass3d/BRepClass3d_SClassifier.cxx

@@ -33,6 +33,8 @@
 #include <math_RealRandom.hxx>
 #include <BRepTopAdaptor_FClass2d.hxx>
 
+#include <vector>
+
 static
   Standard_Boolean FaceNormal (const TopoDS_Face& aF,
                                const Standard_Real U,
@@ -98,76 +100,80 @@ void BRepClass3d_SClassifier::PerformInfinitePoint(BRepClass3d_SolidExplorer& aS
   gp_Pnt aPoint;
   gp_Dir aDN;
 
-  math_RealRandom RandomGenerator(0.1, 0.9);
+  math_RealRandom aRandomGenerator(0.1, 0.9);
   myFace.Nullify();
   myState=2;
 
-  aSE.InitShell();
-  if (aSE.MoreShell())
+  // Collect faces in sequence to iterate
+  std::vector<TopoDS_Face> aFaces;
+  for (aSE.InitShell(); aSE.MoreShell(); aSE.NextShell())
   {
-    aSE.InitFace();
-    if (aSE.MoreFace())
+    for (aSE.InitFace(); aSE.MoreFace(); aSE.NextFace())
     {
-      TopoDS_Face aF = aSE.CurrentFace();
+      aFaces.push_back (aSE.CurrentFace());
+    }
+  }
+
+  // iteratively try up to 10 probing points from each face
+  const int NB_MAX_POINTS_PER_FACE = 10;
+  for (int itry = 0; itry < NB_MAX_POINTS_PER_FACE; itry++)
+  {
+    for (std::vector<TopoDS_Face>::iterator iFace = aFaces.begin(); iFace != aFaces.end(); ++iFace)
+    {
+      TopoDS_Face aF = *iFace;
+
       TopAbs_State aState = TopAbs_OUT;
       IntCurveSurface_TransitionOnCurve aTransition = IntCurveSurface_Tangent;
-      TopoDS_Face MinFace = aF;
-      for (;;)
-      {
-        aParam = RandomGenerator.Next();
-	bFound = aSE.FindAPointInTheFace(aF, aPoint, aU, aV, aParam);
-	if (!bFound)
-	  return;
 
-        if (!FaceNormal(aF, aU, aV, aDN))
-          continue;
-        gp_Lin aLin(aPoint, -aDN);
-        Standard_Real parmin = RealLast();
-        for (aSE.InitShell();aSE.MoreShell();aSE.NextShell()) { 
-          if (aSE.RejectShell(aLin) == Standard_False) { 
-            for (aSE.InitFace();aSE.MoreFace(); aSE.NextFace()) {
-              if (aSE.RejectFace(aLin) == Standard_False) { 
-                TopoDS_Shape aLocalShape = aSE.CurrentFace();
-                TopoDS_Face CurFace = TopoDS::Face(aLocalShape);
-                IntCurvesFace_Intersector& Intersector3d = aSE.Intersector(CurFace);
-                Intersector3d.Perform(aLin,-RealLast(),parmin); 
-                
-                if(Intersector3d.IsDone()) {
-                  if(Intersector3d.NbPnt()) { 
-                    Standard_Integer imin = 1;
-                    for (Standard_Integer i = 2; i <= Intersector3d.NbPnt(); i++)
-                      if (Intersector3d.WParameter(i) < Intersector3d.WParameter(imin))
-                        imin = i;
-                    parmin = Intersector3d.WParameter(imin);
-                    aState = Intersector3d.State(imin);
-                    aTransition = Intersector3d.Transition(imin);
-                    MinFace = CurFace;
-                  }
+      aParam = 0.1 + 0.8 * aRandomGenerator.Next(); // random number in range [0.1, 0.9]
+      bFound = aSE.FindAPointInTheFace(aF, aPoint, aU, aV, aParam);
+      if (!bFound || !FaceNormal(aF, aU, aV, aDN))
+        continue;
+
+      gp_Lin aLin(aPoint, -aDN);
+      Standard_Real parmin = RealLast();
+      for (aSE.InitShell();aSE.MoreShell();aSE.NextShell()) { 
+        if (aSE.RejectShell(aLin) == Standard_False) { 
+          for (aSE.InitFace();aSE.MoreFace(); aSE.NextFace()) {
+            if (aSE.RejectFace(aLin) == Standard_False) { 
+              TopoDS_Shape aLocalShape = aSE.CurrentFace();
+              TopoDS_Face CurFace = TopoDS::Face(aLocalShape);
+              IntCurvesFace_Intersector& Intersector3d = aSE.Intersector(CurFace);
+              Intersector3d.Perform(aLin,-RealLast(),parmin); 
+
+              if(Intersector3d.IsDone()) {
+                if(Intersector3d.NbPnt()) { 
+                  Standard_Integer imin = 1;
+                  for (Standard_Integer i = 2; i <= Intersector3d.NbPnt(); i++)
+                    if (Intersector3d.WParameter(i) < Intersector3d.WParameter(imin))
+                      imin = i;
+                  parmin = Intersector3d.WParameter(imin);
+                  aState = Intersector3d.State(imin);
+                  aTransition = Intersector3d.Transition(imin);
                 }
               }
             }
           }
-          else
-            myState = 1;
-        } //end of loop on the whole solid
-        
-        if (aState == TopAbs_IN)
-        {
-          if (aTransition == IntCurveSurface_Out) { 
-            //-- The line is going from inside the solid to outside 
-            //-- the solid.
-            myState = 3; //-- IN --
-            return;
-          }
-          else if (aTransition == IntCurveSurface_In) { 
-            myState = 4; //-- OUT --
-            return;
-          }
         }
-        aF = MinFace;
+        else
+          myState = 1;
+      } //end of loop on the whole solid
+        
+      if (aState == TopAbs_IN)
+      {
+        if (aTransition == IntCurveSurface_Out) { 
+          //-- The line is going from inside the solid to outside 
+          //-- the solid.
+          myState = 3; //-- IN --
+          return;
+        }
+        else if (aTransition == IntCurveSurface_In) { 
+          myState = 4; //-- OUT --
+          return;
+        }
       }
-    } //if (aSE.MoreFace())
-  } //if (aSE.MoreShell())
+    } // iteration by faces
+  } // iteration by points
 }
 
 //=======================================================================

src/GCPnts/GCPnts_TangentialDeflection.cxx

@@ -21,7 +21,80 @@
 #include <gp_Vec2d.hxx>
 #include <gp_XYZ.hxx>
 #include <Precision.hxx>
+#include <TColStd_Array1OfReal.hxx>
+#include <math_Function.hxx>
+#include <gp_Lin.hxx>
+#include <gp_Lin2d.hxx>
+//
+class GCPnts_DistFunction : public math_Function
+{
+public:
+  Standard_EXPORT GCPnts_DistFunction(const Adaptor3d_Curve& theCurve,
+                                      const Standard_Real U1, const Standard_Real U2)
+  : myCurve(theCurve),
+    myU1(U1), myU2(U2)
+  {
+    gp_Pnt P1 = theCurve.Value(U1), P2 = theCurve.Value(U2);
+    myLin = gp_Lin(P1, P2.XYZ() - P1.XYZ());
+  }
 
+  //
+  Standard_EXPORT GCPnts_DistFunction(const GCPnts_DistFunction& theOther);
+
+
+  Standard_EXPORT virtual Standard_Boolean Value (const Standard_Real X,
+                                                        Standard_Real& F)
+  {
+    if (X < myU1 || X > myU2)
+      return Standard_False;
+    //
+    F = -myLin.SquareDistance(myCurve.Value(X));
+    return Standard_True;
+  }
+
+private:
+  GCPnts_DistFunction & operator = (const GCPnts_DistFunction & theOther);
+
+  const Adaptor3d_Curve& myCurve;
+  gp_Lin myLin;
+  Standard_Real myU1;
+  Standard_Real myU2;
+};
+//
+class GCPnts_DistFunction2d : public math_Function
+{
+public:
+  Standard_EXPORT GCPnts_DistFunction2d(const Adaptor2d_Curve2d& theCurve,
+                                        const Standard_Real U1, const Standard_Real U2)
+  : myCurve(theCurve),
+    myU1(U1), myU2(U2)
+  {
+    gp_Pnt2d P2d1 = theCurve.Value(U1), P2d2 = theCurve.Value(U2);
+    myLin = gp_Lin2d(P2d1, P2d2.XY() - P2d1.XY());
+  }
+  //
+  Standard_EXPORT GCPnts_DistFunction2d(const GCPnts_DistFunction2d& theOther);
+
+  Standard_EXPORT virtual Standard_Boolean Value (const Standard_Real X,
+                                                        Standard_Real& F)
+  {
+    if (X < myU1 || X > myU2)
+      return Standard_False;
+    //
+    gp_Pnt2d aP2d = myCurve.Value(X);
+    F = -myLin.SquareDistance(aP2d);
+    return Standard_True;
+  }
+
+private:
+  GCPnts_DistFunction2d & operator = (const GCPnts_DistFunction2d & theOther);
+
+  const Adaptor2d_Curve2d& myCurve;
+  gp_Lin2d myLin;
+  Standard_Real myU1;
+  Standard_Real myU2;
+};
+//
 inline static void D0 (const Adaptor3d_Curve& C, const Standard_Real U, gp_Pnt& P)
 {
   C.D0 (U, P);
@@ -58,6 +131,37 @@ static void D2 (const Adaptor2d_Curve2d& C, const Standard_Real U,
   VV2.SetCoord (X, Y, 0.0);
 }
 
+static Standard_Real EstimAngl(const gp_Pnt& P1, const gp_Pnt& Pm, const gp_Pnt& P2)
+{
+  gp_Vec V1(P1, Pm), V2(Pm, P2);
+  Standard_Real L = V1.Magnitude() * V2.Magnitude();
+  //
+  if(L > gp::Resolution())
+  {
+    return V1.CrossMagnitude(V2)/L;
+  }
+  else
+  {
+    return 0.;
+  }
+}
+// Return number of interval of continuity on which theParam is located.
+// Last parameter is used to increase search speed.
+static Standard_Integer getIntervalIdx(const Standard_Real theParam, 
+                                       TColStd_Array1OfReal& theIntervs,
+                                       const Standard_Integer thePreviousIdx)
+{
+  Standard_Integer anIdx;
+  for(anIdx = thePreviousIdx; anIdx < theIntervs.Upper(); anIdx++)
+  {
+    if (theParam >= theIntervs(anIdx) && 
+        theParam <= theIntervs(anIdx + 1)) // Inside of anIdx interval.
+    {
+      break;
+    }
+  }
+  return anIdx;
+}
 
 //=======================================================================
 //function : CPnts_TangentialDeflection
@@ -114,10 +218,13 @@ Standard_Integer GCPnts_TangentialDeflection::AddPoint
 #define TheCurve Adaptor3d_Curve
 #define Handle_TheBezierCurve   Handle(Geom_BezierCurve)
 #define Handle_TheBSplineCurve  Handle(Geom_BSplineCurve)
+#define TheMaxCurvLinDist GCPnts_DistFunction
 #include <GCPnts_TangentialDeflection.gxx>
 #undef Handle_TheBezierCurve
 #undef Handle_TheBSplineCurve
 #undef TheCurve
+#undef TheMaxCurvLinDist
+
 
 
 #include <Geom2d_BezierCurve.hxx>
@@ -126,10 +233,13 @@ Standard_Integer GCPnts_TangentialDeflection::AddPoint
 #define TheCurve Adaptor2d_Curve2d
 #define Handle_TheBezierCurve   Handle(Geom2d_BezierCurve)
 #define Handle_TheBSplineCurve  Handle(Geom2d_BSplineCurve)
+#define TheMaxCurvLinDist GCPnts_DistFunction2d
 #include <GCPnts_TangentialDeflection.gxx>
 #undef Handle_TheBezierCurve
 #undef Handle_TheBSplineCurve
 #undef TheCurve
+#undef TheMaxCurvLinDist
+
 
 
 

src/GCPnts/GCPnts_TangentialDeflection.gxx

@@ -21,8 +21,64 @@
 #include <gp_Pnt.hxx>
 #include <gp_Vec.hxx>
 #include <gp.hxx>
-
+#include <TColStd_Array1OfReal.hxx>
+#include <math_BrentMinimum.hxx>
+//
 #define Us3 0.3333333333333333333333333333
+//=======================================================================
+//function : EstimDefl
+//purpose  : Estimation of maximal deflection for interval [U1, U2]
+//           
+//=======================================================================
+static void EstimDefl (const TheCurve& C,
+                       const Standard_Real U1, const Standard_Real U2,
+                       const Standard_Real Tol, 
+                       Standard_Real& MaxDefl, Standard_Real& UMax)
+{
+  //
+  TheMaxCurvLinDist aFunc(C, U1, U2);
+ //
+  math_BrentMinimum anOptLoc(Tol);
+  anOptLoc.Perform(aFunc, U1, (U1+U2)/2., U2);
+  if(anOptLoc.IsDone())
+  {
+    MaxDefl = Sqrt(-anOptLoc.Minimum());
+    UMax = anOptLoc.Location();
+    return;
+  }
+  //
+  Standard_Real Du = (C.LastParameter() - C.FirstParameter());
+  Standard_Integer aNbInt = Max(16, RealToInt(64 * (U2 - U1) / Du));
+  Standard_Real du = (U2 - U1) / aNbInt;
+  Standard_Real umax = (U1 + U2) / 2.;
+  Standard_Real dmax = 0., d = 0.;
+  aFunc.Value(umax, dmax);
+  //
+  Standard_Real u;
+  for(u = U1 + du; u <= U2 - du; u += du)
+  {
+    aFunc.Value(u, d);
+    if(d < dmax) //d, dmax are negative values
+    {
+      dmax = d;
+      umax = u;
+    }
+  }
+  //
+  //math_BrentMinimum anOptLoc(Tol);
+  anOptLoc.Perform(aFunc, umax - du, umax, umax + du);
+  if(anOptLoc.IsDone())
+  {
+    MaxDefl = Sqrt(-anOptLoc.Minimum());
+    UMax = anOptLoc.Location();
+  }
+  else
+  {
+    MaxDefl = Sqrt(-dmax);
+    UMax = umax;
+  }
+
+}
 
 void GCPnts_TangentialDeflection::EvaluateDu (
   const TheCurve& C,
@@ -244,7 +300,7 @@ void GCPnts_TangentialDeflection::PerformCircular (const TheCurve& C)
 void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
 
 {
-  Standard_Integer i;       
+  Standard_Integer i, j;       
   gp_XYZ V1, V2;
   gp_Pnt MiddlePoint, CurrentPoint, LastPoint;   
   Standard_Real Du, Dusave, MiddleU, L1, L2;
@@ -264,62 +320,82 @@ void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
   parameters.Append (U1);
   points    .Append (CurrentPoint);
 
+  // Used to detect "isLine" current bspline and in Du computation in general handling.
+  TheCurve* CPtr = const_cast<TheCurve*>(&C);
+  Standard_Integer NbInterv = CPtr->NbIntervals(GeomAbs_CN);
+  TColStd_Array1OfReal Intervs(1, NbInterv+1);
+  CPtr->Intervals(Intervs, GeomAbs_CN);
+
   if (NotDone) {
     //C'est soit une droite, soit une singularite :
-    V1 = LastPoint.XYZ ();
-    V1.Subtract (CurrentPoint.XYZ());
+    V1 = (LastPoint.XYZ() - CurrentPoint.XYZ());
     L1 = V1.Modulus ();
-    if (L1 > LTol) {
+    if (L1 > LTol)
+    {
       //Si c'est une droite on verifie en calculant minNbPoints :
       Standard_Boolean IsLine   = Standard_True;
-      Standard_Integer NbPoints = 3;
-      if (minNbPnts > 3) NbPoints = minNbPnts;
-      Du = (lastu-firstu)/NbPoints;
-      MiddleU = firstu + Du;
-      for (i = 2; i < NbPoints; i++) {
-        D0 (C, MiddleU, MiddlePoint);
-        V2 = MiddlePoint.XYZ();
-        V2.Subtract (CurrentPoint.XYZ());
-        L2 = V2.Modulus ();
-        if (L2 > LTol) {
-          if (((V2.CrossMagnitude (V1))/(L1*L2)) >= ATol) {
-            //C'etait une singularite
-            IsLine = Standard_False;
-            break;
-          }
-          if (minNbPnts > 2) {
-            parameters.Append (MiddleU);
-            points    .Append (MiddlePoint);
+      Standard_Integer NbPoints = (minNbPnts > 3) ? minNbPnts : 3;
+      ////
+      Standard_Real param = 0.;
+      for (i = 1; i <= NbInterv && IsLine; ++i)
+      {
+        // Avoid usage intervals out of [firstu, lastu].
+        if ((Intervs(i+1) < firstu) ||
+            (Intervs(i)   > lastu))
+        {
+          continue;
+        }
+        // Fix border points in applicable intervals, to avoid be out of target interval.
+        if ((Intervs(i)   < firstu) &&
+            (Intervs(i+1) > firstu))
+        {
+          Intervs(i) = firstu;
+        }
+        if ((Intervs(i)   < lastu) &&
+            (Intervs(i+1) > lastu))
+        {
+          Intervs(i + 1) = lastu;
+        }
+
+        Standard_Real delta = (Intervs(i+1) - Intervs(i))/NbPoints;
+        for (j = 1; j <= NbPoints && IsLine; ++j)
+        {
+          param = Intervs(i) + j*delta;
+          D0 (C, param, MiddlePoint);
+          V2 = (MiddlePoint.XYZ() - CurrentPoint.XYZ());
+          L2 = V2.Modulus ();
+          if (L2 > LTol)
+          {
+            const Standard_Real aAngle = V2.CrossMagnitude(V1)/(L1*L2);
+            IsLine = (aAngle < ATol);
           }
         }
-        MiddleU += Du;
       }
-      if (IsLine) {
-        //C'etait une droite (plusieurs poles alignes), Calcul termine :
-        parameters.Append (lastu);
-        points    .Append (LastPoint);
+
+      if (IsLine)
+      {
+        parameters.Clear();
+        points    .Clear();
+
+        PerformLinear(C);
         return;
       }
-      else {
+      else
+      {
         //c'etait une singularite on continue :
-        Standard_Integer pointsLength=points.Length ();
-        for (i = 2; i <= pointsLength; i++) {
-          points    .Remove (i);
-          parameters.Remove (i);
-          pointsLength--;
-        }
-        Du = Dusave;
+        //Du = Dusave;
+        EvaluateDu (C, param, MiddlePoint, Du, NotDone);
       }
     }
-    else  {
-      
+    else
+    {
       Du = (lastu-firstu)/2.1;
       MiddleU = firstu + Du;
       D0 (C, MiddleU, MiddlePoint);
-      V1 = MiddlePoint.XYZ ();
-      V1.Subtract (CurrentPoint.XYZ());
+      V1 = (MiddlePoint.XYZ() - CurrentPoint.XYZ());
       L1 = V1.Modulus ();
-      if (L1 < LTol) {
+      if (L1 < LTol)
+      {
         // L1 < LTol C'est une courbe de longueur nulle, calcul termine :
         // on renvoi un segment de 2 points   (protection)
         parameters.Append (lastu);
@@ -345,12 +421,14 @@ void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
   Standard_Boolean MorePoints = Standard_True;
   Standard_Real U2            = firstu;   
   Standard_Real AngleMax      = angularDeflection * 0.5;  //car on prend le point milieu
-
+  Standard_Integer aIdx[2] = {Intervs.Lower(), Intervs.Lower()}; // Indexes of intervals of U1 and U2, used to handle non-uniform case.
+  Standard_Boolean isNeedToCheck = Standard_False;
   gp_Pnt aPrevPoint = points.Last();
+  Standard_Integer MaxNbCorr = 100;
 
   while (MorePoints) {
-
-    U2 += Du;                             
+    aIdx[0] = getIntervalIdx(U1, Intervs, aIdx[0]);
+    U2 += Du;
 
     if (U2 >= lastu) {                       //Bout de courbe
       U2 = lastu;
@@ -360,35 +438,62 @@ void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
     }
     else D0 (C, U2, CurrentPoint);           //Point suivant
 
-    Standard_Real Coef, ACoef = 0., FCoef = 0.;
-    Standard_Boolean Correction, TooLarge, TooSmall;
+    Standard_Real Coef = 0.0, ACoef = 0., FCoef = 0.;
+    Standard_Boolean Correction, TooLarge;
     TooLarge   = Standard_False;
-    TooSmall   = Standard_False;
     Correction = Standard_True;
+    Standard_Integer nbcorr = 0;
 
-    while (Correction) {                     //Ajustement Du
+    while (Correction && (++nbcorr <= MaxNbCorr)) {                     //Ajustement Du
+      if (isNeedToCheck)
+      {
+        aIdx[1] = getIntervalIdx(U2, Intervs, aIdx[0]);
+        if (aIdx[1] > aIdx[0]) // Jump to another polynom.
+        {
+          if (Du > (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3) // Set Du to the smallest value and check deflection on it.
+          {
+            Du = (Intervs(aIdx[0] + 1) - Intervs(aIdx[0]) ) * Us3;
+            U2 = U1 + Du;
+            if (U2 > lastu)
+              U2 = lastu;
+            D0 (C, U2, CurrentPoint);
+          }
+        }
+      }
       MiddleU = (U1+U2)*0.5;                 //Verif / au point milieu
       D0 (C, MiddleU, MiddlePoint);
 
-      V1 = CurrentPoint.XYZ ();              //Critere de fleche
-      V1.Subtract (aPrevPoint.XYZ());
-      V2 = MiddlePoint.XYZ ();
-      V2.Subtract (aPrevPoint.XYZ());
+      V1 = (CurrentPoint.XYZ() - aPrevPoint.XYZ()); //Critere de fleche
+      V2 = (MiddlePoint.XYZ()  - aPrevPoint.XYZ());
       L1 = V1.Modulus ();
-      if (L1 > LTol) FCoef = V1.CrossMagnitude(V2)/(L1*curvatureDeflection);
-      else           FCoef = 0.0;
 
-      V1 = CurrentPoint.XYZ ();              //Critere d'angle
-      V1.Subtract (MiddlePoint.XYZ ());
+      FCoef = (L1 > LTol) ? 
+        V1.CrossMagnitude(V2)/(L1*curvatureDeflection) : 0.0;
+
+      V1 = (CurrentPoint.XYZ() - MiddlePoint.XYZ()); //Critere d'angle
       L1 = V1.Modulus ();
       L2 = V2.Modulus ();
-      Standard_Real angg = V1.CrossMagnitude(V2)/(L1*L2);
-      if (L1 > LTol && L2 > LTol) ACoef = angg/AngleMax;
-      else                        ACoef = 0.0;
+      if (L1 > LTol && L2 > LTol)
+      {
+        Standard_Real angg = V1.CrossMagnitude(V2) / (L1 * L2);
+        ACoef = angg / AngleMax;
+      }
+      else
+        ACoef = 0.0;
 
-      if (ACoef >= FCoef) Coef = ACoef;      //On retient le plus penalisant
-      else                Coef = FCoef;
+      //On retient le plus penalisant
+      Coef = Max(ACoef, FCoef);
 
+      if (isNeedToCheck && Coef < 0.55)
+      {
+        isNeedToCheck = Standard_False;
+        Du = Dusave;
+        U2 = U1 + Du;
+        if (U2 > lastu)
+          U2 = lastu;
+        D0 (C, U2, CurrentPoint);
+        continue;
+      }
 
       if (Coef <= 1.0) {
         if (Abs (lastu-U2) < uTol) {
@@ -398,39 +503,35 @@ void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
           Correction = Standard_False;
         }
         else {
-          if (Coef >= 0.55 || TooLarge) { 
+          if (Coef >= 0.55 || TooLarge || nbcorr == MaxNbCorr) { 
             parameters.Append (U2);
             points    .Append (CurrentPoint);
             aPrevPoint = CurrentPoint;
             Correction = Standard_False;
-          }
-          else if (TooSmall) {
-            Correction = Standard_False;
-            aPrevPoint = CurrentPoint;
+            isNeedToCheck = Standard_True;
           }
           else {
-            TooSmall = Standard_True;
             //Standard_Real UUU2 = U2;
             Du += Min((U2-U1)*(1.-Coef), Du*Us3);
 
             U2 = U1 + Du;
-            //if (U2 >= lastu) U2 = UUU2;
-            if (U2 >= lastu) {
-              parameters.Append (lastu);
-              points    .Append (LastPoint);
-              MorePoints = Standard_False;
-              Correction = Standard_False;
-            }
-            else D0 (C, U2, CurrentPoint);
+            if (U2 > lastu)
+              U2 = lastu;
+            D0 (C, U2, CurrentPoint);
           }
         }
       }
       else {
 
         if (Coef >= 1.5) {
-          U2 = MiddleU;
+          if (!aPrevPoint.IsEqual(points.Last(), Precision::Confusion()))
+          {
+            parameters.Append (U1);
+            points    .Append (aPrevPoint);
+          }
+          U2 = 0.9*MiddleU + 0.1*U2;
           Du  = U2-U1;
-          CurrentPoint = MiddlePoint;
+          D0 (C, U2, CurrentPoint);
         }
         else {
           Du*=0.9;
@@ -485,7 +586,6 @@ void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
 //    points.Remove (i+1);
 //    i--;
 //  }
-  
   if (i >= 2) {
     MiddleU = parameters (i-1);
     MiddleU = (lastu + MiddleU)*0.5;
@@ -513,4 +613,39 @@ void GCPnts_TangentialDeflection::PerformCurve (const TheCurve& C)
       i++;
     }
   }
+  //Additional check for intervals
+  Standard_Real LTol2 = LTol * LTol;
+  Standard_Integer MaxNbp = 10 * Nbp;
+  for(i = 1; i < Nbp; ++i)
+  {
+    U1 = parameters(i);
+    U2 = parameters(i + 1);
+    // Check maximal deflection on interval;
+    Standard_Real dmax = 0.;
+    Standard_Real umax = 0.;
+    Standard_Real amax = 0.;
+    EstimDefl(C, U1, U2, uTol, dmax, umax);
+    const gp_Pnt& P1 = points(i);
+    const gp_Pnt& P2 = points(i+1);
+    D0(C, umax, MiddlePoint);
+    amax = EstimAngl(P1, MiddlePoint, P2);
+    if(dmax > curvatureDeflection || amax > AngleMax)
+    {
+      if(umax - U1 > uTol && U2 - umax > uTol)
+      {
+        if (P1.SquareDistance(MiddlePoint) > LTol2 && P2.SquareDistance(MiddlePoint) > LTol2)
+        {
+          parameters.InsertAfter(i, umax);
+          points.InsertAfter(i, MiddlePoint);
+          ++Nbp;
+          --i; //To compensate ++i in loop header: i must point to first part of splitted interval
+          if(Nbp > MaxNbp)
+          {
+            break;
+          }
+        }
+      }
+    }
+  }
+  // 
 }

src/GeometryTest/GeometryTest_CurveCommands.cxx

@@ -925,7 +925,7 @@ static Standard_Integer crvpoints (Draw_Interpretor& di, Standard_Integer /*n*/,
 
   //check deviation
   ComputeDeviation(C,aPnts,dmax,ufmax,ulmax,imax);
-  di << "Max defl: " << dmax << " " << ufmax << " " << ulmax << " " << i << "\n"; 
+  di << "Max defl: " << dmax << " " << ufmax << " " << ulmax << " " << imax << "\n"; 
 
   return 0;
 } 
@@ -980,7 +980,7 @@ static Standard_Integer crvtpoints (Draw_Interpretor& di, Standard_Integer n, co
 
   //check deviation
   ComputeDeviation(C,aPnts,dmax,ufmax,ulmax,imax);
-  di << "Max defl: " << dmax << " " << ufmax << " " << ulmax << " " << i << "\n"; 
+  di << "Max defl: " << dmax << " " << ufmax << " " << ulmax << " " << imax << "\n"; 
 
   return 0;
 } 

src/IntCurvesFace/IntCurvesFace_Intersector.cxx

@@ -16,42 +16,34 @@
 
 #define OPTIMISATION 1 
 
-
 #include <IntCurvesFace_Intersector.ixx>
-
-#include <IntCurveSurface_ThePolyhedronToolOfHInter.hxx>
-#include <Bnd_BoundSortBox.hxx>
-
-#include <IntCurveSurface_IntersectionPoint.hxx>
-#include <gp_Lin.hxx>
-#include <TopoDS_Face.hxx>
-#include <TopAbs.hxx>
-
-
-#include <IntCurveSurface_HInter.hxx>
-#include <BRepAdaptor_HSurface.hxx>
-#include <Geom_Line.hxx>
-#include <gp_Pnt2d.hxx>
-#include <BRepClass_FaceClassifier.hxx>
-
-#include <GeomAdaptor_Curve.hxx>
-
-#include <GeomAdaptor_HCurve.hxx>
-#include <BRepAdaptor_HSurface.hxx>
-
-
-
-#include <Adaptor3d_HSurfaceTool.hxx>
-#include <IntCurveSurface_TheHCurveTool.hxx>
 #include <Adaptor3d_HCurve.hxx>
+#include <Adaptor3d_HSurfaceTool.hxx>
+#include <Bnd_BoundSortBox.hxx>
 #include <Bnd_Box.hxx>
-#include <Intf_Tool.hxx>
-#include <IntCurveSurface_ThePolyhedronOfHInter.hxx>
-#include <IntCurveSurface_ThePolygonOfHInter.hxx>
+#include <BRepAdaptor_HSurface.hxx>
+#include <BRepClass_FaceClassifier.hxx>
+#include <BRepTopAdaptor_TopolTool.hxx>
+#include <Geom_Line.hxx>
+#include <GeomAdaptor_Curve.hxx>
+#include <GeomAdaptor_HCurve.hxx>
+#include <gp_Lin.hxx>
+#include <gp_Pnt.hxx>
+#include <gp_Pnt2d.hxx>
+#include <IntCurvesFace_Intersector.hxx>
+#include <IntCurveSurface_HInter.hxx>
+#include <IntCurveSurface_IntersectionPoint.hxx>
 #include <IntCurveSurface_SequenceOfPnt.hxx>
-
-
-
+#include <IntCurveSurface_TheHCurveTool.hxx>
+#include <IntCurveSurface_ThePolygonOfHInter.hxx>
+#include <IntCurveSurface_ThePolyhedronOfHInter.hxx>
+#include <IntCurveSurface_ThePolyhedronToolOfHInter.hxx>
+#include <Intf_Tool.hxx>
+#include <TopAbs.hxx>
+#include <TopoDS_Face.hxx>
+#include <BRep_Tool.hxx>
+#include <TopoDS.hxx>
+#include <Geom2dAPI_ProjectPointOnCurve.hxx>
 //=======================================================================
 //function : SurfaceType
 //purpose  : 
@@ -65,7 +57,7 @@ GeomAbs_SurfaceType IntCurvesFace_Intersector::SurfaceType() const
 //purpose  : 
 //=======================================================================
 IntCurvesFace_Intersector::IntCurvesFace_Intersector(const TopoDS_Face& Face,
-						     const Standard_Real aTol)
+                                                     const Standard_Real aTol)
 : 
   Tol(aTol),
   done(Standard_False),
@@ -74,8 +66,9 @@ IntCurvesFace_Intersector::IntCurvesFace_Intersector(const TopoDS_Face& Face,
   PtrOnBndBounding(NULL)
 { 
   BRepAdaptor_Surface surface;
+  const Standard_Boolean aRestr = Standard_True;
   face = Face;
-  surface.Initialize(Face,Standard_True);
+  surface.Initialize(Face, aRestr);
   Hsurface = new BRepAdaptor_HSurface(surface);
   myTopolTool = new BRepTopAdaptor_TopolTool(Hsurface);
   
@@ -91,47 +84,35 @@ IntCurvesFace_Intersector::IntCurvesFace_Intersector(const TopoDS_Face& Face,
     U1 = Hsurface->LastUParameter();
     V0 = Hsurface->FirstVParameter();
     V1 = Hsurface->LastVParameter();
-    //modified by NIZNHY-PKV Fri Apr 06 07:30:47 2012f
-    Standard_Boolean bFlag;
-    //
-    {
-      Standard_Real dU, dV, dA, dB, aTresh; 
-      bFlag=Standard_True;
-      //
-      aTresh=100.;
-      dU=U1-U0;
-      dV=V1-V0;
-      dA=dU;
-      dB=dV;
-      if (dV>dU) {
-	dA=dV;
-	dB=dU;
-      }
-      //
-      if (dB < Precision::PConfusion() || dA > dB * aTresh) {
-	bFlag=!bFlag;
-      }
-    }
-    //
-    if (bFlag) {
-      nbsu = myTopolTool->NbSamplesU();
-      nbsv = myTopolTool->NbSamplesV();
-      if(nbsu>40) nbsu = 40;
-      if(nbsv>40) nbsv = 40;
-      PtrOnPolyhedron = (IntCurveSurface_ThePolyhedronOfHInter *) 
-	new IntCurveSurface_ThePolyhedronOfHInter(Hsurface,nbsu,nbsv,U0,V0,U1,V1);
-    }
-    //
-    /*
+
+    Standard_Real aURes = Hsurface->UResolution(1.0);
+    Standard_Real aVRes = Hsurface->VResolution(1.0);
+
+    // Checking correlation between number of samples and length of the face along each axis
+    const Standard_Real aTresh = 100.0;
+    const Standard_Integer aMinSamples = 10;
+    const Standard_Integer aMaxSamples = 40;
+    const Standard_Integer aMaxSamples2 = aMaxSamples * aMaxSamples;
+    Standard_Real dU = (U1 - U0) / aURes;
+    Standard_Real dV = (V1 - V0) / aVRes;
     nbsu = myTopolTool->NbSamplesU();
     nbsv = myTopolTool->NbSamplesV();
-    if(nbsu>40) nbsu = 40;
-    if(nbsv>40) nbsv = 40;
-    PtrOnPolyhedron = (IntCurveSurface_ThePolyhedronOfHInter *) 
-      new IntCurveSurface_ThePolyhedronOfHInter(Hsurface,nbsu,nbsv,U0,V0,U1,V1);
-      */
-    
-    //modified by NIZNHY-PKV Fri Apr 06 07:30:49 2012t
+    if (nbsu > aMaxSamples) nbsu = aMaxSamples;
+    if (nbsv > aMaxSamples) nbsv = aMaxSamples;
+
+    if (Max(dU, dV) > Min(dU, dV) * aTresh)
+    {
+      nbsu = (Standard_Integer)(Sqrt(dU / dV) * aMaxSamples);
+      if (nbsu < aMinSamples) nbsu = aMinSamples;
+      nbsv = aMaxSamples2 / nbsu;
+      if (nbsv < aMinSamples)
+      {
+        nbsv = aMinSamples;
+        nbsu = aMaxSamples2 / aMinSamples;
+      }
+    }
+    PtrOnPolyhedron = (IntCurveSurface_ThePolyhedronOfHInter *)
+        new IntCurveSurface_ThePolyhedronOfHInter(Hsurface,nbsu,nbsv,U0,V0,U1,V1);
   }
 }
 //=======================================================================
@@ -142,65 +123,125 @@ void IntCurvesFace_Intersector::InternalCall(const IntCurveSurface_HInter &HICS,
 					     const Standard_Real parinf,
 					     const Standard_Real parsup) 
 {
-  if(HICS.IsDone()) {
+  if(HICS.IsDone() && HICS.NbPoints() > 0) {
+    //Calculate tolerance for 2d classifier
+    Standard_Real mintol3d = BRep_Tool::Tolerance(face);
+    Standard_Real maxtol3d = mintol3d;
+    Standard_Real mintol2d = Tol, maxtol2d = Tol;
+    TopExp_Explorer anExp(face, TopAbs_EDGE);
+    for(; anExp.More(); anExp.Next())
+    {
+      Standard_Real curtol = BRep_Tool::Tolerance(TopoDS::Edge(anExp.Current()));
+      mintol3d = Min(mintol3d, curtol);
+      maxtol3d = Max(maxtol3d, curtol);
+    }
+    Standard_Real minres = Max(Hsurface->UResolution(mintol3d), Hsurface->VResolution(mintol3d));
+    Standard_Real maxres = Max(Hsurface->UResolution(maxtol3d), Hsurface->VResolution(maxtol3d));
+    mintol2d = Max(minres, Tol); 
+    maxtol2d = Max(maxres, Tol);
+    //
+    Handle(BRepTopAdaptor_TopolTool) anAdditionalTool;
     for(Standard_Integer index=HICS.NbPoints(); index>=1; index--) {  
       const IntCurveSurface_IntersectionPoint& HICSPointindex = HICS.Point(index);
       gp_Pnt2d Puv(HICSPointindex.U(),HICSPointindex.V());
-      
-      TopAbs_State currentstate = myTopolTool->Classify(Puv,Tol);
+
+      //TopAbs_State currentstate = myTopolTool->Classify(Puv,Tol);
+      TopAbs_State currentstate = myTopolTool->Classify(Puv, mintol2d);
+      if(currentstate == TopAbs_OUT && maxtol2d > mintol2d) {
+        if(anAdditionalTool.IsNull())
+        {
+          anAdditionalTool = new BRepTopAdaptor_TopolTool(Hsurface);
+        }
+        currentstate = anAdditionalTool->Classify(Puv,maxtol2d);
+        if(currentstate == TopAbs_ON)
+        {
+          //Find out nearest edge and it's tolerance
+          anExp.Init(face, TopAbs_EDGE);
+          Standard_Real mindist = RealLast();
+          mintol2d = -1.;
+          for(; anExp.More(); anExp.Next())
+          {
+            TopoDS_Edge anE = TopoDS::Edge(anExp.Current());
+            Standard_Real curtol = BRep_Tool::Tolerance(anE);
+            Standard_Real tol2d = Max(Hsurface->UResolution(curtol), Hsurface->VResolution(curtol));
+            tol2d = Max(tol2d, Tol);
+            Standard_Real f, l;
+            Handle(Geom2d_Curve) aPC = BRep_Tool::CurveOnSurface(anE, face, f, l);
+            Geom2dAPI_ProjectPointOnCurve aProj(Puv, aPC, f, l);
+            if(aProj.NbPoints() > 0)
+            {
+              Standard_Real d = aProj.LowerDistance();
+              if(d <= tol2d)
+              {
+                //Nearest edge is found, state is really ON
+                break;
+              }
+              if(d < mindist && d <= maxtol2d)
+              {
+                mindist = d;
+                mintol2d = tol2d;
+              }
+            }
+          }
+          if(mintol2d > 0. && mindist > mintol2d)
+          {
+            currentstate = TopAbs_OUT; 
+          }
+        }
+      }
       if(currentstate==TopAbs_IN || currentstate==TopAbs_ON) { 
-	Standard_Real HICSW = HICSPointindex.W();
-	if(HICSW >= parinf && HICSW <= parsup ) { 
-	  Standard_Real U          = HICSPointindex.U();
-	  Standard_Real V          = HICSPointindex.V();
-	  Standard_Real W          = HICSW; 
-	  IntCurveSurface_TransitionOnCurve transition = HICSPointindex.Transition();
-	  gp_Pnt pnt        = HICSPointindex.Pnt();
-	  //	  state      = currentstate;
-	  //  Modified by skv - Wed Sep  3 16:14:10 2003 OCC578 Begin
-	  Standard_Integer anIntState = (currentstate == TopAbs_IN) ? 0 : 1;
-	  //  Modified by skv - Wed Sep  3 16:14:11 2003 OCC578 End
+        Standard_Real HICSW = HICSPointindex.W();
+        if(HICSW >= parinf && HICSW <= parsup ) { 
+          Standard_Real U          = HICSPointindex.U();
+          Standard_Real V          = HICSPointindex.V();
+          Standard_Real W          = HICSW; 
+          IntCurveSurface_TransitionOnCurve transition = HICSPointindex.Transition();
+          gp_Pnt pnt        = HICSPointindex.Pnt();
+          //	  state      = currentstate;
+          //  Modified by skv - Wed Sep  3 16:14:10 2003 OCC578 Begin
+          Standard_Integer anIntState = (currentstate == TopAbs_IN) ? 0 : 1;
+          //  Modified by skv - Wed Sep  3 16:14:11 2003 OCC578 End
 
-	  if(transition != IntCurveSurface_Tangent && face.Orientation()==TopAbs_REVERSED) { 
-	    if(transition == IntCurveSurface_In) 
-	      transition = IntCurveSurface_Out;
-	    else 
-	      transition = IntCurveSurface_In;
-	  }
-	  //----- Insertion du point 
-	  if(nbpnt==0) { 
-	    IntCurveSurface_IntersectionPoint PPP(pnt,U,V,W,transition);
-	    SeqPnt.Append(PPP);
-	    //  Modified by skv - Wed Sep  3 16:14:10 2003 OCC578 Begin
-	    mySeqState.Append(anIntState);
-	    //  Modified by skv - Wed Sep  3 16:14:11 2003 OCC578 End
-	  }
-	  else { 
-	    Standard_Integer i = 1;
-	    Standard_Integer b = nbpnt+1;                    
-	    while(i<=nbpnt) {
-	      const IntCurveSurface_IntersectionPoint& Pnti=SeqPnt.Value(i);
-	      Standard_Real wi = Pnti.W();
-	      if(wi >= W) { b=i; i=nbpnt; }
-	      i++;
-	    }
-	    IntCurveSurface_IntersectionPoint PPP(pnt,U,V,W,transition);
-	    //  Modified by skv - Wed Sep  3 16:14:10 2003 OCC578 Begin
-// 	    if(b>nbpnt)          { SeqPnt.Append(PPP); } 
-// 	    else if(b>0)             { SeqPnt.InsertBefore(b,PPP); } 
-	    if(b>nbpnt) {
-	      SeqPnt.Append(PPP);
-	      mySeqState.Append(anIntState);
-	    } else if(b>0) {
-	      SeqPnt.InsertBefore(b,PPP);
-	      mySeqState.InsertBefore(b, anIntState);
-	    }
-	    //  Modified by skv - Wed Sep  3 16:14:11 2003 OCC578 End
-	  }
+          if(transition != IntCurveSurface_Tangent && face.Orientation()==TopAbs_REVERSED) { 
+            if(transition == IntCurveSurface_In) 
+              transition = IntCurveSurface_Out;
+            else 
+              transition = IntCurveSurface_In;
+          }
+          //----- Insertion du point 
+          if(nbpnt==0) { 
+            IntCurveSurface_IntersectionPoint PPP(pnt,U,V,W,transition);
+            SeqPnt.Append(PPP);
+            //  Modified by skv - Wed Sep  3 16:14:10 2003 OCC578 Begin
+            mySeqState.Append(anIntState);
+            //  Modified by skv - Wed Sep  3 16:14:11 2003 OCC578 End
+          }
+          else { 
+            Standard_Integer i = 1;
+            Standard_Integer b = nbpnt+1;                    
+            while(i<=nbpnt) {
+              const IntCurveSurface_IntersectionPoint& Pnti=SeqPnt.Value(i);
+              Standard_Real wi = Pnti.W();
+              if(wi >= W) { b=i; i=nbpnt; }
+              i++;
+            }
+            IntCurveSurface_IntersectionPoint PPP(pnt,U,V,W,transition);
+            //  Modified by skv - Wed Sep  3 16:14:10 2003 OCC578 Begin
+            // 	    if(b>nbpnt)          { SeqPnt.Append(PPP); } 
+            // 	    else if(b>0)             { SeqPnt.InsertBefore(b,PPP); } 
+            if(b>nbpnt) {
+              SeqPnt.Append(PPP);
+              mySeqState.Append(anIntState);
+            } else if(b>0) {
+              SeqPnt.InsertBefore(b,PPP);
+              mySeqState.InsertBefore(b, anIntState);
+            }
+            //  Modified by skv - Wed Sep  3 16:14:11 2003 OCC578 End
+          }
 
- 
-	  nbpnt++;
-	} 
+
+          nbpnt++;
+        } 
       } //-- classifier state is IN or ON
     } //-- Loop on Intersection points.
   } //-- HICS.IsDone()

src/IntTools/IntTools_EdgeEdge.cdl

@@ -125,27 +125,39 @@ is
     ---Purpose:
     -- Computes Line/Line intersection.  
 
+    FindSolutions(me:out;  
+        theRanges1 : out SequenceOfRanges from IntTools;
+        theRanges2 : out SequenceOfRanges from IntTools;
+        bSplit2    : out Boolean from Standard)
+    is protected;
+    ---Purpose:
+    -- Intermediate function
+ 
     FindSolutions(me:out;  
         theR1, theR2 : Range from IntTools;
+        theBox2      : Box from Bnd;
         theRanges1   : out SequenceOfRanges from IntTools;
         theRanges2   : out SequenceOfRanges from IntTools)
-    is protected; 
+    is protected;
     ---Purpose:
     -- Looking for the exact intersection ranges
     
     MergeSolutions(me:out; 
-        theRanges1, theRanges2 : SequenceOfRanges from IntTools)
+        theRanges1 : SequenceOfRanges from IntTools;
+        theRanges2 : SequenceOfRanges from IntTools;
+        bSplit2    : Boolean from Standard)
     is protected; 
     ---Purpose:
     -- Merges found solutions
 
     FindParameters(myclass; 
-        theBAC  : Curve from BRepAdaptor;  
-        aT1,aT2 : Real from Standard;  
-        theRes  : Real from Standard;
-        thePTol : Real from Standard;
-        theCBox : Box from Bnd; 
-        aTB1,aTB2 : out Real from Standard)
+        theBAC      : Curve from BRepAdaptor;
+        aT1, aT2    : Real  from Standard;  
+        theRes      : Real  from Standard;
+        thePTol     : Real  from Standard; 
+        theResCoeff : Real  from Standard;
+        theCBox     : Box   from Bnd; 
+        aTB1, aTB2  : out Real from Standard)
     returns Boolean from Standard 
     is protected; 
     ---Purpose:
@@ -213,7 +225,10 @@ fields
     myTol    : Real  from Standard is protected;   
     
     myRes1   : Real  from Standard is protected;  
-    myRes2   : Real  from Standard is protected;  
+    myRes2   : Real  from Standard is protected;   
+     
+    myResCoeff1 : Real  from Standard is protected;
+    myResCoeff2 : Real  from Standard is protected;
      
     myPTol1  : Real  from Standard is protected; 
     myPTol2  : Real  from Standard is protected; 

src/IntTools/IntTools_EdgeEdge.cxx

@@ -24,14 +24,19 @@
 #include <Bnd_Box.hxx>
 #include <BndLib_Add3dCurve.hxx>
 
+#include <Geom_Circle.hxx>
+#include <Geom_Ellipse.hxx>
+#include <Geom_BezierCurve.hxx>
+#include <Geom_BSplineCurve.hxx>
+
 #include <GeomAPI_ProjectPointOnCurve.hxx>
 
 #include <BRep_Tool.hxx>
 #include <BRepAdaptor_Curve.hxx>
 
-#include <IntTools_Tools.hxx>
 #include <IntTools_CommonPrt.hxx>
 
+#include <BOPCol_MapOfInteger.hxx>
 
 static 
   void BndBuildBox(const BRepAdaptor_Curve& theBAC,
@@ -43,11 +48,11 @@ static
   Standard_Real PointBoxDistance(const Bnd_Box& aB,
                                  const gp_Pnt& aP);
 static 
-  void SplitRangeOnSegments(const Standard_Real aT1, 
-                            const Standard_Real aT2,
-                            const Standard_Real theResolution,
-                            const Standard_Integer theNbSeg,
-                            IntTools_SequenceOfRanges& theSegments);
+  Standard_Integer SplitRangeOnSegments(const Standard_Real aT1, 
+                                        const Standard_Real aT2,
+                                        const Standard_Real theResolution,
+                                        const Standard_Integer theNbSeg,
+                                        IntTools_SequenceOfRanges& theSegments);
 static
  Standard_Integer DistPC(const Standard_Real aT1, 
                          const Handle(Geom_Curve)& theC1,
@@ -78,6 +83,23 @@ static
                               Standard_Real& aT1max,
                               Standard_Real& aT2max,
                               const Standard_Boolean bMaxDist = Standard_True);
+static
+  Standard_Real ResolutionCoeff(const BRepAdaptor_Curve& theBAC,
+                                const IntTools_Range& theRange);
+static
+  Standard_Real Resolution(const Handle(Geom_Curve)& theCurve,
+                           const GeomAbs_CurveType theCurveType,
+                           const Standard_Real theResCoeff,
+                           const Standard_Real theR3D);
+static
+  Standard_Real CurveDeflection(const BRepAdaptor_Curve& theBAC,
+                                const IntTools_Range& theRange);
+static 
+  Standard_Integer IsClosed(const Handle(Geom_Curve)& theCurve,
+                            const Standard_Real aT1,
+                            const Standard_Real aT2,
+                            const Standard_Real theTol,
+                            const Standard_Real theRes);
 static 
   Standard_Integer TypeToInteger(const GeomAbs_CurveType theCType);
 
@@ -112,36 +134,11 @@ void IntTools_EdgeEdge::Prepare()
   if (iCT1 == iCT2) {
     if (iCT1 != 0) {
       //compute deflection
-      Standard_Integer i;
-      Standard_Real aDt, aT, aT1, aT2;
-      gp_Vec aV1, aV2;
-      gp_Pnt aP;
+      Standard_Real aC1, aC2;
       //
-      Standard_Real aC1(10.), aC2(10.);
-      for (i = 0; i < 2; ++i) {
-        Standard_Real &aC = !i ? aC2 : aC1;
-        aC = 0;
-        IntTools_Range aR = !i ? myRange2 : myRange1;
-        const BRepAdaptor_Curve& aBAC = !i ? myCurve2 : myCurve1;
-        aR.Range(aT1, aT2);
-        aDt = (aT2 - aT1) / 10.;
-        aT = aT1;
-        aBAC.D1(aT, aP, aV1);
-        while (aT < aT2) {
-          aT += aDt;
-          aBAC.D1(aT, aP, aV2);
-          if (aV1.Magnitude() > gp::Resolution() &&
-              aV2.Magnitude() > gp::Resolution()) {
-            gp_Dir aD1(aV1), aD2(aV2);
-            aC += aD1.Angle(aD2);
-          }
-          aV1 = aV2;
-        }
-        //
-        if (aC < Precision::Confusion()) {
-          break;
-        }
-      }
+      aC2 = CurveDeflection(myCurve2, myRange2);
+      aC1 = (aC2 > Precision::Confusion()) ? 
+        CurveDeflection(myCurve1, myRange1) : 1.;
       //
       if (aC1 < aC2) {
         --iCT1;
@@ -169,15 +166,18 @@ void IntTools_EdgeEdge::Prepare()
   myTol2 = myCurve2.Tolerance();
   myTol = myTol1 + myTol2;
   //
-  myRes1 = myCurve1.Resolution(myTol);
-  myRes2 = myCurve2.Resolution(myTol);
-  //
   if (iCT1 != 0 || iCT2 != 0) {
     Standard_Real f, l, aTM;
     //
     myGeom1 = BRep_Tool::Curve(myEdge1, f, l);
     myGeom2 = BRep_Tool::Curve(myEdge2, f, l);
     //
+    myResCoeff1 = ResolutionCoeff(myCurve1, myRange1);
+    myResCoeff2 = ResolutionCoeff(myCurve2, myRange2);
+    //
+    myRes1 = Resolution(myCurve1.Curve().Curve(), myCurve1.GetType(), myResCoeff1, myTol1);
+    myRes2 = Resolution(myCurve2.Curve().Curve(), myCurve2.GetType(), myResCoeff2, myTol2);
+    //
     myPTol1 = 5.e-13;
     aTM = Max(fabs(myRange1.First()), fabs(myRange1.Last()));
     if (aTM > 999.) {
@@ -217,10 +217,70 @@ void IntTools_EdgeEdge::Perform()
   IntTools_SequenceOfRanges aRanges1, aRanges2;
   //
   //3.2. Find ranges containig solutions
-  FindSolutions(myRange1, myRange2, aRanges1, aRanges2);
+  Standard_Boolean bSplit2;
+  FindSolutions(aRanges1, aRanges2, bSplit2);
   //
   //4. Merge solutions and save common parts
-  MergeSolutions(aRanges1, aRanges2);
+  MergeSolutions(aRanges1, aRanges2, bSplit2);
+}
+
+//=======================================================================
+//function : FindSolutions
+//purpose  : 
+//=======================================================================
+void IntTools_EdgeEdge::FindSolutions(IntTools_SequenceOfRanges& theRanges1,
+                                      IntTools_SequenceOfRanges& theRanges2,
+                                      Standard_Boolean& bSplit2)
+{
+  Standard_Boolean bIsClosed2;
+  Standard_Real aT11, aT12, aT21, aT22;
+  Bnd_Box aB2;
+  //
+  bSplit2 = Standard_False;
+  myRange1.Range(aT11, aT12);
+  myRange2.Range(aT21, aT22);
+  //
+  bIsClosed2 = IsClosed(myGeom2, aT21, aT22, myTol2, myRes2);
+  //
+  if (bIsClosed2) {
+    Bnd_Box aB1;
+    BndBuildBox(myCurve1, aT11, aT12, myTol1, aB1);
+    //
+    gp_Pnt aP = myGeom2->Value(aT21);
+    bIsClosed2 = !aB1.IsOut(aP);
+  }
+  //
+  if (!bIsClosed2) {
+    BndBuildBox(myCurve2, aT21, aT22, myTol2, aB2);
+    FindSolutions(myRange1, myRange2, aB2, theRanges1, theRanges2);
+    return;
+  }
+  //
+  if (!CheckCoincidence(aT11, aT12, aT21, aT22, myTol, myRes1)) {
+    theRanges1.Append(myRange1);
+    theRanges2.Append(myRange2);
+    return;
+  }
+  //
+  Standard_Integer i, j, aNb1, aNb2;
+  IntTools_SequenceOfRanges aSegments1, aSegments2;
+  //
+  aNb1 = IsClosed(myGeom1, aT11, aT12, myTol1, myRes1) ? 2 : 1;
+  aNb2 = 2;
+  //
+  aNb1 = SplitRangeOnSegments(aT11, aT12, myRes1, aNb1, aSegments1);
+  aNb2 = SplitRangeOnSegments(aT21, aT22, myRes2, aNb2, aSegments2);
+  //
+  for (i = 1; i <= aNb1; ++i) {
+    const IntTools_Range& aR1 = aSegments1(i);
+    for (j = 1; j <= aNb2; ++j) {
+      const IntTools_Range& aR2 = aSegments2(j);
+      BndBuildBox(myCurve2, aR2.First(), aR2.Last(), myTol2, aB2);
+      FindSolutions(aR1, aR2, aB2, theRanges1, theRanges2);
+    }
+  }
+  //
+  bSplit2 = aNb2 > 1;
 }
 
 //=======================================================================
@@ -229,82 +289,84 @@ void IntTools_EdgeEdge::Perform()
 //=======================================================================
 void IntTools_EdgeEdge::FindSolutions(const IntTools_Range& theR1,
                                       const IntTools_Range& theR2,
-                                      IntTools_SequenceOfRanges& theRanges1, 
+                                      const Bnd_Box& theBox2,
+                                      IntTools_SequenceOfRanges& theRanges1,
                                       IntTools_SequenceOfRanges& theRanges2)
 {
   Standard_Boolean bOut, bStop, bThin;
   Standard_Real aT11, aT12, aT21, aT22;
   Standard_Real aTB11, aTB12, aTB21, aTB22;
-  Standard_Real aTol, aSmallStep1, aSmallStep2;
+  Standard_Real aSmallStep1, aSmallStep2;
   Standard_Integer iCom;
   Bnd_Box aB1, aB2;
   //
   theR1.Range(aT11, aT12);
   theR2.Range(aT21, aT22);
   //
-  BndBuildBox(myCurve1, aT11, aT12, myTol1, aB1);
-  BndBuildBox(myCurve2, aT21, aT22, myTol2, aB2);
-  if (aB1.IsOut(aB2)) {
-    //no intersection;
-    return;
-  }
+  aB2 = theBox2;
   //
-  bOut  = Standard_False;
   bThin = Standard_False;
   bStop = Standard_False;
-  aTol  = 2*myTol;
   iCom  = 1;
   //
   do {
-    bOut = !FindParameters(myCurve2, aT21, aT22, myRes2, myPTol2, aB1, aTB21, aTB22);
+    aTB11 = aT11;
+    aTB12 = aT12;
+    aTB21 = aT21;
+    aTB22 = aT22;
+    //
+    //1. Build box for first edge and find parameters 
+    //   of the second one in that box
+    BndBuildBox(myCurve1, aT11, aT12, myTol1, aB1);
+    bOut = aB1.IsOut(aB2);
     if (bOut) {
       break;
     }
     //
-    bThin = (aTB22 - aTB21) < myRes2;
-    if (bThin) {
-      bOut = !FindParameters(myCurve1, aT11, aT12, myRes1, myPTol1, aB1, aTB11, aTB12);
+    bThin = ((aT12 - aT11) < myRes1) ||
+      (aB1.IsXThin(myTol) && aB1.IsYThin(myTol) && aB1.IsZThin(myTol));
+    //
+    bOut = !FindParameters(myCurve2, aTB21, aTB22, myRes2, myPTol2, 
+                           myResCoeff2, aB1, aT21, aT22);
+    if (bOut || bThin) {
       break;
     }
     //
-    BndBuildBox(myCurve2, aTB21, aTB22, myTol2, aB2);
-    //
-    bOut = !FindParameters(myCurve1, aT11, aT12, myRes1, myPTol1, aB2, aTB11, aTB12);
+    //2. Build box for second edge and find parameters 
+    //   of the first one in that box
+    BndBuildBox(myCurve2, aT21, aT22, myTol2, aB2);
+    bOut = aB1.IsOut(aB2);
     if (bOut) {
       break;
     }
     //
-    bThin = ((aTB12 - aTB11) < myRes1) ||
-      (aB2.IsXThin(aTol) && aB2.IsYThin(aTol) && aB2.IsZThin(aTol));
+    bThin = ((aT22 - aT21) < myRes2) ||
+      (aB2.IsXThin(myTol) && aB2.IsYThin(myTol) && aB2.IsZThin(myTol));
     //
-    if (!bThin) {
-      aSmallStep1 = (aT12 - aT11) / 250.;
-      aSmallStep2 = (aT22 - aT21) / 250.;
-      //
-      if (aSmallStep1 < myRes1) {
-        aSmallStep1 = myRes1;
-      }
-      if (aSmallStep2 < myRes2) {
-        aSmallStep2 = myRes2;
-      }
-      //
-      if (((aTB11 - aT11) < aSmallStep1) && ((aT12 - aTB12) < aSmallStep1) &&
-          ((aTB21 - aT21) < aSmallStep2) && ((aT22 - aTB22) < aSmallStep2)) {
-        bStop = Standard_True;
-      } else {
-        BndBuildBox(myCurve1, aTB11, aTB12, myTol1, aB1);
-        bOut = aB1.IsOut(aB2);
-        if (bOut) {
-          break;
-        }
-      }
+    bOut = !FindParameters(myCurve1, aTB11, aTB12, myRes1, myPTol1, 
+                           myResCoeff1, aB2, aT11, aT12);
+    //
+    if (bOut || bThin) {
+      break;
     }
     //
-    aT11 = aTB11;
-    aT12 = aTB12;
-    aT21 = aTB21;
-    aT22 = aTB22;
-  } while (!bThin && !bStop);
+    //3. Check if it makes sense to continue
+    aSmallStep1 = (aTB12 - aTB11) / 250.;
+    aSmallStep2 = (aTB22 - aTB21) / 250.;
+    //
+    if (aSmallStep1 < myRes1) {
+      aSmallStep1 = myRes1;
+    }
+    if (aSmallStep2 < myRes2) {
+      aSmallStep2 = myRes2;
+    }
+    //
+    if (((aT11 - aTB11) < aSmallStep1) && ((aTB12 - aT12) < aSmallStep1) &&
+        ((aT21 - aTB21) < aSmallStep2) && ((aTB22 - aT22) < aSmallStep2)) {
+      bStop = Standard_True;
+    }
+    //
+  } while (!bStop);
   //
   if (bOut) {
     //no intersection;
@@ -316,23 +378,37 @@ void IntTools_EdgeEdge::FindSolutions(const IntTools_Range& theR1,
     iCom = CheckCoincidence(aT11, aT12, aT21, aT22, myTol, myRes1);
     if (!iCom) {
       bThin = Standard_True;
-    } 
+    }
   }
   //
   if (bThin) {
     if (iCom != 0) {
       //check intermediate points
-      Standard_Real aT1, aT2, aDist;
-      gp_Pnt aP1, aP2;
-      //
-      aT1 = IntTools_Tools::IntermediatePoint(aT11, aT12);
-      aT2 = IntTools_Tools::IntermediatePoint(aT21, aT22);
+      Standard_Boolean bSol;
+      Standard_Real aT1;
+      gp_Pnt aP1;
+      GeomAPI_ProjectPointOnCurve aProjPC;
       //
+      aT1 = (aT11 + aT12) * .5;
       myGeom1->D0(aT1, aP1);
-      myGeom2->D0(aT2, aP2);
       //
-      aDist = aP1.Distance(aP2);
-      if (aDist > myTol) {
+      aProjPC.Init(myGeom2, aT21, aT22);
+      aProjPC.Perform(aP1);
+      //
+      if (aProjPC.NbPoints()) {
+        bSol = aProjPC.LowerDistance() <= myTol;
+      }
+      else {
+        Standard_Real aT2;
+        gp_Pnt aP2;
+        //
+        aT2 = (aT21 + aT22) * .5;
+        myGeom2->D0(aT2, aP2);
+        //
+        bSol = aP1.IsEqual(aP2, myTol);
+      }
+      //
+      if (!bSol) {
         return;
       }
     }
@@ -353,12 +429,12 @@ void IntTools_EdgeEdge::FindSolutions(const IntTools_Range& theR1,
   IntTools_SequenceOfRanges aSegments1;
   //
   IntTools_Range aR2(aT21, aT22);
+  BndBuildBox(myCurve2, aT21, aT22, myTol2, aB2);
   //
-  SplitRangeOnSegments(aT11, aT12, myRes1, 3, aSegments1);
-  aNb1 = aSegments1.Length();
+  aNb1 = SplitRangeOnSegments(aT11, aT12, myRes1, 3, aSegments1);
   for (i = 1; i <= aNb1; ++i) {
     const IntTools_Range& aR1 = aSegments1(i);
-    FindSolutions(aR1, aR2, theRanges1, theRanges2);
+    FindSolutions(aR1, aR2, aB2, theRanges1, theRanges2);
   }
 }
 
@@ -371,6 +447,7 @@ Standard_Boolean IntTools_EdgeEdge::FindParameters(const BRepAdaptor_Curve& theB
                                                    const Standard_Real aT2, 
                                                    const Standard_Real theRes,
                                                    const Standard_Real thePTol,
+                                                   const Standard_Real theResCoeff,
                                                    const Bnd_Box& theCBox,
                                                    Standard_Real& aTB1, 
                                                    Standard_Real& aTB2)
@@ -391,6 +468,9 @@ Standard_Boolean IntTools_EdgeEdge::FindParameters(const BRepAdaptor_Curve& theB
   aCBx = theCBox;
   aCBx.Enlarge(aTol);
   //
+  const Handle(Geom_Curve)& aCurve = theBAC.Curve().Curve();
+  const GeomAbs_CurveType aCurveType = theBAC.GetType();
+  //
   for (i = 0; i < 2; ++i) {
     aTB = !i ? aT1 : aT2;
     aT = !i ? aT2 : aTB1;
@@ -403,10 +483,10 @@ Standard_Boolean IntTools_EdgeEdge::FindParameters(const BRepAdaptor_Curve& theB
       aDist = PointBoxDistance(theCBox, aP);
       if (aDist > aTol) {
         if (fabs(aDist - aDistP) < aDistTol) {
-          aDt = theBAC.Resolution((++k)*aDist);
+          aDt = Resolution(aCurve, aCurveType, theResCoeff, (++k)*aDist);
         } else {
           k = 0;
-          aDt = theBAC.Resolution(aDist);
+          aDt = Resolution(aCurve, aCurveType, theResCoeff, aDist);
         }
         aTB += aC*aDt;
       } else {
@@ -458,48 +538,94 @@ Standard_Boolean IntTools_EdgeEdge::FindParameters(const BRepAdaptor_Curve& theB
 //purpose  : 
 //=======================================================================
 void IntTools_EdgeEdge::MergeSolutions(const IntTools_SequenceOfRanges& theRanges1, 
-                                       const IntTools_SequenceOfRanges& theRanges2)
+                                       const IntTools_SequenceOfRanges& theRanges2,
+                                       const Standard_Boolean bSplit2)
 {
-  IntTools_Range aRi1, aRi2, aRj1, aRj2;
-  Standard_Integer aNbCP, i, j;
-  TopAbs_ShapeEnum aType;
-  Standard_Real aTi11, aTi12, aTi21, aTi22,
-                aTj11, aTj12, aTj21, aTj22;
+  Standard_Integer aNbCP = theRanges1.Length();
+  if (aNbCP == 0) {
+    return;
+  }
   //
-  aNbCP = theRanges1.Length();
+  IntTools_Range aRi1, aRi2, aRj1, aRj2;
+  Standard_Boolean bCond;
+  Standard_Integer i, j;
+  TopAbs_ShapeEnum aType;
+  Standard_Real aT11, aT12, aT21, aT22;
+  Standard_Real aTi11, aTi12, aTi21, aTi22;
+  Standard_Real aTj11, aTj12, aTj21, aTj22;
+  Standard_Real aRes1, aRes2, dTR1, dTR2;
+  BOPCol_MapOfInteger aMI;
+  //
+  aRes1 = Resolution(myCurve1.Curve().Curve(), 
+                     myCurve1.GetType(), myResCoeff1, myTol);
+  aRes2 = Resolution(myCurve2.Curve().Curve(), 
+                     myCurve2.GetType(), myResCoeff2, myTol);
+  //
+  myRange1.Range(aT11, aT12);
+  myRange2.Range(aT21, aT22);
+  dTR1 = 20*aRes1;
+  dTR2 = 20*aRes2;
   aType = TopAbs_VERTEX;
   //
-  for (i = 1; i <= aNbCP; ) {
+  for (i = 1; i <= aNbCP;) {
+    if (aMI.Contains(i)) {
+      ++i;
+      continue;
+    }
+    //
     aRi1 = theRanges1(i);
     aRi2 = theRanges2(i);
     //
     aRi1.Range(aTi11, aTi12);
     aRi2.Range(aTi21, aTi22);
     //
+    aMI.Add(i);
+    //
     for (j = i+1; j <= aNbCP; ++j) {
+      if (aMI.Contains(j)) {
+        continue;
+      }
+      //
       aRj1 = theRanges1(j);
       aRj2 = theRanges2(j);
       //
       aRj1.Range(aTj11, aTj12);
       aRj2.Range(aTj21, aTj22);
-      if (fabs(aTi12 - aTj11) < 10*myRes1 ||
-          fabs(aTi22 - aTj21) < 10*myRes2) {
+      //
+      bCond = (fabs(aTi12 - aTj11) < dTR1) ||
+        (bSplit2 && (fabs(aTj12 - aTi11) < dTR1));
+      if (bCond && bSplit2) {
+        bCond = (fabs((Max(aTi22, aTj22) - Min(aTi21, aTj21)) - 
+                      ((aTi22 - aTi21) + (aTj22 - aTj21))) < dTR2);
+      }
+      //
+      if (bCond) {
         aTi11 = Min(aTi11, aTj11);
         aTi12 = Max(aTi12, aTj12);
         aTi21 = Min(aTi21, aTj21);
         aTi22 = Max(aTi22, aTj22);
-      } else {
+        aMI.Add(j);
+      }
+      else if (!bSplit2) {
+        i = j;
         break;
       }
     }
-    i = j;
     //
-    if (aTi11 == myRange1.First() && aTi12 == myRange1.Last() &&
-        aTi21 == myRange2.First() && aTi22 == myRange2.Last()) {
+    if (((fabs(aT11 - aTi11) < myRes1) && (fabs(aT12 - aTi12) < myRes1)) ||
+        ((fabs(aT21 - aTi21) < myRes2) && (fabs(aT22 - aTi22) < myRes2))) {
       aType = TopAbs_EDGE;
+      myCommonParts.Clear();
     }
     //
     AddSolution(aTi11, aTi12, aTi21, aTi22, aType);
+    if (aType == TopAbs_EDGE) {
+      break;
+    }
+    //
+    if (bSplit2) {
+      ++i;
+    }
   }
 }
 
@@ -556,35 +682,64 @@ void IntTools_EdgeEdge::FindBestSolution(const Standard_Real aT11,
                                          Standard_Real& aT2)
 {
   Standard_Integer i, aNbS, iErr;
-  Standard_Real aDMin, aD, aCrit;
-  Standard_Real aT1x, aT2x, aT1p, aT2p;
-  GeomAPI_ProjectPointOnCurve aProj;
-  IntTools_SequenceOfRanges aSeg1;
+  Standard_Real aDMin, aD, aRes1, aSolCriteria, aTouchCriteria;
+  Standard_Real aT1A, aT1B, aT1Min, aT2Min;
+  Standard_Real aT1Im, aT2Im, aT1Touch;
+  GeomAPI_ProjectPointOnCurve aProjPC;
+  IntTools_SequenceOfRanges aRanges;
+  Standard_Boolean bTouch;
   //
-  aT1 = IntTools_Tools::IntermediatePoint(aT11, aT12);
-  aT2 = IntTools_Tools::IntermediatePoint(aT21, aT22);
-  //
-  aDMin = 100.;
-  aD = 100.;
-  aCrit = 0.;//1.e-16;
+  aDMin = Precision::Infinite();
+  aSolCriteria   = 5.e-16;
+  aTouchCriteria = 5.e-13;
+  bTouch = Standard_False;
+  aT1Touch = aT11;
   //
+  aRes1 = Resolution(myCurve1.Curve().Curve(), 
+                     myCurve1.GetType(), myResCoeff1, myTol);
   aNbS = 10;
-  SplitRangeOnSegments(aT11, aT12, 3*myRes1, aNbS, aSeg1);
-  aNbS = aSeg1.Length();
+  aNbS = SplitRangeOnSegments(aT11, aT12, 3*aRes1, aNbS, aRanges);
+  //
+  aProjPC.Init(myGeom2, aT21, aT22);
+  //
+  aT1 = (aT11 + aT12) * 0.5;
+  iErr = DistPC(aT1, myGeom1, aSolCriteria, aProjPC, aD, aT2, -1);
+  if (iErr == 1) {
+    aT2 = (aT21 + aT22) * 0.5;
+  }
+  //
+  aT1Im = aT1;
+  aT2Im = aT2;
   //
-  aProj.Init(myGeom2, aT21, aT22);
   for (i = 1; i <= aNbS; ++i) {
-    const IntTools_Range& aR1 = aSeg1(i);
-    aR1.Range(aT1x, aT2x);
+    const IntTools_Range& aR1 = aRanges(i);
+    aR1.Range(aT1A, aT1B);
     //
-    iErr = FindDistPC(aT1x, aT2x, myGeom1, aCrit, myPTol1,
-                      aProj, aD, aT1p, aT2p, Standard_False);
-    if (iErr != 1 && aD < aDMin) {
-      aT1 = aT1p;
-      aT2 = aT2p;
-      aDMin = aD;
-      if (aDMin <= aCrit) {
-        break;
+    aD = myTol;
+    iErr = FindDistPC(aT1A, aT1B, myGeom1, aSolCriteria, myPTol1,
+                      aProjPC, aD, aT1Min, aT2Min, Standard_False);
+    if (iErr != 1) {
+      if (aD < aDMin) {
+        aT1 = aT1Min;
+        aT2 = aT2Min;
+        aDMin = aD;
+      }
+      //
+      if (aD < aTouchCriteria) {
+        if (bTouch) {
+          aT1A = (aT1Touch + aT1Min) * 0.5;
+          iErr = DistPC(aT1A, myGeom1, aTouchCriteria, 
+                        aProjPC, aD, aT2Min, -1);
+          if (aD > aTouchCriteria) {
+            aT1 = aT1Im;
+            aT2 = aT2Im;
+            break;
+          }
+        }
+        else {
+          aT1Touch = aT1Min;
+          bTouch = Standard_True;
+        }
       }
     }
   }
@@ -789,13 +944,14 @@ Standard_Boolean IntTools_EdgeEdge::IsIntersection(const Standard_Real aT11,
     //
     if (((anAngle1 < anAngleCriteria) || ((M_PI - anAngle1) < anAngleCriteria)) ||
         ((anAngle2 < anAngleCriteria) || ((M_PI - anAngle2) < anAngleCriteria))) {
-      GeomAPI_ProjectPointOnCurve aProj;
+      GeomAPI_ProjectPointOnCurve aProjPC;
       Standard_Integer iErr;
-      Standard_Real aD, aT1p, aT2p;
+      Standard_Real aD, aT1Min, aT2Min;
       //
-      aD = 100.;
-      aProj.Init(myGeom2, aT21, aT22);
-      iErr = FindDistPC(aT11, aT12, myGeom1, myTol, myRes1, aProj, aD, aT1p, aT2p, Standard_False);
+      aD = Precision::Infinite();
+      aProjPC.Init(myGeom2, aT21, aT22);
+      iErr = FindDistPC(aT11, aT12, myGeom1, myTol, myRes1, 
+                        aProjPC, aD, aT1Min, aT2Min, Standard_False);
       bRet = (iErr == 2);
     }
   }
@@ -816,7 +972,7 @@ Standard_Integer IntTools_EdgeEdge::CheckCoincidence(const Standard_Real aT11,
   Standard_Integer iErr, aNb, aNb1, i;
   Standard_Real aT1A, aT1B, aT1max, aT2max, aDmax;
   GeomAPI_ProjectPointOnCurve aProjPC;
-  IntTools_SequenceOfRanges aSeg1;
+  IntTools_SequenceOfRanges aRanges;
   //
   iErr  = 0;
   aDmax = -1.;
@@ -824,10 +980,9 @@ Standard_Integer IntTools_EdgeEdge::CheckCoincidence(const Standard_Real aT11,
   //
   // 1. Express evaluation
   aNb = 10; // Number of intervals on the curve #1
-  SplitRangeOnSegments(aT11, aT12, theCurveRes1, aNb, aSeg1);
-  aNb1 = aSeg1.Length();
+  aNb1 = SplitRangeOnSegments(aT11, aT12, theCurveRes1, aNb, aRanges);
   for (i = 1; i < aNb1; ++i) {
-    const IntTools_Range& aR1 = aSeg1(i);
+    const IntTools_Range& aR1 = aRanges(i);
     aR1.Range(aT1A, aT1B);
     //
     iErr = DistPC(aT1B, myGeom1, theCriteria, aProjPC, aDmax, aT2max);
@@ -836,7 +991,7 @@ Standard_Integer IntTools_EdgeEdge::CheckCoincidence(const Standard_Real aT11,
     }
   }
   //
-  // if the ranges in aSeg1 are less than theCurveRes1,
+  // if the ranges in aRanges are less than theCurveRes1,
   // there is no need to do step 2 (deep evaluation)
   if (aNb1 < aNb) {
     return iErr;
@@ -844,7 +999,7 @@ Standard_Integer IntTools_EdgeEdge::CheckCoincidence(const Standard_Real aT11,
   //
   // 2. Deep evaluation
   for (i = 2; i < aNb1; ++i) {
-    const IntTools_Range& aR1 = aSeg1(i);
+    const IntTools_Range& aR1 = aRanges(i);
     aR1.Range(aT1A, aT1B);
     //
     iErr = FindDistPC(aT1A, aT1B, myGeom1, theCriteria, theCurveRes1, 
@@ -886,12 +1041,14 @@ Standard_Integer FindDistPC(const Standard_Real aT1A,
   aB = aT1B;
   //
   // check bounds
-  iErr = DistPC(aA, theC1, theCriteria, theProjPC, aYP, aT2P, aDmax, aT1max, aT2max, iC);
+  iErr = DistPC(aA, theC1, theCriteria, theProjPC, 
+                aYP, aT2P, aDmax, aT1max, aT2max, iC);
   if (iErr == 2) {
     return iErr;
   }
   //
-  iErr = DistPC(aB, theC1, theCriteria, theProjPC, aYL, aT2L, aDmax, aT1max, aT2max, iC);
+  iErr = DistPC(aB, theC1, theCriteria, theProjPC, 
+                aYL, aT2L, aDmax, aT1max, aT2max, iC);
   if (iErr == 2) {
     return iErr;
   }
@@ -899,12 +1056,14 @@ Standard_Integer FindDistPC(const Standard_Real aT1A,
   aXP = aA + (aB - aA)*aGS;
   aXL = aB - (aB - aA)*aGS;
   //
-  iErr = DistPC(aXP, theC1, theCriteria, theProjPC, aYP, aT2P, aDmax, aT1max, aT2max, iC);
+  iErr = DistPC(aXP, theC1, theCriteria, theProjPC, 
+                aYP, aT2P, aDmax, aT1max, aT2max, iC);
   if (iErr) {
     return iErr;
   }
   //
-  iErr = DistPC(aXL, theC1, theCriteria, theProjPC, aYL, aT2L, aDmax, aT1max, aT2max, iC);
+  iErr = DistPC(aXL, theC1, theCriteria, theProjPC, 
+                aYL, aT2L, aDmax, aT1max, aT2max, iC);
   if (iErr) {
     return iErr;
   }
@@ -915,20 +1074,25 @@ Standard_Integer FindDistPC(const Standard_Real aT1A,
       aXL = aXP;
       aYL = aYP;
       aXP = aA + (aB - aA)*aGS;
-      iErr = DistPC(aXP, theC1, theCriteria, theProjPC, aYP, aT2P, aDmax, aT1max, aT2max, iC);
-      if (iErr) {
-        return iErr;
-      }
+      iErr = DistPC(aXP, theC1, theCriteria, theProjPC, 
+                    aYP, aT2P, aDmax, aT1max, aT2max, iC);
     }
     else {
       aB = aXP;
       aXP = aXL;
       aYP = aYL;
       aXL = aB - (aB - aA)*aGS;
-      iErr = DistPC(aXL, theC1, theCriteria, theProjPC, aYL, aT2L, aDmax, aT1max, aT2max, iC);
-      if (iErr) {
-        return iErr;
+      iErr = DistPC(aXL, theC1, theCriteria, theProjPC, 
+                    aYL, aT2L, aDmax, aT1max, aT2max, iC);
+    }
+    //
+    if (iErr) {
+      if ((iErr == 2) && !bMaxDist) {
+        aXP = (aA + aB) * 0.5;
+        DistPC(aXP, theC1, theCriteria, theProjPC, 
+               aYP, aT2P, aDmax, aT1max, aT2max, iC);
       }
+      return iErr;
     }
     //
     if ((aB - aA) < theEps) {
@@ -956,7 +1120,7 @@ Standard_Integer DistPC(const Standard_Real aT1,
   Standard_Integer iErr;
   //
   iErr = DistPC(aT1, theC1, theCriteria, theProjPC, aD, aT2, iC);
-  if (iErr) {
+  if (iErr == 1) {
     return iErr;
   }
   //
@@ -1006,16 +1170,16 @@ Standard_Integer DistPC(const Standard_Real aT1,
 //function : SplitRangeOnSegments
 //purpose  : 
 //=======================================================================
-void SplitRangeOnSegments(const Standard_Real aT1, 
-                          const Standard_Real aT2,
-                          const Standard_Real theResolution,
-                          const Standard_Integer theNbSeg,
-                          IntTools_SequenceOfRanges& theSegments)
+Standard_Integer SplitRangeOnSegments(const Standard_Real aT1, 
+                                      const Standard_Real aT2,
+                                      const Standard_Real theResolution,
+                                      const Standard_Integer theNbSeg,
+                                      IntTools_SequenceOfRanges& theSegments)
 {
   Standard_Real aDiff = aT2 - aT1;
-  if (aDiff < theResolution) {
+  if (aDiff < theResolution || theNbSeg == 1) {
     theSegments.Append(IntTools_Range(aT1, aT2));
-    return;
+    return 1;
   }
   //
   Standard_Real aDt, aT1x, aT2x, aSeg;
@@ -1041,6 +1205,8 @@ void SplitRangeOnSegments(const Standard_Real aT1,
   //
   IntTools_Range aR(aT1x, aT2);
   theSegments.Append(aR);
+  //
+  return aNbSegments;
 }
 
 //=======================================================================
@@ -1104,12 +1270,12 @@ Standard_Integer TypeToInteger(const GeomAbs_CurveType theCType)
   case GeomAbs_Line:
     iRet=0;
     break;
-  case GeomAbs_Circle:
-    iRet=1;
-    break;
-  case GeomAbs_Ellipse:
   case GeomAbs_Hyperbola:
   case GeomAbs_Parabola:
+    iRet=1;
+    break;
+  case GeomAbs_Circle:
+  case GeomAbs_Ellipse:
     iRet=2;
     break;
   case GeomAbs_BezierCurve:
@@ -1123,3 +1289,152 @@ Standard_Integer TypeToInteger(const GeomAbs_CurveType theCType)
   return iRet;
 }
 
+//=======================================================================
+//function : ResolutionCoeff
+//purpose  : 
+//=======================================================================
+Standard_Real ResolutionCoeff(const BRepAdaptor_Curve& theBAC,
+                              const IntTools_Range& theRange)
+{
+  Standard_Real aResCoeff;
+  //
+  const Handle(Geom_Curve)& aCurve = theBAC.Curve().Curve();
+  const GeomAbs_CurveType aCurveType = theBAC.GetType();
+  //
+  switch (aCurveType) {
+  case GeomAbs_Circle :
+    aResCoeff = 1. / (2 * (*((Handle(Geom_Circle)*)&aCurve))->Circ().Radius());
+    break;
+  case GeomAbs_Ellipse :
+    aResCoeff =  1. / (*((Handle(Geom_Ellipse)*)&aCurve))->MajorRadius();
+    break;
+  case GeomAbs_Hyperbola :
+  case GeomAbs_Parabola : 
+  case GeomAbs_OtherCurve :{
+    Standard_Real k, kMin, aDist, aDt, aT1, aT2, aT;
+    Standard_Integer aNbP, i;
+    gp_Pnt aP1, aP2;
+    //
+    aNbP = 30;
+    theRange.Range(aT1, aT2);
+    aDt = (aT2 - aT1) / aNbP;
+    aT = aT1;
+    kMin = 10.;
+    //
+    theBAC.D0(aT1, aP1);
+    for (i = 1; i <= aNbP; ++i) {
+      aT += aDt;
+      theBAC.D0(aT, aP2);
+      aDist = aP1.Distance(aP2);
+      k = aDt / aDist;
+      if (k < kMin) {
+        kMin = k;
+      }
+      aP1 = aP2;
+    }
+    //
+    aResCoeff = kMin;
+    break;
+  }
+  default:
+    aResCoeff = 0.;
+    break;
+  }
+  //
+  return aResCoeff;
+}
+
+//=======================================================================
+//function : Resolution
+//purpose  : 
+//=======================================================================
+Standard_Real Resolution(const Handle(Geom_Curve)& theCurve,
+                         const GeomAbs_CurveType theCurveType,
+                         const Standard_Real theResCoeff,
+                         const Standard_Real theR3D)
+{
+  Standard_Real aRes;
+  //
+  switch (theCurveType) {
+  case GeomAbs_Line :
+    aRes = theR3D;
+    break;
+  case GeomAbs_Circle: {
+    Standard_Real aDt = theResCoeff * theR3D;
+    aRes = (aDt <= 1.) ? 2*ASin(aDt) : 2*M_PI;
+    break;
+  }
+  case GeomAbs_BezierCurve:
+    (*((Handle(Geom_BezierCurve)*)&theCurve))->Resolution(theR3D, aRes);
+    break;
+  case GeomAbs_BSplineCurve:
+    (*((Handle(Geom_BSplineCurve)*)&theCurve))->Resolution(theR3D, aRes);
+    break;
+  default:
+    aRes = theResCoeff * theR3D;
+    break;
+  }
+  //
+  return aRes;
+}
+
+//=======================================================================
+//function : CurveDeflection
+//purpose  : 
+//=======================================================================
+Standard_Real CurveDeflection(const BRepAdaptor_Curve& theBAC,
+                              const IntTools_Range& theRange)
+{
+  Standard_Real aDt, aT, aT1, aT2, aDefl;
+  Standard_Integer i, aNbP;
+  gp_Vec aV1, aV2;
+  gp_Pnt aP;
+  //
+  aDefl = 0;
+  aNbP = 10;
+  theRange.Range(aT1, aT2);
+  aDt = (aT2 - aT1) / aNbP;
+  aT = aT1;
+  //
+  theBAC.D1(aT1, aP, aV1);
+  for (i = 1; i <= aNbP; ++i) {
+    aT += aDt;
+    theBAC.D1(aT, aP, aV2);
+    if (aV1.Magnitude() > gp::Resolution() &&
+        aV2.Magnitude() > gp::Resolution()) {
+      gp_Dir aD1(aV1), aD2(aV2);
+      aDefl += aD1.Angle(aD2);
+    }
+    aV1 = aV2;
+  }
+  //
+  return aDefl;
+}
+
+//=======================================================================
+//function : IsClosed
+//purpose  : 
+//=======================================================================
+Standard_Integer IsClosed(const Handle(Geom_Curve)& theCurve,
+                          const Standard_Real aT1,
+                          const Standard_Real aT2,
+                          const Standard_Real theTol,
+                          const Standard_Real theRes)
+{
+  Standard_Boolean bClosed;
+  Standard_Real aD;
+  gp_Pnt aP1, aP2;
+  //
+  bClosed = Standard_False;
+  if (Abs(aT1 - aT2) < theRes) {
+    return bClosed;
+  }
+  //
+  theCurve->D0(aT1, aP1);
+  theCurve->D0(aT2, aP2);
+  //
+  aD = aP1.Distance(aP2);
+  bClosed = aD < theTol;
+  //
+  return bClosed;
+}

src/IntTools/IntTools_EdgeEdge.lxx

@@ -26,6 +26,8 @@ inline IntTools_EdgeEdge::IntTools_EdgeEdge()
   myTol(0.),
   myRes1(0.),
   myRes2(0.),
+  myResCoeff1(0.),
+  myResCoeff2(0.),
   myPTol1(0.),
   myPTol2(0.),
   myRange1(0., 0.),
@@ -48,6 +50,8 @@ inline IntTools_EdgeEdge::IntTools_EdgeEdge(const TopoDS_Edge&  theEdge1,
   myTol(0.),
   myRes1(0.),
   myRes2(0.),
+  myResCoeff1(0.),
+  myResCoeff2(0.),
   myPTol1(0.),
   myPTol2(0.),
   myRange1(0., 0.),
@@ -74,6 +78,8 @@ inline IntTools_EdgeEdge::IntTools_EdgeEdge(const TopoDS_Edge&  theEdge1,
   myTol(0.),
   myRes1(0.),
   myRes2(0.),
+  myResCoeff1(0.),
+  myResCoeff2(0.),
   myPTol1(0.),
   myPTol2(0.),
   myRange1(aT11, aT12),

src/ProjLib/ProjLib_PrjResolve.cxx

@@ -92,13 +92,13 @@ ProjLib_PrjResolve::ProjLib_PrjResolve(const Adaptor3d_Curve& C,const Adaptor3d_
     ExtraU = Tol2d.X();
     ExtraV = Tol2d.Y();
 //  }
-  if (Abs(mySolution.X()-Inf.X()) < Tol2d.X()) mySolution.SetX(Inf.X());
-  if (Abs(mySolution.X()-Sup.X()) < Tol2d.X()) mySolution.SetX(Sup.X()); 
-  if (Abs(mySolution.Y()-Inf.Y()) < Tol2d.Y()) mySolution.SetY(Inf.Y());
-  if (Abs(mySolution.Y()-Sup.Y()) < Tol2d.Y()) mySolution.SetY(Sup.Y()); 
-  if (mySolution.X() < Inf.X() - ExtraU || 
+  if (mySolution.X() > Inf.X() - Tol2d.X() && mySolution.X() < Inf.X()) mySolution.SetX(Inf.X());
+  if (mySolution.X() > Sup.X() && mySolution.X() < Sup.X() + Tol2d.X()) mySolution.SetX(Sup.X()); 
+  if (mySolution.Y() > Inf.Y() - Tol2d.Y() && mySolution.Y() < Inf.Y()) mySolution.SetY(Inf.Y());
+  if (mySolution.Y() > Sup.Y() && mySolution.Y() < Sup.Y() + Tol2d.Y()) mySolution.SetY(Sup.Y()); 
+  if (mySolution.X() < Inf.X() - ExtraU ||
       mySolution.X() > Sup.X() + ExtraU ||
-      mySolution.Y() < Inf.Y() - ExtraV || 
+      mySolution.Y() < Inf.Y() - ExtraV ||
       mySolution.Y() > Sup.Y() + ExtraV) myDone = Standard_False;
   else if (FuncTol > 0) {
     math_Vector X(1,2,0.), FVal(1,2,0.);

src/ShapeFix/ShapeFix_Shell.cxx

@@ -223,19 +223,19 @@ static  Standard_Boolean GetShells(TopTools_SequenceOfShape& Lface,
         continue;
       
       if((edge.Orientation() == TopAbs_FORWARD && dire.Contains(edge))
-	 || (edge.Orientation() == TopAbs_REVERSED && reve.Contains(edge))) 
-	nbbe++;
+        || (edge.Orientation() == TopAbs_REVERSED && reve.Contains(edge))) 
+        nbbe++;
       else if((edge.Orientation() == TopAbs_FORWARD && reve.Contains(edge))
-	      || (edge.Orientation() == TopAbs_REVERSED && dire.Contains(edge)))   
-	nbe++;
+        || (edge.Orientation() == TopAbs_REVERSED && dire.Contains(edge)))   
+        nbe++;
       
       if(dire.Contains(edge)) dire.Remove(edge);
       else 
-	if(reve.Contains(edge)) reve.Remove(edge);
-	else {
-	  if(edge.Orientation() == TopAbs_FORWARD) dtemp.Add(edge);
-	  if(edge.Orientation() == TopAbs_REVERSED) rtemp.Add(edge);
-	}
+        if(reve.Contains(edge)) reve.Remove(edge);
+        else {
+          if(edge.Orientation() == TopAbs_FORWARD) dtemp.Add(edge);
+          if(edge.Orientation() == TopAbs_REVERSED) rtemp.Add(edge);
+        }
     }
     if(!nbbe && !nbe && dtemp.IsEmpty() && rtemp.IsEmpty()) 
       continue;
@@ -251,21 +251,21 @@ static  Standard_Boolean GetShells(TopTools_SequenceOfShape& Lface,
     
     // Addition of face to shell. In the dependance of orientation faces in the shell 
     //  added face can be reversed.
-    
+
     if((nbe != 0 || nbbe != 0) || j == 1) {
       if(nbbe != 0) {
-	F1.Reverse();
-	for(TopTools_MapIteratorOfMapOfShape ite(dtemp); ite.More(); ite.Next()) 
-	  reve.Add(ite.Key());
-	for(TopTools_MapIteratorOfMapOfShape ite1(rtemp); ite1.More(); ite1.Next())
-	  dire.Add(ite1.Key());
-	done = Standard_True;
+        F1.Reverse();
+        for(TopTools_MapIteratorOfMapOfShape ite(dtemp); ite.More(); ite.Next()) 
+          reve.Add(ite.Key());
+        for(TopTools_MapIteratorOfMapOfShape ite1(rtemp); ite1.More(); ite1.Next())
+          dire.Add(ite1.Key());
+        done = Standard_True;
       }
       else {
-	for(TopTools_MapIteratorOfMapOfShape ite(dtemp); ite.More(); ite.Next()) 
-	  dire.Add(ite.Key());
-	for(TopTools_MapIteratorOfMapOfShape ite1(rtemp); ite1.More(); ite1.Next())
-	  reve.Add(ite1.Key());
+        for(TopTools_MapIteratorOfMapOfShape ite(dtemp); ite.More(); ite.Next()) 
+          dire.Add(ite.Key());
+        for(TopTools_MapIteratorOfMapOfShape ite1(rtemp); ite1.More(); ite1.Next())
+          reve.Add(ite1.Key());
       }
       j++;
       B.Add(nshell,F1);
@@ -273,7 +273,7 @@ static  Standard_Boolean GetShells(TopTools_SequenceOfShape& Lface,
       Lface.Remove(i);
       
       // if closed shell is obtained it adds to sequence of shells and new shell begin to construct.
-      if(isMultiConnex && BRep_Tool::IsClosed(nshell)) {
+      if(isMultiConnex && BRep_Tool::IsClosed (nshell)) {
         aSeqShells.Append(nshell);
         TopoDS_Shell nshellnext;
         B.MakeShell(nshellnext);
@@ -862,8 +862,16 @@ Standard_Boolean ShapeFix_Shell::FixFaceOrientation(const TopoDS_Shell& shell,co
   myShell = shell;
   myShape = shell;
   Standard_Integer aNumMultShell =0;
-  for (TopoDS_Iterator iter(shell); iter.More(); iter.Next()) 
-    Lface.Append(iter.Value());
+  Standard_Integer nbF = 0;
+  TopTools_MapOfShape aMapAdded;
+  for (TopoDS_Iterator iter(shell); iter.More(); iter.Next(),nbF++) 
+  {
+    if(aMapAdded.Add(iter.Value()))
+      Lface.Append(iter.Value());
+  }
+  if(Lface.Length() < nbF)
+    done = Standard_True;
+
   TopTools_IndexedDataMapOfShapeListOfShape aMapEdgeFaces;
   TopExp::MapShapesAndAncestors(myShell,TopAbs_EDGE,TopAbs_FACE,aMapEdgeFaces);
   TopTools_MapOfShape aMapMultiConnectEdges;

src/ShapeFix/ShapeFix_Solid.cxx

@@ -20,7 +20,7 @@
 #include <BRepClass3d_SolidClassifier.hxx>
 #include <Precision.hxx>
 #include <TopoDS_Shape.hxx>
-#include <ShapeBuild_ReShape.hxx> 
+#include <ShapeBuild_ReShape.hxx>
 #include <TopoDS_Iterator.hxx>
 #include <TopoDS.hxx>
 #include <ShapeExtend.hxx>
@@ -244,9 +244,9 @@ static void CollectSolids(const TopTools_SequenceOfShape& aSeqShells ,
     }
   }
   for(TopTools_MapIteratorOfMapOfShape aIterHoles(aMapHoles);aIterHoles.More(); aIterHoles.Next())
-    aMapShellHoles.UnBind(aIterHoles.Key());
-    
-}
+    aMapShellHoles.UnBind (aIterHoles.Key());
+
+    }
 //=======================================================================
 //function : CreateSolids
 //purpose  : 
@@ -343,13 +343,30 @@ static Standard_Boolean CreateSolids(const TopoDS_Shape aShape,TopTools_IndexedM
     BRep_Builder aB;
     aB.MakeCompSolid(aCompSolid);
     isDone = (aShape.ShapeType() != TopAbs_COMPSOLID || isDone);
+    Standard_Integer nbSol = 0;
+
     for(TopTools_ListIteratorOfListOfShape lItSh(lshells);lItSh.More(); lItSh.Next()) {
       if(ShellSolid.Contains(lItSh.Value())) {
-        for(TopExp_Explorer aExpSol(ShellSolid.FindFromKey(lItSh.Value()),TopAbs_SOLID);aExpSol.More(); aExpSol.Next())
+        const TopoDS_Shape& aShape = ShellSolid.FindFromKey(lItSh.Value());
+        TopExp_Explorer aExpSol(aShape, TopAbs_SOLID);
+       
+        for(;aExpSol.More(); aExpSol.Next())
+        {
           aB.Add(aCompSolid,aExpSol.Current());
-        ShellSolid.ChangeFromKey(lItSh.Value()) = aCompSolid;
+          nbSol++;
+        }
+      
       }
     }
+    if(nbSol >1)
+    {
+      for(TopTools_ListIteratorOfListOfShape lItSh1(lshells);lItSh1.More(); lItSh1.Next()) 
+      {
+        if(ShellSolid.Contains(lItSh1.Value())) 
+          ShellSolid.ChangeFromKey(lItSh1.Value()) = aCompSolid;
+      }
+    }
+    
   }
   for(Standard_Integer kk =1 ; kk <= ShellSolid.Extent();kk++)
     if(!aMapSolids.Contains(ShellSolid.FindFromIndex(kk)))
@@ -433,22 +450,24 @@ Standard_Boolean ShapeFix_Solid::Perform(const Handle(Message_ProgressIndicator)
       
     if(isClosed || myCreateOpenSolidMode) {
       if(BRep_Tool::IsClosed(tmpShape)) {
-        TopoDS_Iterator itersh(tmpShape);
-        TopoDS_Shell aShell;
-        if(itersh.More() && itersh.Value().ShapeType() == TopAbs_SHELL)
-          aShell = TopoDS::Shell(itersh.Value());
-        if(!aShell.IsNull()) {
-          TopoDS_Solid aSol = SolidFromShell(aShell);
-          if(ShapeExtend::DecodeStatus(myStatus,ShapeExtend_DONE2)) {
-            SendWarning (Message_Msg ("FixAdvSolid.FixOrientation.MSG20"));// Orientaion of shell was corrected.
-            Context()->Replace(tmpShape,aSol);
-            tmpShape = aSol;
-          }
+      TopoDS_Iterator itersh(tmpShape);
+      TopoDS_Shell aShell;
+      if(itersh.More() && itersh.Value().ShapeType() == TopAbs_SHELL)
+        aShell = TopoDS::Shell(itersh.Value());
+      if(!aShell.IsNull()) {
+        TopoDS_Solid aSol = SolidFromShell(aShell);
+        if(ShapeExtend::DecodeStatus(myStatus,ShapeExtend_DONE2)) {
+          SendWarning (Message_Msg ("FixAdvSolid.FixOrientation.MSG20"));// Orientaion of shell was corrected.
+          Context()->Replace(tmpShape,aSol);
+          tmpShape = aSol;
         }
       }
+      }
       mySolid  = TopoDS::Solid(tmpShape);
     }
     else {
+      status = Standard_True;
+      myStatus |= ShapeExtend::EncodeStatus ( ShapeExtend_DONE3 );
       TopoDS_Iterator aIt(tmpShape,Standard_False);
       Context()->Replace(tmpShape,aIt.Value());
       SendFail (Message_Msg ("FixAdvSolid.FixShell.MSG10")); // Solid can not be created from open shell. 

src/math/FILES

@@ -4,3 +4,10 @@ math_Memory.hxx
 math_Recipes.hxx
 math_GaussPoints.hxx
 math_Kronrod.cxx
+math_Vector.hxx
+math_Vector.cxx
+math_IntegerVector.hxx
+math_IntegerVector.cxx
+math_SingleTab.hxx
+math_GlobOptMin.hxx
+math_GlobOptMin.cxx

src/math/math.cdl

@@ -41,8 +41,8 @@ is
     exception SingularMatrix inherits Failure;
     
 
-    class Vector;
-    class IntegerVector;
+    imported Vector;
+    imported IntegerVector;
     class Matrix;
     deferred class Function;
     deferred class FunctionWithDerivative;
@@ -51,6 +51,7 @@ is
     deferred class MultipleVarFunctionWithHessian;
     deferred class FunctionSet;
     deferred class FunctionSetWithDerivatives;
+    imported GlobOptMin;
     class IntegerRandom;
     class Gauss;
     class GaussLeastSquare;
@@ -93,13 +94,7 @@ is
 	    	    	    	    	       Array1OfValueAndWeight from math,
 	    	    	    	    	       CompareOfValueAndWeight from math);
 
-  generic class SingleTab;
-  generic class DoubleTab;
-
-    --- Instantiate classes:
-    class SingleTabOfReal instantiates SingleTab(Real);
-    class SingleTabOfInteger instantiates SingleTab(Integer);
-    class DoubleTabOfReal instantiates DoubleTab(Real);
+	class DoubleTab;
 
     --- Gauss Points
 

src/math/math_DoubleTab.cdl

@@ -14,17 +14,17 @@
 -- Alternatively, this file may be used under the terms of Open CASCADE
 -- commercial license or contractual agreement.
 
-generic class DoubleTab from math (Item as any)
+class DoubleTab from math
 uses Address from Standard
 is
 
     Create(LowerRow, UpperRow, LowerCol, UpperCol: Integer)
     	returns DoubleTab;
 
-    Create(Tab : Item; LowerRow, UpperRow, LowerCol, UpperCol: Integer)
+    Create(Tab : Address; LowerRow, UpperRow, LowerCol, UpperCol: Integer)
     	returns DoubleTab;
 
-    Init(me : in out; InitValue: Item) is static;
+    Init(me : in out; InitValue: Real) is static;
 
     Create(Other: DoubleTab)
     	returns DoubleTab;
@@ -48,7 +48,7 @@ is
     	---C++: alias operator()
     	---C++: return &
     	---C++: inline
-      returns Item
+      returns Real
       is static;		
 	
 	
@@ -62,7 +62,7 @@ fields
 
 Addr        : Address;
 AddrBuf     : Address[32];
-Buf         : Item[512];
+Buf         : Real[512];
 isAddrAllocated: Boolean;
 isAllocated : Boolean;
 LowR        : Integer;

src/math/math_DoubleTab.cxx

@@ -0,0 +1,142 @@
+// Copyright (c) 1997-1999 Matra Datavision
+// Copyright (c) 1999-2014 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+// Lpa, le 7/02/92
+#include <math_DoubleTab.ixx>
+
+#include <math_Memory.hxx>
+#include <Standard_OutOfRange.hxx>
+#include <Standard_Failure.hxx>
+#include <Standard_Integer.hxx>
+
+// macro to get size of C array
+#define CARRAY_LENGTH(arr) (int)(sizeof(arr)/sizeof(arr[0]))
+
+void math_DoubleTab::Allocate()
+{
+  Standard_Integer RowNumber = UppR - LowR + 1;
+  Standard_Integer ColNumber = UppC - LowC + 1;
+
+  Standard_Real** TheAddr = !isAddrAllocated? (Standard_Real**)&AddrBuf : 
+    (Standard_Real**) Standard::Allocate(RowNumber * sizeof(Standard_Real*));
+  Standard_Real* Address;
+  if(isAllocated) 
+    Address = (Standard_Real*) Standard::Allocate(RowNumber * ColNumber * sizeof(Standard_Real));
+  else
+    Address = (Standard_Real*) Addr;
+  Address -= LowC;
+  
+  for (Standard_Integer Index = 0; Index < RowNumber; Index++) {
+    TheAddr[Index] = Address;
+    Address += ColNumber;
+  }
+  
+  TheAddr -= LowR;
+  Addr = (Standard_Address) TheAddr;
+}
+
+math_DoubleTab::math_DoubleTab(const Standard_Integer LowerRow,
+			       const Standard_Integer UpperRow,
+			       const Standard_Integer LowerCol,
+			       const Standard_Integer UpperCol) :
+  Addr(Buf),
+  isAddrAllocated(UpperRow - LowerRow + 1 > CARRAY_LENGTH(AddrBuf)),
+  isAllocated((UpperRow - LowerRow + 1) * (UpperCol - LowerCol + 1) > CARRAY_LENGTH(Buf)),
+  LowR(LowerRow),
+  UppR(UpperRow),
+  LowC(LowerCol),
+  UppC(UpperCol)
+{
+  Allocate();
+}
+
+math_DoubleTab::math_DoubleTab(const Standard_Address Tab,
+			       const Standard_Integer LowerRow,
+			       const Standard_Integer UpperRow,
+			       const Standard_Integer LowerCol,
+			       const Standard_Integer UpperCol) :
+  Addr(Tab),
+  isAddrAllocated(UpperRow - LowerRow + 1 > CARRAY_LENGTH(AddrBuf)),
+  isAllocated(Standard_False),
+  LowR(LowerRow),
+  UppR(UpperRow),
+  LowC(LowerCol),
+  UppC(UpperCol)
+{
+  Allocate();
+}
+
+void math_DoubleTab::Init(const Standard_Real InitValue) 
+{
+  for (Standard_Integer i = LowR; i <= UppR; i++) {
+    for (Standard_Integer j = LowC; j <= UppC; j++) {
+      ((Standard_Real**) Addr)[i][j] = InitValue;
+    }
+  }
+}
+
+math_DoubleTab::math_DoubleTab(const math_DoubleTab& Other) :
+  Addr(Buf),
+  isAddrAllocated(Other.UppR - Other.LowR + 1 > CARRAY_LENGTH(AddrBuf)),
+  isAllocated((Other.UppR - Other.LowR + 1) * 
+              (Other.UppC - Other.LowC + 1) > CARRAY_LENGTH(Buf)),
+  LowR(Other.LowR),
+  UppR(Other.UppR),
+  LowC(Other.LowC),
+  UppC(Other.UppC)
+{
+  Allocate();
+
+  Standard_Address target = (Standard_Address) &Value(LowR,LowC);
+  Standard_Address source = (Standard_Address) &Other.Value(LowR,LowC);
+
+  memmove(target,source,
+	  (int)((UppR - LowR + 1) * (UppC - LowC + 1) * sizeof(Standard_Real)));
+
+}
+
+void math_DoubleTab::Free()
+{
+  // free the data
+  if(isAllocated) {
+    Standard_Address it = (Standard_Address)&Value(LowR,LowC);
+    Standard::Free(it);
+  }
+  // free the pointers
+  if(isAddrAllocated) {
+    Standard_Address it = (Standard_Address)(((Standard_Real**)Addr) + LowR);
+    Standard::Free (it);
+  }
+  Addr = 0;
+}
+
+void math_DoubleTab::SetLowerRow(const Standard_Integer LowerRow)
+{
+  Standard_Real** TheAddr = (Standard_Real**)Addr;
+  Addr = (Standard_Address) (TheAddr + LowR - LowerRow);
+  UppR = UppR - LowR + LowerRow;
+  LowR = LowerRow;
+}
+
+void math_DoubleTab::SetLowerCol(const Standard_Integer LowerCol)
+{
+  Standard_Real** TheAddr = (Standard_Real**) Addr;
+  for (Standard_Integer Index = LowR; Index <= UppR; Index++) {
+    TheAddr[Index] = TheAddr[Index] + LowC - LowerCol;
+  }
+
+  UppC = UppC - LowC + LowerCol;
+  LowC = LowerCol;
+}
+

src/math/math_DoubleTab.lxx

@@ -17,10 +17,10 @@
 #include <math_Memory.hxx>
 #include <Standard_OutOfRange.hxx>
 
-inline Item& math_DoubleTab::Value (const Standard_Integer RowIndex,
+inline Standard_Real& math_DoubleTab::Value (const Standard_Integer RowIndex,
 				    const Standard_Integer ColIndex) const
 {
-  return ((Item**)Addr)[RowIndex][ColIndex];
+  return ((Standard_Real**)Addr)[RowIndex][ColIndex];
 }
 
 
@@ -29,7 +29,7 @@ inline void math_DoubleTab::Copy(math_DoubleTab& Other)const
 {
   memmove((void*)(& Other.Value(Other.LowR,Other.LowC)),
 	  (void*) (& Value(LowR,LowC)),
-	  (int)((UppR - LowR + 1) * (UppC - LowC + 1) * sizeof(Item)));
+	  (int)((UppR - LowR + 1) * (UppC - LowC + 1) * sizeof(Standard_Real)));
 }
 
 

src/math/math_GlobOptMin.cxx

@@ -0,0 +1,543 @@
+// Created on: 2014-01-20
+// Created by: Alexaner Malyshev
+// Copyright (c) 2014-2014 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement
+
+#include <math_GlobOptMin.hxx>
+
+#include <math_BFGS.hxx>
+#include <math_Matrix.hxx>
+#include <math_MultipleVarFunctionWithGradient.hxx>
+#include <math_MultipleVarFunctionWithHessian.hxx>
+#include <math_NewtonMinimum.hxx>
+#include <math_Powell.hxx>
+#include <Standard_Integer.hxx>
+#include <Standard_Real.hxx>
+#include <Precision.hxx>
+
+
+//=======================================================================
+//function : math_GlobOptMin
+//purpose  : Constructor
+//=======================================================================
+math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc,
+                                 const math_Vector& theA,
+                                 const math_Vector& theB,
+                                 const Standard_Real theC,
+                                 const Standard_Real theDiscretizationTol,
+                                 const Standard_Real theSameTol)
+: myN(theFunc->NbVariables()),
+  myA(1, myN),
+  myB(1, myN),
+  myGlobA(1, myN),
+  myGlobB(1, myN),
+  myX(1, myN),
+  myTmp(1, myN),
+  myV(1, myN),
+  myMaxV(1, myN),
+  myExpandCoeff(1, myN)
+{
+  Standard_Integer i;
+
+  myFunc = theFunc;
+  myC = theC;
+  myZ = -1;
+  mySolCount = 0;
+
+  for(i = 1; i <= myN; i++)
+  {
+    myGlobA(i) = theA(i);
+    myGlobB(i) = theB(i);
+
+    myA(i) = theA(i);
+    myB(i) = theB(i);
+  }
+
+  for(i = 1; i <= myN; i++)
+  {
+    myMaxV(i) = (myB(i) - myA(i)) / 3.0;
+  }
+
+  myExpandCoeff(1) = 1.0;
+  for(i = 2; i <= myN; i++)
+  {
+    myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
+  }
+
+  myTol = theDiscretizationTol;
+  mySameTol = theSameTol;
+
+  myDone = Standard_False;
+}
+
+//=======================================================================
+//function : SetGlobalParams
+//purpose  : Set params without memory allocation.
+//=======================================================================
+void math_GlobOptMin::SetGlobalParams(math_MultipleVarFunction* theFunc,
+                                      const math_Vector& theA,
+                                      const math_Vector& theB,
+                                      const Standard_Real theC,
+                                      const Standard_Real theDiscretizationTol,
+                                      const Standard_Real theSameTol)
+{
+  Standard_Integer i;
+
+  myFunc = theFunc;
+  myC = theC;
+  myZ = -1;
+  mySolCount = 0;
+
+  for(i = 1; i <= myN; i++)
+  {
+    myGlobA(i) = theA(i);
+    myGlobB(i) = theB(i);
+
+    myA(i) = theA(i);
+    myB(i) = theB(i);
+  }
+
+  for(i = 1; i <= myN; i++)
+  {
+    myMaxV(i) = (myB(i) - myA(i)) / 3.0;
+  }
+
+  myExpandCoeff(1) = 1.0;
+  for(i = 2; i <= myN; i++)
+  {
+    myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
+  }
+
+  myTol = theDiscretizationTol;
+  mySameTol = theSameTol;
+
+  myDone = Standard_False;
+}
+
+//=======================================================================
+//function : SetLocalParams
+//purpose  : Set params without memory allocation.
+//=======================================================================
+void math_GlobOptMin::SetLocalParams(const math_Vector& theLocalA,
+                                     const math_Vector& theLocalB)
+{
+  Standard_Integer i;
+
+  myZ = -1;
+  mySolCount = 0;
+
+  for(i = 1; i <= myN; i++)
+  {
+    myA(i) = theLocalA(i);
+    myB(i) = theLocalB(i);
+  }
+
+  for(i = 1; i <= myN; i++)
+  {
+    myMaxV(i) = (myB(i) - myA(i)) / 3.0;
+  }
+
+  myExpandCoeff(1) = 1.0;
+  for(i = 2; i <= myN; i++)
+  {
+    myExpandCoeff(i) = (myB(i) - myA(i)) / (myB(i - 1) - myA(i - 1));
+  }
+
+  myDone = Standard_False;
+}
+
+//=======================================================================
+//function : SetTol
+//purpose  : Set algorithm tolerances.
+//=======================================================================
+void math_GlobOptMin::SetTol(const Standard_Real theDiscretizationTol,
+                             const Standard_Real theSameTol)
+{
+  myTol = theDiscretizationTol;
+  mySameTol = theSameTol;
+}
+
+//=======================================================================
+//function : GetTol
+//purpose  : Get algorithm tolerances.
+//=======================================================================
+void math_GlobOptMin::GetTol(Standard_Real& theDiscretizationTol,
+                             Standard_Real& theSameTol)
+{
+  theDiscretizationTol = myTol;
+  theSameTol = mySameTol;
+}
+
+//=======================================================================
+//function : ~math_GlobOptMin
+//purpose  : 
+//=======================================================================
+math_GlobOptMin::~math_GlobOptMin()
+{
+}
+
+//=======================================================================
+//function : Perform
+//purpose  : Compute Global extremum point
+//=======================================================================
+// In this algo indexes started from 1, not from 0.
+void math_GlobOptMin::Perform()
+{
+  Standard_Integer i;
+
+  // Compute parameters range
+  Standard_Real minLength = RealLast();
+  Standard_Real maxLength = RealFirst();
+  for(i = 1; i <= myN; i++)
+  {
+    Standard_Real currentLength = myB(i) - myA(i);
+    if (currentLength < minLength)
+      minLength = currentLength;
+    if (currentLength > maxLength)
+      maxLength = currentLength;
+  }
+
+  if (minLength < Precision::PConfusion())
+  {
+    #ifdef OCCT_DEBUG
+    cout << "math_GlobOptMin::Perform(): Degenerated parameters space" << endl;
+    #endif
+
+    return;
+  }
+
+  // Compute initial values for myF, myY, myC.
+  computeInitialValues();
+
+  myE1 = minLength * myTol;
+  myE2 = maxLength * myTol;
+    if (myC > 1.0)
+      myE3 = - maxLength * myTol / 4.0;
+    else
+      myE3 = - maxLength * myTol * myC / 4.0;
+
+  computeGlobalExtremum(myN);
+
+  myDone = Standard_True;
+}
+
+//=======================================================================
+//function : computeLocalExtremum
+//purpose  :
+//=======================================================================
+Standard_Boolean math_GlobOptMin::computeLocalExtremum(const math_Vector& thePnt,
+                                                       Standard_Real& theVal,
+                                                       math_Vector& theOutPnt)
+{
+  Standard_Integer i;
+
+  //Newton method
+  if (dynamic_cast<math_MultipleVarFunctionWithHessian*>(myFunc))
+  {
+    math_MultipleVarFunctionWithHessian* myTmp = 
+      dynamic_cast<math_MultipleVarFunctionWithHessian*> (myFunc);
+    
+    math_NewtonMinimum newtonMinimum(*myTmp);
+    newtonMinimum.Perform(*myTmp, thePnt);
+
+    if (newtonMinimum.IsDone())
+    {
+      newtonMinimum.Location(theOutPnt);
+      theVal = newtonMinimum.Minimum();
+    }
+    else return Standard_False;
+  } else
+
+  // BFGS method used.
+  if (dynamic_cast<math_MultipleVarFunctionWithGradient*>(myFunc))
+  {
+    math_MultipleVarFunctionWithGradient* myTmp = 
+      dynamic_cast<math_MultipleVarFunctionWithGradient*> (myFunc);
+    math_BFGS bfgs(*myTmp);
+    bfgs.Perform(*myTmp, thePnt);
+    if (bfgs.IsDone())
+    {
+      bfgs.Location(theOutPnt);
+      theVal = bfgs.Minimum();
+    }
+    else return Standard_False;
+  } else
+
+  // Powell method used.
+  if (dynamic_cast<math_MultipleVarFunction*>(myFunc))
+  {
+    math_Matrix m(1, myN, 1, myN, 0.0);
+    for(i = 1; i <= myN; i++)
+      m(1, 1) = 1.0;
+
+    math_Powell powell(*myFunc, 1e-10);
+    powell.Perform(*myFunc, thePnt, m);
+
+    if (powell.IsDone())
+    {
+      powell.Location(theOutPnt);
+      theVal = powell.Minimum();
+    }
+    else return Standard_False;
+  }
+
+  if (isInside(theOutPnt))
+    return Standard_True;
+  else
+    return Standard_False;
+}
+
+//=======================================================================
+//function : computeInitialValues
+//purpose  : 
+//=======================================================================
+void math_GlobOptMin::computeInitialValues()
+{
+  Standard_Integer i;
+  math_Vector aCurrPnt(1, myN);
+  math_Vector aBestPnt(1, myN);
+  math_Vector aParamStep(1, myN);
+  Standard_Real aCurrVal = RealLast();
+  Standard_Real aBestVal = RealLast();
+
+  // Check functional value in midpoint, low and upp point border and
+  // in each point try to perform local optimization.
+  aBestPnt = (myA + myB) * 0.5;
+  myFunc->Value(aBestPnt, aBestVal);
+
+  for(i = 1; i <= 3; i++)
+  {
+    aCurrPnt = myA + (myB - myA) * (i - 1) / 2.0;
+
+    if(computeLocalExtremum(aCurrPnt, aCurrVal, aCurrPnt))
+    {
+      // Local Extremum finds better solution than current point.
+      if (aCurrVal < aBestVal)
+      {
+        aBestVal = aCurrVal;
+        aBestPnt = aCurrPnt;
+      }
+    }
+  }
+
+  myF = aBestVal;
+  myY.Clear();
+  for(i = 1; i <= myN; i++)
+    myY.Append(aBestPnt(i));
+  mySolCount++;
+
+  // Lipschitz const approximation
+  Standard_Real aLipConst = 0.0, aPrevValDiag, aPrevValProj;
+  Standard_Integer aPntNb = 13;
+  myFunc->Value(myA, aPrevValDiag);
+  aPrevValProj = aPrevValDiag;
+  Standard_Real aStep = (myB - myA).Norm() / aPntNb;
+  aParamStep = (myB - myA) / aPntNb;
+  for(i = 1; i <= aPntNb; i++)
+  {
+    aCurrPnt = myA + aParamStep * i;
+
+    // Walk over diagonal.
+    myFunc->Value(aCurrPnt, aCurrVal);
+    aLipConst = Max (Abs(aCurrVal - aPrevValDiag), aLipConst);
+    aPrevValDiag = aCurrVal;
+
+    // Walk over diag in projected space aPnt(1) = myA(1) = const.
+    aCurrPnt(1) = myA(1);
+    myFunc->Value(aCurrPnt, aCurrVal);
+    aLipConst = Max (Abs(aCurrVal - aPrevValProj), aLipConst);
+    aPrevValProj = aCurrVal;
+  }
+
+  aLipConst *= Sqrt(myN) / aStep;
+
+  if (aLipConst < myC * 0.1)
+  {
+    myC = Max(aLipConst * 0.1, 0.01);
+  }
+  else if (aLipConst > myC * 10)
+  {
+    myC = Min(myC * 2, 30.0);
+  }
+}
+
+//=======================================================================
+//function : ComputeGlobalExtremum
+//purpose  :
+//=======================================================================
+void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
+{
+  Standard_Integer i;
+  Standard_Real d; // Functional in moved point.
+  Standard_Real val = RealLast(); // Local extrema computed in moved point.
+  Standard_Real aStepBestValue = RealLast();
+  Standard_Real aRealStep = 0.0;
+  math_Vector aStepBestPoint(1, myN);
+  Standard_Boolean isInside = Standard_False;
+  Standard_Boolean isReached = Standard_False;
+  Standard_Real r;
+
+  for(myX(j) = myA(j) + myE1; 
+     (myX(j) < myB(j) + myE1) && (!isReached);
+      myX(j) += myV(j))
+  {
+    if (myX(j) > myB(j))
+    {
+      myX(j) = myB(j);
+      isReached = Standard_True;
+    }
+
+    if (j == 1)
+    {
+      isInside = Standard_False;
+      myFunc->Value(myX, d);
+      r = (d + myZ * myC * aRealStep - myF) * myZ;
+      if(r > myE3)
+      {
+        isInside = computeLocalExtremum(myX, val, myTmp);
+      }
+      aStepBestValue = (isInside && (val < d))? val : d;
+      aStepBestPoint = (isInside && (val < d))? myTmp : myX;
+
+      // Solutions are close to each other.
+      if (Abs(aStepBestValue - myF) < mySameTol * 0.01)
+      {
+        if (!isStored(aStepBestPoint))
+        {
+          if ((aStepBestValue - myF) * myZ > 0.0)
+            myF = aStepBestValue;
+          for(i = 1; i <= myN; i++)
+            myY.Append(aStepBestPoint(i));
+          mySolCount++;
+        }
+      }
+
+      // New best solution.
+      if ((aStepBestValue - myF) * myZ > mySameTol * 0.01)
+      {
+        mySolCount = 0;
+        myF = aStepBestValue;
+        myY.Clear();
+        for(i = 1; i <= myN; i++)
+          myY.Append(aStepBestPoint(i));
+        mySolCount++;
+      }
+
+      aRealStep = myE2 + Abs(myF - d) / myC;
+      myV(1) = Min(aRealStep, myMaxV(1));
+    }
+    else
+    {
+      myV(j) = RealLast() / 2.0;
+      computeGlobalExtremum(j - 1);
+
+      // Nullify steps on lower dimensions.
+      for(i = 1; i < j; i++)
+        myV(i) = 0.0;
+    }
+    // Compute step in (j + 1) dimension according to scale.
+    if (j < myN)
+    {
+      Standard_Real aUpperDimStep =  myV(j) * myExpandCoeff(j + 1);
+      if (myV(j + 1) > aUpperDimStep)
+      {
+        if (aUpperDimStep > myMaxV(j + 1)) // Case of too big step.
+          myV(j + 1) = myMaxV(j + 1); 
+        else
+          myV(j + 1) = aUpperDimStep;
+      }
+    }
+  }
+}
+
+//=======================================================================
+//function : IsInside
+//purpose  :
+//=======================================================================
+Standard_Boolean math_GlobOptMin::isInside(const math_Vector& thePnt)
+{
+  Standard_Integer i;
+
+ for(i = 1; i <= myN; i++)
+ {
+   if (thePnt(i) < myGlobA(i) || thePnt(i) > myGlobB(i))
+     return Standard_False;
+ }
+
+ return Standard_True;
+}
+//=======================================================================
+//function : IsStored
+//purpose  :
+//=======================================================================
+Standard_Boolean math_GlobOptMin::isStored(const math_Vector& thePnt)
+{
+  Standard_Integer i,j;
+  Standard_Boolean isSame = Standard_True;
+
+    for(i = 0; i < mySolCount; i++)
+    {
+      isSame = Standard_True;
+      for(j = 1; j <= myN; j++)
+      {
+      if ((Abs(thePnt(j) - myY(i * myN + j))) > (myB(j) -  myA(j)) * mySameTol)
+        {
+          isSame = Standard_False;
+          break;
+        }
+      }
+      if (isSame == Standard_True)
+        return Standard_True;
+
+      }
+  return Standard_False;
+}
+
+//=======================================================================
+//function : NbExtrema
+//purpose  :
+//=======================================================================
+Standard_Integer math_GlobOptMin::NbExtrema()
+{
+  return mySolCount;
+}
+
+//=======================================================================
+//function : GetF
+//purpose  :
+//=======================================================================
+Standard_Real math_GlobOptMin::GetF()
+{
+  return myF;
+}
+
+//=======================================================================
+//function : IsDone
+//purpose  :
+//=======================================================================
+Standard_Boolean math_GlobOptMin::isDone()
+{
+  return myDone;
+}
+
+//=======================================================================
+//function : Points
+//purpose  :
+//=======================================================================
+void math_GlobOptMin::Points(const Standard_Integer theIndex, math_Vector& theSol)
+{
+  Standard_Integer j;
+
+  for(j = 1; j <= myN; j++)
+    theSol(j) = myY((theIndex - 1) * myN + j);
+}

src/math/math_GlobOptMin.hxx

@@ -0,0 +1,123 @@
+// Created on: 2014-01-20
+// Created by: Alexaner Malyshev
+// Copyright (c) 2014-2014 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#ifndef _math_GlobOptMin_HeaderFile
+#define _math_GlobOptMin_HeaderFile
+
+#include <math_MultipleVarFunction.hxx>
+#include <NCollection_Sequence.hxx>
+#include <Standard_Type.hxx>
+
+//! This class represents Evtushenko's algorithm of global optimization based on nonuniform mesh.<br>
+//! Article: Yu. Evtushenko. Numerical methods for finding global extreme (case of a non-uniform mesh). <br>
+//! U.S.S.R. Comput. Maths. Math. Phys., Vol. 11, N 6, pp. 38-54.
+
+class math_GlobOptMin
+{
+public:
+
+  Standard_EXPORT math_GlobOptMin(math_MultipleVarFunction* theFunc,
+                                 const math_Vector& theA,
+                                 const math_Vector& theB,
+                                 const Standard_Real theC = 9,
+                                 const Standard_Real theDiscretizationTol = 1.0e-2,
+                                 const Standard_Real theSameTol = 1.0e-7);
+
+  Standard_EXPORT void SetGlobalParams(math_MultipleVarFunction* theFunc,
+                                       const math_Vector& theA,
+                                       const math_Vector& theB,
+                                       const Standard_Real theC = 9,
+                                       const Standard_Real theDiscretizationTol = 1.0e-2,
+                                       const Standard_Real theSameTol = 1.0e-7);
+
+  Standard_EXPORT void SetLocalParams(const math_Vector& theLocalA,
+                                      const math_Vector& theLocalB);
+
+  Standard_EXPORT void SetTol(const Standard_Real theDiscretizationTol,
+                              const Standard_Real theSameTol);
+
+  Standard_EXPORT void GetTol(Standard_Real& theDiscretizationTol,
+                              Standard_Real& theSameTol);
+
+  Standard_EXPORT ~math_GlobOptMin();
+
+  Standard_EXPORT void Perform();
+
+  //! Get best functional value.
+  Standard_EXPORT Standard_Real GetF();
+
+  //! Return count of global extremas.
+  Standard_EXPORT Standard_Integer NbExtrema();
+
+  //! Return solution i, 1 <= i <= NbExtrema.
+  Standard_EXPORT void Points(const Standard_Integer theIndex, math_Vector& theSol);
+
+  Standard_Boolean isDone();
+
+private:
+
+  math_GlobOptMin & operator = (const math_GlobOptMin & theOther);
+
+  Standard_Boolean computeLocalExtremum(const math_Vector& thePnt, Standard_Real& theVal, math_Vector& theOutPnt);
+
+  void computeGlobalExtremum(Standard_Integer theIndex);
+
+  //! Computes starting value / approximation:
+  // myF - initial best value.
+  // myY - initial best point.
+  // myC - approximation of Lipschitz constant.
+  // to imporve convergence speed.
+  void computeInitialValues();
+
+  //! Check that myA <= pnt <= myB
+  Standard_Boolean isInside(const math_Vector& thePnt);
+
+  Standard_Boolean isStored(const math_Vector &thePnt);
+
+  // Input.
+  math_MultipleVarFunction* myFunc;
+  Standard_Integer myN;
+  math_Vector myA; // Left border on current C2 interval.
+  math_Vector myB; // Right border on current C2 interval.
+  math_Vector myGlobA; // Global left border.
+  math_Vector myGlobB; // Global right border.
+  Standard_Real myTol; // Discretization tolerance, default 1.0e-2.
+  Standard_Real mySameTol; // points with ||p1 - p2|| < mySameTol is equal,
+                           // function values |val1 - val2| * 0.01 < mySameTol is equal,
+                           // default value is 1.0e-7.
+  Standard_Real myC; //Lipschitz constant, default 9
+
+  // Output.
+  Standard_Boolean myDone;
+  NCollection_Sequence<Standard_Real> myY;// Current solutions.
+  Standard_Integer mySolCount; // Count of solutions.
+
+  // Algorithm data.
+  Standard_Real myZ;
+  Standard_Real myE1; // Border coeff.
+  Standard_Real myE2; // Minimum step size.
+  Standard_Real myE3; // Local extrema starting parameter.
+
+  math_Vector myX; // Current modified solution.
+  math_Vector myTmp; // Current modified solution.
+  math_Vector myV; // Steps array.
+  math_Vector myMaxV; // Max Steps array.
+  math_Vector myExpandCoeff; // Define expand coefficient between neighboring indiced dimensions.
+  Standard_Real myF; // Current value of Global optimum.
+};
+
+const Handle(Standard_Type)& TYPE(math_GlobOptMin);
+
+#endif

src/math/math_IntegerVector.cxx

@@ -12,97 +12,91 @@
 // Alternatively, this file may be used under the terms of Open CASCADE
 // commercial license or contractual agreement.
 
-//#ifndef DEB
 #define No_Standard_RangeError
 #define No_Standard_OutOfRange
 #define No_Standard_DimensionError
-//#endif
 
-#include <math_IntegerVector.ixx>
+#include <math_IntegerVector.hxx>
 
 #include <Standard_DimensionError.hxx>
 #include <Standard_RangeError.hxx>
 
-
-
-math_IntegerVector::math_IntegerVector(const Standard_Integer First,
-				       const Standard_Integer Last): 
-				       FirstIndex(First),
-				       LastIndex(Last),
-                                       Array(First, Last) {
-
-  Standard_RangeError_Raise_if(First > Last, " ");
-}
-
-math_IntegerVector::math_IntegerVector(const Standard_Integer First,
-				       const Standard_Integer Last,
-				       const Standard_Integer InitialValue): 
-				       FirstIndex(First),
-				       LastIndex(Last),
-                                       Array(First, Last) {
-
-  Standard_RangeError_Raise_if(First > Last, " ");
-  Array.Init(InitialValue);
-}
-
-
-math_IntegerVector::math_IntegerVector(const Standard_Address Tab,
-				       const Standard_Integer First,
-				       const Standard_Integer Last) :
-				       FirstIndex(First),
-				       LastIndex(Last),
-                                       Array(*((const Standard_Integer *)Tab), 
-					     First, Last)
+math_IntegerVector::math_IntegerVector(const Standard_Integer theFirst, const Standard_Integer theLast) :
+  FirstIndex(theFirst),
+  LastIndex(theLast),
+  Array(theFirst, theLast)
 {
-  Standard_RangeError_Raise_if(First > Last, " ");
+  Standard_RangeError_Raise_if(theFirst > theLast, " ");
 }
 
-
-void math_IntegerVector::Init(const Standard_Integer InitialValue)
+math_IntegerVector::math_IntegerVector(const Standard_Integer theFirst,
+                                       const Standard_Integer theLast,
+                                       const Standard_Integer theInitialValue) :
+  FirstIndex(theFirst),
+  LastIndex(theLast),
+  Array(theFirst, theLast)
 {
-  Array.Init(InitialValue);
+  Standard_RangeError_Raise_if(theFirst > theLast, " ");
+  Array.Init(theInitialValue);
 }
 
-
-math_IntegerVector::math_IntegerVector(const math_IntegerVector& Other):
- 
-                                       FirstIndex(Other.FirstIndex),
-                                       LastIndex(Other.LastIndex),
-                                       Array(Other.Array) {}
-
-
-
-void math_IntegerVector::SetFirst(const Standard_Integer First) {
-  
-  Array.SetLower(First);
-  LastIndex = LastIndex - FirstIndex + First;
-  FirstIndex = First;
+math_IntegerVector::math_IntegerVector(const Standard_Address theTab,
+                                       const Standard_Integer theFirst,
+                                       const Standard_Integer theLast) :
+  FirstIndex(theFirst),
+  LastIndex(theLast),
+  Array(theTab, theFirst, theLast)
+{
+  Standard_RangeError_Raise_if(theFirst > theLast, " ");
 }
 
+void math_IntegerVector::Init(const Standard_Integer theInitialValue)
+{
+  Array.Init(theInitialValue);
+}
 
-Standard_Real math_IntegerVector::Norm() const {
+math_IntegerVector::math_IntegerVector(const math_IntegerVector& theOther) :
+  FirstIndex(theOther.FirstIndex),
+  LastIndex(theOther.LastIndex),
+  Array(theOther.Array)
+{
+}
+
+void math_IntegerVector::SetFirst(const Standard_Integer theFirst)
+{
+  Array.SetLower(theFirst);
+  LastIndex = LastIndex - FirstIndex + theFirst;
+  FirstIndex = theFirst;
+}
+
+Standard_Real math_IntegerVector::Norm() const
+{
   Standard_Real Result = 0;
-  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++) {
+  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++)
+  {
     Result = Result + Array(Index) * Array(Index);
   }
   return Sqrt(Result);
 }
 
-     
-Standard_Real math_IntegerVector::Norm2() const {
+Standard_Real math_IntegerVector::Norm2() const
+{
   Standard_Real Result = 0;
-  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++) {
+  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++)
+  {
     Result = Result + Array(Index) * Array(Index);
   }
   return Result;
 }
 
-
-Standard_Integer math_IntegerVector::Max() const {
+Standard_Integer math_IntegerVector::Max() const
+{
   Standard_Integer I=0;
   Standard_Real X = RealFirst();
-  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++) {
-    if(Array(Index) > X) {
+  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++)
+  {
+    if(Array(Index) > X)
+    {
       X = Array(Index);
       I = Index;
     }
@@ -110,12 +104,14 @@ Standard_Integer math_IntegerVector::Max() const {
   return I;
 }
 
-
-Standard_Integer math_IntegerVector::Min() const {
+Standard_Integer math_IntegerVector::Min() const
+{
   Standard_Integer I=0;
   Standard_Real X = RealLast();
-  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++) {
-    if(Array(Index) < X) {
+  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++)
+  {
+    if(Array(Index) < X)
+    {
       X = Array(Index);
       I = Index;
     }
@@ -123,241 +119,231 @@ Standard_Integer math_IntegerVector::Min() const {
   return I;
 }
 
-
-void math_IntegerVector::Invert() {
-  
+void math_IntegerVector::Invert()
+{
   Standard_Integer J;
   Standard_Integer Temp;
-  
-  for(Standard_Integer Index = FirstIndex; 
-      Index <= FirstIndex + Length() / 2 ; Index++) {
-    J = LastIndex + FirstIndex - Index;
-    Temp = Array(Index);
-    Array(Index) = Array(J);
-    Array(J) = Temp;
+
+  for(Standard_Integer Index = FirstIndex; Index <= FirstIndex + Length() / 2 ; Index++)
+  {
+      J = LastIndex + FirstIndex - Index;
+      Temp = Array(Index);
+      Array(Index) = Array(J);
+      Array(J) = Temp;
   }
 }
 
-
-math_IntegerVector math_IntegerVector::Inverse() const {
-  
+math_IntegerVector math_IntegerVector::Inverse() const
+{
   math_IntegerVector Result = *this;
   Result.Invert();
   return Result;
 }
 
+void math_IntegerVector::Set(const Standard_Integer theI1,
+                             const Standard_Integer theI2,
+                             const math_IntegerVector &theV)
+{
+  Standard_DimensionError_Raise_if((theI1 < FirstIndex) || (theI2 > LastIndex) ||
+    (theI1 > theI2) || (theI2 - theI1 + 1 != theV.Length()), " ");
 
-void math_IntegerVector::Set(const Standard_Integer I1, 
-			     const Standard_Integer I2, 
-			     const math_IntegerVector &V) {
-  
-  Standard_DimensionError_Raise_if((I1 < FirstIndex) || 
-				   (I2 > LastIndex) || 
-				   (I1 > I2)         ||
-				   (I2 - I1 + 1 != V.Length()), " ");
-  
-  Standard_Integer I = V.Lower();
-  for(Standard_Integer Index = I1; Index <= I2; Index++) {
-    Array(Index) = V.Array(I);
+  Standard_Integer I = theV.Lower();
+  for(Standard_Integer Index = theI1; Index <= theI2; Index++)
+  {
+    Array(Index) = theV.Array(I);
     I++;
   }
 }
 
-
-void math_IntegerVector::Multiply(const Standard_Integer Right) {
-  
-  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++) {
-    Array(Index) = Array(Index) * Right;
+void math_IntegerVector::Multiply(const Standard_Integer theRight)
+{
+  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++)
+  {
+    Array(Index) = Array(Index) * theRight;
   }
 }
 
+void math_IntegerVector::Add(const math_IntegerVector& theRight)
+{
+  Standard_DimensionError_Raise_if(Length() != theRight.Length(), " ");
 
-void math_IntegerVector::Add(const math_IntegerVector& Right) {
-  
-  Standard_DimensionError_Raise_if(Length() != Right.Length(), " ");
-
-  Standard_Integer I = Right.FirstIndex;
-  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++) {
-    Array(Index) = Array(Index) + Right.Array(I);
+  Standard_Integer I = theRight.FirstIndex;
+  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++)
+  {
+    Array(Index) = Array(Index) + theRight.Array(I);
     I++;
   }
-}    
+}
 
-
-void math_IntegerVector::Subtract(const math_IntegerVector& Right) {
-  
-  Standard_DimensionError_Raise_if(Length() != Right.Length(), " ");
-  Standard_Integer I = Right.FirstIndex;
-  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++) {
-    Array(Index) = Array(Index) - Right.Array(I);
+void math_IntegerVector::Subtract(const math_IntegerVector& theRight)
+{
+  Standard_DimensionError_Raise_if(Length() != theRight.Length(), " ");
+  Standard_Integer I = theRight.FirstIndex;
+  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++)
+  {
+    Array(Index) = Array(Index) - theRight.Array(I);
     I++;
   }
-}    
+}
 
+math_IntegerVector math_IntegerVector::Slice(const Standard_Integer theI1,
+                                             const Standard_Integer theI2) const
+{
+  Standard_DimensionError_Raise_if((theI1 < FirstIndex) || (theI1 > LastIndex) ||
+    (theI2 < FirstIndex) || (theI2 > LastIndex), " ");
 
-math_IntegerVector math_IntegerVector::Slice(const Standard_Integer I1, 
-					     const Standard_Integer I2) const {
-					       
-  Standard_DimensionError_Raise_if((I1 < FirstIndex) || (I1 > LastIndex) ||
-				   (I2 < FirstIndex) || (I2 > LastIndex),
-				   " ");
-
-  if(I2 >= I1) {
-    math_IntegerVector Result(I1, I2);
-    for(Standard_Integer Index = I1; Index <= I2; Index++) {
-      Result.Array(Index) = Array(Index);
-    }        
-    return Result;  
-  }
-  else {
-    math_IntegerVector Result(I2, I1);
-    for(Standard_Integer Index = I1; Index >= I2; Index--) {
+  if(theI2 >= theI1)
+  {
+    math_IntegerVector Result(theI1, theI2);
+    for(Standard_Integer Index = theI1; Index <= theI2; Index++)
+    {
       Result.Array(Index) = Array(Index);
     }
     return Result;
-  }   
+  }
+  else
+  {
+    math_IntegerVector Result(theI2, theI1);
+    for(Standard_Integer Index = theI1; Index >= theI2; Index--)
+    {
+      Result.Array(Index) = Array(Index);
+    }
+    return Result;
+  }
 }
 
-Standard_Integer math_IntegerVector::Multiplied (const math_IntegerVector& Right) const {
-  
+Standard_Integer math_IntegerVector::Multiplied (const math_IntegerVector& theRight) const
+{
   Standard_Integer Result = 0;
-  
-  Standard_DimensionError_Raise_if(Length() != Right.Length(), " ");
 
-  Standard_Integer I = Right.FirstIndex;
-  for(Standard_Integer Index = 0; Index < Length(); Index++) {
-    Result = Result + Array(Index) * Right.Array(I);
+  Standard_DimensionError_Raise_if(Length() != theRight.Length(), " ");
+
+  Standard_Integer I = theRight.FirstIndex;
+  for(Standard_Integer Index = 0; Index < Length(); Index++)
+  {
+    Result = Result + Array(Index) * theRight.Array(I);
     I++;
   }
   return Result;
-}    
+}
 
-
-math_IntegerVector math_IntegerVector::Multiplied (const Standard_Integer Right)const  {
-       
+math_IntegerVector math_IntegerVector::Multiplied (const Standard_Integer theRight)const
+{
   math_IntegerVector Result(FirstIndex, LastIndex);
-  
-  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++) {
-    Result.Array(Index) = Array(Index) * Right;
+
+  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++)
+  {
+    Result.Array(Index) = Array(Index) * theRight;
   }
   return Result;
-}    
-     
-math_IntegerVector math_IntegerVector::TMultiplied (const Standard_Integer Right)const  {
-       
+}
+
+math_IntegerVector math_IntegerVector::TMultiplied (const Standard_Integer theRight) const
+{
   math_IntegerVector Result(FirstIndex, LastIndex);
-  
-  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++) {
-    Result.Array(Index) = Array(Index) * Right;
+
+  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++)
+  {
+    Result.Array(Index) = Array(Index) * theRight;
   }
   return Result;
-}    
-     
+}
+
+math_IntegerVector math_IntegerVector::Added (const math_IntegerVector& theRight) const
+{
+  Standard_DimensionError_Raise_if(Length() != theRight.Length(), " ");
 
-math_IntegerVector math_IntegerVector::Added (const math_IntegerVector& Right) const {
-       
-  Standard_DimensionError_Raise_if(Length() != Right.Length(), " ");
-  
   math_IntegerVector Result(FirstIndex, LastIndex);
-  
-  Standard_Integer I = Right.FirstIndex;
-  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++) {
-    Result.Array(Index) = Array(Index) + Right.Array(I);
+
+  Standard_Integer I = theRight.FirstIndex;
+  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++)
+  {
+    Result.Array(Index) = Array(Index) + theRight.Array(I);
     I++;
   }
   return Result;
-}    
+}
 
-
-     
-math_IntegerVector math_IntegerVector::Opposite() {
-   
+math_IntegerVector math_IntegerVector::Opposite()
+{
   math_IntegerVector Result(FirstIndex, LastIndex);
-  
-  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++) {
+
+  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++)
+  {
     Result.Array(Index) = - Array(Index);
   }
   return Result;
-}    
-     
+}
+
+math_IntegerVector math_IntegerVector::Subtracted (const math_IntegerVector& theRight) const
+{
+  Standard_DimensionError_Raise_if(Length() != theRight.Length(), " ");
 
-math_IntegerVector math_IntegerVector::Subtracted (const math_IntegerVector& Right) const {
-  
-  Standard_DimensionError_Raise_if(Length() != Right.Length(), " ");
-  
   math_IntegerVector Result(FirstIndex, LastIndex);
-  
-  Standard_Integer I = Right.FirstIndex;
-  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++) {
-    Result.Array(Index) = Array(Index) - Right.Array(I);
+
+  Standard_Integer I = theRight.FirstIndex;
+  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++)
+  {
+    Result.Array(Index) = Array(Index) - theRight.Array(I);
     I++;
   }
   return Result;
-}    
-     
-void math_IntegerVector::Add (const math_IntegerVector& Left, 
-			      const math_IntegerVector& Right) {
-  
-  Standard_DimensionError_Raise_if((Length() != Right.Length()) ||
-				   (Right.Length() != Left.Length()), " ");
+}
 
-  Standard_Integer I = Left.FirstIndex;
-  Standard_Integer J = Right.FirstIndex;
-  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++) {
-    Array(Index) = Left.Array(I) + Right.Array(J);
-    I++;
-    J++;
-  }
-}    
-
-
-void math_IntegerVector::Subtract (const math_IntegerVector& Left, 
-				    const math_IntegerVector& Right) {
-  
-  Standard_DimensionError_Raise_if((Length() != Right.Length()) ||
-				   (Right.Length() != Left.Length()), " ");
-  
-  Standard_Integer I = Left.FirstIndex;
-  Standard_Integer J = Right.FirstIndex;
-  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++) {
-    Array(Index) = Left.Array(I) - Right.Array(J);
-    I++;
-    J++;
-  }
-}    
-
-
-void math_IntegerVector::Multiply(const Standard_Integer Left,
-				  const math_IntegerVector& Right) 
+void math_IntegerVector::Add (const math_IntegerVector& theLeft, const math_IntegerVector& theRight)
 {
-  Standard_DimensionError_Raise_if((Length() != Right.Length()),
-				   " ");
-  for(Standard_Integer I = FirstIndex; I <= LastIndex; I++) {
-    Array(I) = Left * Right.Array(I);
+  Standard_DimensionError_Raise_if((Length() != theRight.Length()) ||
+    (theRight.Length() != theLeft.Length()), " ");
+
+  Standard_Integer I = theLeft.FirstIndex;
+  Standard_Integer J = theRight.FirstIndex;
+  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++)
+  {
+    Array(Index) = theLeft.Array(I) + theRight.Array(J);
+    I++;
+    J++;
   }
 }
 
+void math_IntegerVector::Subtract (const math_IntegerVector& theLeft,
+                                   const math_IntegerVector& theRight)
+{
+  Standard_DimensionError_Raise_if((Length() != theRight.Length()) ||
+    (theRight.Length() != theLeft.Length()), " ");
 
-math_IntegerVector& math_IntegerVector::Initialized (const math_IntegerVector& Other) {
+  Standard_Integer I = theLeft.FirstIndex;
+  Standard_Integer J = theRight.FirstIndex;
+  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++)
+  {
+    Array(Index) = theLeft.Array(I) - theRight.Array(J);
+    I++;
+    J++;
+  }
+}
 
-  Standard_DimensionError_Raise_if(Length() != Other.Length(), " ");
+void math_IntegerVector::Multiply(const Standard_Integer theLeft, const math_IntegerVector& theRight)
+{
+  Standard_DimensionError_Raise_if((Length() != theRight.Length()), " ");
+  for(Standard_Integer I = FirstIndex; I <= LastIndex; I++)
+  {
+    Array(I) = theLeft * theRight.Array(I);
+  }
+}
 
-  (Other.Array).Copy(Array);
+math_IntegerVector& math_IntegerVector::Initialized(const math_IntegerVector& theOther)
+{
+  Standard_DimensionError_Raise_if(Length() != theOther.Length(), " ");
+
+  (theOther.Array).Copy(Array);
   return *this;
 }
 
-
-
-void math_IntegerVector::Dump(Standard_OStream& o) const
+void math_IntegerVector::Dump(Standard_OStream& theO) const
 {
-   o << "math_IntegerVector of Range = " << Length() << "\n";
-   for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++) {
-     o << "math_IntegerVector(" << Index << ") = " << Array(Index) << "\n";
-   }
+  theO << "math_IntegerVector of Range = " << Length() << "\n";
+  for(Standard_Integer Index = FirstIndex; Index <= LastIndex; Index++)
+  {
+    theO << "math_IntegerVector(" << Index << ") = " << Array(Index) << "\n";
+  }
 }
 
-
-
-
-
-

src/math/math_IntegerVector.hxx

@@ -0,0 +1,261 @@
+// Copyright (c) 1997-1999 Matra Datavision
+// Copyright (c) 1999-2014 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#ifndef _math_IntegerVector_HeaderFile
+#define _math_IntegerVector_HeaderFile
+
+#include <math_SingleTab.hxx>
+
+class Standard_DimensionError;
+class Standard_DivideByZero;
+class Standard_RangeError;
+class math_Matrix;
+
+//! This class implements the real IntegerVector abstract data type.
+//! IntegerVectors can have an arbitrary range which must be define at
+//! the declaration and cannot be changed after this declaration.
+//! Example:
+//! @code
+//!    math_IntegerVector V1(-3, 5); // an IntegerVector with range [-3..5]
+//! @endcode
+//!
+//! IntegerVector is copied through assignement :
+//! @code
+//!    math_IntegerVector V2( 1, 9);
+//!    ....
+//!    V2 = V1;
+//!    V1(1) = 2.0; // the IntegerVector V2 will not be modified.
+//! @endcode
+//!
+//! The Exception RangeError is raised when trying to access outside
+//! the range of an IntegerVector :
+//! @code
+//!    V1(11) = 0 // --> will raise RangeError;
+//! @endcode
+//!
+//! The Exception DimensionError is raised when the dimensions of two
+//! IntegerVectors are not compatible :
+//! @code
+//!    math_IntegerVector V3(1, 2);
+//!    V3 = V1;    // --> will raise DimensionError;
+//!    V1.Add(V3)  // --> will raise DimensionError;
+//! @endcode
+class math_IntegerVector
+{
+public:
+
+  DEFINE_STANDARD_ALLOC
+
+  //! contructs an IntegerVector in the range [Lower..Upper]
+  Standard_EXPORT math_IntegerVector(const Standard_Integer theFirst, const Standard_Integer theLast);
+
+  //! contructs an IntegerVector in the range [Lower..Upper]
+  //! with all the elements set to theInitialValue.
+  Standard_EXPORT math_IntegerVector(const Standard_Integer theFirst, const Standard_Integer theLast, const Standard_Integer theInitialValue);
+
+  //! Initialize an IntegerVector with all the elements
+  //! set to theInitialValue.
+  Standard_EXPORT void Init(const Standard_Integer theInitialValue);
+
+  //! constructs an IntegerVector in the range [Lower..Upper]
+  //! which share the "c array" theTab.
+  Standard_EXPORT math_IntegerVector(const Standard_Address theTab, const Standard_Integer theFirst, const Standard_Integer theLast);
+
+  //! constructs a copy for initialization.
+  //! An exception is raised if the lengths of the IntegerVectors
+  //! are different.
+  Standard_EXPORT math_IntegerVector(const math_IntegerVector& theOther);
+
+  //! returns the length of an IntegerVector
+  inline Standard_Integer Length() const
+  {
+    return LastIndex - FirstIndex +1;
+  }
+
+  //! returns the value of the Lower index of an IntegerVector.
+  inline Standard_Integer Lower() const
+  {
+    return FirstIndex;
+  }
+
+  //! returns the value of the Upper index of an IntegerVector.
+  inline Standard_Integer Upper() const
+  {
+    return LastIndex;
+  }
+
+  //! returns the value of the norm of an IntegerVector.
+  Standard_EXPORT Standard_Real Norm() const;
+
+  //! returns the value of the square of the norm of an IntegerVector.
+  Standard_EXPORT Standard_Real Norm2() const;
+
+  //! returns the value of the Index of the maximum element of an IntegerVector.
+  Standard_EXPORT Standard_Integer Max() const;
+
+  //! returns the value of the Index of the minimum element of an IntegerVector.
+  Standard_EXPORT Standard_Integer Min() const;
+
+  //! inverses an IntegerVector.
+  Standard_EXPORT void Invert();
+
+  //! returns the inverse IntegerVector of an IntegerVector.
+  Standard_EXPORT math_IntegerVector Inverse() const;
+
+  //! sets an IntegerVector from "theI1" to "theI2" to the IntegerVector "theV";
+  //! An exception is raised if "theI1" is less than "LowerIndex" or "theI2" is greater than "UpperIndex" or "theI1" is greater than "theI2".
+  //! An exception is raised if "theI2-theI1+1" is different from the Length of "theV".
+  Standard_EXPORT void Set(const Standard_Integer theI1, const Standard_Integer theI2, const math_IntegerVector& theV);
+
+  //! slices the values of the IntegerVector between "theI1" and "theI2":
+  //! Example: [2, 1, 2, 3, 4, 5] becomes [2, 4, 3, 2, 1, 5] between 2 and 5.
+  //! An exception is raised if "theI1" is less than "LowerIndex" or "theI2" is greater than "UpperIndex".
+  Standard_EXPORT math_IntegerVector Slice(const Standard_Integer theI1, const Standard_Integer theI2) const;
+
+  //! returns the product of an IntegerVector by an integer value.
+  Standard_EXPORT void Multiply(const Standard_Integer theRight);
+
+  void operator *=(const Standard_Integer theRight)
+  {
+    Multiply(theRight);
+  }
+
+  //! returns the product of an IntegerVector by an integer value.
+  Standard_EXPORT math_IntegerVector Multiplied(const Standard_Integer theRight) const;
+
+  math_IntegerVector operator*(const Standard_Integer theRight) const
+  {
+    return Multiplied(theRight);
+  }
+
+  //! returns the product of a vector and a real value.
+  Standard_EXPORT math_IntegerVector TMultiplied(const Standard_Integer theRight) const;
+
+  friend inline math_IntegerVector operator* (const Standard_Integer theLeft, const math_IntegerVector& theRight)
+  {
+    return theRight.Multiplied(theLeft);
+  }
+
+  //! adds the IntegerVector "theRight" to an IntegerVector.
+  //! An exception is raised if the IntegerVectors have not the same length.
+  //! An exception is raised if the lengths are not equal.
+  Standard_EXPORT void Add(const math_IntegerVector& theRight);
+
+  void operator +=(const math_IntegerVector& theRight)
+  {
+    Add(theRight);
+  }
+
+  //! adds the IntegerVector "theRight" to an IntegerVector.
+  //! An exception is raised if the IntegerVectors have not the same length.
+  //! An exception is raised if the lengths are not equal.
+  Standard_EXPORT math_IntegerVector Added(const math_IntegerVector& theRight) const;
+
+  math_IntegerVector operator+(const math_IntegerVector& theRight) const
+  {
+    return Added(theRight);
+  }
+
+  //! sets an IntegerVector to the sum of the IntegerVector
+  //! "theLeft" and the IntegerVector "theRight".
+  //! An exception is raised if the lengths are different.
+  Standard_EXPORT void Add(const math_IntegerVector& theLeft, const math_IntegerVector& theRight);
+
+  //! sets an IntegerVector to the substraction of "theRight" from "theLeft".
+  //! An exception is raised if the IntegerVectors have not the same length.
+  Standard_EXPORT void Subtract(const math_IntegerVector& theLeft, const math_IntegerVector& theRight);
+
+  //! accesses (in read or write mode) the value of index theNum of an IntegerVector.
+  inline Standard_Integer& Value(const Standard_Integer theNum) const
+  {
+    Standard_RangeError_Raise_if(theNum < FirstIndex || theNum > LastIndex, " ");
+    return Array(theNum);
+  }
+
+  Standard_EXPORT Standard_Integer& operator()(const Standard_Integer theNum) const
+  {
+    return Value(theNum);
+  }
+
+  //! Initialises an IntegerVector by copying "theOther".
+  //! An exception is raised if the Lengths are different.
+  Standard_EXPORT math_IntegerVector& Initialized(const math_IntegerVector& theOther);
+
+  math_IntegerVector& operator=(const math_IntegerVector& theOther) 
+  {
+    return Initialized(theOther);
+  }
+
+  //! returns the inner product of 2 IntegerVectors.
+  //! An exception is raised if the lengths are not equal.
+  Standard_EXPORT Standard_Integer Multiplied(const math_IntegerVector& theRight) const;
+
+  Standard_Integer operator*(const math_IntegerVector& theRight) const
+  {
+    return Multiplied(theRight);
+  }
+
+  //! returns the opposite of an IntegerVector.
+  Standard_EXPORT math_IntegerVector Opposite();
+
+  math_IntegerVector operator-()
+  {
+    return Opposite();
+  }
+
+  //! returns the subtraction of "theRight" from "me".
+  //! An exception is raised if the IntegerVectors have not the same length.
+  Standard_EXPORT void Subtract(const math_IntegerVector& theRight);
+
+  void operator-=(const math_IntegerVector& theRight)
+  {
+    Subtract(theRight);
+  }
+
+  //! returns the subtraction of "theRight" from "me".
+  //! An exception is raised if the IntegerVectors have not the same length.
+  Standard_EXPORT math_IntegerVector Subtracted(const math_IntegerVector& theRight) const;
+
+  math_IntegerVector operator-(const math_IntegerVector& theRight) const
+  {
+    return Subtracted(theRight);
+  }
+
+  //! returns the multiplication of an integer by an IntegerVector.
+  Standard_EXPORT void Multiply(const Standard_Integer theLeft,const math_IntegerVector& theRight);
+
+  //! Prints on the stream theO information on the current state of the object.
+  //! Is used to redefine the operator <<.
+  Standard_EXPORT void Dump(Standard_OStream& theO) const;
+
+  friend inline Standard_OStream& operator<<(Standard_OStream& theO, const math_IntegerVector& theVec)
+  {
+    theVec.Dump(theO);
+    return theO;
+  }
+
+protected:
+
+  //! is used internally to set the Lower value of the IntegerVector.
+  void SetFirst(const Standard_Integer theFirst);
+
+private:
+
+  Standard_Integer FirstIndex;
+  Standard_Integer LastIndex;
+  math_SingleTab<Standard_Integer> Array;
+};
+
+#endif
+

src/math/math_Matrix.cdl

@@ -51,7 +51,7 @@ class Matrix from math
 
 
 uses Vector from math,
-     DoubleTabOfReal from math,
+     DoubleTab from math,
      OStream from Standard
 
 raises DimensionError from Standard, RangeError from Standard, 
@@ -562,7 +562,7 @@ LowerRowIndex:      Integer;
 UpperRowIndex:      Integer;
 LowerColIndex:      Integer;
 UpperColIndex:      Integer;
-Array:              DoubleTabOfReal;
+Array:              DoubleTab;
 
 friends
 class Vector from math

src/math/math_Matrix.cxx

@@ -85,9 +85,7 @@ math_Matrix::math_Matrix (const Standard_Address Tab,
 			  UpperRowIndex(UpperRow),
 			  LowerColIndex(LowerCol),
 			  UpperColIndex(UpperCol),
-			  Array(*((const Standard_Real *)Tab), 
-				LowerRow, UpperRow,
-				LowerCol, UpperCol) 
+			  Array(Tab, LowerRow, UpperRow, LowerCol, UpperCol) 
 { 
   
   Standard_RangeError_Raise_if((LowerRow > UpperRow) ||

src/math/math_MultipleVarFunction.cdl

@@ -20,6 +20,8 @@ deferred class MultipleVarFunction from math
 uses Vector from math
 
 is
+    Delete(me:out) is virtual;
+    ---C++: alias "Standard_EXPORT virtual ~math_MultipleVarFunction(){Delete();}"
 
     NbVariables(me)
     	---Purpose:

src/math/math_MultipleVarFunction.cxx

@@ -16,3 +16,6 @@
 #include <math_MultipleVarFunction.ixx>
 
 Standard_Integer math_MultipleVarFunction::GetStateNumber() { return 0; }
+
+void math_MultipleVarFunction::Delete()
+{}

src/math/math_MultipleVarFunctionWithGradient.cdl

@@ -24,7 +24,7 @@ uses Vector from math
 
 is
 
-    Delete(me:out) is virtual;
+    Delete(me:out) is redefined virtual;
     ---C++: alias "Standard_EXPORT virtual ~math_MultipleVarFunctionWithGradient(){Delete();}"
     
     NbVariables(me)

src/math/math_MultipleVarFunctionWithHessian.cdl

@@ -24,6 +24,9 @@ uses Matrix from math,
 
 is
 
+    Delete(me:out) is redefined virtual;
+    ---C++: alias "Standard_EXPORT virtual ~math_MultipleVarFunctionWithHessian(){Delete();}"
+
     NbVariables(me)
     	---Purpose: returns the number of variables of the function.
 

src/math/math_MultipleVarFunctionWithHessian.cxx

@@ -15,3 +15,6 @@
 // commercial license or contractual agreement.
 
 #include <math_MultipleVarFunctionWithHessian.ixx>
+
+void math_MultipleVarFunctionWithHessian::Delete()
+{}

src/math/math_NewtonMinimum.cxx

@@ -143,12 +143,19 @@ void math_NewtonMinimum::Perform(math_MultipleVarFunctionWithHessian& F,
     }
    
     LU.Solve(TheGradient, TheStep);
-    *suivant = *precedent - TheStep;
+    Standard_Boolean hasProblem = Standard_False;
+    do
+    {
+      *suivant = *precedent - TheStep;
 
+      //  Gestion de la convergence
+      hasProblem = !(F.Value(*suivant, TheMinimum));
 
-    //  Gestion de la convergence
-
-    F.Value(*suivant, TheMinimum);
+      if (hasProblem)
+      {
+        TheStep /= 2.0;
+      }
+    } while (hasProblem);
 
     if (IsConverged()) { NbConv++; }
     else               { NbConv=0; }

src/math/math_SingleTab.hxx

@@ -0,0 +1,115 @@
+// Copyright (c) 1997-1999 Matra Datavision
+// Copyright (c) 1999-2014 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#ifndef _math_SingleTab_HeaderFile
+#define _math_SingleTab_HeaderFile
+
+#include <math_Memory.hxx>
+#include <Standard_OutOfRange.hxx>
+#include <Standard_Failure.hxx>
+
+static const Standard_Integer aLengthOfBuf = 512;
+
+template<class T> class math_SingleTab
+{
+public:
+
+  DEFINE_STANDARD_ALLOC
+
+  math_SingleTab(const Standard_Integer LowerIndex, const Standard_Integer UpperIndex) :
+      Addr(Buf),
+      isAllocated(UpperIndex - LowerIndex + 1 > aLengthOfBuf),
+      First(LowerIndex), Last(UpperIndex)
+  {
+    T* TheAddr = !isAllocated? Buf :
+      (T*) Standard::Allocate((Last-First+1) * sizeof(T));
+    Addr = (Standard_Address) (TheAddr - First);
+  }
+
+  math_SingleTab(const Standard_Address Tab, const Standard_Integer LowerIndex, const Standard_Integer UpperIndex) :
+      Addr((void*)((const T*)Tab - LowerIndex)), 
+      isAllocated(Standard_False),
+      First(LowerIndex), Last(UpperIndex)
+  {
+  }
+
+  void Init(const T InitValue)
+  {
+    for(Standard_Integer i = First; i<= Last; i++)
+    {
+      ((T*)Addr)[i] = InitValue;
+    }
+  }
+
+  math_SingleTab(const math_SingleTab& Other) :
+      isAllocated(Other.Last - Other.First + 1 > aLengthOfBuf),
+      First(Other.First),
+      Last(Other.Last)
+  {
+    T* TheAddr = !isAllocated? Buf : (T*) Standard::Allocate((Last-First+1) * sizeof(T));
+    Addr = (Standard_Address) (TheAddr - First);
+    T* TheOtherAddr = (T*) Other.Addr;
+    memmove((void*) TheAddr, (const void*) (TheOtherAddr + First), (size_t)(Last - First + 1) * sizeof(T));
+  }
+
+  inline void Copy(math_SingleTab& Other) const
+  {
+    memmove((void*) (((T*)Other.Addr) + Other.First),
+      (const void*) (((T*)Addr) + First),
+      (size_t)(Last - First + 1) * sizeof(T));
+  }
+
+  void SetLower(const Standard_Integer LowerIndex)
+  {
+    T* TheAddr = (T*) Addr;
+    Addr = (Standard_Address) (TheAddr + First - LowerIndex);
+    Last = Last - First + LowerIndex;
+    First = LowerIndex;
+  }
+
+  inline T& Value(const Standard_Integer Index) const
+  {
+    return ((T*)Addr)[Index];
+  }
+
+  T& operator()(const Standard_Integer Index) const
+  {
+    return Value(Index);
+  }
+  
+  void Free()
+  {
+    if(isAllocated)
+    {
+      Standard_Address it = (Standard_Address)&((T*)Addr)[First];
+      Standard::Free(it);
+      Addr = 0;
+    }
+  }
+
+  ~math_SingleTab()
+  {
+    Free();
+  }
+
+private:
+  
+  Standard_Address Addr;
+  T Buf[aLengthOfBuf];
+  Standard_Boolean isAllocated;
+  Standard_Integer First;
+  Standard_Integer Last;
+};
+
+#endif

src/math/math_Vector.cxx

@@ -14,7 +14,7 @@
 
 #include <stdio.h>
 
-#include <math_Vector.ixx>
+#include <math_Vector.hxx>
 #include <math_Matrix.hxx>
 
 #include <Standard_DimensionError.hxx>
@@ -22,84 +22,85 @@
 #include <Standard_RangeError.hxx>
 #include <Standard_NullValue.hxx>
 
-
-math_Vector::math_Vector(const Standard_Integer Lower,
-			 const Standard_Integer Upper):
-			 
-			 LowerIndex(Lower),
-			 UpperIndex(Upper),
-			 Array(Lower,Upper) {
-			   Standard_RangeError_Raise_if(Lower > Upper, "");		       
-			 }
-
-math_Vector::math_Vector(const Standard_Integer Lower,
-			 const Standard_Integer Upper,
-			 const Standard_Real    InitialValue):
-			 
-			 LowerIndex(Lower),
-			 UpperIndex(Upper),
-			 Array(Lower,Upper) 
+math_Vector::math_Vector(const Standard_Integer theLower, const Standard_Integer theUpper) :
+  LowerIndex(theLower),
+  UpperIndex(theUpper),
+  Array(theLower,theUpper)
 {
-  Standard_RangeError_Raise_if(Lower > Upper, "");		       
-  Array.Init(InitialValue);
+  Standard_RangeError_Raise_if(theLower > theUpper, "");
 }
 
-math_Vector::math_Vector(const Standard_Address Tab,
-			 const Standard_Integer Lower,
-			 const Standard_Integer Upper) :
-			 
-			 LowerIndex(Lower),
-			 UpperIndex(Upper),
-			 Array(*((const Standard_Real *)Tab), Lower,Upper) 
+math_Vector::math_Vector(const Standard_Integer theLower,
+                         const Standard_Integer theUpper,
+                         const Standard_Real    theInitialValue):
+  LowerIndex(theLower),
+  UpperIndex(theUpper),
+  Array(theLower,theUpper)
 {
-  Standard_RangeError_Raise_if((Lower > Upper) , "");		       
+  Standard_RangeError_Raise_if(theLower > theUpper, "");
+  Array.Init(theInitialValue);
 }
 
-void math_Vector::Init(const Standard_Real InitialValue) {
-  Array.Init(InitialValue);
+math_Vector::math_Vector(const Standard_Address theTab,
+                         const Standard_Integer theLower,
+                         const Standard_Integer theUpper) :
+  LowerIndex(theLower),
+  UpperIndex(theUpper),
+  Array(theTab, theLower,theUpper)
+{
+  Standard_RangeError_Raise_if((theLower > theUpper) , "");
 }
 
-math_Vector::math_Vector(const math_Vector& Other): 
-
-LowerIndex(Other.LowerIndex),
-UpperIndex(Other.UpperIndex),
-Array(Other.Array) {}
-
-
-void math_Vector::SetLower(const Standard_Integer Lower) {
-
-  Array.SetLower(Lower);
-  UpperIndex = UpperIndex - LowerIndex + Lower;
-  LowerIndex = Lower;
+void math_Vector::Init(const Standard_Real theInitialValue)
+{
+  Array.Init(theInitialValue);
 }
 
-Standard_Real math_Vector::Norm() const {
+math_Vector::math_Vector(const math_Vector& theOther) :
+  LowerIndex(theOther.LowerIndex),
+  UpperIndex(theOther.UpperIndex),
+  Array(theOther.Array)
+{
+}
 
+void math_Vector::SetLower(const Standard_Integer theLower)
+{
+  Array.SetLower(theLower);
+  UpperIndex = UpperIndex - LowerIndex + theLower;
+  LowerIndex = theLower;
+}
+
+Standard_Real math_Vector::Norm() const
+{
   Standard_Real Result = 0;
 
-  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) {
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
     Result = Result + Array(Index) * Array(Index);
   }
   return Sqrt(Result);
 }
 
-Standard_Real math_Vector::Norm2() const {
-
+Standard_Real math_Vector::Norm2() const
+{
   Standard_Real Result = 0;
 
-  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) {
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
     Result = Result + Array(Index) * Array(Index);
   }
   return Result;
 }
 
-Standard_Integer math_Vector::Max() const {
-
+Standard_Integer math_Vector::Max() const
+{
   Standard_Integer I=0;
   Standard_Real X = RealFirst();
 
-  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) {
-    if(Array(Index) > X) {
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
+    if(Array(Index) > X)
+    {
       X = Array(Index);
       I = Index;
     }
@@ -107,13 +108,15 @@ Standard_Integer math_Vector::Max() const {
   return I;
 }
 
-Standard_Integer math_Vector::Min() const {
-
+Standard_Integer math_Vector::Min() const
+{
   Standard_Integer I=0;
   Standard_Real X = RealLast();
 
-  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) {
-    if(Array(Index) < X) {
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
+    if(Array(Index) < X)
+    {
       X = Array(Index);
       I = Index;
     }
@@ -121,45 +124,45 @@ Standard_Integer math_Vector::Min() const {
   return I;
 }
 
-void math_Vector::Set(const Standard_Integer I1, 
-		      const Standard_Integer I2, 
-		      const math_Vector &V) {
+void math_Vector::Set(const Standard_Integer theI1,
+                      const Standard_Integer theI2,
+                      const math_Vector &theV)
+{
+  Standard_RangeError_Raise_if((theI1 < LowerIndex) || (theI2 > UpperIndex) ||
+    (theI1 > theI2) || (theI2 - theI1 + 1 != theV.Length()), "");
 
-  Standard_RangeError_Raise_if((I1 < LowerIndex) || 
-			       (I2 > UpperIndex) || 
-			       (I1 > I2)         ||
-			       (I2 - I1 + 1 != V.Length()), "");
-  
-  Standard_Integer I = V.Lower();
-  for(Standard_Integer Index = I1; Index <= I2; Index++) {
-    Array(Index) = V.Array(I);
+  Standard_Integer I = theV.Lower();
+  for(Standard_Integer Index = theI1; Index <= theI2; Index++)
+  {
+    Array(Index) = theV.Array(I);
     I++;
   }
 }
 
-void math_Vector::Normalize() {
-
+void math_Vector::Normalize()
+{
   Standard_Real Result = Norm();
   Standard_NullValue_Raise_if((Result <= RealEpsilon()), "");
-  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) {
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
     Array(Index) = Array(Index) / Result;
   }
 }
 
-math_Vector math_Vector::Normalized() const {
-
+math_Vector math_Vector::Normalized() const
+{
   math_Vector Result = *this;
 
   Result.Normalize();
   return Result;
 }
 
-void math_Vector::Invert() {
+void math_Vector::Invert()
+{
   Standard_Integer J;
   Standard_Real Temp;
-  for(Standard_Integer Index = LowerIndex; 
-//      Index <= LowerIndex + (Length()) >> 1 ; Index++) {
-      Index <= (LowerIndex + Length()) >> 1 ; Index++) {
+  for(Standard_Integer Index = LowerIndex; Index <= (LowerIndex + Length()) >> 1 ; Index++)
+  {
     J = UpperIndex + LowerIndex - Index;
     Temp = Array(Index);
     Array(Index) = Array(J);
@@ -167,316 +170,310 @@ void math_Vector::Invert() {
   }
 }
 
-math_Vector math_Vector::Inverse() const {
+math_Vector math_Vector::Inverse() const
+{
   math_Vector Result = *this;
   Result.Invert();
   return Result;
 }
 
-math_Vector math_Vector::Multiplied(const Standard_Real Right) const{
-
+math_Vector math_Vector::Multiplied(const Standard_Real theRight) const
+{
   math_Vector Result (LowerIndex, UpperIndex);
 
-  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) {
-    Result.Array(Index) = Array(Index) * Right;
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
+    Result.Array(Index) = Array(Index) * theRight;
   }
   return Result;
 }
 
-math_Vector math_Vector::TMultiplied(const Standard_Real Right) const{
-
+math_Vector math_Vector::TMultiplied(const Standard_Real theRight) const
+{
   math_Vector Result (LowerIndex, UpperIndex);
 
-  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) {
-    Result.Array(Index) = Array(Index) * Right;
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
+    Result.Array(Index) = Array(Index) * theRight;
   }
   return Result;
 }
 
-
-void math_Vector::Multiply(const Standard_Real Right) {
-
-  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) {
-    Array(Index) = Array(Index) * Right;
+void math_Vector::Multiply(const Standard_Real theRight)
+{
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
+    Array(Index) = Array(Index) * theRight;
   }
 }
 
+void math_Vector::Divide(const Standard_Real theRight)
+{
+  Standard_DivideByZero_Raise_if(Abs(theRight) <= RealEpsilon(), "");
 
-void math_Vector::Divide(const Standard_Real Right) {
-
-  Standard_DivideByZero_Raise_if(Abs(Right) <= RealEpsilon(), "");
-
-  for(Standard_Integer Index =LowerIndex; Index <=UpperIndex; Index++) {
-    Array(Index) = Array(Index) / Right;
+  for(Standard_Integer Index =LowerIndex; Index <=UpperIndex; Index++)
+  {
+    Array(Index) = Array(Index) / theRight;
   }
 }
 
-
-math_Vector math_Vector::Divided (const Standard_Real Right) const {
-  
-  Standard_DivideByZero_Raise_if(Abs(Right) <= RealEpsilon(), "");
-  math_Vector temp = Multiplied(1./Right);
+math_Vector math_Vector::Divided (const Standard_Real theRight) const
+{
+  Standard_DivideByZero_Raise_if(Abs(theRight) <= RealEpsilon(), "");
+  math_Vector temp = Multiplied(1./theRight);
   return temp;
 }
 
-
-void math_Vector::Add(const math_Vector& Right) {
-
-  Standard_DimensionError_Raise_if(Length() != Right.Length(), "");
-
-  Standard_Integer I = Right.LowerIndex;
-  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) {
-    Array(Index) = Array(Index) + Right.Array(I);
-    I++;
-  }
-}    
-
-math_Vector math_Vector::Added(const math_Vector& Right) const{
-
-  Standard_DimensionError_Raise_if(Length() != Right.Length(), "");
-
-  math_Vector Result(LowerIndex, UpperIndex);
-  
-  Standard_Integer I = Right.LowerIndex;
-  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) {
-    Result.Array(Index) = Array(Index) + Right.Array(I);
-    I++;
-  }
-  return Result;
-}    
-
-
-
-void math_Vector::Subtract(const math_Vector& Right) {
-
-  Standard_DimensionError_Raise_if(Length() != Right.Length(), "");
-
-  Standard_Integer I = Right.LowerIndex;
-  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) {
-    Array(Index) = Array(Index) - Right.Array(I);
-    I++;
-  }
-}    
-
-
-math_Vector math_Vector::Subtracted (const math_Vector& Right) const {
-
-  Standard_DimensionError_Raise_if(Length() != Right.Length(), "");
-  
-  math_Vector Result(LowerIndex, UpperIndex);
-  
-  Standard_Integer I = Right.LowerIndex;
-  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) {
-    Result.Array(Index) = Array(Index) - Right.Array(I);
-    I++;
-  }
-  return Result;
-}    
-
-
-math_Vector math_Vector::Slice(const Standard_Integer I1, 
-			       const Standard_Integer I2) const 
+void math_Vector::Add(const math_Vector& theRight)
 {
-  
-  Standard_RangeError_Raise_if((I1 < LowerIndex) ||
-			       (I1 > UpperIndex) ||
-			       (I2 < LowerIndex) ||
-			       (I2 > UpperIndex) , "");
-  
+  Standard_DimensionError_Raise_if(Length() != theRight.Length(), "");
 
-  if(I2 >= I1) {
-    math_Vector Result(I1, I2);
-    for(Standard_Integer Index = I1; Index <= I2; Index++) {
-      Result.Array(Index) = Array(Index);
-    }        
-    return Result;  
+  Standard_Integer I = theRight.LowerIndex;
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
+    Array(Index) = Array(Index) + theRight.Array(I);
+    I++;
   }
-  else {
-    math_Vector Result(I2, I1);
-    for(Standard_Integer Index = I1; Index >= I2; Index--) {
+}
+
+math_Vector math_Vector::Added(const math_Vector& theRight) const
+{
+  Standard_DimensionError_Raise_if(Length() != theRight.Length(), "");
+
+  math_Vector Result(LowerIndex, UpperIndex);
+
+  Standard_Integer I = theRight.LowerIndex;
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
+    Result.Array(Index) = Array(Index) + theRight.Array(I);
+    I++;
+  }
+  return Result;
+}
+
+void math_Vector::Subtract(const math_Vector& theRight)
+{
+  Standard_DimensionError_Raise_if(Length() != theRight.Length(), "");
+
+  Standard_Integer I = theRight.LowerIndex;
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
+    Array(Index) = Array(Index) - theRight.Array(I);
+    I++;
+  }
+}
+
+math_Vector math_Vector::Subtracted (const math_Vector& theRight) const
+{
+  Standard_DimensionError_Raise_if(Length() != theRight.Length(), "");
+
+  math_Vector Result(LowerIndex, UpperIndex);
+
+  Standard_Integer I = theRight.LowerIndex;
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
+    Result.Array(Index) = Array(Index) - theRight.Array(I);
+    I++;
+  }
+  return Result;
+}
+
+math_Vector math_Vector::Slice(const Standard_Integer theI1, const Standard_Integer theI2) const
+{
+  Standard_RangeError_Raise_if((theI1 < LowerIndex) || (theI1 > UpperIndex) || (theI2 < LowerIndex) || (theI2 > UpperIndex) , "");
+
+  if(theI2 >= theI1)
+  {
+    math_Vector Result(theI1, theI2);
+    for(Standard_Integer Index = theI1; Index <= theI2; Index++)
+    {
       Result.Array(Index) = Array(Index);
     }
     return Result;
-  }   
+  }
+  else
+  {
+    math_Vector Result(theI2, theI1);
+    for(Standard_Integer Index = theI1; Index >= theI2; Index--)
+    {
+      Result.Array(Index) = Array(Index);
+    }
+    return Result;
+  }
 }
 
+void math_Vector::Add (const math_Vector& theLeft, const math_Vector& theRight)
+{
+  Standard_DimensionError_Raise_if((Length() != theRight.Length()) || (theRight.Length() != theLeft.Length()), "");
 
-void math_Vector::Add (const math_Vector& Left, const math_Vector& Right) {
-
-  Standard_DimensionError_Raise_if((Length() != Right.Length()) ||
-				   (Right.Length() != Left.Length()), "");
-
-  
-  Standard_Integer I = Left.LowerIndex;
-  Standard_Integer J = Right.LowerIndex;
-  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) {
-    Array(Index) = Left.Array(I) + Right.Array(J);
+  Standard_Integer I = theLeft.LowerIndex;
+  Standard_Integer J = theRight.LowerIndex;
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
+    Array(Index) = theLeft.Array(I) + theRight.Array(J);
     I++;
     J++;
   }
-}    
+}
 
-void math_Vector::Subtract (const math_Vector& Left, 
-			    const math_Vector& Right) {
+void math_Vector::Subtract (const math_Vector& theLeft, const math_Vector& theRight)
+{
+  Standard_DimensionError_Raise_if((Length() != theRight.Length()) || (theRight.Length() != theLeft.Length()), "");
 
-  Standard_DimensionError_Raise_if((Length() != Right.Length()) ||
-				   (Right.Length() != Left.Length()), "");
-  
-  Standard_Integer I = Left.LowerIndex;
-  Standard_Integer J = Right.LowerIndex;
-  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) {
-    Array(Index) = Left.Array(I) - Right.Array(J);
+  Standard_Integer I = theLeft.LowerIndex;
+  Standard_Integer J = theRight.LowerIndex;
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
+    Array(Index) = theLeft.Array(I) - theRight.Array(J);
     I++;
     J++;
   }
-}    
+}
 
-void math_Vector::Multiply(const math_Matrix& Left, 
-			   const math_Vector& Right) {
-
-  Standard_DimensionError_Raise_if((Length() != Left.RowNumber()) ||
-				   (Left.ColNumber() != Right.Length()),
-				   "");
+void math_Vector::Multiply(const math_Matrix& theLeft, const math_Vector& theRight)
+{
+  Standard_DimensionError_Raise_if((Length() != theLeft.RowNumber()) ||
+    (theLeft.ColNumber() != theRight.Length()), "");
 
   Standard_Integer Index = LowerIndex;
-  for(Standard_Integer I = Left.LowerRowIndex; I <= Left.UpperRowIndex; I++) {
+  for(Standard_Integer I = theLeft.LowerRowIndex; I <= theLeft.UpperRowIndex; I++)
+  {
     Array(Index) = 0.0;
-    Standard_Integer K = Right.LowerIndex;
-    for(Standard_Integer J = Left.LowerColIndex; J <= Left.UpperColIndex; J++) {
-      Array(Index) = Array(Index) + Left.Array(I, J) * Right.Array(K);
+    Standard_Integer K = theRight.LowerIndex;
+    for(Standard_Integer J = theLeft.LowerColIndex; J <= theLeft.UpperColIndex; J++)
+    {
+      Array(Index) = Array(Index) + theLeft.Array(I, J) * theRight.Array(K);
       K++;
     }
     Index++;
   }
-}    
+}
 
-void math_Vector::Multiply(const math_Vector& Left,
-			   const math_Matrix& Right) {
-
-  Standard_DimensionError_Raise_if((Length() != Right.ColNumber()) ||
-				   (Left.Length() != Right.RowNumber()),
-				   "");
+void math_Vector::Multiply(const math_Vector& theLeft, const math_Matrix& theRight)
+{
+  Standard_DimensionError_Raise_if((Length() != theRight.ColNumber()) ||
+    (theLeft.Length() != theRight.RowNumber()), "");
 
   Standard_Integer Index = LowerIndex;
-  for(Standard_Integer J = Right.LowerColIndex; J <= Right.UpperColIndex; J++) {
+  for(Standard_Integer J = theRight.LowerColIndex; J <= theRight.UpperColIndex; J++)
+  {
     Array(Index) = 0.0;
-    Standard_Integer K = Left.LowerIndex;
-    for(Standard_Integer I = Right.LowerRowIndex; I <= Right.UpperRowIndex; I++) {
-      Array(Index) = Array(Index) + Left.Array(K) * Right.Array(I, J);
+    Standard_Integer K = theLeft.LowerIndex;
+    for(Standard_Integer I = theRight.LowerRowIndex; I <= theRight.UpperRowIndex; I++)
+    {
+      Array(Index) = Array(Index) + theLeft.Array(K) * theRight.Array(I, J);
       K++;
     }
     Index++;
   }
-}    
+}
 
-void math_Vector::TMultiply(const math_Matrix& TLeft,
-			    const math_Vector&  Right) {
-
-  Standard_DimensionError_Raise_if((Length() != TLeft.ColNumber()) ||
-				   (TLeft.RowNumber() != Right.Length()),
-				   "");
+void math_Vector::TMultiply(const math_Matrix& theTLeft, const math_Vector&  theRight)
+{
+  Standard_DimensionError_Raise_if((Length() != theTLeft.ColNumber()) ||
+    (theTLeft.RowNumber() != theRight.Length()), "");
 
   Standard_Integer Index = LowerIndex;
-  for(Standard_Integer I = TLeft.LowerColIndex; I <= TLeft.UpperColIndex; I++) {
+  for(Standard_Integer I = theTLeft.LowerColIndex; I <= theTLeft.UpperColIndex; I++)
+  {
     Array(Index) = 0.0;
-    Standard_Integer K = Right.LowerIndex;
-    for(Standard_Integer J = TLeft.LowerRowIndex; J <= TLeft.UpperRowIndex; J++) {
-      Array(Index) = Array(Index) + TLeft.Array(J, I) * Right.Array(K);
+    Standard_Integer K = theRight.LowerIndex;
+    for(Standard_Integer J = theTLeft.LowerRowIndex; J <= theTLeft.UpperRowIndex; J++)
+    {
+      Array(Index) = Array(Index) + theTLeft.Array(J, I) * theRight.Array(K);
       K++;
     }
     Index++;
   }
-}    
+}
 
-void math_Vector::TMultiply(const math_Vector&  Left,
-			    const math_Matrix& TRight) {
-
-  Standard_DimensionError_Raise_if((Length() != TRight.RowNumber()) ||
-				   (Left.Length() != TRight.ColNumber()),
-				   "");
+void math_Vector::TMultiply(const math_Vector&  theLeft, const math_Matrix& theTRight)
+{
+  Standard_DimensionError_Raise_if((Length() != theTRight.RowNumber()) ||
+    (theLeft.Length() != theTRight.ColNumber()), "");
 
   Standard_Integer Index = LowerIndex;
-  for(Standard_Integer J = TRight.LowerRowIndex; J <= TRight.UpperRowIndex; J++) {
+  for(Standard_Integer J = theTRight.LowerRowIndex; J <= theTRight.UpperRowIndex; J++)
+  {
     Array(Index) = 0.0;
-    Standard_Integer K = Left.LowerIndex;
-    for(Standard_Integer I = TRight.LowerColIndex; 
-	I <= TRight.UpperColIndex; I++) {
-      Array(Index) = Array(Index) + Left.Array(K) * TRight.Array(J, I);
-      K++;
+    Standard_Integer K = theLeft.LowerIndex;
+    for(Standard_Integer I = theTRight.LowerColIndex;
+      I <= theTRight.UpperColIndex; I++)
+    {
+        Array(Index) = Array(Index) + theLeft.Array(K) * theTRight.Array(J, I);
+        K++;
     }
     Index++;
   }
-}    
+}
 
-
-
-
-Standard_Real math_Vector::Multiplied(const math_Vector& Right) const{
+Standard_Real math_Vector::Multiplied(const math_Vector& theRight) const
+{
   Standard_Real Result = 0;
 
-  Standard_DimensionError_Raise_if(Length() != Right.Length(), "");
+  Standard_DimensionError_Raise_if(Length() != theRight.Length(), "");
 
-  Standard_Integer I = Right.LowerIndex;
-  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) {
-    Result = Result + Array(Index) * Right.Array(I);
+  Standard_Integer I = theRight.LowerIndex;
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
+    Result = Result + Array(Index) * theRight.Array(I);
     I++;
   }
   return Result;
-}    
+}
 
-math_Vector math_Vector::Opposite() {
+math_Vector math_Vector::Opposite()
+{
   math_Vector Result(LowerIndex, UpperIndex);
 
-  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) {
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
     Result.Array(Index) = - Array(Index);
   }
   return Result;
-}    
+}
 
-math_Vector math_Vector::Multiplied(const math_Matrix& Right)const {    
-  Standard_DimensionError_Raise_if(Length() != Right.RowNumber(), "");
+math_Vector math_Vector::Multiplied(const math_Matrix& theRight)const
+{
+  Standard_DimensionError_Raise_if(Length() != theRight.RowNumber(), "");
 
-  math_Vector Result(Right.LowerColIndex, Right.UpperColIndex);
-  for(Standard_Integer J2 = Right.LowerColIndex; 
-      J2 <= Right.UpperColIndex; J2++) {
-    Array(J2) = 0.0;
-    Standard_Integer I2 = Right.LowerRowIndex;
-    for(Standard_Integer I = LowerIndex; I <= UpperIndex; I++) {
-      Result.Array(J2) = Result.Array(J2) + Array(I) * 
-	Right.Array(I2, J2);
-      I2++;
-    }
+  math_Vector Result(theRight.LowerColIndex, theRight.UpperColIndex);
+  for(Standard_Integer J2 = theRight.LowerColIndex; J2 <= theRight.UpperColIndex; J2++)
+  {
+      Array(J2) = 0.0;
+      Standard_Integer theI2 = theRight.LowerRowIndex;
+      for(Standard_Integer I = LowerIndex; I <= UpperIndex; I++)
+      {
+        Result.Array(J2) = Result.Array(J2) + Array(I) * theRight.Array(theI2, J2);
+        theI2++;
+      }
   }
   return Result;
-}    
+}
 
-
-void math_Vector::Multiply(const Standard_Real Left,
-			   const math_Vector& Right) 
+void math_Vector::Multiply(const Standard_Real theLeft, const math_Vector& theRight)
 {
-  Standard_DimensionError_Raise_if((Length() != Right.Length()),
-				   "");
-  for(Standard_Integer I = LowerIndex; I <= UpperIndex; I++) {
-    Array(I) = Left * Right.Array(I);
+  Standard_DimensionError_Raise_if((Length() != theRight.Length()), "");
+  for(Standard_Integer I = LowerIndex; I <= UpperIndex; I++)
+  {
+    Array(I) = theLeft * theRight.Array(I);
   }
 }
 
+math_Vector& math_Vector::Initialized(const math_Vector& theOther)
+{
+  Standard_DimensionError_Raise_if(Length() != theOther.Length(), "");
 
-math_Vector& math_Vector::Initialized(const math_Vector& Other) {
-
-  Standard_DimensionError_Raise_if(Length() != Other.Length(), "");
-  
-  (Other.Array).Copy(Array);
+  (theOther.Array).Copy(Array);
   return *this;
 }
 
-
-
-void math_Vector::Dump(Standard_OStream& o) const
+void math_Vector::Dump(Standard_OStream& theO) const
 {
-  o << "math_Vector of Length = " << Length() << "\n";
-  for(Standard_Integer Index = LowerIndex; 
-      Index <= UpperIndex; Index++) {
-    o << "math_Vector(" << Index << ") = " << Array(Index) << "\n";
+  theO << "math_Vector of Length = " << Length() << "\n";
+  for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++)
+  {
+    theO << "math_Vector(" << Index << ") = " << Array(Index) << "\n";
   }
 }
+

src/math/math_Vector.hxx

@@ -0,0 +1,326 @@
+// Copyright (c) 1997-1999 Matra Datavision
+// Copyright (c) 1999-2014 OPEN CASCADE SAS
+//
+// This file is part of Open CASCADE Technology software library.
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License version 2.1 as published
+// by the Free Software Foundation, with special exception defined in the file
+// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
+// distribution for complete text of the license and disclaimer of any warranty.
+//
+// Alternatively, this file may be used under the terms of Open CASCADE
+// commercial license or contractual agreement.
+
+#ifndef _math_Vector_HeaderFile
+#define _math_Vector_HeaderFile
+
+#include <math_SingleTab.hxx>
+
+class Standard_DimensionError;
+class Standard_DivideByZero;
+class Standard_RangeError;
+class Standard_NullValue;
+class math_Matrix;
+
+//! This class implements the real vector abstract data type.
+//! Vectors can have an arbitrary range which must be defined at
+//! the declaration and cannot be changed after this declaration.
+//! @code
+//!    math_Vector V1(-3, 5); // a vector with range [-3..5]
+//! @endcode
+//!
+//! Vector are copied through assignement :
+//! @code
+//!    math_Vector V2( 1, 9);
+//!    ....
+//!    V2 = V1;
+//!    V1(1) = 2.0; // the vector V2 will not be modified.
+//! @endcode
+//!
+//! The Exception RangeError is raised when trying to access outside
+//! the range of a vector :
+//! @code
+//!    V1(11) = 0.0 // --> will raise RangeError;
+//! @endcode
+//!
+//! The Exception DimensionError is raised when the dimensions of two
+//! vectors are not compatible :
+//! @code
+//!    math_Vector V3(1, 2);
+//!    V3 = V1;    // --> will raise DimensionError;
+//!    V1.Add(V3)  // --> will raise DimensionError;
+//! @endcode
+class math_Vector
+{
+public:
+
+  DEFINE_STANDARD_ALLOC
+
+  //! Contructs a non-initialized vector in the range [theLower..theUpper]
+  //! "theLower" and "theUpper" are the indexes of the lower and upper bounds of the constructed vector.
+  Standard_EXPORT math_Vector(const Standard_Integer theLower, const Standard_Integer theUpper);
+
+  //! Contructs a vector in the range [theLower..theUpper]
+  //! whose values are all initialized with the value "theInitialValue"
+  Standard_EXPORT math_Vector(const Standard_Integer theLower, const Standard_Integer theUpper, const Standard_Real theInitialValue);
+
+  //! Constructs a vector in the range [theLower..theUpper]
+  //! with the "c array" theTab.
+  Standard_EXPORT math_Vector(const Standard_Address theTab, const Standard_Integer theLower, const Standard_Integer theUpper);
+
+  //! Initialize all the elements of a vector with "theInitialValue".
+  Standard_EXPORT void Init(const Standard_Real theInitialValue);
+
+  //! Constructs a copy for initialization.
+  //! An exception is raised if the lengths of the vectors are different.
+  Standard_EXPORT math_Vector(const math_Vector& theOther);
+
+  //! Returns the length of a vector
+  inline Standard_Integer Length() const
+  {
+    return UpperIndex - LowerIndex +1;
+  }
+
+  //! Returns the value of the theLower index of a vector.
+  inline Standard_Integer Lower() const
+  {
+    return LowerIndex;
+  }
+
+  //! Returns the value of the theUpper index of a vector.
+  inline Standard_Integer Upper() const
+  {
+    return UpperIndex;
+  }
+
+  //! Returns the value or the square  of the norm of this vector.
+  Standard_EXPORT Standard_Real Norm() const;
+
+  //! Returns the value of the square of the norm of a vector.
+  Standard_EXPORT Standard_Real Norm2() const;
+
+  //! Returns the value of the "Index" of the maximum element of a vector.
+  Standard_EXPORT Standard_Integer Max() const;
+
+  //! Returns the value of the "Index" of the minimum element  of a vector.
+  Standard_EXPORT Standard_Integer Min() const;
+
+  //! Normalizes this vector (the norm of the result
+  //! is equal to 1.0) and assigns the result to this vector
+  //! Exceptions
+  //! Standard_NullValue if this vector is null (i.e. if its norm is
+  //! less than or equal to Standard_Real::RealEpsilon().
+  Standard_EXPORT void Normalize();
+
+  //! Normalizes this vector (the norm of the result
+  //! is equal to 1.0) and creates a new vector
+  //! Exceptions
+  //! Standard_NullValue if this vector is null (i.e. if its norm is
+  //! less than or equal to Standard_Real::RealEpsilon().
+  Standard_EXPORT math_Vector Normalized() const;
+
+  //! Inverts this vector and assigns the result to this vector.
+  Standard_EXPORT void Invert();
+
+  //! Inverts this vector and creates a new vector.
+  Standard_EXPORT math_Vector Inverse() const;
+
+  //! sets a vector from "theI1" to "theI2" to the vector "theV";
+  //! An exception is raised if "theI1" is less than "LowerIndex" or "theI2" is greater than "UpperIndex" or "theI1" is greater than "theI2".
+  //! An exception is raised if "theI2-theI1+1" is different from the "Length" of "theV".
+  Standard_EXPORT void Set(const Standard_Integer theI1, const Standard_Integer theI2, const math_Vector& theV);
+
+  //!Creates a new vector by inverting the values of this vector
+  //! between indexes "theI1" and "theI2".
+  //! If the values of this vector were (1., 2., 3., 4.,5., 6.),
+  //! by slicing it between indexes 2 and 5 the values
+  //! of the resulting vector are (1., 5., 4., 3., 2., 6.)
+  Standard_EXPORT math_Vector Slice(const Standard_Integer theI1, const Standard_Integer theI2) const;
+
+  //! returns the product of a vector and a real value.
+  Standard_EXPORT void Multiply(const Standard_Real theRight);
+
+  void operator *=(const Standard_Real theRight)
+  {
+    Multiply(theRight);
+  }
+
+  //! returns the product of a vector and a real value.
+  Standard_EXPORT math_Vector Multiplied(const Standard_Real theRight) const;
+
+  math_Vector operator*(const Standard_Real theRight) const
+  {
+    return Multiplied(theRight);
+  }
+
+  //! returns the product of a vector and a real value.
+  Standard_EXPORT math_Vector TMultiplied(const Standard_Real theRight) const;
+
+  friend inline math_Vector operator* (const Standard_Real theLeft, const math_Vector& theRight) 
+  {
+    return theRight.Multiplied(theLeft);
+  }
+
+  //! divides a vector by the value "theRight".
+  //! An exception is raised if "theRight" = 0.
+  Standard_EXPORT void Divide(const Standard_Real theRight);
+
+  void operator /=(const Standard_Real theRight) 
+  {
+    Divide(theRight);
+  }
+
+  //! divides a vector by the value "theRight".
+  //! An exception is raised if "theRight" = 0.
+  Standard_EXPORT math_Vector Divided(const Standard_Real theRight) const;
+
+  math_Vector operator/(const Standard_Real theRight) const
+  {
+    return Divided(theRight);
+  }
+
+  //! adds the vector "theRight" to a vector.
+  //! An exception is raised if the vectors have not the same length.
+  //! Warning
+  //! In order to avoid time-consuming copying of vectors, it
+  //! is preferable to use operator += or the function Add whenever possible.
+  Standard_EXPORT void Add(const math_Vector& theRight);
+
+  void operator +=(const math_Vector& theRight) 
+  {
+    Add(theRight);
+  }
+
+  //! adds the vector theRight to a vector.
+  //! An exception is raised if the vectors have not the same length.
+  //! An exception is raised if the lengths are not equal.
+  Standard_EXPORT math_Vector Added(const math_Vector& theRight) const;
+
+  math_Vector operator+(const math_Vector& theRight) const
+  {
+    return Added(theRight);
+  }
+
+  //! sets a vector to the product of the vector "theLeft"
+  //! with the matrix "theRight".
+  Standard_EXPORT void Multiply(const math_Vector& theLeft, const math_Matrix& theRight);
+
+  //!sets a vector to the product of the matrix "theLeft"
+  //! with the vector "theRight".
+  Standard_EXPORT void Multiply(const math_Matrix& theLeft, const math_Vector& theRight);
+
+  //! sets a vector to the product of the transpose
+  //! of the matrix "theTLeft" by the vector "theRight".
+  Standard_EXPORT void TMultiply(const math_Matrix& theTLeft, const math_Vector& theRight);
+
+  //! sets a vector to the product of the vector
+  //! "theLeft" by the transpose of the matrix "theTRight".
+  Standard_EXPORT void TMultiply(const math_Vector& theLeft, const math_Matrix& theTRight);
+
+  //! sets a vector to the sum of the vector "theLeft"
+  //! and the vector "theRight".
+  //! An exception is raised if the lengths are different.
+  Standard_EXPORT void Add(const math_Vector& theLeft, const math_Vector& theRight);
+
+  //! sets a vector to the Subtraction of the
+  //! vector theRight from the vector theLeft.
+  //! An exception is raised if the vectors have not the same length.
+  //! Warning
+  //! In order to avoid time-consuming copying of vectors, it
+  //! is preferable to use operator -= or the function
+  //! Subtract whenever possible.
+  Standard_EXPORT void Subtract(const math_Vector& theLeft,const math_Vector& theRight);
+
+  //! accesses (in read or write mode) the value of index "theNum" of a vector.
+  inline Standard_Real& Value(const Standard_Integer theNum) const
+  {
+    Standard_RangeError_Raise_if(theNum < LowerIndex || theNum > UpperIndex, " ");
+    return Array(theNum);
+  }
+
+  Standard_Real& operator()(const Standard_Integer theNum) const
+  {
+    return Value(theNum);
+  }
+
+  //! Initialises a vector by copying "theOther".
+  //! An exception is raised if the Lengths are differents.
+  Standard_EXPORT math_Vector& Initialized(const math_Vector& theOther);
+
+  math_Vector& operator=(const math_Vector& theOther)
+  {
+    return Initialized(theOther);
+  }
+
+  //! returns the inner product of 2 vectors.
+  //! An exception is raised if the lengths are not equal.
+  Standard_EXPORT Standard_Real Multiplied(const math_Vector& theRight) const;
+  Standard_Real operator*(const math_Vector& theRight) const
+  {
+    return Multiplied(theRight);
+  }
+
+  //! returns the product of a vector by a matrix.
+  Standard_EXPORT math_Vector Multiplied(const math_Matrix& theRight) const;
+
+  math_Vector operator*(const math_Matrix& theRight) const
+  {
+    return Multiplied(theRight);
+  }
+
+  //! returns the opposite of a vector.
+  Standard_EXPORT math_Vector Opposite();
+
+  math_Vector operator-()
+  {
+    return Opposite();
+  }
+
+  //! returns the subtraction of "theRight" from "me".
+  //! An exception is raised if the vectors have not the same length.
+  Standard_EXPORT void Subtract(const math_Vector& theRight);
+
+  void operator-=(const math_Vector& theRight)
+  {
+    Subtract(theRight);
+  }
+
+  //! returns the subtraction of "theRight" from "me".
+  //! An exception is raised if the vectors have not the same length.
+  Standard_EXPORT math_Vector Subtracted(const math_Vector& theRight) const;
+
+  math_Vector operator-(const math_Vector& theRight) const
+  {
+    return Subtracted(theRight);
+  }
+
+  //! returns the multiplication of a real by a vector.
+  //! "me" = "theLeft" * "theRight"
+  Standard_EXPORT void Multiply(const Standard_Real theLeft,const math_Vector& theRight);
+
+  //! Prints information on the current state of the object.
+  //! Is used to redefine the operator <<.
+  Standard_EXPORT void Dump(Standard_OStream& theO) const;
+
+  friend inline Standard_OStream& operator<<(Standard_OStream& theO, const math_Vector& theVec)
+  {
+    theVec.Dump(theO);
+    return theO;
+  }
+
+  friend class math_Matrix;
+
+protected:
+
+  //! Is used internally to set the "theLower" value of the vector.
+  void SetLower(const Standard_Integer theLower);
+
+private:
+
+  Standard_Integer LowerIndex;
+  Standard_Integer UpperIndex;
+  math_SingleTab<Standard_Real> Array;
+};
+
+#endif