// File:	FEmTool_Curve.cxx
// Created:	Tue Nov  4 10:45:06 1997
// Author:	Roman BORISOV / Igor FEOKTISTOV
//		<rbv@redfox.nnov.matra-dtv.fr>


#define No_Standard_RangeError
#define No_Standard_OutOfRange

#include <FEmTool_Curve.ixx>
#include <PLib.hxx>
#include <PLib_JacobiPolynomial.hxx>
#include <PLib_HermitJacobi.hxx>

//=======================================================================
//function : FEmTool_Curve
//purpose  :
//=======================================================================
FEmTool_Curve::FEmTool_Curve(const Standard_Integer Dimension,
			     const Standard_Integer NbElements,
			     const Handle(PLib_Base)& TheBase, 
			     const Standard_Real Tolerance) : 
       myNbElements(NbElements), myDimension(Dimension), 
       myTolerance(Tolerance), 
       myBase(TheBase), myDegree(1, myNbElements), 
       myCoeff(1, myDimension*myNbElements*(myBase->WorkDegree() + 1)), 
       myPoly(1, myDimension*myNbElements*(myBase->WorkDegree() + 1)), 
       myDeri(1, myDimension*myNbElements*(myBase->WorkDegree())), 
       myDsecn(1, myDimension*myNbElements*(myBase->WorkDegree() - 1)), 
       HasPoly(1, myNbElements), HasDeri(1, myNbElements), 
       HasSecn(1, myNbElements), myLength(1, myNbElements),
       myIndex(0)
{
  myKnots = new TColStd_HArray1OfReal(1, myNbElements + 1);
  myDegree.Init(myBase->WorkDegree());
  HasPoly.Init(0);
  HasDeri.Init(0);
  HasSecn.Init(0);
  myLength.Init(-1);
}

 TColStd_Array1OfReal& FEmTool_Curve::Knots() const
{
  return myKnots->ChangeArray1();
}



//=======================================================================
//function : SetElement
//purpose  :
//=======================================================================
 void FEmTool_Curve::SetElement(const Standard_Integer IndexOfElement,
				const TColStd_Array2OfReal& Coeffs) 
{
  Standard_Integer i, j, degBase, deg;
  if (IndexOfElement > myNbElements || IndexOfElement < 1) Standard_OutOfRange::Raise(); 
  degBase = myBase->WorkDegree();
  deg = myDegree(IndexOfElement);
  Standard_Integer iBase = (IndexOfElement - 1)*(degBase + 1)*myDimension, 
                   i1 = iBase-myDimension, i2 = Coeffs.LowerRow() - 1, 
                   j1 = Coeffs.LowerCol() - 1;
  for(i = 1; i <= deg + 1; i++) {
    i1 += myDimension; i2++; 
    for(j = 1; j <= myDimension; j++)
      myCoeff(i1 + j) = Coeffs(i2, j1 + j);
  }

  Standard_Real stenor = (myKnots->Value(IndexOfElement + 1) - myKnots->Value(IndexOfElement)) / 2.,
                mfact;
  Handle(PLib_HermitJacobi) myHermitJacobi = (*((Handle(PLib_HermitJacobi)*)&myBase));

  i1 = iBase;
  i2 = iBase + (myHermitJacobi->NivConstr() + 1) * myDimension;
  for(i = 1; i <= myHermitJacobi->NivConstr(); i++) {
    i1 += myDimension; i2 += myDimension; 
    mfact = Pow(stenor, i);
    for(j = 1; j <= myDimension; j++) {
      myCoeff(i1 + j) *= mfact;
      myCoeff(i2 + j) *= mfact;
    }
  }

  HasPoly(IndexOfElement) = HasDeri(IndexOfElement) = HasSecn(IndexOfElement) = 0;
  myLength(IndexOfElement) = - 1;
}


//=======================================================================
//function : GetElement
//purpose  :
//=======================================================================
 void FEmTool_Curve::GetElement(const Standard_Integer IndexOfElement, TColStd_Array2OfReal& Coeffs) 
{
  Standard_Integer i, j, degBase, deg;
  if (IndexOfElement > myNbElements || IndexOfElement < 1) Standard_OutOfRange::Raise(); 
  degBase = myBase->WorkDegree();
  deg = myDegree(IndexOfElement);
  Standard_Integer iBase = (IndexOfElement - 1)*(degBase + 1)*myDimension, 
                   i1 = iBase-myDimension, i2 = Coeffs.LowerRow() - 1, 
                   j1 = Coeffs.LowerCol() - 1;
  for(i = 1; i <= deg + 1; i++) {
    i1 += myDimension; i2++; 
    for(j = 1; j <= myDimension; j++)
      Coeffs(i2, j1 + j) = myCoeff.Value(i1 + j);
  }

  Standard_Real stenor = 2. / (myKnots->Value(IndexOfElement + 1) - myKnots->Value(IndexOfElement)),
                mfact;

  Handle(PLib_HermitJacobi) myHermitJacobi = (*((Handle(PLib_HermitJacobi)*)&myBase));

  i2 = Coeffs.LowerRow();
  Standard_Integer i3 = i2 + myHermitJacobi->NivConstr() + 1;

  for(i = 1; i <= myHermitJacobi->NivConstr(); i++) {
    mfact = Pow(stenor, i);
    for(j = j1+1; j <= myDimension; j++) {
      Coeffs(i2 + i, j) *= mfact;
      Coeffs(i3 + i, j) *= mfact;
    }
  }
}


//=======================================================================
//function : GetPolynom
//purpose  :
//=======================================================================
 void FEmTool_Curve::GetPolynom(TColStd_Array1OfReal& Coeffs) 

{
  Standard_Integer IndexOfElement, i, di = Coeffs.Lower() - myPoly.Lower();


  for(IndexOfElement = 1; IndexOfElement <= myNbElements; IndexOfElement++)
    if (!HasPoly.Value(IndexOfElement)) Update(IndexOfElement, 0);

  for(i = myPoly.Lower(); i <= myPoly.Upper(); i++)
    Coeffs(di + i) = myPoly(i);

}


//=======================================================================
//function : D0
//purpose  :
//=======================================================================
 void FEmTool_Curve::D0(const Standard_Real U, TColStd_Array1OfReal& Pnt) 
{
  Standard_Integer deg;
  Standard_Real S;

  if (!myIndex || (U<Uf) || (U>Ul) ||
      (myKnots->Value(myIndex)!=Uf) ||
      (myKnots->Value(myIndex+1)!=Ul)) { 
    // Search the span
    if(U <= myKnots->Value(2)) myIndex = 1;
    else {
      for(myIndex = 2; myIndex <= myNbElements; myIndex++)
	if(U >= myKnots->Value(myIndex) && U <= myKnots->Value(myIndex+1)) break;
      if (myIndex > myNbElements) myIndex = myNbElements;
    }
    Uf = myKnots->Value(myIndex);
    Ul = myKnots->Value(myIndex+1);
    Denom = 1. /(Ul-Uf);
    USum = Uf+Ul;
    myPtr = (myIndex - 1)*(myBase->WorkDegree()+ 1)*myDimension + 1;
  }

  deg = myDegree(myIndex);
  if (!HasPoly.Value(myIndex)) Update(myIndex, 0);

//Parameter normalization: S [-1, 1]
  S = (2*U - USum)*Denom;
  PLib::NoDerivativeEvalPolynomial(S, deg, myDimension, deg*myDimension, 
				   myPoly(myPtr), Pnt(Pnt.Lower()));
}


//=======================================================================
//function : D1
//purpose  :
//=======================================================================
 void FEmTool_Curve::D1(const Standard_Real U, TColStd_Array1OfReal& Vec) 
{
  Standard_Integer deg, i;
  Standard_Real S;

  if (!myIndex || (U<Uf) || (U>Ul) ||
      (myKnots->Value(myIndex)!=Uf) ||
      (myKnots->Value(myIndex+1)!=Ul)) { 
    // Search the span
    if(U <= myKnots->Value(2)) myIndex = 1;
    else {
      for(myIndex = 2; myIndex <= myNbElements; myIndex++)
	if(U >= myKnots->Value(myIndex) && U <= myKnots->Value(myIndex+1)) break;
      if (myIndex > myNbElements) myIndex = myNbElements;
    }
    Uf = myKnots->Value(myIndex);
    Ul = myKnots->Value(myIndex+1);
    Denom = 1. /(Ul-Uf);
    USum = Uf+Ul;
    myPtr = (myIndex - 1)*(myBase->WorkDegree()+ 1)*myDimension + 1;
  }

  deg = myDegree(myIndex);
  if (!HasDeri.Value(myIndex)) Update(myIndex, 1);
  
//Parameter normalization: S [-1, 1]
  S = (2*U - USum) * Denom;
  PLib::NoDerivativeEvalPolynomial(S, deg-1, myDimension, (deg-1)*myDimension, 
				   myDeri(1+(myIndex - 1)*myBase->WorkDegree()*
					  myDimension), 
				   Vec(Vec.Lower()));
  
  S = 2*Denom;
  for(i=Vec.Lower(); i<=Vec.Upper(); i++) Vec(i) *= S;
}



//=======================================================================
//function : D2
//purpose  :
//=======================================================================
 void FEmTool_Curve::D2(const Standard_Real U, TColStd_Array1OfReal& Vec) 
{
  Standard_Integer deg, i;
  Standard_Real S;

  if (!myIndex || (U<Uf) || (U>Ul) ||
      (myKnots->Value(myIndex)!=Uf) ||
      (myKnots->Value(myIndex+1)!=Ul)) { 
    // Search the span
    if(U <= myKnots->Value(2)) myIndex = 1;
    else {
      for(myIndex = 2; myIndex <= myNbElements; myIndex++)
	if(U >= myKnots->Value(myIndex) && U <= myKnots->Value(myIndex+1)) break;
      if (myIndex > myNbElements) myIndex = myNbElements;
    }
    Uf = myKnots->Value(myIndex);
    Ul = myKnots->Value(myIndex+1);
    Denom = 1. /(Ul-Uf);
    USum = Uf+Ul;
    myPtr = (myIndex - 1)*(myBase->WorkDegree()+ 1)*myDimension + 1;
  }

  deg = myDegree(myIndex);
  if (!HasSecn.Value(myIndex)) Update(myIndex, 2); 

//Parameter normalization: S [-1, 1]
  S = (2*U - USum)*Denom;
  PLib::NoDerivativeEvalPolynomial(S, deg-2, myDimension, (deg-2)*myDimension, 
				   myDsecn(1+(myIndex - 1)*
					   (myBase->WorkDegree() - 1)*myDimension), 
				   Vec(Vec.Lower()));

  S = 4*Denom*Denom;
  for(i=Vec.Lower(); i<=Vec.Upper(); i++) Vec(i) *= S;
}



//=======================================================================
//function : Length
//purpose  :
//=======================================================================
 void FEmTool_Curve::Length(const Standard_Real FirstU,
			    const Standard_Real LastU,
			    Standard_Real& Length) 
{
  Standard_Integer Low, High, deg, degBase, i, Ptr;
  if(FirstU > LastU) Standard_OutOfRange::Raise("FEmTool_Curve::Length");  
  
  if(myKnots->Value(1) > FirstU) Low = 1;
  else
    for(Low = 1; Low <= myNbElements; Low++)
      if(FirstU >= myKnots->Value(Low) && FirstU <= myKnots->Value(Low+1)) break;
  if(Low > myNbElements) Low = myNbElements;
  
  if(myKnots->Value(1) > LastU) High = 1;
  else
    for(High = Low; High <= myNbElements; High++)
      if(LastU >= myKnots->Value(High) && LastU <= myKnots->Value(High+1)) break;
  if(myKnots->Value(myNbElements + 1) < LastU) High = myNbElements;
  
  Standard_Real Li;
  degBase = myBase->WorkDegree();
  Length = 0;

  Standard_Real FirstS, LastS;
  FirstS = (2*FirstU - myKnots->Value(Low) - myKnots->Value(Low + 1))/
    (myKnots->Value(Low + 1) - myKnots->Value(Low));
  LastS = (2*LastU - myKnots->Value(High) - myKnots->Value(High + 1))/
    (myKnots->Value(High + 1) - myKnots->Value(High));
  
  if(Low == High) {
    
    Ptr = (Low - 1)*(degBase + 1)*myDimension + 1;
    deg = myDegree(Low);
    
    if(!HasPoly(Low)) Update(Low, 0);
    PLib::EvalLength(deg, myDimension, myPoly(Ptr), 
		     FirstS, LastS, Length);
    return;
  }
  
  deg = myDegree(Low);
  Ptr = (Low - 1)*(degBase + 1)*myDimension + 1;

  if(!HasPoly(Low))  Update(Low, 0);
  if(FirstS < -1.) {
    PLib::EvalLength(deg, myDimension, myPoly(Ptr), FirstS, -1., Li);
    Length += Li;
    if(myLength(Low) < 0.) { 
      PLib::EvalLength(deg, myDimension, myPoly(Ptr), -1., 1.,Li);
      myLength(Low) = Li;
    }
    Length += myLength(Low);
  }
  else {
    PLib::EvalLength(deg, myDimension, myPoly(Ptr), FirstS, 1., Li);
    Length += Li;
  }

  deg = myDegree(High);
  Ptr = (High - 1)*(degBase + 1)*myDimension + 1;

  if(!HasPoly(High)) Update(High, 0);
  if(LastS > 1.) {
    PLib::EvalLength(deg, myDimension, myPoly(Ptr), 1., LastS, Li);
    Length += Li;
    if(myLength(High) < 0.) { 
      PLib::EvalLength(deg, myDimension, myPoly(Ptr), -1., 1., Li);
      myLength(High) = Li;
    }
    Length += myLength(High);
  }
  else {
    PLib::EvalLength(deg, myDimension, myPoly(Ptr), -1., LastS, Li);
    Length += Li;
  }

  
  for(i = Low + 1; i < High; i++) {
    if (myLength.Value(i) < 0) {
      Ptr = (i - 1)*(degBase + 1)*myDimension + 1;
      deg = myDegree(i);
      if(!HasPoly(i)) Update(i, 0);
      PLib::EvalLength(deg, myDimension, myPoly(Ptr), -1., 1., Li);
      myLength(i) = Li; 
    }
    Length += myLength(i);
  }
}

 Standard_Integer FEmTool_Curve::NbElements() const
{
  return myNbElements;
}

 Standard_Integer FEmTool_Curve::Dimension() const
{
  return myDimension;
}

 Handle(PLib_Base) FEmTool_Curve::Base() const
{
  return myBase;
}

 Standard_Integer FEmTool_Curve::Degree(const Standard_Integer IndexOfElement) const
{
  return myDegree.Value(IndexOfElement);
}


//=======================================================================
//function : SetDegree
//purpose  :
//=======================================================================
 void FEmTool_Curve::SetDegree(const Standard_Integer IndexOfElement,
			       const Standard_Integer Degree)
{
  if (Degree <= myBase->WorkDegree()) {
    myDegree(IndexOfElement) = Degree;
    HasPoly(IndexOfElement) = HasDeri(IndexOfElement) = HasSecn(IndexOfElement) = 0;
    myLength(IndexOfElement) = -1;
  }
  else if(Degree > myBase->WorkDegree()) Standard_OutOfRange::Raise("FEmTool_Curve::SetDegree");
}



//=======================================================================
//function : ReduceDegree
//purpose  :
//=======================================================================
 void FEmTool_Curve::ReduceDegree(const Standard_Integer IndexOfElement,
				  const Standard_Real Tol,
				  Standard_Integer& NewDegree,
				  Standard_Real& MaxError)
{
  Standard_Integer deg = myDegree(IndexOfElement);
  
  Standard_Integer Ptr = (IndexOfElement - 1)*
    (myBase->WorkDegree() + 1)*myDimension + 1;
  
  myBase->ReduceDegree(myDimension, deg, Tol, myCoeff.ChangeValue(Ptr), NewDegree, MaxError);
  Handle(PLib_HermitJacobi) myHermitJacobi = (*((Handle(PLib_HermitJacobi)*)&myBase));
  
  NewDegree = Max(NewDegree, 2 * myHermitJacobi->NivConstr() + 1);
  
  if (NewDegree < deg) {
    myDegree(IndexOfElement) = NewDegree;
    HasPoly(IndexOfElement) = HasDeri(IndexOfElement) = HasSecn(IndexOfElement) = 0;
    myLength(IndexOfElement) = -1;
  }
}

//=======================================================================
//function : Update
//purpose  :
//=======================================================================
 void FEmTool_Curve::Update(const Standard_Integer Index,
			    const Standard_Integer Order)
{
  Standard_Integer  degBase = myBase->WorkDegree(), deg = myDegree(Index);

  if (!HasPoly(Index)) {
    Standard_Integer Ptr = (Index - 1)*(degBase + 1)*myDimension + 1;

    TColStd_Array1OfReal Coeff(myPoly.ChangeValue(Ptr), 0, myDimension*(deg + 1) - 1);
    TColStd_Array1OfReal BaseCoeff(myCoeff.ChangeValue(Ptr), 0, myDimension*(deg + 1) - 1);

    myBase->ToCoefficients(myDimension, deg, BaseCoeff, Coeff);
    HasPoly(Index) = 1;
  }
 
  if (Order >=1) {
    Standard_Integer i1 = (Index - 1)*(degBase)*myDimension - myDimension, 
                     i2 = (Index - 1)*(degBase + 1)*myDimension;
    Standard_Integer i,j;
    if (!HasDeri.Value(Index)) {
      for(i = 1; i <= deg; i++) {
	i1 += myDimension; i2 += myDimension;
	for(j = 1; j<= myDimension; j++)
	  myDeri(i1 + j) = i*myPoly.Value(i2 + j);
      }
      HasDeri(Index) = 1;
    }
    if ((Order>=2) &&!HasSecn.Value(Index)) {
      i1 = (Index - 1)*(degBase-1)*myDimension - myDimension, 
      i2 = (Index - 1)*(degBase)*myDimension;
      for(i = 1; i < deg; i++) {
	i1 += myDimension; i2 += myDimension;
	for(j = 1; j<= myDimension; j++)
	  myDsecn(i1 + j) = i*myDeri.Value(i2 + j);
      }
      HasSecn(Index) = 1;
    }  
  }
}