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

// File:	Approx_CurvlinFunc.cxx
// Created:	Tue May 12 14:03:09 1998
// Author:	Roman BORISOV
//		<rbv@sgi38>


#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>

static Standard_Real myPrevS, myPrevU;
#ifdef DEB
#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)
{
  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)
{  
  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)
{  
  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;

  myPrevS = myFirstS;
  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 DEB
  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();

  myPrevS = S;
  myPrevU = U;

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