From ac1a15a827a00652b284de9fb3239478e0cdaae3 Mon Sep 17 00:00:00 2001 From: emv Date: Thu, 23 Oct 2014 10:03:01 +0400 Subject: [PATCH] 0025410: Tool for extended check of validity of the curve on the surface Patch for 6.7.1 version of OCCT. --- src/BOPAlgo/BOPAlgo.cdl | 1 + src/BOPAlgo/BOPAlgo_ArgumentAnalyzer.cdl | 42 +- src/BOPAlgo/BOPAlgo_ArgumentAnalyzer.cxx | 57 ++ src/BOPAlgo/BOPAlgo_ArgumentAnalyzer.lxx | 4 + src/BOPAlgo/BOPAlgo_CheckResult.cdl | 42 +- src/BOPAlgo/BOPAlgo_CheckResult.cxx | 48 +- src/BOPTest/BOPTest_CheckCommands.cxx | 266 ++++++++- src/BOPTools/BOPTools_AlgoTools.cdl | 25 + src/BOPTools/BOPTools_AlgoTools_1.cxx | 288 ++++++++++ src/math/FILES | 7 + src/math/math.cdl | 13 +- src/math/math_DoubleTab.cdl | 10 +- src/math/math_DoubleTab.cxx | 142 +++++ src/math/math_DoubleTab.lxx | 6 +- src/math/math_GlobOptMin.cxx | 488 ++++++++++++++++ src/math/math_GlobOptMin.hxx | 123 ++++ src/math/math_IntegerVector.cxx | 428 +++++++------- src/math/math_IntegerVector.hxx | 261 +++++++++ src/math/math_Matrix.cdl | 4 +- src/math/math_Matrix.cxx | 4 +- src/math/math_MultipleVarFunction.cdl | 2 + src/math/math_MultipleVarFunction.cxx | 3 + .../math_MultipleVarFunctionWithGradient.cdl | 2 +- .../math_MultipleVarFunctionWithHessian.cdl | 3 + .../math_MultipleVarFunctionWithHessian.cxx | 3 + src/math/math_NewtonMinimum.cxx | 19 +- src/math/math_SingleTab.hxx | 115 ++++ src/math/math_Vector.cxx | 529 +++++++++--------- src/math/math_Vector.hxx | 326 +++++++++++ 29 files changed, 2715 insertions(+), 546 deletions(-) create mode 100644 src/math/math_DoubleTab.cxx create mode 100644 src/math/math_GlobOptMin.cxx create mode 100644 src/math/math_GlobOptMin.hxx create mode 100644 src/math/math_IntegerVector.hxx create mode 100644 src/math/math_SingleTab.hxx create mode 100644 src/math/math_Vector.hxx diff --git a/src/BOPAlgo/BOPAlgo.cdl b/src/BOPAlgo/BOPAlgo.cdl index 4b3be38a7b..dadb51980f 100644 --- a/src/BOPAlgo/BOPAlgo.cdl +++ b/src/BOPAlgo/BOPAlgo.cdl @@ -52,6 +52,7 @@ is IncompatibilityOfFace, OperationAborted, GeomAbs_C0, + InvalidCurveOnSurface, NotValid end CheckStatus; diff --git a/src/BOPAlgo/BOPAlgo_ArgumentAnalyzer.cdl b/src/BOPAlgo/BOPAlgo_ArgumentAnalyzer.cdl index 47fa9f4943..0c15f28392 100644 --- a/src/BOPAlgo/BOPAlgo_ArgumentAnalyzer.cdl +++ b/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; diff --git a/src/BOPAlgo/BOPAlgo_ArgumentAnalyzer.cxx b/src/BOPAlgo/BOPAlgo_ArgumentAnalyzer.cxx index 3326863a2d..c782f2ccd9 100644 --- a/src/BOPAlgo/BOPAlgo_ArgumentAnalyzer.cxx +++ b/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: diff --git a/src/BOPAlgo/BOPAlgo_ArgumentAnalyzer.lxx b/src/BOPAlgo/BOPAlgo_ArgumentAnalyzer.lxx index e248a057cc..d5102f0302 100644 --- a/src/BOPAlgo/BOPAlgo_ArgumentAnalyzer.lxx +++ b/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; diff --git a/src/BOPAlgo/BOPAlgo_CheckResult.cdl b/src/BOPAlgo/BOPAlgo_CheckResult.cdl index c733ec953c..541a91bffa 100644 --- a/src/BOPAlgo/BOPAlgo_CheckResult.cdl +++ b/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; diff --git a/src/BOPAlgo/BOPAlgo_CheckResult.cxx b/src/BOPAlgo/BOPAlgo_CheckResult.cxx index 2e0fd35013..9a87b82ca0 100644 --- a/src/BOPAlgo/BOPAlgo_CheckResult.cxx +++ b/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; +} diff --git a/src/BOPTest/BOPTest_CheckCommands.cxx b/src/BOPTest/BOPTest_CheckCommands.cxx index c93f5f8c0a..8b30917b50 100644 --- a/src/BOPTest/BOPTest_CheckCommands.cxx +++ b/src/BOPTest/BOPTest_CheckCommands.cxx @@ -15,6 +15,11 @@ #include #include +#include + +#include +#include + #include #include @@ -50,11 +55,20 @@ #include #include +#include + 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**); + // //======================================================================= @@ -77,6 +91,12 @@ 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); } //======================================================================= @@ -293,7 +313,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 +332,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 +355,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 << " -" << "\n"; di << " F (fuse)" << "\n"; di << " O (common)" << "\n"; @@ -344,6 +375,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 +459,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 +532,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 +585,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 +594,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 +706,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 +740,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 +821,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 +913,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 +934,167 @@ 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, aMaxPar, aTolE; + // + 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; + } + // + aTolE = BRep_Tool::Tolerance(anEdge); + 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, aD, aDMax; + TopExp_Explorer aExpF, aExpE; + char buf[200], aFName[10], anEName[10]; + NCollection_DataMap 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)) { + aD = aDMETol.Find(aE); + if (aDMax > aD) { + aDMETol.UnBind(aE); + aDMETol.Bind(aE, 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; +} diff --git a/src/BOPTools/BOPTools_AlgoTools.cdl b/src/BOPTools/BOPTools_AlgoTools.cdl index 541cd52f9a..8f2868cbdd 100644 --- a/src/BOPTools/BOPTools_AlgoTools.cdl +++ b/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,26 @@ is returns Boolean from Standard; ---Purpose: Returns true if the solid is inverted + ComputeTolerance(myclass; + theCurve3D : Curve from Geom; + theCurve2D : Curve from Geom2d; + theSurf : Surface from Geom; + theFirst : Real from Standard; + theLast : Real from Standard; + 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 + end AlgoTools; diff --git a/src/BOPTools/BOPTools_AlgoTools_1.cxx b/src/BOPTools/BOPTools_AlgoTools_1.cxx index 97d2efb51d..3efbc1bbdd 100644 --- a/src/BOPTools/BOPTools_AlgoTools_1.cxx +++ b/src/BOPTools/BOPTools_AlgoTools_1.cxx @@ -111,6 +111,118 @@ static static void UpdateVertices(const TopoDS_Edge& aE); + +//======================================================================= +//class : BOPTools_CheckCurveOnSurface +//purpose : it is used to check the curve on the surface +//======================================================================= +#include +#include +#include +#include + +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) { + 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) { + 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; +}; +//======================================================================= + //======================================================================= // Function : CorrectTolerances // purpose : @@ -812,3 +924,179 @@ void CheckEdge (const TopoDS_Edge& Ed, const Standard_Real aMaxTol) } } } + +//======================================================================= +// Function : ComputeTolerance +// purpose : +//======================================================================= +Standard_Boolean BOPTools_AlgoTools::ComputeTolerance + (const Handle(Geom_Curve)& theCurve3D, + const Handle(Geom2d_Curve)& theCurve2D, + const Handle(Geom_Surface)& theSurf, + const Standard_Real theFirst, + const Standard_Real theLast, + Standard_Real& theMaxDist, + Standard_Real& theMaxPar) +{ + if (theCurve3D.IsNull() || + theCurve2D.IsNull() || + theSurf.IsNull()) { + return Standard_False; + } + // + try { + const Standard_Integer aNbParticles = 100; + // + BOPTools_CheckCurveOnSurface aFunc(theCurve3D, theCurve2D, theSurf); + // + math_Vector anOutputParam(1, 1); + anOutputParam(1) = theFirst; + // + theMaxDist = 0.; + theMaxPar = anOutputParam(1); + aFunc.Value(anOutputParam, theMaxDist); + // + Standard_Real aValue = 0.; + math_Vector aFirstV(1, 1), aLastV(1, 1), aStepV(1, 1); + aFirstV(1) = theFirst; + aLastV(1) = theLast; + aStepV = (aLastV - aFirstV)/(100.0*aNbParticles); + // + math_GlobOptMin aFinder(&aFunc, aFirstV, aLastV); + aFinder.SetTol(1.0e-2, 1.0e-3); + aFinder.Perform(); + // + Standard_Integer i, aNbExtr = aFinder.NbExtrema(); + for(i = 1; i <= aNbExtr; i++) { + aFinder.Points(i, anOutputParam); + aFunc.Value(anOutputParam, aValue); + aValue = sqrt(Abs(aValue)); + // + if(aValue > theMaxDist) { + theMaxDist = aValue; + theMaxPar = anOutputParam(1); + } + } + } + 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, aFirst, aLast, 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, aFirst, aLast, 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); + // + if(BOPTools_AlgoTools::ComputeTolerance + (aC, aC2d, aPlane, aFirst, aLast, aD, aT)) { + bRet = Standard_True; + if (aD > theMaxDist) { + theMaxDist = aD; + theParameter = aT; + } + } + // + return bRet; +} \ No newline at end of file diff --git a/src/math/FILES b/src/math/FILES index 7988203b0e..dc3df00e7e 100755 --- a/src/math/FILES +++ b/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 diff --git a/src/math/math.cdl b/src/math/math.cdl index 20eb358cfb..ea190ffe9d 100644 --- a/src/math/math.cdl +++ b/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 diff --git a/src/math/math_DoubleTab.cdl b/src/math/math_DoubleTab.cdl index a71a003c15..f534559ad9 100644 --- a/src/math/math_DoubleTab.cdl +++ b/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; diff --git a/src/math/math_DoubleTab.cxx b/src/math/math_DoubleTab.cxx new file mode 100644 index 0000000000..670f4b14bc --- /dev/null +++ b/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 + +#include +#include +#include +#include + +// 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; +} + diff --git a/src/math/math_DoubleTab.lxx b/src/math/math_DoubleTab.lxx index 89782feea4..1a467bc568 100644 --- a/src/math/math_DoubleTab.lxx +++ b/src/math/math_DoubleTab.lxx @@ -17,10 +17,10 @@ #include #include -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))); } diff --git a/src/math/math_GlobOptMin.cxx b/src/math/math_GlobOptMin.cxx new file mode 100644 index 0000000000..2ad9b165e6 --- /dev/null +++ b/src/math/math_GlobOptMin.cxx @@ -0,0 +1,488 @@ +// 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 + +#include +#include +#include +#include +#include +#include +#include +#include + +const Handle(Standard_Type)& STANDARD_TYPE(math_GlobOptMin) +{ + static Handle(Standard_Type) _atype = new Standard_Type ("math_GlobOptMin", sizeof (math_GlobOptMin)); + return _atype; +} + +//======================================================================= +//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) +{ + 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; + } + + 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); + } + + 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; + } + + 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 initial values for myF, myY, myC. + computeInitialValues(); + + // 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; + } + + 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(myFunc)) + { + math_MultipleVarFunctionWithHessian* myTmp = + dynamic_cast (myFunc); + + math_NewtonMinimum newtonMinimum(*myTmp, thePnt); + if (newtonMinimum.IsDone()) + { + newtonMinimum.Location(theOutPnt); + theVal = newtonMinimum.Minimum(); + } + else return Standard_False; + } else + + // BFGS method used. + if (dynamic_cast(myFunc)) + { + math_MultipleVarFunctionWithGradient* myTmp = + dynamic_cast (myFunc); + math_BFGS bfgs(*myTmp, thePnt); + if (bfgs.IsDone()) + { + bfgs.Location(theOutPnt); + theVal = bfgs.Minimum(); + } + else return Standard_False; + } else + + // Powell method used. + if (dynamic_cast(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, thePnt, m, 1e-10); + + 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); + + 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, aPrevVal; + Standard_Integer aPntNb = 13; + myFunc->Value(myA, aPrevVal); + Standard_Real aStep = (myB - myA).Norm() / aPntNb; + for(i = 1; i <= aPntNb; i++) + { + aCurrPnt = myA + (myB - myA) * i / (aPntNb - 1); + myFunc->Value(aCurrPnt, aCurrVal); + + if(Abs(aCurrVal - aPrevVal) / aStep > aLipConst) + aLipConst = Abs(aCurrVal - aPrevVal) / aStep; + + aPrevVal = aCurrVal; + } + aLipConst *= Sqrt(myN); + + 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 stepBestValue = RealLast(); + Standard_Real realStep = RealLast(); + math_Vector stepBestPoint(1, myN); + Standard_Boolean isInside = Standard_False; + Standard_Real r; + + + for(myX(j) = myA(j) + myE1; myX(j) < myB(j) + myE1; myX(j) += myV(j)) + { + if (myX(j) > myB(j)) + myX(j) = myB(j); + + if (j == 1) + { + isInside = Standard_False; + myFunc->Value(myX, d); + r = (d - myF) * myZ; + if(r > myE3) + { + isInside = computeLocalExtremum(myX, val, myTmp); + } + stepBestValue = (isInside && (val < d))? val : d; + stepBestPoint = (isInside && (val < d))? myTmp : myX; + + // Solutions are close to each other. + if (Abs(stepBestValue - myF) < mySameTol * 0.01) + { + if (!isStored(stepBestPoint)) + { + if ((stepBestValue - myF) * myZ > 0.0) + myF = stepBestValue; + for(i = 1; i <= myN; i++) + myY.Append(stepBestPoint(i)); + mySolCount++; + } + } + + // New best solution. + if ((stepBestValue - myF) * myZ > mySameTol * 0.01) + { + mySolCount = 0; + myF = stepBestValue; + myY.Clear(); + for(i = 1; i <= myN; i++) + myY.Append(stepBestPoint(i)); + mySolCount++; + } + + realStep = myE2 + Abs(myF - d) / myC; + myV(1) = Min(realStep, myMaxV(1)); + } + else + { + myV(j) = RealLast() / 2.0; + computeGlobalExtremum(j - 1); + } + if ((j < myN) && (myV(j + 1) > realStep)) + { + if (realStep > myMaxV(j + 1)) // Case of too big step. + myV(j + 1) = myMaxV(j + 1); + else + myV(j + 1) = realStep; + } + } +} + +//======================================================================= +//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); +} diff --git a/src/math/math_GlobOptMin.hxx b/src/math/math_GlobOptMin.hxx new file mode 100644 index 0000000000..ac5a4db9e9 --- /dev/null +++ b/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 +#include +#include + +//! This class represents Evtushenko's algorithm of global optimization based on nonuniform mesh.
+//! Article: Yu. Evtushenko. Numerical methods for finding global extreme (case of a non-uniform mesh).
+//! 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 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. + + Standard_Real myF; // Current value of Global optimum. +}; + +const Handle(Standard_Type)& TYPE(math_GlobOptMin); + +#endif diff --git a/src/math/math_IntegerVector.cxx b/src/math/math_IntegerVector.cxx index 7ef36981ae..ae3b192965 100644 --- a/src/math/math_IntegerVector.cxx +++ b/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 +#include #include #include - - -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 theFirst, const Standard_Integer theLast) : + FirstIndex(theFirst), + LastIndex(theLast), + Array(theFirst, theLast) +{ + Standard_RangeError_Raise_if(theFirst > theLast, " "); } -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_Integer theFirst, + const Standard_Integer theLast, + const Standard_Integer theInitialValue) : + FirstIndex(theFirst), + LastIndex(theLast), + Array(theFirst, theLast) +{ + Standard_RangeError_Raise_if(theFirst > theLast, " "); + Array.Init(theInitialValue); } - -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_Address theTab, + const Standard_Integer theFirst, + const Standard_Integer theLast) : + FirstIndex(theFirst), + LastIndex(theLast), + Array(theTab, theFirst, theLast) { - Standard_RangeError_Raise_if(First > Last, " "); + Standard_RangeError_Raise_if(theFirst > theLast, " "); } - -void math_IntegerVector::Init(const Standard_Integer InitialValue) +void math_IntegerVector::Init(const Standard_Integer theInitialValue) { - Array.Init(InitialValue); + 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 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 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 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++) { +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), " "); + + if(theI2 >= theI1) + { + math_IntegerVector Result(theI1, theI2); + for(Standard_Integer Index = theI1; Index <= theI2; Index++) + { Result.Array(Index) = Array(Index); - } - return Result; + } + return Result; } - else { - math_IntegerVector Result(I2, I1); - for(Standard_Integer Index = I1; Index >= I2; Index--) { + 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); +} + +void math_IntegerVector::Add (const math_IntegerVector& theLeft, const math_IntegerVector& theRight) +{ + 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& 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); +} + +void math_IntegerVector::Subtract (const math_IntegerVector& theLeft, + const math_IntegerVector& theRight) +{ + 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::Multiply(const Standard_Integer Left, - const math_IntegerVector& Right) +void math_IntegerVector::Multiply(const Standard_Integer 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()), " "); + for(Standard_Integer I = FirstIndex; I <= LastIndex; I++) + { + Array(I) = theLeft * theRight.Array(I); } } +math_IntegerVector& math_IntegerVector::Initialized(const math_IntegerVector& theOther) +{ + Standard_DimensionError_Raise_if(Length() != theOther.Length(), " "); -math_IntegerVector& math_IntegerVector::Initialized (const math_IntegerVector& Other) { - - Standard_DimensionError_Raise_if(Length() != Other.Length(), " "); - - (Other.Array).Copy(Array); + (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"; + } } - - - - - diff --git a/src/math/math_IntegerVector.hxx b/src/math/math_IntegerVector.hxx new file mode 100644 index 0000000000..2e9563c6a2 --- /dev/null +++ b/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 + +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 Array; +}; + +#endif + diff --git a/src/math/math_Matrix.cdl b/src/math/math_Matrix.cdl index e6ffd3e4c9..95949206e3 100644 --- a/src/math/math_Matrix.cdl +++ b/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 diff --git a/src/math/math_Matrix.cxx b/src/math/math_Matrix.cxx index 36228ec32b..7383359baa 100644 --- a/src/math/math_Matrix.cxx +++ b/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) || diff --git a/src/math/math_MultipleVarFunction.cdl b/src/math/math_MultipleVarFunction.cdl index d77983fdaf..cb0d481d38 100644 --- a/src/math/math_MultipleVarFunction.cdl +++ b/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: diff --git a/src/math/math_MultipleVarFunction.cxx b/src/math/math_MultipleVarFunction.cxx index fa5c832313..b1cb1877c3 100644 --- a/src/math/math_MultipleVarFunction.cxx +++ b/src/math/math_MultipleVarFunction.cxx @@ -16,3 +16,6 @@ #include Standard_Integer math_MultipleVarFunction::GetStateNumber() { return 0; } + +void math_MultipleVarFunction::Delete() +{} diff --git a/src/math/math_MultipleVarFunctionWithGradient.cdl b/src/math/math_MultipleVarFunctionWithGradient.cdl index c2bbc398c1..8804416238 100644 --- a/src/math/math_MultipleVarFunctionWithGradient.cdl +++ b/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) diff --git a/src/math/math_MultipleVarFunctionWithHessian.cdl b/src/math/math_MultipleVarFunctionWithHessian.cdl index 4f387e5ed0..b20035f2ec 100644 --- a/src/math/math_MultipleVarFunctionWithHessian.cdl +++ b/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. diff --git a/src/math/math_MultipleVarFunctionWithHessian.cxx b/src/math/math_MultipleVarFunctionWithHessian.cxx index 972f474112..2c0044e33a 100644 --- a/src/math/math_MultipleVarFunctionWithHessian.cxx +++ b/src/math/math_MultipleVarFunctionWithHessian.cxx @@ -15,3 +15,6 @@ // commercial license or contractual agreement. #include + +void math_MultipleVarFunctionWithHessian::Delete() +{} diff --git a/src/math/math_NewtonMinimum.cxx b/src/math/math_NewtonMinimum.cxx index f337945220..7788c6915f 100644 --- a/src/math/math_NewtonMinimum.cxx +++ b/src/math/math_NewtonMinimum.cxx @@ -143,12 +143,19 @@ void math_NewtonMinimum::Perform(math_MultipleVarFunctionWithHessian& F, } LU.Solve(TheGradient, TheStep); - *suivant = *precedent - TheStep; - - - // Gestion de la convergence - - F.Value(*suivant, TheMinimum); + Standard_Boolean hasProblem = Standard_False; + do + { + *suivant = *precedent - TheStep; + + // Gestion de la convergence + hasProblem = !(F.Value(*suivant, TheMinimum)); + + if (hasProblem) + { + TheStep /= 2.0; + } + } while (hasProblem); if (IsConverged()) { NbConv++; } else { NbConv=0; } diff --git a/src/math/math_SingleTab.hxx b/src/math/math_SingleTab.hxx new file mode 100644 index 0000000000..9b6ffd5a82 --- /dev/null +++ b/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 +#include +#include + +static const Standard_Integer aLengthOfBuf = 512; + +template 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 diff --git a/src/math/math_Vector.cxx b/src/math/math_Vector.cxx index acf11c8f69..799b8e5765 100644 --- a/src/math/math_Vector.cxx +++ b/src/math/math_Vector.cxx @@ -14,7 +14,7 @@ #include -#include +#include #include #include @@ -22,84 +22,85 @@ #include #include - -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) { +void math_Vector::Init(const Standard_Real theInitialValue) +{ + Array.Init(theInitialValue); +} - Array.SetLower(Lower); - UpperIndex = UpperIndex - LowerIndex + Lower; - LowerIndex = Lower; +math_Vector::math_Vector(const math_Vector& theOther) : + LowerIndex(theOther.LowerIndex), + UpperIndex(theOther.UpperIndex), + Array(theOther.Array) +{ } -Standard_Real math_Vector::Norm() const { +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& theRight) +{ + Standard_DimensionError_Raise_if(Length() != theRight.Length(), ""); -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); + 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::Added(const math_Vector& Right) const{ +} - Standard_DimensionError_Raise_if(Length() != Right.Length(), ""); +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 = Right.LowerIndex; - for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) { - Result.Array(Index) = Array(Index) + Right.Array(I); + + 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& Right) { +} - Standard_DimensionError_Raise_if(Length() != Right.Length(), ""); +void math_Vector::Subtract(const math_Vector& theRight) +{ + Standard_DimensionError_Raise_if(Length() != theRight.Length(), ""); - Standard_Integer I = Right.LowerIndex; - for(Standard_Integer Index = LowerIndex; Index <= UpperIndex; Index++) { - Array(Index) = Array(Index) - Right.Array(I); + 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& Right) const { +math_Vector math_Vector::Subtracted (const math_Vector& theRight) const +{ + Standard_DimensionError_Raise_if(Length() != theRight.Length(), ""); - 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); + + 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 I1, - const Standard_Integer I2) const +math_Vector math_Vector::Slice(const Standard_Integer theI1, const Standard_Integer theI2) const { - - Standard_RangeError_Raise_if((I1 < LowerIndex) || - (I1 > UpperIndex) || - (I2 < LowerIndex) || - (I2 > UpperIndex) , ""); - + Standard_RangeError_Raise_if((theI1 < LowerIndex) || (theI1 > UpperIndex) || (theI2 < LowerIndex) || (theI2 > UpperIndex) , ""); - if(I2 >= I1) { - math_Vector Result(I1, I2); - for(Standard_Integer Index = I1; Index <= I2; Index++) { + if(theI2 >= theI1) + { + math_Vector Result(theI1, theI2); + for(Standard_Integer Index = theI1; Index <= theI2; Index++) + { Result.Array(Index) = Array(Index); - } - return Result; + } + return Result; } - else { - math_Vector Result(I2, I1); - for(Standard_Integer Index = I1; Index >= I2; Index--) { + 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) { - - 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); +} + +void math_Vector::Subtract (const math_Vector& theLeft, const math_Vector& theRight) +{ + Standard_DimensionError_Raise_if((Length() != theRight.Length()) || (theRight.Length() != theLeft.Length()), ""); + + 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 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 math_Vector::Multiplied(const math_Matrix& theRight)const +{ + Standard_DimensionError_Raise_if(Length() != theRight.RowNumber(), ""); + + 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"; } } + diff --git a/src/math/math_Vector.hxx b/src/math/math_Vector.hxx new file mode 100644 index 0000000000..e3b1b54fe5 --- /dev/null +++ b/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 + +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 Array; +}; + +#endif -- 2.39.5