| 1 | // Created on: 1998-05-12 |
| 2 | // Created by: Roman BORISOV |
| 3 | // Copyright (c) 1998-1999 Matra Datavision |
| 4 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
| 5 | // |
| 6 | // This file is part of Open CASCADE Technology software library. |
| 7 | // |
| 8 | // This library is free software; you can redistribute it and / or modify it |
| 9 | // under the terms of the GNU Lesser General Public version 2.1 as published |
| 10 | // by the Free Software Foundation, with special exception defined in the file |
| 11 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT |
| 12 | // distribution for complete text of the license and disclaimer of any warranty. |
| 13 | // |
| 14 | // Alternatively, this file may be used under the terms of Open CASCADE |
| 15 | // commercial license or contractual agreement. |
| 16 | |
| 17 | #include <Approx_CurvlinFunc.ixx> |
| 18 | #include <Adaptor3d_CurveOnSurface.hxx> |
| 19 | #include <Adaptor3d_HCurveOnSurface.hxx> |
| 20 | #include <TColStd_SequenceOfReal.hxx> |
| 21 | #include <GeomLib.hxx> |
| 22 | #include <GCPnts_AbscissaPoint.hxx> |
| 23 | #include <Precision.hxx> |
| 24 | |
| 25 | #ifdef __OCC_DEBUG_CHRONO |
| 26 | #include <OSD_Timer.hxx> |
| 27 | static OSD_Chronometer chr_uparam; |
| 28 | Standard_EXPORT Standard_Integer uparam_count; |
| 29 | Standard_EXPORT Standard_Real t_uparam; |
| 30 | |
| 31 | //Standard_IMPORT extern void InitChron(OSD_Chronometer& ch); |
| 32 | Standard_IMPORT void InitChron(OSD_Chronometer& ch); |
| 33 | //Standard_IMPORT extern void ResultChron( OSD_Chronometer & ch, Standard_Real & time); |
| 34 | Standard_IMPORT void ResultChron( OSD_Chronometer & ch, Standard_Real & time); |
| 35 | #endif |
| 36 | |
| 37 | static Standard_Real cubic(const Standard_Real X, const Standard_Real *Xi, const Standard_Real *Yi) |
| 38 | { |
| 39 | Standard_Real I1, I2, I3, I21, I22, I31, Result; |
| 40 | |
| 41 | I1 = (Yi[0] - Yi[1])/(Xi[0] - Xi[1]); |
| 42 | I2 = (Yi[1] - Yi[2])/(Xi[1] - Xi[2]); |
| 43 | I3 = (Yi[2] - Yi[3])/(Xi[2] - Xi[3]); |
| 44 | |
| 45 | I21 = (I1 - I2)/(Xi[0] - Xi[2]); |
| 46 | I22 = (I2 - I3)/(Xi[1] - Xi[3]); |
| 47 | |
| 48 | I31 = (I21 - I22)/(Xi[0] - Xi[3]); |
| 49 | |
| 50 | Result = Yi[0] + (X - Xi[0])*(I1 + (X - Xi[1])*(I21 + (X - Xi[2])*I31)); |
| 51 | |
| 52 | return Result; |
| 53 | } |
| 54 | |
| 55 | //static void findfourpoints(const Standard_Real S, |
| 56 | static void findfourpoints(const Standard_Real , |
| 57 | Standard_Integer NInterval, |
| 58 | const Handle(TColStd_HArray1OfReal)& Si, |
| 59 | Handle(TColStd_HArray1OfReal)& Ui, |
| 60 | const Standard_Real prevS, |
| 61 | const Standard_Real prevU, Standard_Real *Xi, |
| 62 | Standard_Real *Yi) |
| 63 | { |
| 64 | Standard_Integer i, j; |
| 65 | Standard_Integer NbInt = Si->Length() - 1; |
| 66 | if (NbInt < 3) Standard_ConstructionError::Raise("Approx_CurvlinFunc::GetUParameter"); |
| 67 | |
| 68 | if(NInterval < 1) NInterval = 1; |
| 69 | else if(NInterval > NbInt - 2) NInterval = NbInt - 2; |
| 70 | |
| 71 | for(i = 0; i < 4; i++) { |
| 72 | Xi[i] = Si->Value(NInterval - 1 + i); |
| 73 | Yi[i] = Ui->Value(NInterval - 1 + i); |
| 74 | } |
| 75 | // try to insert (S, U) |
| 76 | for(i = 0; i < 3; i++) { |
| 77 | if(Xi[i] < prevS && prevS < Xi[i+1]) { |
| 78 | for(j = 0; j < i; j++) { |
| 79 | Xi[j] = Xi[j+1]; |
| 80 | Yi[j] = Yi[j+1]; |
| 81 | } |
| 82 | Xi[i] = prevS; |
| 83 | Yi[i] = prevU; |
| 84 | break; |
| 85 | } |
| 86 | } |
| 87 | } |
| 88 | |
| 89 | /*static Standard_Real curvature(const Standard_Real U, const Adaptor3d_Curve& C) |
| 90 | { |
| 91 | Standard_Real k, tau, mod1, mod2, OMEGA; |
| 92 | gp_Pnt P; |
| 93 | gp_Vec D1, D2, D3; |
| 94 | C.D3(U, P, D1, D2, D3); |
| 95 | mod1 = D1.Magnitude(); |
| 96 | mod2 = D1.Crossed(D2).Magnitude(); |
| 97 | k = mod2/(mod1*mod1*mod1); |
| 98 | tau = D1.Dot(D2.Crossed(D3)); |
| 99 | tau /= mod2*mod2; |
| 100 | OMEGA = Sqrt(k*k + tau*tau); |
| 101 | |
| 102 | return OMEGA; |
| 103 | } |
| 104 | */ |
| 105 | |
| 106 | Approx_CurvlinFunc::Approx_CurvlinFunc(const Handle(Adaptor3d_HCurve)& C, const Standard_Real Tol) : myC3D(C), |
| 107 | myCase(1), |
| 108 | myFirstS(0), |
| 109 | myLastS(1), |
| 110 | myTolLen(Tol), |
| 111 | myPrevS (0.0), |
| 112 | myPrevU (0.0) |
| 113 | { |
| 114 | Init(); |
| 115 | } |
| 116 | |
| 117 | Approx_CurvlinFunc::Approx_CurvlinFunc(const Handle(Adaptor2d_HCurve2d)& C2D, const Handle(Adaptor3d_HSurface)& S, const Standard_Real Tol) : |
| 118 | myC2D1(C2D), |
| 119 | mySurf1(S), |
| 120 | myCase(2), |
| 121 | myFirstS(0), |
| 122 | myLastS(1), |
| 123 | myTolLen(Tol), |
| 124 | myPrevS (0.0), |
| 125 | myPrevU (0.0) |
| 126 | { |
| 127 | Init(); |
| 128 | } |
| 129 | |
| 130 | Approx_CurvlinFunc::Approx_CurvlinFunc(const Handle(Adaptor2d_HCurve2d)& C2D1, const Handle(Adaptor2d_HCurve2d)& C2D2, const Handle(Adaptor3d_HSurface)& S1, const Handle(Adaptor3d_HSurface)& S2, const Standard_Real Tol) : |
| 131 | myC2D1(C2D1), |
| 132 | myC2D2(C2D2), |
| 133 | mySurf1(S1), |
| 134 | mySurf2(S2), |
| 135 | myCase(3), |
| 136 | myFirstS(0), |
| 137 | myLastS(1), |
| 138 | myTolLen(Tol), |
| 139 | myPrevS (0.0), |
| 140 | myPrevU (0.0) |
| 141 | { |
| 142 | Init(); |
| 143 | } |
| 144 | |
| 145 | void Approx_CurvlinFunc::Init() |
| 146 | { |
| 147 | Adaptor3d_CurveOnSurface CurOnSur; |
| 148 | |
| 149 | switch(myCase) { |
| 150 | case 1: |
| 151 | Init(myC3D->GetCurve(), mySi_1, myUi_1); |
| 152 | myFirstU1 = myC3D->FirstParameter(); |
| 153 | myLastU1 = myC3D->LastParameter(); |
| 154 | myFirstU2 = myLastU2 = 0; |
| 155 | break; |
| 156 | case 2: |
| 157 | CurOnSur.Load(myC2D1); |
| 158 | CurOnSur.Load(mySurf1); |
| 159 | Init(CurOnSur, mySi_1, myUi_1); |
| 160 | myFirstU1 = CurOnSur.FirstParameter(); |
| 161 | myLastU1 = CurOnSur.LastParameter(); |
| 162 | myFirstU2 = myLastU2 = 0; |
| 163 | break; |
| 164 | case 3: |
| 165 | CurOnSur.Load(myC2D1); |
| 166 | CurOnSur.Load(mySurf1); |
| 167 | Init(CurOnSur, mySi_1, myUi_1); |
| 168 | myFirstU1 = CurOnSur.FirstParameter(); |
| 169 | myLastU1 = CurOnSur.LastParameter(); |
| 170 | CurOnSur.Load(myC2D2); |
| 171 | CurOnSur.Load(mySurf2); |
| 172 | Init(CurOnSur, mySi_2, myUi_2); |
| 173 | myFirstU2 = CurOnSur.FirstParameter(); |
| 174 | myLastU2 = CurOnSur.LastParameter(); |
| 175 | } |
| 176 | |
| 177 | Length(); |
| 178 | } |
| 179 | |
| 180 | |
| 181 | //======================================================================= |
| 182 | //function : Init |
| 183 | //purpose : Init the values |
| 184 | //history : 23/10/1998 PMN : Cut at curve's discontinuities |
| 185 | //======================================================================= |
| 186 | void Approx_CurvlinFunc::Init(Adaptor3d_Curve& C, Handle(TColStd_HArray1OfReal)& Si, |
| 187 | Handle(TColStd_HArray1OfReal)& Ui) const |
| 188 | { |
| 189 | Standard_Real Step, FirstU, LastU; |
| 190 | Standard_Integer i, j, k, NbInt, NbIntC3; |
| 191 | FirstU = C.FirstParameter(); |
| 192 | LastU = C.LastParameter(); |
| 193 | |
| 194 | NbInt = 10; |
| 195 | NbIntC3 = C.NbIntervals(GeomAbs_C3); |
| 196 | TColStd_Array1OfReal Disc(1, NbIntC3+1); |
| 197 | |
| 198 | if (NbIntC3 >1) { |
| 199 | C.Intervals(Disc, GeomAbs_C3); |
| 200 | } |
| 201 | else { |
| 202 | Disc(1) = FirstU; |
| 203 | Disc(2) = LastU; |
| 204 | } |
| 205 | |
| 206 | Ui = new TColStd_HArray1OfReal (0,NbIntC3*NbInt); |
| 207 | Si = new TColStd_HArray1OfReal (0,NbIntC3*NbInt); |
| 208 | |
| 209 | Ui->SetValue(0, FirstU); |
| 210 | Si->SetValue(0, 0); |
| 211 | |
| 212 | for(j = 1, i=1; j<=NbIntC3; j++) { |
| 213 | Step = (Disc(j+1) - Disc(j))/NbInt; |
| 214 | for(k = 1; k <= NbInt; k++, i++) { |
| 215 | Ui->ChangeValue(i) = Ui->Value(i-1) + Step; |
| 216 | Si->ChangeValue(i) = Si->Value(i-1) + Length(C, Ui->Value(i-1), Ui->Value(i)); |
| 217 | } |
| 218 | } |
| 219 | |
| 220 | Standard_Real Len = Si->Value(Si->Upper()); |
| 221 | for(i = Si->Lower(); i<= Si->Upper(); i++) |
| 222 | Si->ChangeValue(i) /= Len; |
| 223 | |
| 224 | // TODO - fields should be mutable |
| 225 | const_cast<Approx_CurvlinFunc*>(this)->myPrevS = myFirstS; |
| 226 | const_cast<Approx_CurvlinFunc*>(this)->myPrevU = FirstU; |
| 227 | } |
| 228 | |
| 229 | void Approx_CurvlinFunc::SetTol(const Standard_Real Tol) |
| 230 | { |
| 231 | myTolLen = Tol; |
| 232 | } |
| 233 | |
| 234 | Standard_Real Approx_CurvlinFunc::FirstParameter() const |
| 235 | { |
| 236 | return myFirstS; |
| 237 | } |
| 238 | |
| 239 | Standard_Real Approx_CurvlinFunc::LastParameter() const |
| 240 | { |
| 241 | return myLastS; |
| 242 | } |
| 243 | |
| 244 | Standard_Integer Approx_CurvlinFunc::NbIntervals(const GeomAbs_Shape S) const |
| 245 | { |
| 246 | Adaptor3d_CurveOnSurface CurOnSur; |
| 247 | |
| 248 | switch(myCase) { |
| 249 | case 1: |
| 250 | return myC3D->NbIntervals(S); |
| 251 | case 2: |
| 252 | CurOnSur.Load(myC2D1); |
| 253 | CurOnSur.Load(mySurf1); |
| 254 | return CurOnSur.NbIntervals(S); |
| 255 | case 3: |
| 256 | Standard_Integer NbInt; |
| 257 | CurOnSur.Load(myC2D1); |
| 258 | CurOnSur.Load(mySurf1); |
| 259 | NbInt = CurOnSur.NbIntervals(S); |
| 260 | TColStd_Array1OfReal T1(1, NbInt+1); |
| 261 | CurOnSur.Intervals(T1, S); |
| 262 | CurOnSur.Load(myC2D2); |
| 263 | CurOnSur.Load(mySurf2); |
| 264 | NbInt = CurOnSur.NbIntervals(S); |
| 265 | TColStd_Array1OfReal T2(1, NbInt+1); |
| 266 | CurOnSur.Intervals(T2, S); |
| 267 | |
| 268 | TColStd_SequenceOfReal Fusion; |
| 269 | GeomLib::FuseIntervals(T1, T2, Fusion); |
| 270 | return Fusion.Length() - 1; |
| 271 | } |
| 272 | |
| 273 | //POP pour WNT |
| 274 | return 1; |
| 275 | } |
| 276 | |
| 277 | void Approx_CurvlinFunc::Intervals(TColStd_Array1OfReal& T, const GeomAbs_Shape S) const |
| 278 | { |
| 279 | Adaptor3d_CurveOnSurface CurOnSur; |
| 280 | Standard_Integer i; |
| 281 | |
| 282 | switch(myCase) { |
| 283 | case 1: |
| 284 | myC3D->Intervals(T, S); |
| 285 | break; |
| 286 | case 2: |
| 287 | CurOnSur.Load(myC2D1); |
| 288 | CurOnSur.Load(mySurf1); |
| 289 | CurOnSur.Intervals(T, S); |
| 290 | break; |
| 291 | case 3: |
| 292 | Standard_Integer NbInt; |
| 293 | CurOnSur.Load(myC2D1); |
| 294 | CurOnSur.Load(mySurf1); |
| 295 | NbInt = CurOnSur.NbIntervals(S); |
| 296 | TColStd_Array1OfReal T1(1, NbInt+1); |
| 297 | CurOnSur.Intervals(T1, S); |
| 298 | CurOnSur.Load(myC2D2); |
| 299 | CurOnSur.Load(mySurf2); |
| 300 | NbInt = CurOnSur.NbIntervals(S); |
| 301 | TColStd_Array1OfReal T2(1, NbInt+1); |
| 302 | CurOnSur.Intervals(T2, S); |
| 303 | |
| 304 | TColStd_SequenceOfReal Fusion; |
| 305 | GeomLib::FuseIntervals(T1, T2, Fusion); |
| 306 | |
| 307 | for (i = 1; i <= Fusion.Length(); i++) |
| 308 | T.ChangeValue(i) = Fusion.Value(i); |
| 309 | } |
| 310 | |
| 311 | for(i = 1; i <= T.Length(); i++) |
| 312 | T.ChangeValue(i) = GetSParameter(T.Value(i)); |
| 313 | } |
| 314 | |
| 315 | void Approx_CurvlinFunc::Trim(const Standard_Real First, const Standard_Real Last, const Standard_Real Tol) |
| 316 | { |
| 317 | if (First < 0 || Last >1) Standard_OutOfRange::Raise("Approx_CurvlinFunc::Trim"); |
| 318 | if ((Last - First) < Tol) return; |
| 319 | |
| 320 | Standard_Real FirstU, LastU; |
| 321 | Adaptor3d_CurveOnSurface CurOnSur; |
| 322 | Handle(Adaptor3d_HCurve) HCurOnSur; |
| 323 | |
| 324 | switch(myCase) { |
| 325 | case 1: |
| 326 | myC3D = myC3D->Trim(myFirstU1, myLastU1, Tol); |
| 327 | FirstU = GetUParameter(myC3D->GetCurve(), First, 1); |
| 328 | LastU = GetUParameter(myC3D->GetCurve(), Last, 1); |
| 329 | myC3D = myC3D->Trim(FirstU, LastU, Tol); |
| 330 | break; |
| 331 | case 3: |
| 332 | CurOnSur.Load(myC2D2); |
| 333 | CurOnSur.Load(mySurf2); |
| 334 | HCurOnSur = CurOnSur.Trim(myFirstU2, myLastU2, Tol); |
| 335 | myC2D2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve(); |
| 336 | mySurf2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface(); |
| 337 | CurOnSur.Load(myC2D2); |
| 338 | CurOnSur.Load(mySurf2); |
| 339 | |
| 340 | FirstU = GetUParameter(CurOnSur, First, 1); |
| 341 | LastU = GetUParameter(CurOnSur, Last, 1); |
| 342 | HCurOnSur = CurOnSur.Trim(FirstU, LastU, Tol); |
| 343 | myC2D2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve(); |
| 344 | mySurf2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface(); |
| 345 | |
| 346 | case 2: |
| 347 | CurOnSur.Load(myC2D1); |
| 348 | CurOnSur.Load(mySurf1); |
| 349 | HCurOnSur = CurOnSur.Trim(myFirstU1, myLastU1, Tol); |
| 350 | myC2D1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve(); |
| 351 | mySurf1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface(); |
| 352 | CurOnSur.Load(myC2D1); |
| 353 | CurOnSur.Load(mySurf1); |
| 354 | |
| 355 | FirstU = GetUParameter(CurOnSur, First, 1); |
| 356 | LastU = GetUParameter(CurOnSur, Last, 1); |
| 357 | HCurOnSur = CurOnSur.Trim(FirstU, LastU, Tol); |
| 358 | myC2D1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve(); |
| 359 | mySurf1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface(); |
| 360 | } |
| 361 | myFirstS = First; |
| 362 | myLastS = Last; |
| 363 | } |
| 364 | |
| 365 | void Approx_CurvlinFunc::Length() |
| 366 | { |
| 367 | Adaptor3d_CurveOnSurface CurOnSur; |
| 368 | Standard_Real FirstU, LastU; |
| 369 | |
| 370 | switch(myCase){ |
| 371 | case 1: |
| 372 | FirstU = myC3D->FirstParameter(); |
| 373 | LastU = myC3D->LastParameter(); |
| 374 | myLength = Length(myC3D->GetCurve(), FirstU, LastU); |
| 375 | myLength1 = myLength2 = 0; |
| 376 | break; |
| 377 | case 2: |
| 378 | CurOnSur.Load(myC2D1); |
| 379 | CurOnSur.Load(mySurf1); |
| 380 | FirstU = CurOnSur.FirstParameter(); |
| 381 | LastU = CurOnSur.LastParameter(); |
| 382 | myLength = Length(CurOnSur, FirstU, LastU); |
| 383 | myLength1 = myLength2 = 0; |
| 384 | break; |
| 385 | case 3: |
| 386 | CurOnSur.Load(myC2D1); |
| 387 | CurOnSur.Load(mySurf1); |
| 388 | FirstU = CurOnSur.FirstParameter(); |
| 389 | LastU = CurOnSur.LastParameter(); |
| 390 | myLength1 = Length(CurOnSur, FirstU, LastU); |
| 391 | CurOnSur.Load(myC2D2); |
| 392 | CurOnSur.Load(mySurf2); |
| 393 | FirstU = CurOnSur.FirstParameter(); |
| 394 | LastU = CurOnSur.LastParameter(); |
| 395 | myLength2 = Length(CurOnSur, FirstU, LastU); |
| 396 | myLength = (myLength1 + myLength2)/2; |
| 397 | } |
| 398 | } |
| 399 | |
| 400 | |
| 401 | Standard_Real Approx_CurvlinFunc::Length(Adaptor3d_Curve& C, const Standard_Real FirstU, const Standard_Real LastU) const |
| 402 | { |
| 403 | Standard_Real Length; |
| 404 | |
| 405 | Length = GCPnts_AbscissaPoint::Length(C, FirstU, LastU, myTolLen); |
| 406 | return Length; |
| 407 | } |
| 408 | |
| 409 | |
| 410 | Standard_Real Approx_CurvlinFunc::GetLength() const |
| 411 | { |
| 412 | return myLength; |
| 413 | } |
| 414 | |
| 415 | Standard_Real Approx_CurvlinFunc::GetSParameter(const Standard_Real U) const |
| 416 | { |
| 417 | Standard_Real S=0, S1, S2; |
| 418 | Adaptor3d_CurveOnSurface CurOnSur; |
| 419 | |
| 420 | switch (myCase) { |
| 421 | case 1: |
| 422 | S = GetSParameter(myC3D->GetCurve(), U, myLength); |
| 423 | break; |
| 424 | case 2: |
| 425 | CurOnSur.Load(myC2D1); |
| 426 | CurOnSur.Load(mySurf1); |
| 427 | S = GetSParameter(CurOnSur, U, myLength); |
| 428 | break; |
| 429 | case 3: |
| 430 | CurOnSur.Load(myC2D1); |
| 431 | CurOnSur.Load(mySurf1); |
| 432 | S1 = GetSParameter(CurOnSur, U, myLength1); |
| 433 | CurOnSur.Load(myC2D2); |
| 434 | CurOnSur.Load(mySurf2); |
| 435 | S2 = GetSParameter(CurOnSur, U, myLength2); |
| 436 | S = (S1 + S2)/2; |
| 437 | } |
| 438 | return S; |
| 439 | } |
| 440 | |
| 441 | |
| 442 | |
| 443 | Standard_Real Approx_CurvlinFunc::GetUParameter(Adaptor3d_Curve& C, |
| 444 | const Standard_Real S, |
| 445 | const Standard_Integer NumberOfCurve) const |
| 446 | { |
| 447 | Standard_Real deltaS, base, U, Length; |
| 448 | Standard_Integer NbInt, NInterval, i; |
| 449 | Handle(TColStd_HArray1OfReal) InitUArray, InitSArray; |
| 450 | #ifdef __OCC_DEBUG_CHRONO |
| 451 | InitChron(chr_uparam); |
| 452 | #endif |
| 453 | if(S < 0 || S > 1) Standard_ConstructionError::Raise("Approx_CurvlinFunc::GetUParameter"); |
| 454 | |
| 455 | if(NumberOfCurve == 1) { |
| 456 | InitUArray = myUi_1; |
| 457 | InitSArray = mySi_1; |
| 458 | if(myCase == 3) |
| 459 | Length = myLength1; |
| 460 | else |
| 461 | Length = myLength; |
| 462 | } |
| 463 | else { |
| 464 | InitUArray = myUi_2; |
| 465 | InitSArray = mySi_2; |
| 466 | Length = myLength2; |
| 467 | } |
| 468 | |
| 469 | NbInt = InitUArray->Length() - 1; |
| 470 | |
| 471 | if(S == 1) NInterval = NbInt - 1; |
| 472 | else { |
| 473 | for(i = 0; i < NbInt; i++) { |
| 474 | if((InitSArray->Value(i) <= S && S < InitSArray->Value(i+1))) |
| 475 | break; |
| 476 | } |
| 477 | NInterval = i; |
| 478 | } |
| 479 | if(S==InitSArray->Value(NInterval)) { |
| 480 | return InitUArray->Value(NInterval); |
| 481 | } |
| 482 | if(S==InitSArray->Value(NInterval+1)) { |
| 483 | return InitUArray->Value(NInterval+1); |
| 484 | } |
| 485 | |
| 486 | base = InitUArray->Value(NInterval); |
| 487 | deltaS = (S - InitSArray->Value(NInterval))*Length; |
| 488 | |
| 489 | // to find an initial point |
| 490 | Standard_Real Xi[4], Yi[4], UGuess; |
| 491 | findfourpoints(S, NInterval, InitSArray, InitUArray, myPrevS, myPrevU, Xi, Yi); |
| 492 | UGuess = cubic(S , Xi, Yi); |
| 493 | |
| 494 | U = GCPnts_AbscissaPoint(C, deltaS, base, UGuess, myTolLen).Parameter(); |
| 495 | |
| 496 | // TODO - fields should be mutable |
| 497 | const_cast<Approx_CurvlinFunc*>(this)->myPrevS = S; |
| 498 | const_cast<Approx_CurvlinFunc*>(this)->myPrevU = U; |
| 499 | |
| 500 | #ifdef __OCC_DEBUG_CHRONO |
| 501 | ResultChron(chr_uparam, t_uparam); |
| 502 | uparam_count++; |
| 503 | #endif |
| 504 | |
| 505 | return U; |
| 506 | } |
| 507 | |
| 508 | Standard_Real Approx_CurvlinFunc::GetSParameter(Adaptor3d_Curve& C, const Standard_Real U, const Standard_Real Len) const |
| 509 | { |
| 510 | Standard_Real S, Origin; |
| 511 | |
| 512 | Origin = C.FirstParameter(); |
| 513 | S = myFirstS + Length(C, Origin, U)/Len; |
| 514 | return S; |
| 515 | } |
| 516 | |
| 517 | Standard_Boolean Approx_CurvlinFunc::EvalCase1(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result) const |
| 518 | { |
| 519 | if(myCase != 1) Standard_ConstructionError::Raise("Approx_CurvlinFunc::EvalCase1"); |
| 520 | |
| 521 | gp_Pnt C; |
| 522 | gp_Vec dC_dU, dC_dS, d2C_dU2, d2C_dS2; |
| 523 | Standard_Real U, Mag, dU_dS, d2U_dS2; |
| 524 | |
| 525 | U = GetUParameter(myC3D->GetCurve(), S, 1); |
| 526 | |
| 527 | switch(Order) { |
| 528 | |
| 529 | case 0: |
| 530 | myC3D->D0(U, C); |
| 531 | |
| 532 | Result(0) = C.X(); |
| 533 | Result(1) = C.Y(); |
| 534 | Result(2) = C.Z(); |
| 535 | break; |
| 536 | |
| 537 | case 1: |
| 538 | myC3D->D1(U, C, dC_dU); |
| 539 | Mag = dC_dU.Magnitude(); |
| 540 | dU_dS = myLength/Mag; |
| 541 | dC_dS = dC_dU*dU_dS; |
| 542 | |
| 543 | Result(0) = dC_dS.X(); |
| 544 | Result(1) = dC_dS.Y(); |
| 545 | Result(2) = dC_dS.Z(); |
| 546 | break; |
| 547 | |
| 548 | case 2: |
| 549 | myC3D->D2(U, C, dC_dU, d2C_dU2); |
| 550 | Mag = dC_dU.Magnitude(); |
| 551 | dU_dS = myLength/Mag; |
| 552 | d2U_dS2 = -myLength*dC_dU.Dot(d2C_dU2)*dU_dS/(Mag*Mag*Mag); |
| 553 | d2C_dS2 = d2C_dU2*dU_dS*dU_dS + dC_dU*d2U_dS2; |
| 554 | |
| 555 | Result(0) = d2C_dS2.X(); |
| 556 | Result(1) = d2C_dS2.Y(); |
| 557 | Result(2) = d2C_dS2.Z(); |
| 558 | break; |
| 559 | |
| 560 | default: Result(0) = Result(1) = Result(2) = 0; |
| 561 | return Standard_False; |
| 562 | } |
| 563 | return Standard_True; |
| 564 | } |
| 565 | |
| 566 | Standard_Boolean Approx_CurvlinFunc::EvalCase2(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result) const |
| 567 | { |
| 568 | if(myCase != 2) Standard_ConstructionError::Raise("Approx_CurvlinFunc::EvalCase2"); |
| 569 | |
| 570 | Standard_Boolean Done; |
| 571 | |
| 572 | Done = EvalCurOnSur(S, Order, Result, 1); |
| 573 | |
| 574 | return Done; |
| 575 | } |
| 576 | |
| 577 | Standard_Boolean Approx_CurvlinFunc::EvalCase3(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result) |
| 578 | { |
| 579 | if(myCase != 3) Standard_ConstructionError::Raise("Approx_CurvlinFunc::EvalCase3"); |
| 580 | |
| 581 | TColStd_Array1OfReal tmpRes1(0, 4), tmpRes2(0, 4); |
| 582 | Standard_Boolean Done; |
| 583 | |
| 584 | Done = EvalCurOnSur(S, Order, tmpRes1, 1); |
| 585 | |
| 586 | Done = EvalCurOnSur(S, Order, tmpRes2, 2) && Done; |
| 587 | |
| 588 | Result(0) = tmpRes1(0); |
| 589 | Result(1) = tmpRes1(1); |
| 590 | Result(2) = tmpRes2(0); |
| 591 | Result(3) = tmpRes2(1); |
| 592 | Result(4) = 0.5*(tmpRes1(2) + tmpRes2(2)); |
| 593 | Result(5) = 0.5*(tmpRes1(3) + tmpRes2(3)); |
| 594 | Result(6) = 0.5*(tmpRes1(4) + tmpRes2(4)); |
| 595 | |
| 596 | return Done; |
| 597 | } |
| 598 | |
| 599 | Standard_Boolean Approx_CurvlinFunc::EvalCurOnSur(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result, const Standard_Integer NumberOfCurve) const |
| 600 | { |
| 601 | Handle(Adaptor2d_HCurve2d) Cur2D; |
| 602 | Handle(Adaptor3d_HSurface) Surf; |
| 603 | Standard_Real U=0, Length=0; |
| 604 | |
| 605 | if (NumberOfCurve == 1) { |
| 606 | Cur2D = myC2D1; |
| 607 | Surf = mySurf1; |
| 608 | Adaptor3d_CurveOnSurface CurOnSur(myC2D1, mySurf1); |
| 609 | U = GetUParameter(CurOnSur, S, 1); |
| 610 | if(myCase == 3) Length = myLength1; |
| 611 | else Length = myLength; |
| 612 | } |
| 613 | else if (NumberOfCurve == 2) { |
| 614 | Cur2D = myC2D2; |
| 615 | Surf = mySurf2; |
| 616 | Adaptor3d_CurveOnSurface CurOnSur(myC2D2, mySurf2); |
| 617 | U = GetUParameter(CurOnSur, S, 2); |
| 618 | Length = myLength2; |
| 619 | } |
| 620 | else |
| 621 | Standard_ConstructionError::Raise("Approx_CurvlinFunc::EvalCurOnSur"); |
| 622 | |
| 623 | Standard_Real Mag, dU_dS, d2U_dS2, dV_dU, dW_dU, dV_dS, dW_dS, d2V_dS2, d2W_dS2, d2V_dU2, d2W_dU2; |
| 624 | gp_Pnt2d C2D; |
| 625 | gp_Pnt C; |
| 626 | gp_Vec2d dC2D_dU, d2C2D_dU2; |
| 627 | gp_Vec dC_dU, d2C_dU2, dC_dS, d2C_dS2, dS_dV, dS_dW, d2S_dV2, d2S_dW2, d2S_dVdW; |
| 628 | |
| 629 | switch(Order) { |
| 630 | case 0: |
| 631 | Cur2D->D0(U, C2D); |
| 632 | Surf->D0(C2D.X(), C2D.Y(), C); |
| 633 | |
| 634 | Result(0) = C2D.X(); |
| 635 | Result(1) = C2D.Y(); |
| 636 | Result(2) = C.X(); |
| 637 | Result(3) = C.Y(); |
| 638 | Result(4) = C.Z(); |
| 639 | break; |
| 640 | |
| 641 | case 1: |
| 642 | Cur2D->D1(U, C2D, dC2D_dU); |
| 643 | dV_dU = dC2D_dU.X(); |
| 644 | dW_dU = dC2D_dU.Y(); |
| 645 | Surf->D1(C2D.X(), C2D.Y(), C, dS_dV, dS_dW); |
| 646 | dC_dU = dS_dV*dV_dU + dS_dW*dW_dU; |
| 647 | Mag = dC_dU.Magnitude(); |
| 648 | dU_dS = Length/Mag; |
| 649 | |
| 650 | dV_dS = dV_dU*dU_dS; |
| 651 | dW_dS = dW_dU*dU_dS; |
| 652 | dC_dS = dC_dU*dU_dS; |
| 653 | |
| 654 | Result(0) = dV_dS; |
| 655 | Result(1) = dW_dS; |
| 656 | Result(2) = dC_dS.X(); |
| 657 | Result(3) = dC_dS.Y(); |
| 658 | Result(4) = dC_dS.Z(); |
| 659 | break; |
| 660 | |
| 661 | case 2: |
| 662 | Cur2D->D2(U, C2D, dC2D_dU, d2C2D_dU2); |
| 663 | dV_dU = dC2D_dU.X(); |
| 664 | dW_dU = dC2D_dU.Y(); |
| 665 | d2V_dU2 = d2C2D_dU2.X(); |
| 666 | d2W_dU2 = d2C2D_dU2.Y(); |
| 667 | Surf->D2(C2D.X(), C2D.Y(), C, dS_dV, dS_dW, d2S_dV2, d2S_dW2, d2S_dVdW); |
| 668 | dC_dU = dS_dV*dV_dU + dS_dW*dW_dU; |
| 669 | d2C_dU2 = (d2S_dV2*dV_dU + d2S_dVdW*dW_dU)*dV_dU + dS_dV*d2V_dU2 + |
| 670 | (d2S_dVdW*dV_dU + d2S_dW2*dW_dU)*dW_dU + dS_dW*d2W_dU2; |
| 671 | Mag = dC_dU.Magnitude(); |
| 672 | dU_dS = Length/Mag; |
| 673 | d2U_dS2 = -Length*dC_dU.Dot(d2C_dU2)*dU_dS/(Mag*Mag*Mag); |
| 674 | |
| 675 | dV_dS = dV_dU * dU_dS; |
| 676 | dW_dS = dW_dU * dU_dS; |
| 677 | d2V_dS2 = d2V_dU2*dU_dS*dU_dS + dV_dU*d2U_dS2; |
| 678 | d2W_dS2 = d2W_dU2*dU_dS*dU_dS + dW_dU*d2U_dS2; |
| 679 | |
| 680 | d2U_dS2 = -dC_dU.Dot(d2C_dU2)*dU_dS/(Mag*Mag); |
| 681 | d2C_dS2 = (d2S_dV2 * dV_dS + d2S_dVdW * dW_dS) * dV_dS + dS_dV * d2V_dS2 + |
| 682 | (d2S_dW2 * dW_dS + d2S_dVdW * dV_dS) * dW_dS + dS_dW * d2W_dS2; |
| 683 | |
| 684 | Result(0) = d2V_dS2; |
| 685 | Result(1) = d2W_dS2; |
| 686 | Result(2) = d2C_dS2.X(); |
| 687 | Result(3) = d2C_dS2.Y(); |
| 688 | Result(4) = d2C_dS2.Z(); |
| 689 | break; |
| 690 | |
| 691 | default: Result(0) = Result(1) = Result(2) = Result(3) = Result(4) = 0; |
| 692 | return Standard_False; |
| 693 | } |
| 694 | return Standard_True; |
| 695 | } |