[occt.git] / src / Approx / Approx_CurvlinFunc.cxx

// Created on: 1998-05-12
// Created by: Roman BORISOV
// Copyright (c) 1998-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.

#include <Approx_CurvlinFunc.ixx>
#include <Adaptor3d_CurveOnSurface.hxx>
#include <Adaptor3d_HCurveOnSurface.hxx>
#include <TColStd_SequenceOfReal.hxx>
#include <GeomLib.hxx>
#include <GCPnts_AbscissaPoint.hxx>
#include <Precision.hxx>

#ifdef __OCC_DEBUG_CHRONO
#include <OSD_Timer.hxx>
static OSD_Chronometer chr_uparam;
Standard_EXPORT Standard_Integer uparam_count;
Standard_EXPORT Standard_Real t_uparam;

//Standard_IMPORT extern void InitChron(OSD_Chronometer& ch);
Standard_IMPORT void InitChron(OSD_Chronometer& ch);
//Standard_IMPORT extern void ResultChron( OSD_Chronometer & ch, Standard_Real & time);
Standard_IMPORT void ResultChron( OSD_Chronometer & ch, Standard_Real & time);
#endif

static Standard_Real cubic(const Standard_Real X, const Standard_Real *Xi, const Standard_Real *Yi)
{
  Standard_Real I1, I2, I3, I21, I22, I31, Result;

  I1 = (Yi[0] - Yi[1])/(Xi[0] - Xi[1]);
  I2 = (Yi[1] - Yi[2])/(Xi[1] - Xi[2]);
  I3 = (Yi[2] - Yi[3])/(Xi[2] - Xi[3]);

  I21 = (I1 - I2)/(Xi[0] - Xi[2]);
  I22 = (I2 - I3)/(Xi[1] - Xi[3]);
  
  I31 = (I21 - I22)/(Xi[0] - Xi[3]);  

  Result = Yi[0] + (X - Xi[0])*(I1 + (X - Xi[1])*(I21 + (X - Xi[2])*I31));

  return Result;
}

//static void findfourpoints(const Standard_Real S, 
static void findfourpoints(const Standard_Real , 
			   Standard_Integer NInterval, 
			   const Handle(TColStd_HArray1OfReal)& Si, 
			   Handle(TColStd_HArray1OfReal)& Ui, 
			   const Standard_Real prevS, 
			   const Standard_Real prevU, Standard_Real *Xi, 
			   Standard_Real *Yi)
{
  Standard_Integer i, j;
  Standard_Integer NbInt = Si->Length() - 1;
  if (NbInt < 3) Standard_ConstructionError::Raise("Approx_CurvlinFunc::GetUParameter");

  if(NInterval < 1) NInterval = 1;
  else if(NInterval > NbInt - 2) NInterval = NbInt - 2;

  for(i = 0; i < 4; i++) {
    Xi[i] = Si->Value(NInterval - 1 + i);
    Yi[i] = Ui->Value(NInterval - 1 + i);
  }
// try to insert (S, U)  
  for(i = 0; i < 3; i++) {
    if(Xi[i] < prevS && prevS < Xi[i+1]) {
      for(j = 0; j < i; j++) {
	Xi[j] = Xi[j+1];
	Yi[j] = Yi[j+1];
      }
      Xi[i] = prevS;
      Yi[i] = prevU;
      break;
    }
  }
}

/*static Standard_Real curvature(const Standard_Real U, const Adaptor3d_Curve& C)
{
  Standard_Real k, tau, mod1, mod2, OMEGA;
  gp_Pnt P;
  gp_Vec D1, D2, D3;
  C.D3(U, P, D1, D2, D3);
  mod1 = D1.Magnitude();
  mod2 = D1.Crossed(D2).Magnitude();
  k = mod2/(mod1*mod1*mod1);
  tau = D1.Dot(D2.Crossed(D3));
  tau /= mod2*mod2;
  OMEGA = Sqrt(k*k + tau*tau);

  return OMEGA;
}
*/

Approx_CurvlinFunc::Approx_CurvlinFunc(const Handle(Adaptor3d_HCurve)& C, const Standard_Real Tol) : myC3D(C), 
                    myCase(1), 
                    myFirstS(0), 
                    myLastS(1), 
                    myTolLen(Tol),
                    myPrevS (0.0),
                    myPrevU (0.0)
{
  Init();
}

Approx_CurvlinFunc::Approx_CurvlinFunc(const Handle(Adaptor2d_HCurve2d)& C2D, const Handle(Adaptor3d_HSurface)& S, const Standard_Real Tol) : 
                    myC2D1(C2D), 
                    mySurf1(S), 
                    myCase(2), 
                    myFirstS(0), 
                    myLastS(1), 
                    myTolLen(Tol),
                    myPrevS (0.0),
                    myPrevU (0.0)
{  
  Init();
}

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) : 
                    myC2D1(C2D1), 
                    myC2D2(C2D2), 
                    mySurf1(S1), 
                    mySurf2(S2), 
                    myCase(3), 
                    myFirstS(0), 
                    myLastS(1), 
                    myTolLen(Tol),
                    myPrevS (0.0),
                    myPrevU (0.0)
{  
  Init();
}

void Approx_CurvlinFunc::Init()
{
  Adaptor3d_CurveOnSurface CurOnSur;
  
  switch(myCase) {
  case 1:
    Init(myC3D->GetCurve(), mySi_1, myUi_1);
    myFirstU1 = myC3D->FirstParameter();
    myLastU1 = myC3D->LastParameter();
    myFirstU2 = myLastU2 = 0;
    break;
  case 2:
    CurOnSur.Load(myC2D1);
    CurOnSur.Load(mySurf1);
    Init(CurOnSur, mySi_1, myUi_1);
    myFirstU1 = CurOnSur.FirstParameter();
    myLastU1 = CurOnSur.LastParameter();
    myFirstU2 = myLastU2 = 0;
    break;
  case 3:
    CurOnSur.Load(myC2D1);
    CurOnSur.Load(mySurf1);
    Init(CurOnSur, mySi_1, myUi_1);
    myFirstU1 = CurOnSur.FirstParameter();
    myLastU1 = CurOnSur.LastParameter();
    CurOnSur.Load(myC2D2);
    CurOnSur.Load(mySurf2);
    Init(CurOnSur, mySi_2, myUi_2);
    myFirstU2 = CurOnSur.FirstParameter();
    myLastU2 = CurOnSur.LastParameter();
  }

  Length();
}


//=======================================================================
//function : Init
//purpose  : Init the values
//history  : 23/10/1998 PMN : Cut at curve's discontinuities
//=======================================================================
void Approx_CurvlinFunc::Init(Adaptor3d_Curve& C, Handle(TColStd_HArray1OfReal)& Si, 
			      Handle(TColStd_HArray1OfReal)& Ui) const
{
  Standard_Real Step, FirstU, LastU;
  Standard_Integer i, j, k, NbInt, NbIntC3;
  FirstU = C.FirstParameter();
  LastU  = C.LastParameter();

  NbInt = 10;
  NbIntC3 = C.NbIntervals(GeomAbs_C3);
  TColStd_Array1OfReal Disc(1, NbIntC3+1);

  if (NbIntC3 >1) {
     C.Intervals(Disc, GeomAbs_C3);
  }
  else {
    Disc(1) = FirstU;
    Disc(2) = LastU;
  }

  Ui = new TColStd_HArray1OfReal (0,NbIntC3*NbInt);
  Si = new TColStd_HArray1OfReal (0,NbIntC3*NbInt);

  Ui->SetValue(0, FirstU);
  Si->SetValue(0, 0);

  for(j = 1, i=1; j<=NbIntC3; j++) {
    Step = (Disc(j+1) - Disc(j))/NbInt;
    for(k = 1; k <= NbInt; k++, i++) {
      Ui->ChangeValue(i) = Ui->Value(i-1) + Step;
      Si->ChangeValue(i) = Si->Value(i-1) + Length(C, Ui->Value(i-1), Ui->Value(i));
    }
  }

  Standard_Real Len = Si->Value(Si->Upper());
  for(i = Si->Lower(); i<= Si->Upper(); i++)
    Si->ChangeValue(i) /= Len;

  // TODO - fields should be mutable
  const_cast<Approx_CurvlinFunc*>(this)->myPrevS = myFirstS;
  const_cast<Approx_CurvlinFunc*>(this)->myPrevU = FirstU;
}

void  Approx_CurvlinFunc::SetTol(const Standard_Real Tol)
{
  myTolLen = Tol;
}

Standard_Real Approx_CurvlinFunc::FirstParameter() const
{
  return myFirstS;
}

Standard_Real Approx_CurvlinFunc::LastParameter() const
{
  return myLastS;
}

Standard_Integer Approx_CurvlinFunc::NbIntervals(const GeomAbs_Shape S) const
{
  Adaptor3d_CurveOnSurface CurOnSur;

  switch(myCase) {
  case 1:
    return myC3D->NbIntervals(S);
  case 2:
    CurOnSur.Load(myC2D1);
    CurOnSur.Load(mySurf1);
    return CurOnSur.NbIntervals(S);
  case 3:
    Standard_Integer NbInt;
    CurOnSur.Load(myC2D1);
    CurOnSur.Load(mySurf1);
    NbInt = CurOnSur.NbIntervals(S);
    TColStd_Array1OfReal T1(1, NbInt+1);
    CurOnSur.Intervals(T1, S);
    CurOnSur.Load(myC2D2);
    CurOnSur.Load(mySurf2);
    NbInt = CurOnSur.NbIntervals(S);
    TColStd_Array1OfReal T2(1, NbInt+1);
    CurOnSur.Intervals(T2, S);

    TColStd_SequenceOfReal Fusion;
    GeomLib::FuseIntervals(T1, T2, Fusion);
    return Fusion.Length() - 1;
  }

  //POP pour WNT
  return 1;
}

void Approx_CurvlinFunc::Intervals(TColStd_Array1OfReal& T, const GeomAbs_Shape S) const
{
  Adaptor3d_CurveOnSurface CurOnSur;
  Standard_Integer i;
    
  switch(myCase) {
  case 1:
    myC3D->Intervals(T, S);
    break;
  case 2:
    CurOnSur.Load(myC2D1);
    CurOnSur.Load(mySurf1);
    CurOnSur.Intervals(T, S);
    break;
  case 3:
    Standard_Integer NbInt;
    CurOnSur.Load(myC2D1);
    CurOnSur.Load(mySurf1);
    NbInt = CurOnSur.NbIntervals(S);
    TColStd_Array1OfReal T1(1, NbInt+1);
    CurOnSur.Intervals(T1, S);
    CurOnSur.Load(myC2D2);
    CurOnSur.Load(mySurf2);
    NbInt = CurOnSur.NbIntervals(S);
    TColStd_Array1OfReal T2(1, NbInt+1);
    CurOnSur.Intervals(T2, S);

    TColStd_SequenceOfReal Fusion;
    GeomLib::FuseIntervals(T1, T2, Fusion);

    for (i = 1; i <= Fusion.Length(); i++)
      T.ChangeValue(i) = Fusion.Value(i);
  }

  for(i = 1; i <= T.Length(); i++)
    T.ChangeValue(i) = GetSParameter(T.Value(i));
}

void Approx_CurvlinFunc::Trim(const Standard_Real First, const Standard_Real Last, const Standard_Real Tol)
{
  if (First < 0 || Last >1) Standard_OutOfRange::Raise("Approx_CurvlinFunc::Trim");
  if ((Last - First) < Tol) return;

  Standard_Real FirstU, LastU;
  Adaptor3d_CurveOnSurface CurOnSur;
  Handle(Adaptor3d_HCurve) HCurOnSur;

  switch(myCase) {
  case 1:
    myC3D = myC3D->Trim(myFirstU1, myLastU1, Tol);
    FirstU = GetUParameter(myC3D->GetCurve(), First, 1);
    LastU = GetUParameter(myC3D->GetCurve(), Last, 1);
    myC3D = myC3D->Trim(FirstU, LastU, Tol);
    break;
  case 3:
    CurOnSur.Load(myC2D2);
    CurOnSur.Load(mySurf2);
    HCurOnSur = CurOnSur.Trim(myFirstU2, myLastU2, Tol);
    myC2D2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve();
    mySurf2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface();
    CurOnSur.Load(myC2D2);
    CurOnSur.Load(mySurf2);

    FirstU = GetUParameter(CurOnSur, First, 1);
    LastU = GetUParameter(CurOnSur, Last, 1);
    HCurOnSur = CurOnSur.Trim(FirstU, LastU, Tol);
    myC2D2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve();
    mySurf2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface();

  case 2:
    CurOnSur.Load(myC2D1);
    CurOnSur.Load(mySurf1);
    HCurOnSur = CurOnSur.Trim(myFirstU1, myLastU1, Tol);
    myC2D1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve();
    mySurf1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface();
    CurOnSur.Load(myC2D1);
    CurOnSur.Load(mySurf1);

    FirstU = GetUParameter(CurOnSur, First, 1);
    LastU = GetUParameter(CurOnSur, Last, 1);
    HCurOnSur = CurOnSur.Trim(FirstU, LastU, Tol);
    myC2D1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve();
    mySurf1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface();
  }
  myFirstS = First;
  myLastS = Last;  
}

void Approx_CurvlinFunc::Length()
{
  Adaptor3d_CurveOnSurface CurOnSur;
  Standard_Real FirstU, LastU;

  switch(myCase){
  case 1:
    FirstU = myC3D->FirstParameter();
    LastU = myC3D->LastParameter();
    myLength = Length(myC3D->GetCurve(), FirstU, LastU);    
    myLength1 = myLength2 = 0;    
    break;
  case 2:
    CurOnSur.Load(myC2D1);
    CurOnSur.Load(mySurf1);
    FirstU = CurOnSur.FirstParameter();
    LastU = CurOnSur.LastParameter();
    myLength = Length(CurOnSur, FirstU, LastU);
    myLength1 = myLength2 = 0;    
    break;
  case 3:
    CurOnSur.Load(myC2D1);
    CurOnSur.Load(mySurf1);
    FirstU = CurOnSur.FirstParameter();
    LastU = CurOnSur.LastParameter();
    myLength1 = Length(CurOnSur, FirstU, LastU);
    CurOnSur.Load(myC2D2);
    CurOnSur.Load(mySurf2);
    FirstU = CurOnSur.FirstParameter();
    LastU = CurOnSur.LastParameter();
    myLength2 = Length(CurOnSur, FirstU, LastU);
    myLength = (myLength1 + myLength2)/2;
  }
}


Standard_Real Approx_CurvlinFunc::Length(Adaptor3d_Curve& C, const Standard_Real FirstU, const Standard_Real LastU) const
{
  Standard_Real Length;

  Length = GCPnts_AbscissaPoint::Length(C, FirstU, LastU, myTolLen);
  return Length;
}


Standard_Real Approx_CurvlinFunc::GetLength() const
{
  return myLength;
}

Standard_Real Approx_CurvlinFunc::GetSParameter(const Standard_Real U) const
{
  Standard_Real S=0, S1, S2;
  Adaptor3d_CurveOnSurface CurOnSur;

  switch (myCase) {
  case 1:
    S = GetSParameter(myC3D->GetCurve(), U, myLength);
    break;
  case 2:
    CurOnSur.Load(myC2D1);
    CurOnSur.Load(mySurf1);
    S = GetSParameter(CurOnSur, U, myLength);
    break;
  case 3:
    CurOnSur.Load(myC2D1);
    CurOnSur.Load(mySurf1);
    S1 = GetSParameter(CurOnSur, U, myLength1);    
    CurOnSur.Load(myC2D2);
    CurOnSur.Load(mySurf2);
    S2 = GetSParameter(CurOnSur, U, myLength2);    
    S = (S1 + S2)/2;
  }
  return S;
}


Standard_Real Approx_CurvlinFunc::GetUParameter(Adaptor3d_Curve& C, 
						const Standard_Real S, 
						const Standard_Integer NumberOfCurve) const
{
  Standard_Real deltaS, base, U, Length;
  Standard_Integer NbInt, NInterval, i;
  Handle(TColStd_HArray1OfReal) InitUArray, InitSArray;
#ifdef __OCC_DEBUG_CHRONO
  InitChron(chr_uparam);
#endif
  if(S < 0 || S > 1) Standard_ConstructionError::Raise("Approx_CurvlinFunc::GetUParameter");

  if(NumberOfCurve == 1) {
    InitUArray = myUi_1;
    InitSArray = mySi_1;
    if(myCase == 3)
      Length = myLength1;
    else 
      Length = myLength;
  }
  else {
    InitUArray = myUi_2;
    InitSArray = mySi_2;
    Length = myLength2;
  }

  NbInt = InitUArray->Length() - 1;

  if(S == 1) NInterval = NbInt - 1;
  else {
    for(i = 0; i < NbInt; i++) {
      if((InitSArray->Value(i) <= S && S < InitSArray->Value(i+1)))
	break;
    }
    NInterval = i;
  }
  if(S==InitSArray->Value(NInterval)) {
    return InitUArray->Value(NInterval);
  }
  if(S==InitSArray->Value(NInterval+1)) {
    return InitUArray->Value(NInterval+1);
  } 

  base = InitUArray->Value(NInterval);
  deltaS = (S - InitSArray->Value(NInterval))*Length;

// to find an initial point
  Standard_Real Xi[4], Yi[4], UGuess;
  findfourpoints(S, NInterval, InitSArray, InitUArray, myPrevS, myPrevU, Xi, Yi);
  UGuess = cubic(S , Xi, Yi);

  U = GCPnts_AbscissaPoint(C, deltaS, base, UGuess, myTolLen).Parameter();

  // TODO - fields should be mutable
  const_cast<Approx_CurvlinFunc*>(this)->myPrevS = S;
  const_cast<Approx_CurvlinFunc*>(this)->myPrevU = U;

#ifdef __OCC_DEBUG_CHRONO
  ResultChron(chr_uparam, t_uparam);
  uparam_count++;
#endif

  return U;
}

Standard_Real Approx_CurvlinFunc::GetSParameter(Adaptor3d_Curve& C, const Standard_Real U, const Standard_Real Len) const
{
  Standard_Real S, Origin;

  Origin = C.FirstParameter();
  S = myFirstS + Length(C, Origin, U)/Len;    
  return S;
}

Standard_Boolean Approx_CurvlinFunc::EvalCase1(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result) const
{
  if(myCase != 1) Standard_ConstructionError::Raise("Approx_CurvlinFunc::EvalCase1");

  gp_Pnt C;
  gp_Vec dC_dU, dC_dS, d2C_dU2, d2C_dS2;
  Standard_Real U, Mag, dU_dS, d2U_dS2;
  
  U = GetUParameter(myC3D->GetCurve(), S, 1);

  switch(Order) {

  case 0: 
    myC3D->D0(U, C);

    Result(0) = C.X();
    Result(1) = C.Y();
    Result(2) = C.Z();
    break;
    
  case 1:
    myC3D->D1(U, C, dC_dU);
    Mag = dC_dU.Magnitude();
    dU_dS = myLength/Mag;
    dC_dS = dC_dU*dU_dS;

    Result(0) = dC_dS.X();
    Result(1) = dC_dS.Y();
    Result(2) = dC_dS.Z();
    break;

  case 2:
    myC3D->D2(U, C, dC_dU, d2C_dU2);
    Mag = dC_dU.Magnitude();
    dU_dS = myLength/Mag;
    d2U_dS2 = -myLength*dC_dU.Dot(d2C_dU2)*dU_dS/(Mag*Mag*Mag);
    d2C_dS2 = d2C_dU2*dU_dS*dU_dS + dC_dU*d2U_dS2;

    Result(0) = d2C_dS2.X();
    Result(1) = d2C_dS2.Y();
    Result(2) = d2C_dS2.Z();
    break;

  default: Result(0) = Result(1) = Result(2) = 0;
    return Standard_False;
  }
  return Standard_True;
}

Standard_Boolean Approx_CurvlinFunc::EvalCase2(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result) const
{
  if(myCase != 2) Standard_ConstructionError::Raise("Approx_CurvlinFunc::EvalCase2");

  Standard_Boolean Done;

  Done = EvalCurOnSur(S, Order, Result, 1);

  return Done;
}

Standard_Boolean Approx_CurvlinFunc::EvalCase3(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result)
{
  if(myCase != 3) Standard_ConstructionError::Raise("Approx_CurvlinFunc::EvalCase3");
  
  TColStd_Array1OfReal tmpRes1(0, 4), tmpRes2(0, 4);
  Standard_Boolean Done;

  Done = EvalCurOnSur(S, Order, tmpRes1, 1);

  Done = EvalCurOnSur(S, Order, tmpRes2, 2) && Done;

  Result(0) = tmpRes1(0);
  Result(1) = tmpRes1(1);
  Result(2) = tmpRes2(0);
  Result(3) = tmpRes2(1);
  Result(4) = 0.5*(tmpRes1(2) + tmpRes2(2));
  Result(5) = 0.5*(tmpRes1(3) + tmpRes2(3));
  Result(6) = 0.5*(tmpRes1(4) + tmpRes2(4));

  return Done;
}

Standard_Boolean Approx_CurvlinFunc::EvalCurOnSur(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result, const  Standard_Integer NumberOfCurve) const
{
  Handle(Adaptor2d_HCurve2d) Cur2D;
  Handle(Adaptor3d_HSurface) Surf;
  Standard_Real U=0, Length=0;

  if (NumberOfCurve == 1) {
    Cur2D = myC2D1;
    Surf = mySurf1;    
    Adaptor3d_CurveOnSurface CurOnSur(myC2D1, mySurf1);
    U = GetUParameter(CurOnSur, S, 1);
    if(myCase == 3) Length = myLength1;
    else Length = myLength;
  }
  else if (NumberOfCurve == 2) {
    Cur2D = myC2D2;
    Surf = mySurf2;    
    Adaptor3d_CurveOnSurface CurOnSur(myC2D2, mySurf2);
    U = GetUParameter(CurOnSur, S, 2);
    Length = myLength2;
  }
  else 
    Standard_ConstructionError::Raise("Approx_CurvlinFunc::EvalCurOnSur");

  Standard_Real Mag, dU_dS, d2U_dS2, dV_dU, dW_dU, dV_dS, dW_dS, d2V_dS2, d2W_dS2, d2V_dU2, d2W_dU2;
  gp_Pnt2d C2D;
  gp_Pnt C;
  gp_Vec2d dC2D_dU, d2C2D_dU2;
  gp_Vec dC_dU, d2C_dU2, dC_dS, d2C_dS2, dS_dV, dS_dW, d2S_dV2, d2S_dW2, d2S_dVdW;

  switch(Order) {
  case 0:
    Cur2D->D0(U, C2D);
    Surf->D0(C2D.X(), C2D.Y(), C);
    
    Result(0) = C2D.X();
    Result(1) = C2D.Y();
    Result(2) = C.X();
    Result(3) = C.Y();
    Result(4) = C.Z();
    break;

  case 1:
    Cur2D->D1(U, C2D, dC2D_dU);
    dV_dU = dC2D_dU.X();
    dW_dU = dC2D_dU.Y();
    Surf->D1(C2D.X(), C2D.Y(), C, dS_dV, dS_dW);
    dC_dU = dS_dV*dV_dU + dS_dW*dW_dU;
    Mag = dC_dU.Magnitude();
    dU_dS = Length/Mag;

    dV_dS = dV_dU*dU_dS;
    dW_dS = dW_dU*dU_dS;
    dC_dS = dC_dU*dU_dS;

    Result(0) = dV_dS;
    Result(1) = dW_dS;
    Result(2) = dC_dS.X();
    Result(3) = dC_dS.Y();
    Result(4) = dC_dS.Z();    
    break;

  case 2:
    Cur2D->D2(U, C2D, dC2D_dU, d2C2D_dU2);
    dV_dU = dC2D_dU.X();
    dW_dU = dC2D_dU.Y();
    d2V_dU2 = d2C2D_dU2.X();
    d2W_dU2 = d2C2D_dU2.Y();
    Surf->D2(C2D.X(), C2D.Y(), C, dS_dV, dS_dW, d2S_dV2, d2S_dW2, d2S_dVdW);
    dC_dU = dS_dV*dV_dU + dS_dW*dW_dU;
    d2C_dU2 = (d2S_dV2*dV_dU + d2S_dVdW*dW_dU)*dV_dU + dS_dV*d2V_dU2 + 
              (d2S_dVdW*dV_dU + d2S_dW2*dW_dU)*dW_dU + dS_dW*d2W_dU2;
    Mag = dC_dU.Magnitude();
    dU_dS = Length/Mag;
    d2U_dS2 = -Length*dC_dU.Dot(d2C_dU2)*dU_dS/(Mag*Mag*Mag);
    
    dV_dS = dV_dU * dU_dS;
    dW_dS = dW_dU * dU_dS;
    d2V_dS2 = d2V_dU2*dU_dS*dU_dS + dV_dU*d2U_dS2;
    d2W_dS2 = d2W_dU2*dU_dS*dU_dS + dW_dU*d2U_dS2;
    
    d2U_dS2 = -dC_dU.Dot(d2C_dU2)*dU_dS/(Mag*Mag);
    d2C_dS2 = (d2S_dV2 * dV_dS + d2S_dVdW * dW_dS) * dV_dS + dS_dV * d2V_dS2 +
              (d2S_dW2 * dW_dS + d2S_dVdW * dV_dS) * dW_dS + dS_dW * d2W_dS2;
    
    Result(0) = d2V_dS2;
    Result(1) = d2W_dS2;
    Result(2) = d2C_dS2.X();
    Result(3) = d2C_dS2.Y();
    Result(4) = d2C_dS2.Z();      
    break;

  default: Result(0) = Result(1) = Result(2) = Result(3) = Result(4) = 0;
    return Standard_False;
  }
  return Standard_True;
}
Commit	Line	Data
b311480e	1	// Created on: 1998-05-12
	2	// Created by: Roman BORISOV
	3	// Copyright (c) 1998-1999 Matra Datavision
973c2be1	4	// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e	5	//
973c2be1	6	// This file is part of Open CASCADE Technology software library.
b311480e	7	//
d5f74e42	8	// This library is free software; you can redistribute it and/or modify it under
d5f74e42	9	// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1	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.
b311480e	13	//
973c2be1	14	// Alternatively, this file may be used under the terms of Open CASCADE
973c2be1	15	// commercial license or contractual agreement.
7fd59977	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>
b311480e	22	#include <GCPnts_AbscissaPoint.hxx>
7fd59977	23	#include <Precision.hxx>
7fd59977	24
41194117	25	#ifdef __OCC_DEBUG_CHRONO
7fd59977	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/(mod1mod1mod1);
98	tau = D1.Dot(D2.Crossed(D3));
99	tau /= mod2*mod2;
100	OMEGA = Sqrt(kk + tautau);
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),
41194117 K	110	myTolLen(Tol),
	111	myPrevS (0.0),
	112	myPrevU (0.0)
7fd59977	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),
41194117 K	123	myTolLen(Tol),
	124	myPrevS (0.0),
	125	myPrevU (0.0)
7fd59977	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),
41194117 K	138	myTolLen(Tol),
	139	myPrevS (0.0),
	140	myPrevU (0.0)
7fd59977	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
41194117 K	224	// TODO - fields should be mutable
	225	const_cast<Approx_CurvlinFunc*>(this)->myPrevS = myFirstS;
	226	const_cast<Approx_CurvlinFunc*>(this)->myPrevU = FirstU;
7fd59977	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;
41194117	450	#ifdef __OCC_DEBUG_CHRONO
7fd59977	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
41194117 K	496	// TODO - fields should be mutable
	497	const_cast<Approx_CurvlinFunc*>(this)->myPrevS = S;
	498	const_cast<Approx_CurvlinFunc*>(this)->myPrevU = U;
7fd59977	499
41194117	500	#ifdef __OCC_DEBUG_CHRONO
7fd59977	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 = -myLengthdC_dU.Dot(d2C_dU2)dU_dS/(MagMagMag);
	553	d2C_dS2 = d2C_dU2dU_dSdU_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_dVdV_dU + dS_dWdW_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_dVdV_dU + dS_dWdW_dU;
669	d2C_dU2 = (d2S_dV2dV_dU + d2S_dVdWdW_dU)dV_dU + dS_dVd2V_dU2 +
670	(d2S_dVdWdV_dU + d2S_dW2dW_dU)dW_dU + dS_dWd2W_dU2;
671	Mag = dC_dU.Magnitude();
672	dU_dS = Length/Mag;
673	d2U_dS2 = -LengthdC_dU.Dot(d2C_dU2)dU_dS/(MagMagMag);
674
675	dV_dS = dV_dU * dU_dS;
676	dW_dS = dW_dU * dU_dS;
677	d2V_dS2 = d2V_dU2dU_dSdU_dS + dV_dU*d2U_dS2;
678	d2W_dS2 = d2W_dU2dU_dSdU_dS + dW_dU*d2U_dS2;
679
680	d2U_dS2 = -dC_dU.Dot(d2C_dU2)dU_dS/(MagMag);
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	}